UFDC Home  myUFDC Home  Help 



Full Text  
NUMERICAL AND EXACT DENSITY FUNCTIONAL STUDIES OF LIGHT ATOMS IN STRONG MAGNETIC FIELDS By WUMING ZHU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2005 Copyright 2005 by Wuming Zhu Dedicated to Mom, to Dad, and to my wife. ACKNOWLEDGMENTS First of all, I would like to thank Professor Samuel B. Trickey, my research advisor and committee chair, for the guidance he provided throughout the course of my graduate study at the University of Florida. His patience and constant encouragement are truly appreciated. Besides physics, I have also learned a lot from him about life and language skills which also are indispensable for becoming a successful physicist. I would also like to thank Professor HaiPing Cheng, Professor Jeffrey L. Krause, Professor Susan B. Sinnott, and Professor David B. Tanner for serving in my supervisory committee, and for the guidance and advice they have given me. Professor David A. Micha is acknowledged for his help when I was in his class and for being a substitute committee member in my qualifying exam. My gratitude goes to Dr. John Ashley Alford II for many helpful discussions, and to Dr. Chun Zhang, Dr. LinLin Wang, and Dr. MaoHua Du for their academic and personal help. Besides them, many other friends have also enriched my life in Gainesville. They are Dr. Rongliang Liu, Dr. Linlin Qiu, Dr. Xu Du, Dr. Zhihong Chen, and Dr. Lingyin Zhu, who have moved to other places to advance their academic careers, and Guangyu Sun, Haifeng Pi, Minghan Chen, Yongke Sun, and Hui Xiong, who are continuing to make progress in their Ph.D. research. Many thanks go to the incredible staff at the Department of Physics and at QTP. I especially would like to thank Darlene Latimer, Coralu Clements, and Judy Parker for the assistance they provided during my graduate study. Financial support from NSF grants DMR0218957 and DMR9980015 is acknowledged. Lastly, I thank my parents, who will never read this dissertation but can feel as much as I do about it, for their boundless love. I thank my wife for her special patience and understanding during the days I wrote my dissertation, and for all the wonderful things she brings to me. TABLE OF CONTENTS page A C K N O W L ED G M EN T S ............................... ......... ................................................... iv L IST O F T A B L E S ................. .................................................................. .. viii LIST OF FIGURES ......... ......................... ...... ........ ............ xi AB STRA CT.................................................... ............................. xiii CHAPTER 1 BASICS OF DENSITY FUNCTIONAL THEORY AND CURRENT DENSITY FU N C TIO N A L TH EO R Y ........................................ ........ ................. ..................... Intro du action ............................................ ..... ............................ D ensity Functional Theory ........................................ ................................. 3 F foundations for D F T ............................................................. ....................... 4 T he K ohnSham Schem e ................................................................. ..... ............ 7 Current Density Functional Theory (CDFT)............................................................10 B asic Form ulations ................... ............. ........................... 10 VignaleRasoltGeldart (VRG) Functional .................................... ............... 13 Survey on the Applications of CDFT....................................... ....... ............ 15 Other D evelopm ents in CD FT....................................... .......................... 16 2 ATOMS IN UNIFORM MAGNETIC FIELDS THEORY .................................18 Single P article E qu ation s ................................................................. ..................... 18 H artreeFock A pproxim ation ........................................ .........................18 Simple DFT Approximation...................................................... 19 C D F T A pproxim ation ........................................... ........................................ 19 Exchangecorrelation Potentials ...........................................................................20 3 BASIS SET AND BASIS SET OPTIMIZATION...............................................25 Survey of Basis Sets Used in Other Work..... ........................ ..............25 SphericalGTO and AnisotropicGTO Representations.................. .. ............. 27 Spherical G TO Basis Set Expansion ............................. ............................... 27 Anisotropic GTO (AGTO) Basis Set Expansion.....................................29 Connection between GTOs and AGTOs .................................. ............... 30 Primary and Secondary Sequences in AGTO............... ...........................................31 O ptim ized A G T O B asis Sets ........................................................... .....................33 4 ATOMS IN UNIFORM MAGNETIC FIELDS NUMERICAL RESULTS .........48 C om prison w ith D ata in L literature ................................................ .....................48 Magnetic Field Induced Ground State Transitions.............................................51 Atom ic Density Profile as a Function of B ................... ............... .................. 53 Total Atomic Energies and Their Exchange and Correlation Components ..............55 Ionization Energies and Highest Occupied Orbitals for Magnetized Atoms ............61 Current Density Correction and Other Issues...... ............................65 5 HOOKE'S ATOM AS AN INSTRUCTIVE MODEL............... .............. 71 H ooke's Atom in V anishing B Field ............................................... ............... 71 Hooke's Atom in B Field, Analytical Solution ......................................................75 Hooke's Atom in B Field, Numerical Solution .................................. ............... 84 Phase Diagram for Hooke's Atom in B Field................. ........ ......... .......... 89 Electron Density and Paramagnetic Current Density ...........................................91 Construction of KohnSham Orbitals from Densities ............................................. 95 Exact DFT/CDFT Energy Components and Exchangecorrelation Potentials...........97 Comparison of Exact and Approximate Functionals..........................................102 6 SUMM ARY AND CONCLUSION .................................................. .............. 110 APPENDIX A HAMILTONIAN AND MATRIX ELEMENTS IN SPHERICAL GAUSSIAN B A S IS .............................................................................................1 12 B ATOMIC ENERGIES FOR ATOMS He, Li, Be, B, C AND THEIR POSITIVE IONS Li+, Be+, B IN MAGNETIC FIELDS ............................... ....................115 C EXCHANGE AND CORRELATION ENERGIES OF ATOMS He, Li, Be, and POSITIVE IONS Li+, Be IN MAGNETIC FIELDS ...........................................132 D EFFECTIVE POTENTIAL INTEGRALS WITH RESPECT TO LANDAU ORBITALS IN EQUATION (5.30) ..................... .................................... 140 E ENERGY VALUES FOR HOOKE'S ATOM IN MAGNETIC FIELDS ...............143 L IST O F R E F E R E N C E S ...................................................................... ..................... 154 BIOGRAPHICAL SKETCH ............................................................. ............... 160 LIST OF TABLES Table p 31 Basis set effect on the HF energies of the H and C atoms with B = 0.........................35 32 Basis set errors for the ground state energy of the H atom in B = 10 au ...................36 33 Optimized basis set and expansion coefficients for the wavefunction of the H atom in B = 10 au. .....................................................................37 34 Test of basis sets including 1, 2, and 3 sequences on the energies of the H atom in B fi e ld s ........................................................................ 4 2 35 Energies for high angular momentum states of the H atom in B fields....................43 36 Basis sets for the H atom in B fields with accuracy of 1 ,H ................... ..............44 37 Basis set effect on the HF energies of the C atom in B = 10 au............. ..................45 38 Construction of the AGTO basis set for the C atom in B = 10 au.............................47 39 Overlaps between HF orbitals for the C atom in B = 10 au and hydrogenlike sy stem s in the sam e fi eld ............................................................... .....................47 41 Atomic ionization energies in magnetic fields ................................. ..................... 62 42 Eigenvalues for the highest occupied orbitals of magnetized atoms...........................63 43 CDFT corrections to LDA results within VRG approximation ................................68 44 Effect of cutoff function on CDFT corrections for the He atom ls2p.l state in m magnetic field B = 1 au ................................................. ............................... 69 51 Confinement frequencies co for HA that have analytical solutions to eqn. (5.5) ........74 52 Confinement frequencies which have analytical solutions to eqn.(5.12) ..................82 53 Som e solutions to eqn. (5.12) ...................................................................... 84 54 Field strengths for configuration changes for the ground states of HA ....................... 89 55 SCF results for HF and approximate DFT functionals.................... ............... 109 Bl Atomic energies of the He atom in B fields............... ........................................... 115 B2 Atomic energies of the Li ion in B fields .......................................... ............121 B3 Atomic energies of the Li atom in B fields ............................... .... ...........122 B4 Atomic energies of the Be+ ion in B fields ............. ............................................. 124 B5 Atomic energies of the Be atom in B fields ............. ............................................. 125 B6 Atomic energies of the B+ ion in B fields. ...................................... ............126 B7 Atomic energies of the B atom in B fields............. ................... .... ...........128 B8 Atomic energies of the C atom in B fields......... .............................................130 C1 Exchange and correlation energies of the He atom in magnetic fields.....................132 C2 Exchange and correlation energies of the Li ion in magnetic fields .....................135 C3 Exchange and correlation energies of the Li atom in magnetic fields...................... 136 C4 Exchange and correlation energies of the Be+ ion in magnetic fields ......................137 C5 Exchange and correlation energies of the Be atom in magnetic fields.....................138 D1 Expressions for Vs(z) with z1<2 and 0 El Relative motion and spin energies for the HA in B fields (0 = 1/2)........................143 E2 A s in Table El, but for = 1/10. ............................. ..... ............................. 145 E3 Contributions to the total energy for the HA in Zero B field (B = 0, m = 0) ..........147 E4 Contributions to the total energy for the HA in B fields (co = 1, m = 0, singlet) ..147 E5 Contributions to the total energy for the HA in B fields (co = 1/10, m = 0, singlet)148 E6 Contributions to the total energy for the HA in B fields (c = 1/, m = 1, triplet) ...149 E7 Contributions to the total energy for the HA in B fields (c = 1/10, m = 1, triplet)150 E8 Exact and approximate XC energies for the HA in Zero B field (B = 0, m = 0, singlet) .......................... ....... ........ ................ .............. 151 E9 Exact and approximate XC energies for the HA in B fields (co =1/2, m = 0, singlet)151 E10 Exact and approximate XC energies for the HA in B fields (co = 1/10, m = 0, singlet) ............... ........ ........ ....................................... 152 E11 Exact and approximate XC energies for the HA in B fields (wc = 1/2, m = 1, triplet) ............................................................... .... ...... ......... 152 E12 Exact and approximate XC energies for the HA in B fields (co = 1/10, m=l, triplet) ............................................................... .... ...... ......... 153 LIST OF FIGURES Figure pge 31 Exponents of optimized basis sets for the H, He+, Li Be C5 and 0+ in reduced magnetic fields y = 0.1, 1, 10, and 100.................................................40 32 Fitting the parameter b(y =1) using the function (3.26)....................................41 41 UHF total energies for different electronic states of the He atom in B fields ............52 42 Crosssectional view of the HF total electron densities of the He atom ls2 and ls2p.l states as a function of magnetic field strength.................. ............. 54 43 Differences of the HF and DFT total atomic energies of the He atom Is2, ls2po, and s2p. states with respect to the corresponding CI energies as functions of B field strength ......... .............. ...... ....................................... ................... .......... 56 44 Differences of DFT exchange, correlation, and exchangecorrelation energies with HF ones, for the H e atom in B fields. .................................................................... 58 45 Atomic ground state ionization energies with increasing B field.............................64 46 Various quantities for the helium atom ls2p.1 state in B = 1 au........................... 66 51 Confinement strengths subject to analytical solution to eqn. (5.12) .........................83 52 Phase diagram for the HA in B fields......................................................... .. ........ 90 53 Crosssectional view of the electron density and paramagnetic current density for the ground state HA with co = 1/10 in B = 0.346 au............................................94 54 Energy components of HA with B = 0. ................................... .......... ....... ........ 99 55 Comparison of exact and approximate XC functionals for the HA with different confinement frequency wc in vanishing B field (B = 0) .......................................103 56 Comparison of exact and approximate exchange, correlation, and XC energies of the H A w ith co = 1/2 in B fields...................................... ............................ 105 57 Sam e as Fig. 56, except for co = 1/10. ................... ............................................. 106 58 Crosssectional views of the exact and approximate XC potentials for the ground state HA with co = 1/10 in B = 0.346 au............... ....................... ............... 107 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy NUMERICAL AND EXACT DENSITY FUNCTIONAL STUDIES OF LIGHT ATOMS IN STRONG MAGNETIC FIELDS By Wuming Zhu August 2005 Chair: Samuel B. Trickey Major Department: Physics Although current density functional theory (CDFT) was proposed almost two decades ago, rather little progress has been made in development and application of this theory, in contrast to many successful applications that ordinary density functional theory (DFT) has enjoyed. In parallel with early DFT exploration, we have made extensive studies on atomlike systems in an external magnetic field. The objectives are to advance our comparative understanding of the DFT and CDFT descriptions of such systems. A subsidiary objective is to provide extensive data on light atoms in high fields, notably those of astrophysical interest. To address the cylindrical symmetry induced by the external field, an efficient, systematic way to construct high quality basis sets within anisotropic Gaussians is provided. Using such basis sets, we did extensive HartreeFock and DFT calculations on helium through carbon atoms in a wide range of B fields. The applicability and limitations of modern DFT and CDFT functionals for atomic systems in such fields is analyzed. An exactly soluble twoelectron model system, Hooke's atom (HA), is studied in detail. Analogously with known results for zero field, we developed exact analytical solutions for some specific confinement and field strengths. Exact DFT and CDFT quantities for the HA in B fields, specifically exchange and correlation functionals, were obtained and compared with results from approximate functionals. Major qualitative differences were identified. A major overall conclusion of the work is that the vorticity variable, introduced in CDFT to ensure gauge invariance, is rather difficult to handle computationally. The difficulty is severe enough to suggest that it might be profitable to seek an alternative gaugeinvariant formulation of the currentdependence in DFT. CHAPTER 1 BASICS OF DENSITY FUNCTIONAL THEORY AND CURRENT DENSITY FUNCTIONAL THEORY Introduction Ambient and lowtemperature properties of normal bulk materials are largely determined by knowledge of the motion of the nuclei in the field of the electrons. In essence, this is a statement that the BomOppenheimer approximation [1] is widely relevant. For materials drawn from the lighter elements of the periodic table, the electrons even can be treated nonrelativistically [2]. While doing some electronic structure calculations on aquartz [3], and some classical internuclear potential molecular dynamics (MD) simulations on silicalike nanorods [4], a feature of modern computational materials physics became obvious. Very little is done with external magnetic fields. This scarcity seems like a missed opportunity. Even with no external field, within the BomOppenheimer approximation, a non relativistic approach to solution of the Nelectron Schrodinger equation is not a trivial task. For simple systems, e.g. the He atom, highly accurate approximate variational wavefunctions exist [5], but these are too complicated to extend. Much of the work of modern quantum chemistry involves extremely sophisticated sequences of approximations to the exact system wavefunction [6]. The HartreeFock (HF) approximation, which uses a single Slater determinant as the approximation to the many electron wavefunction, usually constitutes the first step toward a more accurate, sophisticated method. Several approaches, such as configuration interaction (CI), many body perturbation theory (MBPT), and coupled cluster (CC), are widely used in practice to improve HF results. It is worthwhile mentioning that such methods are extremely demanding computationally. Their computational cost scales as some high power of the number of electrons, typically 57th power. Thus these methods are only affordable for systems having up to tens of electrons. An external magnetic field which could not be treated perturbatively would make things much worse. The largest system that has been investigated with the full CI method as of today is a fourelectron system, beryllium atom [7, 8]. On the other hand, people always have interests in larger systems and more accurate results than those achievable, no matter how fast and how powerful the computers are; thus theorists continue to conceive all kinds of clever approximations and theories to cope with this problem. Density functional theory (DFT) [9, 10] is an alternative approach to the many electron problem that avoids explicit contact with the Nelectron wavefunction. DFT developed mostly in the materials physics community until the early 1990s when it reappeared in the quantum chemistry community as a result of the success of new approximate functionals. These aspects will be discussed below. Two other aspects are worth emphasizing. DFT has been remarkably successful in predicting and interpreting materials properties. Almost none of those predictions involve an external magnetic field. Particularly in Florida, with the National High Magnetic Field Laboratory, that is striking. Even for very simple atoms, inclusion of an external B field is not easy. Only recently have the calculations on the helium atom in a high field been pushed beyond the HF approximation [11, 12, 13]. Although a version of DFT called current density functional theory (CDFT) [14, 15, 16] exists for external magnetic fields, it has seen very little application or development. As discussed below, there is a lack of good approximate CDFT functionals and a lack of studies on which to try to build such improved functionals. One of the foundations of the success of ordinary DFT has been the availability of exact analytical and highly precise numerical data for atoms for comparison of various functionals and understanding their behavior. The main purpose of this dissertation is to find how the effect of an external magnetic field on electron motion should be incorporated in the DFT functional. In particular, I obtain numerical results on various atomlike systems in an external field, with and without CDFT approximate functionals. In addition, I give exact solutions for a model twoelectron atom in a nonzero external B field, the socalled Hooke's atom (HA), that has provided valuable insight for DFT at B = 0. Density Functional Theory Attempts to avoid calculation of the manyelectron wavefunction began almost simultaneously with the emergence of modem quantum mechanics. In 1927, Thomas and Fermi proposed a model in which the electron kinetic energy is expressed as a functional of the electron density, totally neglecting exchange and correlation effects [17]. The kinetic energy density is assumed to be solely determined by the electron density at that point, and approximated by the kinetic energy density of a noninteracting uniform electron gas having the same density. Later this approach was called the "local density approximation" (LDA) in DFT. The ThomasFermi (TF) model was refined subsequently by Dirac to include exchange effects (antisymmetry of identicalparticle wavefunction) and by von Weizsacker to include spatial gradient corrections for the kinetic energy. The result is called the TFDW model. Though useful, it fails as a candidate for a model of materials behavior. Teller proved that the model will not provide binding even for a simple diatomic molecule [17]. The modern form of DFT is rooted in the 1964 paper of Hohenberg and Kohn [9] which put forth two basic theorems, and the subsequent paper by Kohn and Sham [10], which gave an ingenious scheme for the use of those theorems. A difficulty with the KS scheme is that it lumps all of the subtlety of the manyelectron problem, exchange and correlation, in one approximation. The popularity of DFT depends on the availability of reasonably accurate, tractable approximate functionals. To make the point clear and establish notation, next I give the bare essentials of ordinary DFT. Foundations for DFT The Hamiltonian of an interacting Nelectron system is 1 i 1 ( 1 HV2+Zv(i)+ (1.1) 2 ,1 1=1 2 ',,=1 ~' i where i labels the space coordinate of the ith electron. Hartree atomic units are used throughout. The Schr6dinger equation specifies the map from the external potential v (i) to the ground state manybody wavefunction, and the electron number density can be obtained by integrating out N1 space variables. Schematically, v( (fl, ,..., N > n() (1.2) Hohenberg and Kohn noticed that the inversion of the above maps is also true [9], even though it is not as obvious as above. Because of the key importance of this observation, their proof is included here. For simplicity, they considered the spin independent, nondegenerate ground states. Let the Hamiltonian H, ground state wavefunction Y, density n(F), and energy E associated with the specific external potential v(F), v(F) : H, T, n (), E. (1.3) Similarly define a primed system, v'(f) :  ',n'( ), E'. (1.4) where v(F) v'(F)+ C, and hence Y T Y'. By the variational principle, E=( H ) Interchanging the primed and unprimed systems gives us another inequality. Summation of those two inequalities leads to a contradiction E + E' < E + E' if we assume n'() = n(F) Thus, different potentials must generate different ground state electron densities. Equivalently speaking, the knowledge of the ground state density n(f) uniquely determines the external potential v(r) up to a physically irrelevant addictive constant. This assertion is referred as HohenbergKohn (HK) theorem I. Now the maps in eqn. (1.2) are both bijective, v(f) (f, ,..f.,F) n(f) (1.6) An immediate consequence of HK theorem I is that the ground state electron density n(f) can be chosen as the basic variable to describe the interacting Nelectron system, since it is as good as the manybody wavefunction. Here "as good as" means that the ground state density n(r) contains no more or less information about the system than the wavefunction does. It does not mean the density is a variable as easy as, or as hard as, the manybody wavefunction to handle. Actually, since the density is a 3dimensional physically observable variable, whereas the spatial part of wavefunction is a 3N dimensional variable, the density is a much simpler variable to manipulate and to think about. On the other hand, by switching from the wavefunction to the density, we also lose some tools from quantum mechanics (QM) which we are quite adept at using. For example, in QM, an observable can be calculated by evaluating the expectation value of its corresponding operator. This approach often does not work in DFT. The best we can say is that the observable is a functional of the ground state density. In contrast to the explicit dependence on the wavefunction in the QM formulas for the expectation value, the implicit dependence on the ground state density in DFT is rarely expressible in a form useful for calculation. Most such functionals are not known as of today. Among them, the most exploited and the most successful one is the exchangecorrelation energy functional, which is amenable to approximations for large varieties of systems. Another one being extensively studied but not so successfully is the kinetic energy functional, already mentioned in the paragraph about TFtype models. While DFT is a whole new theory that does not need to resort to the manybody wavefunction, to make use of the RayleighRitz variational principle to find the ground state energy Eo, we retain that concept for a while. By the socalled constrained search scheme independently given by Lieb [18] and Levy [19] in 1982, all the trial wavefunctions are sorted into classes according to the densities n'(F) to which they give rise. The minimization is split into two steps, E, = min Yf' Hi '') = min f v(r)n' (i)dr + FL [n'(r)]} (1.7) where FLL[n'] min(Y T+U Y' (1.8) LL J X'T' + /T) 18 The LiebLevy functional FLL is defined on all the possible densities realizable from some antisymmetric, normalized Nparticle functions, or N representable densities. Both the densities of degenerate states and even excited states are included. One good thing about FLL is that we have a simple criterion for N representable densities: all non negative, integrable densities are N representable. The KohnSham Scheme The HK theorem showed that the ground state energy of a manyelectron system can be obtained by minimizing the energy functional E[n'(f)]. TFtype models constitute a direct approach to attack this problem, in which energy functionals are constructed as explicit approximate forms dependent upon the electron density. However, the accuracy of TFtype models is far from acceptable in most applications, and there are seemingly insurmountable difficulties to improve those models significantly. The reason is that the kinetic energy functionals in TFtype models bring in a large error. To circumvent this difficulty, an ingenious indirect approach to the kinetic energy was invented by Kohn and Sham [10]. A fictitious noninteracting system having the same ground state electron density as the one under study is introduced. Because the kinetic energy of this KS system, Ts, can be calculated exactly, and because T, includes almost all the true kinetic energy T, the dominant part of the error in TFtype models is eliminated. Since then, DFT has become a practical tool for realistic calculations. It is advantageous to decompose the total energy in the following way, E[n] = F, [n(r) + v(r)n(r)dr = T[n] + [n] + Ju()n(i)d = T[n] + J +u(F)n(F)c + + V, [n] + (T[n] T[n]) (1.9) where Eee is the total QM electronelectron interaction energy, J j= n(fn(' dFdF' is 2 i r r' the classical electronelectron repulsion energy, E is the conventional exchange energy, and Vc is the conventional correlation energy. Ts is defined in terms of the non interacting system as usual, T = v22 (1.10) Then EHF is replaced by Ex[n], the exchange energy calculated using the single determinant HF formula but with the same orbitals in Ts, and Ec is defined as all the remaining energy Ek [n]= V=[n]+(T[n] T [I])+(E E E [n]) (1.11) The total energy is finally expressed as E[n] = T [n]+ J+ f u()n(F)dF + E [n] (1.12) where E[n] = E[n] + E[n]. (1.13) Each of the first three terms in eqn. (1.12) usually makes a large contribution to the total energy, but they all can be calculated exactly. The remainder, Ec, is normally a small fraction of the total energy and is more amenable to approximation than the kinetic energy. Even though the equations in this section are all exact, approximations to the E, functional ultimately must be introduced. The variational principle leads to the socalled KohnSham selfconsistent equations, v + (r (2 )()) (1.14) where v (t ) v(0) + vH (0) + v~ (i) (1.15) where )+ (r) (1.15) vH ()= n(r' d (1.16) Sr' v(f) En()] (1.17) dn(r) N and n(F) = (F)2. (1.18) i=1 Again, Hartree atomic units are used throughout. Equations (1.14) (1.18) constitute the basic formulas for KS calculations. The detailed derivation and elaborations can be found in the abundant literature on DFT, for example, references 2022. Since the foundations of DFT were established, there have been many generalizations to this theory. The most important include spin density functional theory (SDFT) [23], DFT for multicomponent systems [24, 25], thermal DFT for finite temperature ensembles [26], DFT of excited states, superconductors [27], relativistic electrons [28], timedependent density functional theory (TDDFT) [29], and current density functional theory (CDFT) for systems with external magnetic fields [1416]. Among them, SDFT is the most welldeveloped and successful one. TDDFT has attracted much attention in recent years and shows great promise. Compared to the thousands of papers published on DFT and SDFT, we have fewer than 80 papers on CDFT in any form. Thus CDFT seems to be the least developed DFT generalization, perhaps surprisingly since there is a great deal of experimental work on systems in external B fields. That disparity is the underlying motivation for this thesis. Current Density Functional Theory (CDFT) Basic Formulations One of the striking features of the very limited CDFT literature is the extremely restricted choice of functionals. A second striking feature is that most of the work using CDFT has been at B= 0, in essence using CDFT either to gain access to magnetic susceptibility [30, 31] or to provide a richer parameterization of the B = 0 ground state than that provided by SDFT [32]. In order to comprehend the challenge it is first necessary to outline the essentials of CDFT. For an interacting Nelectron system under both a scalar potential v(r) and a vector potential A (f), its Hamiltonian reads, S ( ++A(I))2+ ) (1.19) 2 2, , The paramagnetic current density jP (f) is the expectation value of the corresponding operator Ji(r), p V/ ) : (ri, r2, f, ),)l  (v) I  (1,f2,", )) (1.20) 1(1.20) where JJ(F)= I (i)V(7) q()V+ ()] (1.21) in terms of the usual fermion field operators. CDFT is an extension of DFT to include the vector potential A () The original papers [1416] followed the HK argument by contradiction and purported thereby to prove not only that the ground state is uniquely parameterized by the density n (F) and paramagnetic current density jp (f), but also that the vt(f) and A(F) are uniquely determined. Later it was found that a subtlety was overlooked. It is obvious that two Hamiltonians with different scalar external potentials cannot even have a common eigenstate, e.g., the first map in eqn. (1.6) is bijective, but this is not true when a vector potential is introduced. It is possible that two sets of potentials u(F), A(F) and v'(F), A' () could have the same ground state wavefunction. This nonuniqueness was later realized [33]. Fortunately, the HKlike variational principle only needs the oneto one map between ground state wavefunction and densities, without recourse to the system's external potentials being functionals of densities [34]. To avoid the difficulties of representability problems, we follow the LiebLevy constrained search approach. Sort all the trial wavefunctions according to the densities n(F) and jp (F) they would generate. The ground state wavefunction, which generates correct densities, will give the minimum of the total energy. E, = min (' to A 12 (1.22) = min (f) n'(f>)df+ J (f).J ;(F) +F[n (),J (f) ] n'() J(r) (2)f+ i (r),j where FIn',1] min jw' TU f 'f) (1.23) A nonphysical noninteracting KS system is now introduced, which generates correct densities n(r) and j (f). The functional F is customarily decomposed as F n, = T jp] +J + E j] (1.24) The variational principle gives us the selfconsistent equations [I +A(ff) +uvdft (i)) (1.25) where Aef (r) = A(r) + A, (r) (1.26) ff(r) = v[(r+ ()+ (r)] + [2(f) Af()] (1.27) It is easy to see that vucf (r) reduces to v, (r) when we set Axc (r) = 0. Exchange correlation potentials are defined as functional derivatives, SE [n(r, j ()] (.2) Axc (') = (1.28) s j (tr) cft 5E' [""(f),(' n((r) v ()) (1.29) The electron density n(F) can be calculated just as in eqn. (1.18). The paramagnetic current density is constructed from the KS orbitals according to jp(V) = '*()V, () 0' ()V #(rt)] (1.30) 1 S2i (1.30) The total energy expression of the system is Ef = T, + + Ec [n(), j,(F)] + n() v(F) + A + f J ) A(F)d = J + E{f [n(Fr), ()] (i)] f f (i)n(F)r JP (i). A (F)r (1.31) Equation (1.25) can be rewritten as V. 1 2 i 2i (1.32) which is more suitable for application. VignaleRasoltGeldart (VRG) Functional On grounds of gaugeinvariance, Vignale and Rasolt argued that the exchange correlation functional Ef should be expressed as a functional of n(f) and the socalled voticity [1416] (f) Vx (1.33) Following their proposal, if we choose v(f) as the second basic variable in CDFT, Ef [n(i), Ji (l)] = Ef [n(ri), (i)] (1.34) exchangecorrelation potentials can be found from functional derivatives Axc (i) L v I )l (1.35) oftft = ) A () (1.36) n(r) n(r) To make use of the already proven successful DFT functionals, it is useful to separate the exchangecorrelation functional into a currentindependent term and an explicitly currentdependent term, E j [n(ir), v(r)] = Edf [n()] + AE [n(r), v(r)] (1.37) The currentindependent term can be any widely used XC functional, such as LDA, or the generalized gradient approximation (GGA) functional. The currentdependent term is presumably small, and should vanish for zero current system, for example, the ground state of the helium atom. Next we proceed in a slightly different, more general way. By homogeneous scaling of both the n(r) and j (f), Erhard and Gross deduced that the currentdependent exchange functional scales homogeneously as [35] E [ftn[1 ]J= AEcdf[n, j] (1.38) 4. where h is the scaling factor, nz = An and j% = 24 j are scaled charge density and paramagnetic current density, respectively. Assuming that the exchange part dominates the exchangecorrelation energy, a local approximation for the Eft takes the form f [n(i), ()]= Edft[n(j)]+ Jg([n(i), (f)],) i(r) 12 F (1.39) The foregoing expression is derived based on the assumption that v(r) is a basic variable in CDFT. A further (drastic) approximation is to assume g([n(r), (r)], ) = g(n(F)) (1.40) which is done in the VRG approximation. By considering the perturbative energy of a homogeneous electron gas (HEG) in a uniform B field (the question whether the HEG remains uniform after the B field is turned on was not discussed), Vignale and Rasolt[14 16] gave the form for g, kF (n (F)) g(F) =g(n(F))= [ 1] (1.41) 247T2 Zo(n(f)) Here kF is the Fermi momentum, and Z0 are the orbital magnetic susceptibilities for the interacting and noninteracting HEG, respectively. From the tabulated data for 1 < r < 10 in reference 36, Lee, Colwell, and Handy (LCH) obtained a fitted form [31], SLCH = / 0 = (1.0 + 0.028rs)e ... (1.42) where r, =\( (1.43) k accordingly, g, = (LCH 1) (1.44) 24;r Orestes, Marcasso, and Capelle proposed two other fits, both polynomial [37] sOC3 = 0.99560.01254r 0.0002955r2 (1.45) soMc5 = 1.1038 0.4991 /3 + 0.4423 0.06696ur + 0.0008432r2 (1.46) Those fits all give rise to divergence problems in the low density region. A cutoff function needs to be introduced, which will be discussed in the next chapter. Survey on the Applications of CDFT The VRG functional has been applied in the calculation of magnetizabilities [30, 31], nuclear shielding constants [38], and frequencydependent polarizabilities [39, 40] for small molecules, and ionization energies for atoms [37, 41]. In those calculations, the vector potential was treated perturbatively. Fully selfconsistent calculations are still lacking. None of those studies has a conclusive result. The first calculation in Handy's group was plagued with problems arising from an insufficiently large basis set [31]. In their second calculation, they found that the VRG functional would cause divergences and set g(n(r)) = 0 for rs >10. The small VRG contribution was overwhelmed by the limitations of the local density functional [38]. The VRG contribution to the frequency dependent polarizability was also found to be negligible, and several other issues emerge as more important than explicitly including the current density functional [40]. Contrary to the properties of small molecules studied by Handy's group, Orestes et al. found that the current contribution to the atomic ionization energy is nonnegligible, even though use of VRG did not improve the energy systematically [37]. All those investigations are based on the assumption that the VRG functional at least can give the correct order of magnitude of current contributions to the properties under study. Actually this is not guaranteed. The errors from ordinary DFT functionals and from the current part are always intertwined. To see how much the current term contributes, an exact solution is desired. Vignale's group has never done any actual numerical calculation based on the VRG functional. Either for an electronic system [42, 43] or for an electronhole liquid [44, 45], they used Danz and Glasser's approximation [46] for the exchange, and the random phase approximation for the correlation energy, which is known to be problematic in the low density regime. The kinetic energy was approximated from a noninteracting particle model or a TFtype model. Even though correlation effects were included in their formulas, the numerical errors introduced in each part were uncontrollable, and their calculations could only be thought as being very crude at best. This is somewhat inconsistent with invoking CDFT to do a better calculation than the ordinary DFT calculation does. While the fully CDFT calculations on threedimensional (3D) systems are scarce, there are more applications of CDFT to 2D systems. Examples include the 2D Wigner crystal transition [47], quantum dots [48], and quantum rings [49] in a magnetic field. In these cases, the Ecf was interpolated between the zerofield value from the Monte Carlo calculation by Tanatar and Ceperley [50], and the strong field limit [51]. Other Developments in CDFT Some formal properties and virial theorems for CDFT have been derived from density scaling arguments [35, 5254] or density matrix theory [55]. A connection between CDFT and SDFT functionals is also established [56]. Those formal relations could be used as guidance in the construction of CDFT functionals, but as of today, there is no functional derived from them as far as I know. CDFT is also extended to TDDFT in the linear response regime [57], which is called timedependent CDFT (TDCDFT). There is not much connection between TD CDFT [57] and the originally proposed CDFT formulation [1416]. Notice an important change in reference 57, the basic variables are electron density n(F) and physical current density j(f), as opposed to the paramagnetic current density j (f), which is argued in references 1416 to be the basic variable. In TDCDFT, the frequency dependent XC kernel functions are approximated from the HEG [58, 59], and the formalism is used in the calculations of polarizabilities of polymers and optical spectra of group IV semiconductors [60, 61, 62]. TDCDFT has also been extended to weakly disordered systems [63] and solids [64]. CHAPTER 2 ATOMS IN UNIFORM MAGNETIC FIELDS THEORY Single Particle Equations When a uniform external magnetic field B, which we choose along the z direction, is imposed on the central field atom, its symmetry goes over to cylindrical. The Hamiltonian of the system commutes with a rotation operation about the direction of the B field, so the magnetic quantum number m is still a good quantum number. The natural gauge origin for an atomlike system is its center, e.g. the position of its nucleus. In the coulomb gauge, the external vector potential is expressed as 1 A()= Bx (2.1) 2 The total manyelectron Hamiltonian (in Hartree atomic units) then becomes "I 1 _Z B2 B 1+2 1 H= IV + x + y+ n +2n,) + 1 (2.2) 1=1 2 8 2 2, ,1 r where Z is the nuclear charge, M, m,, mi,, are the space coordinate, magnetic quantum number, and spin z component for the ith electron. HartreeFock Approximation In the HartreeFock (HF) approximation, correlation effects among electrons are totally neglected. The simplest case is restricted HartreeFock (RHF), which corresponds to a singledeterminant variational wavefunction of doubly occupied orbitals. We use q, (f)to denote singleparticle orbitals. For spinunrestricted HartreeFock (UHF) theory, different spatial orbitals are assigned to the spinup (a) and spindown (fl) electrons. The resulting singleparticle equation is V2 Z B 2 X 7\ B + vH () + (x2 +2 )+(m, +2m) F (f ) 2 r r L 2 1(2.3) (f&07 HF (iC) HF10 Notice that exchange contributes only for likespin orbitals. Simple DFT Approximation It seems plausible to assume that the Ax contribution to the total energy is small compared to the ordinary DFT Exc. Then the zerothorder approximation to the CDFT exchangecorrelation functional Efj [n, j ] can be taken to be the same form as the XC functional in ordinary DFT, Ex [n], with the current dependence in the XC functional totally neglected. Notice that in this scheme, the interaction between the B field and the orbitals is still partially included. From eqn. (1.28), we see that this approximation amounts to setting the XC vector potential identically zero everywhere, Ac(r) = 0. The corresponding singleparticle KohnSham equation is h, 'q (F) ( = )= (f) (2.4) V2 Z B2 B where / =+vH ( r + 8 2)+x +(m + 2m,) F) (2.5) 2 r 2 The scalar XC potential is defined as in eqn. (1.17). CDFT Approximation In this case, both the density dependence and the current dependence of the XC energy functional Efj [n, jp are included. If we knew the exact form of this functional, this scheme in principal would be an exact theory including all manybody effects. In practice, just as in ordinary DFT, the XC functional must be approximated. Unfortunately very little is known about it. A major theme of this work is to develop systematic knowledge about the exact CDFT functional and the one available generalpurpose approximation, VRG. The CDFT KS equation reads Scf (f) = cd (Fr) (2.6) where kdf =(c(F))+ (r)+ (i.V+V.,) (2.7) and vcu (f) is defined by eqn. (1.29). Notice the last term means V (A,~ c) when the operator is applied to a KS orbital. Exchangecorrelation Potentials For ordinary DFT, both LDA and GGA approximations were implemented. Specifically, the XC functionals include HL [65], VWN [66, 67], PZ [68], PW92 [69], PBE [70], PW91 [71], and BLYP [7275]. jPBE is an extension of the PBE functional that includes a current term [32], but does not treat j, as an independent variable, which means Axc ()= 0. For GGAs, the XC scalar potential is calculated according to df X A[n(F),Vn(F)] OE xGA V V G (2.8) G4n(f) an(F) SVn(Fr) Before considering any specific approximate XC functional in CDFT, we point out several cases for which CDFT should reduce to ordinary DFT. The errors in those DFT calculations are solely introduced by the approximate DFT functionals, not by neglecting the effects of the current. Such systems can provide estimates of the accuracy of DFT functionals. Comparing their residual errors with the errors in corresponding current carrying states can give us some clues about the magnitude of current effects. The ground states of several small atoms have zero angular momentum for sufficiently small external fields. These are the hydrogen atom in an arbitrary field, the helium atom inB < 0.711 au. [76] (1 au. of B field =2.3505 x105 Tesla), the lithium atom inB < 0.1929 au. [77], and the beryllium atom inB < 0.0612 au. [8]. Since their paramagnetic current density j, vanishes everywhere, the proper CDFT and DFT descriptions must coincide. Notice (for future reference) that their density distributions are not necessarily spherically symmetric. This argument also holds for positive ions with four or fewer electrons and any closed shell atom. If we admit the vorticity 9(F) to be one basic variable in CDFT, as proposed by Vignale and Rasolt [1416], there is another kind of system for which the DFT and CDFT descriptions must be identical. As Lee, Handy, and Colwell pointed out [38], for any system that can be described by a single complex wavefunction y/(i), v9() vanishes everywhere. The proof is trivial, v() = Vx i i Vx = VxV In =0 n(F) 2i L// 2i Cases include any single electron system, and the singlet states for twoelectron systems in which the two electrons have the same spatial parts, such as H2 and HeH+ molecules. Notice that the system can have nonvanishing paramagnetic current density, jp (t) 0O. A puzzling implication would seem to be that the choice of parameterization by v9() is not adequate to capture all the physics of imposed B fields. For CDFT calculations, we have mainly investigated the VRG functional already introduced. It is the only explicitly parameterized CDFT functional designed for B > 0 and applicable to 3D systems that we have encountered in the literature: ERG [n(f), i(f)] Jg(n(r))l(rI2dr (2.9) x G 2vdr (2.9) where g(n(F)) and i(f) are defined in (1.41) and (1.33). Substitution of (2.9) into (1.35) gives the expression for the vector XC potential, A(P) = 2 Vx[g(n(r))V()] (2.10) n(7) In actual calculations it generally was necessary to compute the curl in this equation numerically. In CDFT, the scalar potential has two more terms beyond those found in ordinary DFT, namely f (r)= x (r) +dg(n) 2 () (2.11) dn n(i) There are three fits for g (n) to the same set of data tabulated in the range of 1 < r, <10 from random phase approximation (RPA) on the diamagnetic susceptibility of a uniform electron gas [36], namely eqns. (1.42), (1.45) and (1.46). Their derivatives are dgLCH dr, dgLCH dn dn dr, SY e0 042rs 0.042 0.028 1+ 1 (2.12) 3n 24;42 4 r 2_ 9OMC3 i 1(9 Y3 00 0.0002955 (2.13) dn 3n 24;r21 4) r, dgoc r 1 9T 2 53 23 dgn 3n 242 (0.1038r2 +0.3327r3 0.2212r2 +0.0008432) dn 3n 242 4 (2.14) The three fits are very close in the range of 1< r < 10, but differ wildly in other regions due to different chosen functional forms for g(n). They all cause divergence problems in any low density region. Without improving the reliability, precision, and the valid range of the original data set, it seems impossible to improve the quality of the fitted functions. It is desirable to know its behavior in the low density region, especially for finite system calculations, but unfortunately, reference 36 did not give any data for r, >10, nor do we know its asymptotic form. Because dg/dn is required for all r, hence for all n, yet g(n) is undefined for low densities, we must introduce a cutoff function. After some numerical experiment we chose gtf F (c +c2 [sr ac' (2.15) 24)2 where a.,cOf is the cutoff exponent, which determines how fast the function dies out. The two constants cl and c2 are determined by the smooth connection between g(n) and gctoff at the designated cutoff density n cutoff toff (ncutoff LCH OMC(ncu dgcutoff dgLCHIOMC (2.16) cun tonff to In this work, we use ncutof = 0.001ao3, cutoff = 2.0ao', unless other values are explicitly specified. There is an identity about the vector XC potential Axc derived from the VRG functional, 24 Ji, (') Jp ()= n V x [2g(n())(ir)]. j (i)dF = 2g(n(i))i(r). V x ()0 d (2.17) = 2AER [n(f), G (f) Since Ax (i)C j (r) and AERG can be computed independently, this equation can provide a useful check in the code for whether the mesh is adequate and whether numerical accuracy is acceptable. CHAPTER 3 BASIS SET AND BASIS SET OPTIMIZATION Survey of Basis Sets Used in Other Work For numerical calculations, the single particle orbitals in eqn. (2.3), or (2.4), or (2.6), can be represented in several ways. One is straightforward discretization on a mesh. For compatibility with extended system and molecular techniques, however, we here consider basis set expansions. For zero B field, the usual choices are Gaussiantype orbitals (GTO), or, less commonly, Slatertype orbitals (STO). Plane wave basis sets are more commonly seen in calculations on extended systems. Large B fields impose additional demands on the basis set, as discussed below. Here we summarize various basis sets that have been used for direct solution of the fewelectron Schrodinger equation and in variational approaches such as the HF approximation, DFT, etc. For the oneelectron problem, the hydrogen atom in an arbitrary B field, the typical treatment is a mixture of numerical mesh and basis functions. The wavefunction is expanded in spherical harmonics Ym (0, 9) in the low field regime, and in Landau orbitals 4L (p, ) for large B fields. Here r, 0, p are spherical coordinates, and z, p, p are cylindrical coordinates. The radial part (for low B) or the z part (for high B) of the wavefunction is typically represented by numerical values on a onedimensional mesh [78]. In Chapter 5 we will also use this technique for the relative motion part of the Hooke's atom in a B field. Of course, the hydrogen atom has also been solved algebraically, an approach in which the wavefunction takes the form of a polynomial multiplying an exponential. This is by no means a trivial task. To get an accurate description for the wavefunction, the polynomial may have to include thousands of terms, and the recursion relation for the polynomial coefficients is complicated [7982]. The multichannel Landau orbital expansion was also used in DFT calculations on manyelectron atoms [83]. Another approach is the twodimensional finite element method [84]. Dirac exchangeonly or similar functionals were used in those two calculations. In the series of HartreeFock calculations on the atoms hydrogen through neon by Ivanov, and by Ivanov and Schmelcher, the wavefunctions were expressed on twodimensional meshes [8590, 76]. Slatertype orbitals were chosen by Jones, Ortiz, and Ceperley for their HF orbitals to provide the input to quantum Monte Carlo calculations, with the aim to develop XC functionals in the context of CDFT [9193]. Later they found that the STO basis was not sufficient and turned to anisotropic Gaussian type orbitals (AGTO) [94]. Apparently their interests changed since no subsequent publications along this thread were found in the literature. Schmelcher's group also employed AGTOs in their full CI calculations on the helium [1113], lithium [77], and beryllium [8] atoms. At present, AGTOs seem to be the basis set of choice for atomic calculations which span a wide range of field strengths. This basis has the flexibility of adjusting to different field strengths, and the usual advantage of converting the onebody differential eigenvalue problem into a matrix eigenvalue problem. Moreover, the one center coulomb integral can be expressed in a closed form in this basis, though the expression is lengthy [11, 12]. The disadvantage of AGTOs is that one has to optimize their exponents nonlinearly for each value of the B field, which is not an easy task, and a simple, systematic optimization is lacking. We will come back to this issue and prescribe an efficient, systematic procedure. SphericalGTO and AnisotropicGTO Representations As with any finite GTO basis, there is also the improper representation of the nuclear cusp. Given the predominance of GTO basis sets in molecular calculations and the local emphasis on their use in periodic system calculations, this limitation does not seem to be a barrier. Spherical GTOs are most widely used in electronic structure calculations on finite systems without external magnetic field. The periodic system code we use and develop, GTOFF [95], also uses a GTO basis. Several small molecules in high B fields were investigated by Runge and Sabin with relatively small GTO basis sets [96]. To understand the performance of GTOs in nonzero field and make a connection to the code GTOFF, our implementation includes both GTO and AGTO basis sets. The former is, of course, a special case of the latter, in which the exponents in the longitudinal and transverse directions are the same. Spherical GTO Basis Set Expansion The form of spherical Gaussian basis we used is G1 (f)= Nmrlear2Y (0,9) (3.1) where N, is the normalization factor. The KS or Slater orbitals (DFT or HF) are expanded in the G1, (r), w here mI c)= R (r))Y(0, (o )o) (3.2) I a I where R (r) r' a,1 (3.3) (3.3 Notice m is understood as m,, the magnetic quantum number of the ith orbital. For simplicity, the subscript i is omitted when that does not cause confusion. The electron density and its gradient can be evaluated conveniently as n(i) () 2 R, Z (r)Y,(0,(p SI I (3.4) r + ( = 1R,, (r)R,,,(r)Y, (0, (p)r,,(, () vn()R = ( r [R, (r)R (r) + R' (r)R,, ( r)Y, (, )Y,,(, 3p) (3.5) i r 80 The paramagnetic current density is ( r)= 0 mZZ M (r)R, (rfm (0,f m, )= p(r,) (3.6) rsinO and the curl of jp (r) is Vxj (i) = I n +, im Vx ; r2 80 sin0 80 sin8 r sin8 1 ;;M The vorticity is evaluated analytically according to +) V lJp) Vn(_ ) p Vx () v(F)=Vx x v=^ (F)+ n(i) n (i) n(i) (3.8) = r (r, )ri + v (r, 0) For the VRG functional, the vector XC potential is expressed in spherical coordinates as A )= 2 1 [r g(n(r, ))vi (r,)] [ g(n(r,0))v, (r, o)] = ,q (2r,0) The last term in eqn. (2.7) becomes 1 aA (r, ) r sin 9o (3.10) Appendix A includes the matrix elements in this basis for each term in the Hamiltonian. Anisotropic GTO (AGTO) Basis Set Expansion An external B field effectively increases the confinement of the electron motions in the xy plane, and causes an elongation of the electron density distribution along the z axis. It is advantageous to reflect this effect in the basis set by having different decay rates along directions parallel and perpendicular to the B field. AGTOs are devised precisely in this way: X(p,z,p) =Np "z 'e p pe j =1,2,3, (3.11) where = m +2k, with k= 0,1,... m ,2,1,0,1,2,... where with nz = Z, +2/1, = 0,1,. n, = 0,1. and N, is the normalization factor. If we leta, = /,, this basis of course recovers the isotropic Gaussian basis, appropriate for B = 0. The basis sets used in reference 94 were limited only tok, = = 0, which are more restrictive than those used by Becken and Schmelcher [1113] and ours. Singleparticle wavefunctions expanded in AGTOs have the general form 0 (Fr)= C bzX, (p,z,,o) Io) (3.12) Various quantities can be calculated in this basis according to their expressions in cylindrical coordinates, v ,(F) 0 + + im, . V4 (r)=p '+z '+;p = b \p 2a P + (3Z13) 2 =IA P ill a ) yz ) 'P (pXAP]^} n"(F) = ()2 = bV p "zn'e,2/zL =n (pz) (3.14) Vn() an() + an() 0 \ 0, + 8, (3.15) V =n!,ap + = 2 /3 + z (3.15) p Oz p 8z (i) M1 0' I 1[ i) 2 m b Np'z7e JZ,p2_f z2 p (3.16) =p jp(p, z) V x JP (r) = 2.m + z 0 S1 v f 8 P O z Op] (3.17) On .n S z v j ())+ v j ()) p S2 + n i + z .i (3.18) = vP (p, z)p + v, (p, z)2 Again, for the VRG functional we have f)= { [g(n(p, z))v(p, z) [g(n(p,z)(p,z) (3.19) n(p,z) 8z Bp? We follow the scheme in references 1112 for evaluation of matrix elements, in which all the integrals, including Coulomb integrals, are expressed in closed form. Connection between GTOs and AGTOs As pointed out before, a GTO basis is a special case of an AGTO basis. Conversely, a particular AGTO can be expanded in GTOs. X (p, z, (p)= N p' e _,, e  k= k! k =N ( r e ++2k jr2em (sin O)nf (COS O)"+2k (3.20) k=0 k It is easy to see this is a linear combination of Ga' with = n + n + 2k, k = 0, 1, 2, An ordinary contracted Gaussian basis is a fixed linear combination of several primitive Gaussians having same the I and m but different exponents a, Similarly, an AGTO can also be thought as a contracted GTO that contains infinitely many GTOs (in principle) having the same exponent and m but different / values with increment of 2. This establishes the equivalence of the two kinds of orbitals. The relative efficiency of the AGTO basis in cylindrically confined systems is apparent for B 0. Primary and Secondary Sequences in AGTO While the AGTO basis provides extra flexibility, its optimization is more complicated than for a GTO set of comparable size. Kravchenko and Liberman investigated the performance of AGTO basis sets in oneelectron systems, the hydrogen atom and the hydrogen molecular ion, and showed that they could provide accuracy of 106 Hartree or better [97]. Jones, Ortiz, and Ceperley estimated their basis set truncation error for the helium atom in B < 8au. to be less than one milihartree in the total atomic energy of about 2 Hartrees [94]. Eventempered Gaussian (ETG) sequences often are used in zerofield calculations. For a sequence of primitive spherical Gaussians having the same quantum numbers, their exponents are given by a, = ,j = pqJ, j = 1,2,.Nb. (3.21) where p and q are determined by In p = aln(q 1)+ a' (3.22) In(ln q) = b In Nb+ b' and Nb is the basis size. For the hydrogen atom, Schmidt and Ruedenburg [98] recommended the following parameters: a = 0.3243, a'= 3.6920, b = 0.4250, b'= 0.9280. Since the external magnetic field only increases the confinement in the horizontal direction, we may expect eqn. (3.21) to be equally useful for generating longitudinal exponents f/ for the AGTO basis. The choice of a, is more subtle. Jones, Ortiz, and Ceperley [94] used several tempered sequences of the types aj = Pf 2,8, 4,j, 8/,j, ... (3.23) For convenience, we refer to the first sequence (a, = /,) as the primary sequence, and the second, the third, the fourth, ... sequences as the secondary sequences in our discussion. The primary JOC sequence in eqn. (3.23) is obviously as same as the spherical GTOs, for which the transverse and longitudinal exponents are the same. However for the second JOC sequence, the transverse exponents a, 's are twice the longitudinal exponents / 's, and for the third sequence, a, = 4/8 etc. The basis set is the sum of all those sequences. The total number of basis functions is Nb multiplied by the number of different sequences. Reference 94 used 25 sequences in the expansion of HF orbitals. Obviously, when several sequences are included, which is necessary for large B fields, very large basis sets can result. Kravchenko and Liberman [97] chose a, =,j +BAKL,f +1.2BAKL,/ +0.8BAKL, /, +1.4BAKL,/ +0.6BAKL (3.24) where AKL is a value between 0 and 0.3 which minimizes the basis set truncation error compared to more accurate results. Here we still refer to the sequences in eqn. (3.24) as primary and secondary sequences. In each KL sequence, the differences between the transverse and longitudinal exponents are the same for all the basis functions. An improvement of KL basis sets over JOC basis sets is that the former have shorter secondary sequences, which helps to keep basis size within reason. Namely, the second and the third KL sequences have lengths of onehalf of KL primary sequence, e.g., Nb /2 is used in eqns. (3.21) and (3.22) to generate them, and the fourth and fifth KL sequences have lengths of Nb /4. Becken et al. [11] used a seemingly different algorithm to optimize both a, and 8j in the same spirit of minimizing the oneparticle HF energy, H atom or He in a B field, but they did not give enough details for one to repeat their optimization procedure. Optimized AGTO Basis Sets In this section, I give some numerical illustrations of the basis set issues. These examples illustrate the importance, difficulties, and what can be expected from a reasonably welloptimized basis. Our goal is set to reduce basis set error in the total energy of a light atom to below one milihartree. This criterion is based on two considerations. One is the observation by Orestes, Marcasso, and Capelle that the magnitude of current effects in CDFT is of the same order as the accuracy reached by modern DFT functionals [37]. They compared atomic ionization energies from experiment with DFTbased calculations. A typical difference is 0.4eV, or 15 mH. To study the current effect in CDFT, we need to reduce the basis set errors to considerably below this value. Another factor considered is the wellknown standard of chemical accuracy, usually taken to be 1 kcal/mol, or 1.6 mH. It turns out that this goal is much harder to reach for multielectron atoms in a large B field than for the fieldfree case. Two systems I choose for comparison are the hydrogen and carbon atoms. There are extensive tabulations for the magnetized hydrogen atom [78], and even more accurate data from the algebraic method [80] against which to compare. However, the hydrogen atom does not include electronelectron interaction, which is exactly the subject of our interest. For the carbon atom, our comparison mainly will be made with numerical HartreeFock data [90]. Without external field, the correlation energy of the C atom is about 0.15 Hartree [99], two orders of magnitude larger than our goal. This difference also makes the choice of one mH basis set error plausible. Examine the zerofield case first. It is well known that the nonrelativistic energy of the hydrogen atom is exactly 0.5 Hartree. For the carbon atom, the numerical HF data taken from reference 90 are treated as the exact reference. Calculated HF energies in various basis sets are listed in Table 31, together with basis set errors in parentheses. We first tested the widely used 631G basis sets. Those basis sets are obtained from the GAMESS code [100]. In primitive Gaussians, they include up through 4s for the hydrogen, and 10s4p for the carbon atom. As expected, the accuracies in total energy that they deliver increase only slightly after decontraction. A sequence of exponents derived from eqn. (3.21) with length Nb = 8 has a comparable size with the 631G basis for the carbon atom. It gives rather bad results, but recall that a significant deficiency of GTOs is that they cannot describe the nuclear cusp condition. By adding five tighter s orbitals extrapolated from eqn. (3.21) with = 9,10,...,13, the basis set error is reduced by 99%. To further reduce the remaining 1.6 mH error, higher angular momentum orbitals are required. Addition of four d orbitals and removal of the tightest, unnecessary p orbital gives a 13s7p4d basis set, with only 0.8 mH truncation error left. A larger basis set, 20s1 lp6d, similarly constructed from the Nb = 16 sequence by adding 4 tighter s orbitals has error only of 0.05 mH. Table 31 Basis set effect on the HF energies of the H and C atoms with B = 0 (energies in Hartree) Basis Set a Hydrogen atom Carbon atom 631G 0.498233 (0.001767) 37.67784 (0.01312) Decontracted 631G 0.498655 (0.001345) 37.67957 (0.01139) Sequence Nb = 8 0.499974 (0.000026) 37.51166 (0.17930) Nb = 8, plus 5 tighter s 0.499989 (0.000011) 37.68938 (0.00158) 13s7p4d 0.499989 (0.000011) 37.69018 (0.00078) Sequence Nb = 16 0.49999992(0.00000008) 37.68949 (0.00147) 20s1 p6d 0.49999996(0.00000004) 37.69091 (0.00005) oo 0.5 37.69096 b (a) see text for definitions; (b) from reference 90; (c) numbers in parentheses are basis set errors. The situation changes greatly when a substantial external B field is added. Let us first take the example of the H atom ground state in a field B = 10 au. Its energy is known accurately to be 1.747 797 163 714 Hartree [80]. The sequence ofNb = 16 included in Table 31 works remarkably well for the fieldfree energy, but gives 24% error in the B = 10 au field. See Table 32. Adding a sequence of d orbitals that has same length and same exponents as that for the s orbitals, which doubles the basis size, reduces the error by 80%. Further supplementation by g and i orbitals in the same way decreases the error by another order of magnitude. But this is still far from satisfactory. To reduce the basis set error below 1 /H, higher angular momentum basis with I up to 20 must be included. Obviously, this is a very inefficient approach. The basis sets used by Jones, Ortiz, and Ceperley [94] (see eqn. 3.23) converge the total energy more rapidly than these spherical bases. The primary sequence in the JOC basis sets is the same as the spherical basis, but subsequent secondary sequences double the transverse exponents aj's successively. With four sequences the error is less than 1% of the error in a spherical basis set having the same size. Another significant gain can be obtained if we move to the KL basis sets [97] (see eqn. 3.24). Here we choose AKL = 0.18, which is obtained by searching with a step of 0.01 to minimize the basis set truncation error. Including only the primary KL sequence gains about the accuracy of the threesequence JOC basis set. Recall that the subsequent KL secondary sequences have shorter lengths than the primary one (refer to the discussion after eqn. 3.24). Specifically, the second and the third sequences have length of Nb /2 = 8, and the fourth and the fifth sequences have length of Nb /4 = 4. Thus, the basis size will be 16 + 8 + 8 + 4 + 4 = 40 if we include five KL sequences, with accuracy of 1 uH. Table 32 Basis set errors for the ground state energy of the hydrogen atom in B = 10 au. (energies in Hartree) Basis size Spherical JOC a KL b Optimized Eqn. (3.26) 16 0.4198 0.41978728 0.00373820 0.00000060 0.00104451 32 0.081 5 0.027 124 87 0.000 005 39 0.000 00036 0.000 000 50 48 0.021 7 0.001 008 57 64 0.008 1 0.000 075 02 40 0.000 001 12 0.000 000 30 0.000 000 28 (a) JonesOrtizCeperley basis sets, see ref [94] and eqn. (3.23); (b) KravchenkoLiberman basis sets, see ref [97] and eqn. (3.24). AK is chosen to be 0.18. However, this does not mean there is no opportunity left for basis set optimization. Starting from the primary sequence in the KL basis set, we then searched in the parameter space {a,} to minimize the total energy of the H atom. First, the energy gradient in parameter space is calculated, then a walk is made in the steepest descent direction. These two steps are repeated until OE = 0 (3.25) ca] The error left in this optimized basis set is only 0.6 pH, six orders of magnitude smaller than the error of a spherical basis set of the same size! The resulting exponents are listed in Table 33, together with the coefficients used in the wavefunction expansion. Addition of the same secondary KL sequences can further reduce the remaining error by one half. This improvement is not as spectacular as that for the KL basis set because those exponents have already been optimized. It is worth mentioning that, while it is easy to optimize the basis set for the H atom fully, it is very hard to do so for multielectron atoms. We usually only get partially optimized results, but by including secondary sequences, the basis error can be greatly reduced, as demonstrated here. Table 33 Optimized basis set and expansion coefficients for the wavefunction of the hydrogen atom in B = 10 au. j Coefficients a aj j 1 0.000493 1.8886 0.0573 1.8313 2 0.011007 2.8640 0.1247 2.7393 3 0.184818 2.4462 0.2717 2.1745 4 0.372811 2.5541 0.5917 1.9624 5 0.277663 2.6442 1.2890 1.3552 6 0.132857 3.7690 2.8077 0.9613 7 0.050662 6.7855 6.1159 0.6696 8 0.019285 13.9048 13.3221 0.5827 9 0.007139 29.3287 29.0190 0.3097 10 0.002772 64.4702 63.2111 1.2591 11 0.000955 139.3854 137.6904 1.6950 12 0.000418 301.7122 299.9260 1.7862 13 0.000114 655.1166 653.3180 1.7986 14 0.000082 1424.8988 1423.0988 1.8000 15 0.000001 3101.6845 3099.8845 1.8000 16 0.000021 6754.1687 6752.3659 1.8028 From Table 33, we see that the wavefunction is mainly expanded in basis = 3, 4, 5, 6, and a, p, is not a constant as the KL sequences suggest. The smaller the exponent, the larger the difference between the transverse and the longitudinal exponents. This is quite understandable. A smaller exponent means that the electron density extends far from the nucleus, and the magnetic field will overpower the nuclear attraction, thus the distortion from the fieldfree spherical shape will be relatively larger. In the limit of fj 0, which can be equivalently thought of as the large B limit, or zero nuclear charge, the electron wavefunction is a Landau orbital with an exponential parameter B a = aB = The opposite limit, /,* oo, corresponds to B = 0, for which a,= /,. A 4 natural measure of the orbital exponents is aB Now we can make an explicit construction (discussed below) which incorporates all these behaviors, namely B 4 f, 4 , a= + 4 1+4 j2 + 1+[ ] (3.26) S 20 bG () B b (y) B where b(y) =0.16[tan 1()]2 + 0.77 tan 1(7)+ 0.74 (3.27) = pq' j = 1,2, Nb. (3.28) and y = B/Z is the reduced field strength for an effective nuclear charge Z. The parametersp and q are defined in eqn. (3.22). For the innermost electrons, Z is close to the bare nuclear charge; for the outmost electrons in a neutral atom, it is close to unity. Nevertheless, accurate Z values do not need to be provided. Nominal values turn out to be good enough for the input to eqn. (3.26). Secondary sequences are defined similarly as in KL basis set, using a factor of 1.2, 0.8, 1.4, and 0.6 for the second, third, fourth, and fifth sequences, respectively. Next I show how eqn. (3.26) was obtained. Start from the basis set of one sequence Nb = 16. Full basis set optimizations were done for H, He+, Li+, Be+++, C+ and 07+ in reduced fields y = 0.1, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 800, 1000, 2000, and 4000. Results for y = 0.1, 1, 10, and 100 are plotted in Fig. 31. The first observation is that data points for different nuclear charges with the same y are on the same curve. One can show that this must be the case. Suppose the wavefunction for the hydrogen atom in an external B field is WH (, B), + 2 +y2 H (r, B)= EHr H (, B) (3.29) Scaling r Z r leads to V2 1 ZZ2B2 2) 2_ 2 H (Z i, B)= EH, H (Z i, B) 2Z2 Zr 8 V2 Z (Z2B)2 [ 2 + z(X)2( +Y H (Z ,B) =(z2EH (Z i,B) 2 r 8 The Hamiltonian on the left side is the same as that of a hydrogenlike atom with nuclear charge Z in an external fieldB'= Z2B. The scaled hydrogenatom wavefunction is precisely the eigenfunction of this Hamiltonian with energy of Z2EH. Now we expand WH (, B) in the optimized basis set, /IH (, B) = aNj p" zP e'e, p2 #e (3.30) and the scaled wavefunction is 40 VZ (F, ZB)= /H (ZF, B) = aJN. Z" '"%" p z 'e e", (3.31) a, =2 a' 8j aj 8j 8j' 8 where a' = Zaj, f/ = Z2fj. Obviously, a' and B' B B' B Another observation is that the curvatures for different y values are slightly different. To describe the rapid decrease in the small /8 / B region and the slowly decaying long tail, we used a fixed combination of inverse square and inverse square root terms, which proved to be superior to a single reciprocal function. In Fig. 32, the functional forms are compared with data points from optimized basis sets for y= 1. 0.25# 0.2  A 0.15 m CC 0 0.1 E x A + x 0.05 + [] 0 + 0 O + 0 0 0 0 0.5 1 1.5 2 2.5 3 3.5 4 FIG. 31 Exponents of optimized basis sets for the H(+), He+(x), Li++(o), Be +(A), C5+(0) and 07+() in reduced magnetic fields y = 0.1 (blue), 1 (black), and 100 (red). 41 0 optimized for7=l 0.2 Eqn. (3.26) using b=1.286 2 Eqn. (3.26) using b=1.246 C 1/4/(1+4/0.722*l/B) 0.18  i 0.16  0.14  m 0.12  S0.1 o '.\ 0.08 0 .\ 0o.\ 0.06 0.04 .Q >.. ., 0.02 0..... @ " ... o0 0 0.5 1 1.5 2 2.5 3 3.5 4 Pj/B FIG. 32 Fitting the parameter b(y =1) using the function (3.26). Fitted result is 1.286, compare to the calculated value 1.246 from eqn. (3.27). Two curves are shown by dotted black line and solid green line. They are almost indistinguishable in the graph. A reciprocal fitting result B + 4 6, is also shown by a 4 0.722 B dashed blue line. A calculation using the basis set derived from eqn. (3.26) also is included in Table 32. The new primary sequence outperforms the KL primary sequence by a factor of four, but by including only the second and the third sequences, the basis set is almost saturated, compared to other, more slowly converging basis sets. Another advantage of the present basis set is the explicit expressions eqns. (3.26) and (3.27), whereas searching for the best parameter AKL in the KL basis sets is quite time consuming. Table 34 Test of basis sets including 1, 2, and 3 negative signs are omitted) sequences on the energies of the hydrogen atom in B fields (energies in Hartree, State Is 2po 2p_ 3dl 3d2 4f 2 43 5g4 B (au.) 0.1 1 10 100 1000 4000 10000 1 10 1000 1 10 100 1000 10 1000 10 1000 100 10 1000 1 100 1000 Reference 80 0.547526480401 0.831168896733 1.747797163714 3.789804236305 7.662423247755 11.204145206603 0.260006615944 0.382649848306 0.4924950 0.456597058424 1.125422341840 2.634760665299 5.63842108 0.3389561898 0.4869777 0.9082147755 4.80511067 1 sequence 0.5475263 0.83109 1.74675 3.78933 7.66224 11.20372 14.14037 0.25991 0.38263 0.49248 0.45658 1.12521 2.63472 5.63792 0.33890 0.48696 0.90813 4.80432 0.42767 0.78773 4.30738 0.26570 1.75448 3.96338 Reference 76 0.83116892 1.74779718 3.7898043 7.6624234 14.14097 0.260007 0.382650 (a) Numbers before/after slashes are the upper/lower bounds to the energy. 2 sequences 0.5475263 0.8311680 1.747781 3.789790 7.662419 11.204139 14.140959 0.2600055 0.382641 0.4924948 0.456596 1.125415 2.634756 5.638413 0.3389555 0.4869775 0.908212 4.805094 0.427756 0.787768 4.308344 0.266184 1.754848 3.964471 3 sequences 0.547526461 0.83116886 1.74779694 3.78980395 7.66242306 11.20414499 14.14096829 0.26000652 0.38264977 0.49249495 0.45659703 1.12542217 2.63476052 5.63842079 0.33895610 0.48697789 0.90821466 4.80511012 0.42775840 0.78776910 4.30836962 0.26618782 1.75485593 3.96450833 Reference 78a 0.5475265 0.831169 1.747797 3.78905/90250 7.66205/65 0.2600066 0.38264875/5180 0.492495 0.45659705 1.1254225 2.634740/95 5.638405/35 0.33895610/45 0.48697795 0.9082115/235 4.8051095/125 0.4277585 0.7877685/705 4.3083700/05 0.26618770/875 1.754856 3.9645095 The basis set constructed from eqn. (3.26) not only works well for Is electrons, but also for higher excited states. Table 34 includes some test results on the hydrogen atom in a wide range of B fields. The primary sequence was derived from eqns. (3.26) through (3.28) using Nb = 16. Extrapolations toj > Nb orj < 0 were made for extremely tight or diffuse orbitals whenever necessary. The averaged basis set error for the primary sequence is 0.3 mH, which is reduced to 7 /H if the second sequence is added. With three sequences, the accuracy of our basis sets reaches 1 /H level. Notice that energies quoted from reference 76 are slightly lower than the more accurate algebraic results [80]. This implies that Ivanov and Schmelcher's data are not necessarily upper bounds to the energy. We need to keep this in mind when we compare our results with theirs. Table 35 Energies for high angular momentum states of the hydrogen atom in B fields (B in au; energy in Hartree) State B = 1 B= 10 B= 100 B = 1000 5g4 0.2661880 0.7080264 1.7548563 3.9645100 6h_5 0.2421928 0.6499941 1.6252244 3.7061998 7i6 0.2239757 0.6051943 1.5238725 3.5018527 8j7 0.2095131 0.5691841 1.4415788 3.3343126 9k_ 0.1976562 0.5393750 1.3728860 3.1933035 10/19 0.1876974 0.5141408 1.3143222 3.0722218 The accuracy of the previous basis sets can be improved further by increasing Nb and including the fourth and the fifth sequences. Using larger basis sets having five sequences with the primary sequence derived from Nb = 28, we obtained energies of 5.638 421 065, 4.805 110 65, 4.308 370 6, and 3.964 510 0 for the H atom 2p.i, 3d2, 4f3, and 5g4 states in B = 1000, respectively. Their accuracy is at the same level as the best available data in the literature (see Table 34). There are no specific difficulties for the higher angular momentum states in our expansion. The energy values for the excited states with quantum numbers I = m = 4, 5, ..., 9 of the hydrogen atom in B fields are listed in Table 35. Those values will be used in the next step, construction of Table 36. Actually we do not need to use the entire sequences generated from eqns. (3.26) through (3.28). In the expansion of nonzero angular momentum orbitals, very tight basis functions (j close to Nb ) are not necessary, but extrapolation to negative may be required in order to include very diffuse basis functions. Table 36 lists the subsets which have accuracy of 1 uH. The one to five segments in each basis set means the ranges ofj values selected from the primary and the subsequent secondary sequences. Numbers underlined identify the negative values. Nb values used for the five sequences are 16, 8, 8, 4, and 4. Also recall a factor of 1.2, 0.8, 1.4, and 0.6 is used for the second, third, fourth, and fifth sequences, respectively. Numbers in parentheses are the sizes of the basis sets. Table 36 Basis sets for the hydrogen atom in B fields with accuracy of 1 /H (B in au) State B=0 B= B=10 B = 100 B=1000 Is 114(14) 214,13,12(18) 116,26,26(26) 420,28,26(29) 421,410,510(31) 2po 17(9) 08,03,3(14) 010,15(16) 112,25(16) 114,35(17) 2p1 17(9) 19,03,3(14) 111,5,26(17) 314,48,36(21) 416,410,510(26) 3di 34(8) 16,03,03(16) 08,15,23(16) 010,16,2(18) 112,26,3(18) 3d_2 24(7) 07,13,03(15) 18,35,26(16) 312,37,37(20) 415,410,47(23) 4f2 52(8) 15,03,03(15) 07,14,15(17) 010,16,23(19) 112,26,3(18) 4f3 42(7) 04,03,04(14) 17,25,26(16) 311,17,17(23) 314,28,47,3,3(25) 5g_3 40(5) 15,12,01(13) 07,14,1(13) 09,15,23(17) 111,26,34(18) 5g4 40,21(7) 04,03,04(14) 16,25,13(13) 211,17,36(21) 313,38,47,3,3(23) 6h_5 52,41(8) 04,03,03(13) 16,25,13(13) 210,37,25(18) 313,38,47,3,3(23) 7i6 62,53(8) 12,02,03(11) 15,25,13(12) 210,15,36(18) 313,38,47,3,3(23) 8j_7 12,02,03(11) 15,25,13(12) 210,26,3,2,2(17) 313,38,47,3,3(23) 9k8 15,25,12(11) 29,26,3,2,2(16) 312,38,46,3,3(21) 1019 15,25,12(11) 29,26,3,2,2(16) 312,38,46,3,3(21) While the previous optimization scheme is quite impressive for the hydrogen atom in a B field, we want to know whether it also works equivalently well for multielectron atoms. Thus we do another case study, the carbon atom in the same field B = 10 au. Its ground state configuration is s22p 13d 24f_ 5g 4, and HF energy is 44.3872 Hartree from calculations on a numerical grid [90]. The performance of various basis sets is summarized in Table 37. Table 37 Basis set effect on the HF energies of the carbon atom in B = 10 au. Basis set Basis size HF energy (Hartree) Error(Hartree) Spherical(spdfghi) 112 43.6157 0.7715 2sequences, JOC 160 44.1572 0.2300 3 sequences, JOC 240 44.3529 0.0343 sequence, KL 80 44.1629 0.2243 2sequences, KL 120 44.3824 0.0048 3 sequences, KL 160 44.3863 0.0009 sequences, KL 200 44.3867 0.0005 1 sequence, present 50 44.3859 0.0013 2sequences, present 72 44.38704 0.0002 sequences, present 91 44.38714 <0.0001 The spherical basis set includes 16s16p16d16f16g16h16i orbitals. Again it gives a large basis set error. For the KL basis sets, I used the same parameter as before, AKL = 0.18. Its accuracy can also be improved greatly by systematic augmentation of secondary sequences. However, the basis size will be increased considerably. Based on the previous detailed study on the basis set for the hydrogen atom, here we prescribe a procedure to construct the basis set for a multielectron atom in a B field, with the C atom as the example. We first assign the effective nuclear charge Zf for each electron roughly. The approach is by approximate isoelectronic sequences. Since the Is electrons feel the whole strength of the nuclear attraction, we use 6 for them. For the 2p electron, the nucleus is screened by the two Is electrons, so we use 4, and so forth. Next basis functions are generated according to eqns. (3.26) through (3.28) using Nb= 16, 8, 8, 4, 4 for the primary, the second, ..., the fifth sequences. To use Table 36 as guidance in selecting subsets of basis functions from the previously generated sequences, first recall the scaling argument after eqn. (3.29). The Is wavefunction with an effective nuclear charge Zeffis = 6 in a field B = 10 can be scaled from the hydrogen atom is wavefunction in a field B' = B/ZbX, = 0.28. The value ofB' falls in the range of 0 to 1. By inspection of Table 36, we find a sufficient choice of basis set includes the first through the fourteenth elements in the primary, 13 elements in the second, and 12 in the third sequences. But do not forget the scaling factor. Since the C atom Is wavefunction is tighter than the H atom Is wavefunction approximately by Zeffs, = 6 times, the basis function exponents used in the expansion of the C atom Is orbital should be larger than those used for the H atom by Z2 = 36 times (recall eqns. 3.30 and 3.31). Remember the exponents fl,'s consist of a geometrical series (eqn. 3.28). The increase of the exponents amounts to a 21n6 shift of logq Z2 21n, = 4.6 elements in the primary sequence, and a shift of In 2.18 21n 6  = 3.4 elements in the second and the third sequences. Hence, we should pick the In 2.84 519, 47, and 46 elements in the primary, the second, and the third sequences, respectively. Basis function selections for other electron orbitals are similar. The final basis set is 22s19pI6d21f3g, which is summarized in Table 38. Among the total of 91 gaussians in this basis set, 50 are from the primary, 22 from the second, and 19 from the third sequences. From Table 37, we can see the accuracy of this basis set is remarkably higher than that of the others. By using only the primary sequence, the error left is close to 1 mH. Supplementation with the second sequence reduces the error to 0.2 mH. We estimate the error of the 3sequence present basis to be less than 0.1 mH. Table 38 Construction of the AGTO basis set for the carbon atom in B = 10 au. Orbital Zeff B'= B/Z2 H atom basis shifts C atom basis eft is 6 0.28 114, 13, 12 4.6, 3.4, 3.4 519, 47, 46 2p1 4 0.63 19, 03, 3 3.6, 2.6, 2.6 213, 26, 56 3d2 3 1.1 07, 13, 03 2.8, 2.1, 2.1 210, 35, 25 4f3 2 2.5 07, 05, 06 1.8, 1.3, 1.3 29, 16, 17 5g4 1 10 16, 25, 13 0, 0, 0 16, 25, 13 One may wonder why this procedure works so well, or even why it works at all. The main reason is that each electron orbital can be approximated by a hydrogenlike problem fairly closely. For example, the overlap between the Is HF orbital for the carbon atom in B = 10 and the Is orbital for C5 in the same field is 0.998. See Table 39. Actually by adjusting the nuclear charge to 5.494 and 5.572, the overlaps for Is spin down and spin up orbitals with their corresponding hydrogenlike counterparts reach the maxima of 0.9998 and 0.9999, respectively. Other orbitals are similar. Table 39 Overlaps between HF orbitals for the carbon atom in B = 10 au and hydrogen like systems in the same field Orbital Zf overlap Z' overlap Is, [ 6 0.99773 5.494 0.99984 Is, 6 0.99838 5.572 0.99989 2p1 4 0.99967 3.944 0.99968 3d2 3 0.99974 3.145 0.99983 4f3 2 0.99731 2.644 0.99986 5g4 1 0.98224 2.108 0.99981 Now that we have a systematic way to construct reasonably accurate basis sets for atoms in a B field. In the next section, I will use those basis sets for the DFT and CDFT studies. CHAPTER 4 ATOMS IN UNIFORM FIELDS NUMERICAL RESULTS Comparison with Data in Literature I did extensive unrestricted HartreeFock (UHF) and conventional DFT calculations on the atoms He, Li, Be, B, C, and their positive ions Li Be B+ in a large range of B fields with basis sets constructed according to the procedure outlined in the previous chapter. Total energies are compiled in appendix B. Ground states are indicated in orange. Data available from the literature are also included for comparison. The UHF calculations were primarily for purposes of validation. The agreement of our calculations with those from other groups is excellent. For the helium atom, our HF energies are generally slightly lower than those by Jones, Ortiz, and Ceperley [91, 94]. Their earlier calculations used an STO expansion [91] which is labeled as JOCHF1 in Table B1. Later they utilized JOC basis sets within AGTO [94] (also refer eqn. 3.23), which we call JOCHF2. Presumably the small distinction between their data and ours results from better optimized basis sets I generated, as already illustrated in the previous chapter. This observation is supported by the generally closer agreement of their anisotropic basis set results with ours (in contrast to their spherical basis results). One notable exception to the overall agreement is the 1s4f3 state in B = 800 au. Our result is 23.42398 Hartree versus theirs 23.4342. Another is ls3d2 at B = 560 au: 21.59002 versus their 21.5954 Hartree. The reason for these discrepancies is unclear. It may be some peculiarity of the basis for a particular field strength. For other atoms and ions, the data for comparison are mainly from the series of studies by Ivanov and Schmelcher [86, 8890]. Our HF energies generally match or are slightly higher than theirs, and the overall agreement is quite satisfactory. Differences are usually less than 0.1 mH, far surpassing our goal of 1 mH accuracy for the basis set. The remaining differences arise, presumably from our basis set truncation error and their numerical mesh errors. As for any basisset based calculations, we can only use a finite number of basis functions, which will cause basis set truncation error. Since this error in our calculation is always positive (by the variational principle), one can use our data as an upper bound for the HF energies. However, the numerical error in Ivanov and Schmelcher's 2D HF mesh method seems to tend to be negative. Recall the comparison made in Table 34 for the hydrogen atom. Their energies are always lower than the more accurate algebraic result [80]. Another indication is the zero field atomic energies. For example, the HF energy for the beryllium atom is known accurately to be 14.57302316 Hartree [101], which agrees well with our result of 14.57302, but Ivanov and Schmelcher gave a lower value of 14.57336 [88]. They commented that this configuration has large correlation energy and the contribution from the ls22p2 configuration should be considered, but did not specify whether their result was from single determinant or multideterminant calculation. Since multi configuration HF (MCHF) is not our main interest here, our data are solely from single configuration HF calculations. They also noted that the precision of their mesh approach decreases for the is22s2 configuration in a strong field. This can also contribute to the discrepancy. From previous observations, one may speculate that the true HF energies lie between our data and theirs. Furthermore, the accuracy of our data is ready to be improved by invoking a larger basis set, but this is not necessary for the purpose of the present study. There are fewer DFT calculations for atoms in B fields than HF studies. As far as I know, appendix B is the first extensive compilation of magnetized atomic energies based on modern DFT functionals. I chose PW92 [69] for LDA and PBE [70] for GGA calculations. For the fieldfree case, the present results agree well with published data [102]. Since reference 102 only gave spin nonpolarized DFT energies for spherically distributed densities (which is no problem for the helium and beryllium atoms), one needs to use fractional occupation numbers in other atoms for comparison. For example, in the carbon atom, twop electrons are placed in six spin orbitals, p ,p' and p+1 with equal occupation number of 1/3 for each of them. Actually there are more accurate data for the VWN functional [67]. My results differ from those from reference 67 by no more than 5 /H if I choose the VWN functional, either for spinpolarized or spin nonpolarized energies, neutral atoms or their positive ions. Comparison of nonvanishing B field DFT calculations is handicapped by the different magnetic field grids on which different authors present their results, and by the different functionals implemented. The functional due to Jones [103], which is at the level of the local density approximation, was used by Neuhauser, Koonin, and Langanke [104], and by Braun [84]. The simple Dirac exchangeonly functional was used by Relovsky and Ruder [83], and by Braun [84]. When I choose the Dirac exchange functional, good agreement is found with Braun's data. However, no application of densitygradientdependent fiunctionals on atoms in a strong magnetic field is found in the literature so far as I can tell. For the CDFT study, the present work apparently is the first fully selfconsistent calculation on atoms in large B fields based on the VRG functional. The most closely related study is the perturbative implementation of the VRG functional by Orestes, Marcasso, and Capelle [37] on the atomic ionization energies in vanishing B field. For comparison, CI results are available for helium, lithium, and beryllium atoms, and their positive ions in external fields [7, 8, 1113, 77, 105]. Those data will be treated as the most reliable ones in my comparison. Magnetic Field Induced Ground State Transitions The most drastic change of the ground state atoms caused by an external magnetic field is a series of configuration transitions from the fieldfree ground state to the high field limit ground state. It is well known that for the fieldfree case, two competing factors the spherical nucleus attractive potential and the electronelectron interactions  lead to the shell structure of atoms. This structure is perturbed slightly if the external field is relatively small, but when the field is strong enough that the Lorentz forces exerted on the electrons are comparable to nuclear attraction and electron repulsions, the original shell structure is crushed, and the electrons make a new arrangement. Thus, a series of configuration transitions happens as the B field becomes arbitrarily large. In the largefield limit, the ground state is a fully spinpolarized state in which the electrons take the minimum value of spin Zeeman energy. Ivanov and Schmelcher further analyzed the spatial distribution of electrons and described the highfield limit ground state "... with no nodal surfaces crossing the z axis (the field axis), and with nonpositive magnetic quantum numbers decreasing from m = 0 to m = N + 1, where Nis the number of electrons" [76]. In HartreeFock language, the high field configuration is ls2p _3d 24f3 **. In this regime, the magnetic field is the dominant factor, and coulombic forces can be treated as small perturbations. A cylindrical separation of the z part from x and y parts is usually 52 made for the electron state. Its motion in the xy plane is described by a Landau orbital, and the question is reduced to a quasiID problem. This technique is often referred to as the adiabatic approximation, valid only in the limit of very large field. Many early investigations on matter in aB field concentrated on this regime [103, 104, 106]. C .24 '1'., s  2 m ,f 22 4 2.4  s. .. . 32 0 0,5 1 1.5 2 B(au) FIG. 41 UHF total energies for different electronic states of the helium atom in B fields. Curves 1 to 9 represent configurations Is2, ls2s, 1s2po, ls2pi, Is3Lj, s3d2, ls4f2, ls4pfJ3, and ls5g3, respectively. The most difficult part is the region of intermediate B field, for which both cylindrical and spherical expansions are inefficient, and where the ground state configurations can only be determined by explicit, accurate calculations. Figure 41 displays the HF energies of various configurations for the helium atom in B fields as listed in Table Bl, which includes one singlet state and eight triplet states. Curves 1 to 9 represent configurations of Is2, ls2s, ls2po, ls2pi, ls3Lj, ls32, ls4f2, ls4pf3, and Ss5g3, respectively. Each configuration belongs to a spin subspace according to its total spin z component. For convenience, we use the inexact terminology of "local ground state" and "global ground state". For a specific field strength, the configuration which has the lowest energy within a spin subspace is called the local ground state for this spin multiplet. Among them, the one taking the minimum regardless of its spin is called the global ground state. Thus in Fig. 41, the singlet state remains as the global ground state until B reaches 0.71 au., then a configuration transition to ls2pil occurs. This triplet state is the global ground state for B fields larger than 0.71 au. Atoms with more electrons can have more complicated series of configuration transitions. For example, the carbon atom undergoes six transitions with seven electronic configurations involved being the ground states in different regions of B field strengths [90]. This scenario is basically the same if one uses DFT or CI energies instead of HF energies, except the crossing points for different configurations usually change. Atomic Density Profile as a Function of B Within each configuration, the electron density is squeezed toward the z axis with increasing B field. This follows from energy minimization: the electron density shrinks toward the z axis to alleviate the corresponding diamagnetic energy increment (expectation of B2(x2+y2)/8). Figure 42 shows the density profiles for the ls2 and ls2pl states of the helium atom at field strengths B = 0, 0.5, 1, and 10 au. The transverse shrinkage is quite evident. However, this shrinkage increases the electron repulsion energy. A configuration transition therefore will happen at some field strength (for He, B = 0.71 au.), accompanied by a change of quantum numbers, and eventually a spin flip. Note Figs. 42 (a), (b), and (g). The energy increase in the diamagnetic term caused by FIG. 42 Crosssectional view of the HF total electron densities of the helium atom Is2 (panels ad) and s2pl (panels eh) states as a function of magnetic field strength. Each large tick mark is 2 bohr radii. The B field orientation is in the plane of the paper from bottom to top. The density at the outermost contour lines isl0 6ao3, with a factor of 10 increase for each neighboring curves. Panels (a), (b), (g), and (h) are global ground states. See text for details. the electron density expansion in the new spinconfiguration for the global ground state is more than compensated by the accompany energy lowering of the Zeeman and electron electron repulsion terms. In fact, the same net lowering can occur (and sometimes does occur) without change of spin symmetry. Total Atomic Energies and Their Exchange and Correlation Components Atomic total energies of He, Li, Be, B, C, and their positive ions Li Be B+ in a large range of B fields within the HF, DFTLDA, and DFTGGA approximations are compiled in Appendix B. The exchangecorrelation energy Exc, and its exchange and correlation components Ex, Ec corresponding to the total energies in appendix B, are given in appendix C. As is conventional, the HF correlation energy is defined as the difference between the CI total atomic energy and HF total energy tabulated in appendix B. DFT exchange and correlation energies are defined at the selfconsistent electron densities within the corresponding XC energy functionals (LDA or PBE). Keep in mind that exchange and correlation energies are not defined identically in the wavefunction and DFT approaches; recall the discussion in Chap. 1, particularly eqn. (1.11). Because energies of different states in different field strengths vary considerably, we compare their differences from the corresponding CI total energies. Figure 43 shows those differences for the HF and DFT total atomic energies of helium atom ls2, 1 2po, and s2pl states as functions of B field strength. Since the HF calculation includes exchange exactly, the difference for the HF energy is the negative of the conventional correlation energy EcH (the superscript "HF" is for clarify). First we observe that the conventional correlation energies for the states ls2po and Is2p., are extremely small. 56 0 25 02 0 015 a)He l2 0 0 S 01 o 0 05  005 A A 0 0 HF 0 05 PBE 01 ..... 102 10 100 10' 102 B (au) 005 O O o oo oo oo0 o 0a 9o A A 0 0A S005 b) He 1s2po 015 F HF O LDA 0 2 GGA 0 25' 10 2 10 1 10 101 102 B (au) 01 0 0  c) He 1s2p _ 0 0o 005 o o o o o oo o F D functions of B field strength. Blue squares o: HF; Black circles o: DFTLDA; Red triangles A: DFTGGA. 0 1 O LDA A GGA 10 2 101 100 101 102 B(au) FIG. 43 Differences of the HF and DFT total atomic energies of the helium atom 1s2, ls2po, and ls2p.i states with respect to the corresponding CI energies as functions ofB field strength. Blue squares D: HF; Black circles o: DFTLDA; Red triangles A: DFTGGA. This is because the Is and 2p electrons in the atom are well separated, unlike the two Is electrons in the ls2 configuration. The increase of the absolute value of EHF in the large B field regime for the two states Is2 and ls2pl seems to be the result of the compression of electron densities illustrated in Fig. 42. The PBE generalized gradient functional gives extremely good results for both the singlet state Is2, and triplet sates 1s2po and Ss2p.i when the B field magnitude is less than 1 au. Both LDA and GGA approximations fail in the large field regime, B > 10 au. Notice the similar performance of DFT functionals for the two triplet states s2po and Is2p.,. The former one does not carry paramagnetic current density, thus there is no CDFT current correction for this configuration, whereas the later one is a currentcarrying state. This observation implies that the success or failure of these particular LDA and PBE functionals is not because they omit current terms. The success of DFT calculations mainly depends on accurate approximations for the system exchange and correlation energies. As given in detail in Chap. 1, DFT exchange and correlation energies differ subtly from conventional exchange and correlation energies. The DFT quantities refer to the auxiliary KS determinant (and include a kinetic energy contribution) whereas the conventional quantities are defined with respect to the HF determinant. Nevertheless, conventional exchangecorrelation energies often are used as the quantity to approximate in DFT exchangecorrelation functionals, mostly for pragmatic reasons. The difficulty is that exact DFT quantities and KS orbitals are only available for a few, very small systems. One of those is discussed in the next chapter. For most systems, the exact KS orbitals are unknown. However, there 1s2, LDA 1s2, PBE 1s2p LDA 1s2p PBE 002 0 0 02 004 0 06 0 0 0 0 08 ls2, LDA ls2, PBE S1s2po, LDA 01 A s2p, PBE S s2p I, LDA S1s2p PBE 02 0 15 005o 8 8 8 00 005 0 1s2, LDA 1s2, PBE 0 15 0 1s2p1' LDA 02 A 1s2p, PBE 0 25 102 101 FIG. 44 Differences of DFT exchange (top panel), correlation (middle panel), and exchangecorrelation (bottom panel) energies with HF ones, for the helium atom in B fields. are abundant HF and correlated calculations for many finite systems, providing good reference densities and energies. To gain a better understanding of the behavior of DFT exchange and correlation approximations, we make a separate comparison of DFT exchangecorrelation energies with the HF ones in Fig. 44. Note that "DFT exchange" here means the Ex term in a particular functional and not exact DFT exchange calculated from KS orbitals. We can see from Fig. 44 that the LDA approximation shows its typical underestimation of exchange and overestimation of correlation energies. The PBE functional gives good approximations to the exchange and correlation energies separately when the B field is less than 1 au., but it seriously overestimates the exchange when B>10 au., while the correlation energy does not depend on the field strength very much. Since exchange dominates the XC energies, the error in the exchange term overwhelms the correlation term in large B fields. Of course both the LDA and PBE functionals are based on analysis of the fieldfree electron gas, in which the exchangecorrelation hole is centered at the position of the electron, and only the spherically averaged hole density enters. This picture breaks down for an atom in a large B field. Because of the strong confinement from the B field, there is strong angular correlation among electrons. The XC hole is not centered at the electron, and is not isotropic. Moreover, the external B field will effectively elongate XC hole as well as the electron density. If one wishes to improve XC functionals for applications in the large B field regime, those factors need to be considered. Another observation from Fig. 44 is that PBE overestimates the correlation energies for the two triplet states, in which the HF correlation is very small. This presumably is due to the imperfect cancellation of selfinteractions in the functional. The lithium positive ion is a twoelectron system isoelectronic with neutral He. It has approximately the same correlation energy as that of the helium atom in the fieldfree case. For nonvanishing field, recall the scaling argument for the wavefunction of hydrogenlike atoms in a B field. The deformation of the atomic density induced by the B field is measured by its reduced strength = B / Z2, rather than by its absolute value B, where Z is the nuclear charge (refer to eqns. 3.30 and 3.31 and discussion). Of course, the atomic configuration is another important factor. For the same electronic configuration, the helium atom in B = 5 au. and the lithium positive ion in B = 10 au. have about the same y values. Indeed they have about the same HF and PBE correlation energies. On the other hand, an attempt at a similar comparison between the lithium atom and the beryllium ion is obscured by two factors. One is the large correlation energy between the two Is electrons for the doublet states. The effect of the external B field on its correlation energy is hardly discernable in the studied range. For the quadruplet state, notice the tabulated conventional correlation energy for the beryllium ion is much smaller in magnitude than that of the lithium atom, giving rise to the suspicion of systematic errors existing in those data. Also notice that the conventional correlation energy of the lithium atom ls2p.13d2 state in vanishing B field is even larger than that of its ground state ls22s. Since the electrons are wellseparated in the ls2p.13d2 state, its correlation energy is expected to be smaller than that of a more compact state. Even for vanishing B field, large discrepancies on the correlated energies are found in the literature. For example, Al Hujaj and Schmelcher gave 14.6405 Hartree for the ground state of the beryllium atom from a full CI calculation [8], versus 14.66287 Hartree from a frozencore approximation by Guan et al. [7]. The difference is more than 20 mH. This shows it is difficult to systematically extract atomic correlation energies from the literature, especially for nonvanishing field data. The DFT functionals investigated here fail spectacularly in a large B field, mainly from their exchange part. However, the PBE correlation still gives a large portion of the correlation energy even though its performance degrades somewhat with increasing B field. On the other hand, the HF approximation is more robust than DFTbased calculations and includes exchange exactly, but it totally neglects correlation. From those analyses, it seems a better estimation to the total atomic energies in a large B field may be achieved by combining HF exact exchange and PBE correlation energy rather than using solely the HF or DFT approximations. Ionization Energies and Highest Occupied Orbitals for Magnetized Atoms Because of the magneticfieldinduced configuration transitions for both magnetized atoms and their positive ions, atomic ionization energies are not monotonic increasing or decreasing smooth functions of the applied B field. This is already obvious from Koopmans' theorem and the UHF total energies in Fig. 41. Here we use the total energy difference between the neutral atom and its positive ion, AEsCF, for estimation of the ionization energy. For each field strength, the ground state configurations for the atom and its positive ion must be determined first. Table 41 and Fig. 45 show the change of ionization energies of the atoms He, Li, and Be with increasing B field. Results from different methods are close. For the beryllium atom, a frozencore calculation [7] gave a larger ionization energy, by 26 mH, in the nearzerofield region than the one derived from AlHujaj and Schmelcher's data [8]. This difference is mainly caused by the lower ground state atomic energy obtained in the former reference, which has already been mentioned near the end of the previous section. Table 41 Atomic ionization energies in magnetic fields (energy in Hartree) Atom Configurations B(au) HF CI LDA PBE Hea Is2 1s 0 0.8617 0.9034 0.8344 0.8929 0.02 0.8516 0.8933 0.8244 0.8828 0.04 0.8415 0.8831 0.8142 0.8727 0.08 0.8208 0.8625 0.7935 0.8520 0.16 0.7782 0.8199 0.7506 0.8092 0.24 0.7340 0.7756 0.7059 0.7647 0.4 0.6409 0.6824 0.6113 0.6706 0.5 0.5798 0.6212 0.5492 0.6089 ls2p_1  Is 0.8 0.4687 0.4741 0.4199 0.4734 1 0.5187 0.5245 0.4685 0.5225 1.6 0.6438 0.6504 0.5887 0.6452 2 0.7132 0.7201 0.6549 0.7136 5 1.0734 1.0816 0.9978 1.0739 10 1.4394 1.4493 1.3519 1.4527 20 1.9061 1.9190 1.8161 1.9554 50 2.7182 2.7378 2.6627 2.8829 100 3.5161 3.5442 3.5445 3.8593 Li 1s22s 1s2 0 0.1963 0.2006 0.2011 0.2054 0.01 0.2012 0.2056 0.2059 0.2102 0.02 0.2068 0.2136 0.2135 0.2178 0.05 0.2177 0.2254 0.2226 0.2269 0.1 0.2329 0.2403 0.2380 0.2425 1s22p 1  1s2 0.2 0.2587 0.2614 0.2691 0.2729 0.5 0.3699 0.3750 0.3844 0.3870 1 0.5025 0.5111 0.5216 0.5226 2 0.6995 0.7113 0.7226 0.7229 ls2p_13d2 I 1s2p_1 3 0.7525 0.7572 0.7635 5 0.9475 0.9555 0.9558 0.9644 10 1.2877 1.2982 1.3074 1.3219 Be 1s22S2 1S22s 0 0.2956 0.3158 0.3318 0.3306 0.001 0.2951 0.3159 0.3313 0.3302 0.01 0.2905 0.3112 0.3267 0.3255 0.02 0.2852 0.3313 0.3214 0.3203 0.05 0.2683 0.2911 0.3047 0.3035 ls22s2pi_ ls22s 0.1 0.3234 0.3242 0.3304 0.3312 0.2 0.3941 0.3941 0.4010 0.4019 0.3 0.4531 0.4597 0.4603 0.4612 1s22s2p1 * 1s22p 0.4 0.4687 0.4758 0.4677 0.4717 0.5 0.4710 0.4749 0.4713 0.4758 0.6 0.4696 0.4767 0.4718 0.4766 0.8 0.4593 0.4650 0.4663 0.4718 1s22p_13d2 Is22pl 1 0.4559 0.4455 0.4575 0.4636 2 0.6231 0.6217 0.6257 0.6336 ls2p_13d24f3 ls2p_13d2 5 0.8787 0.8772 0.8895 0.9019 10 1.1973 1.1959 1.2223 1.2401 (a) Exact energies are used for the oneelectron system He+. Table 42 Eigenvalues for the highest occupied orbitals of magnetized atoms (energy in Hartree) Atom Configuration B (au) HOMO HF LDA PBE He 1s2 0 Is 0.91795 0.5702 0.5792 0.02 0.90789 0.5601 0.5692 0.04 0.89771 0.5499 0.5589 0.08 0.87699 0.5289 0.5379 0.16 0.83412 0.4848 0.4940 0.24 0.78942 0.4383 0.4476 0.4 0.69501 0.3387 0.3485 0.5 0.63298 0.2728 0.2829 ls2p_1 0.8 2p_1 0.47120 0.3184 0.3184 1 0.52183 0.3529 0.3532 1.6 0.64824 0.4389 0.4408 2 0.71820 0.4867 0.4900 5 1.07974 0.7379 0.7502 10 1.44629 0.9994 1.0230 20 1.91388 1.3424 1.3824 50 2.72841 1.9648 2.0381 100 3.52994 2.6077 2.7192 Li 1s22s 0 2s 0.19636 0.1162 0.1185 0.01 0.20122 0.1211 0.1234 0.02 0.20668 0.1282 0.1306 0.05 0.21778 0.1364 0.1389 0.1 0.23293 0.1487 0.1516 1s22_1 0.2 2p_1 0.25885 0.1728 0.1751 0.5 0.37038 0.2529 0.2549 1 0.50398 0.3472 0.3489 2 0.70293 0.4873 0.4903 1s2p_13d2 5 3d2 0.95259 0.6626 0.6763 10 1.29348 0.9139 0.9355 Be 1s22s2 0 2s 0.30927 0.2058 0.2061 0.01 0.30417 0.2007 0.2010 0.02 0.29888 0.1953 0.1957 0.05 0.28186 0.1779 0.1783 ls22s2pi_ 0.1 2p_ 0.33120 0.1959 0.1961 0.2 0.40159 0.2560 0.2566 0.3 0.46016 0.3046 0.3054 0.4 2s 0.47732 0.3108 0.3163 0.5 0.47908 0.3099 0.3161 0.6 0.47722 0.3064 0.3134 0.8 0.46613 0.2953 0.3037 1s22p_13d2 1 3d2 0.46092 0.3105 0.3180 2 0.62799 0.4283 0.4391 ls2p_13d24f3 5 4f3 0.88345 0.6259 0.6421 10 1.20284 0.8671 0.8908 64 1.5 O O He E Li O Be ' 0 0 0.5  O' FIG. 45 Atomic ground state ionization energies with increasing B field. Data plotted are from CI calculations shown in Table 41. Dotted lines are the guides to the eye. Even though the ionization energies in both the low and intermediate field regions are rather complicated as the result of atomic configuration changes, their behaviors are similar for the strong field limit configurations. This is an indication that the original atomic shell structure has been effectively obliterated by the field. Eigenvalues of the highest occupied orbital are reported in Table 42. In all the cases, the F orbital energies give the closest approximation to the atomic ionization energies, with an average deviation of only 7.6 mH. KS eigenvalues, as usual, significantly underestimate the ionization energy. This is because both LDA and P 101 B(au) FIG. 45 Atomic ground state ionization energies with increasing B field. Data plotted are from CI calculations shown in Table 41. Dotted lines are the guides to the eye. Even though the ionization energies in both the low and intermediate field regions are rather complicated as the result of atomic configuration changes, their behaviors are similar for the strong field limit configurations. This is an indication that the original atomic shell structure has been effectively obliterated by the field. Eigenvalues of the highest occupied orbitals are reported in Table 42. In all the cases, the HF orbital energies give the closest approximation to the atomic ionization energies, with an average deviation of only 7.6 mH. KS eigenvalues, as usual, significantly underestimate the ionization energy. This is because both LDA and PBE functionals give approximate potentials too shallow compared with the exact DFT XC potential. Selfinteraction correction (SIC) could significantly improve these values, but we will not pursue it further here, because our focus is on Exc functionals that are not explicitly orbitally dependent. Current Density Correction and Other Issues Advancement in CDFT, especially in applications, is hindered by the lack of reliable, tractable functionals. In comparison with the vast literature of ordinary DFT approximate XC functionals, the total number of papers about CDFT approximate functional is substantially less than 50. The earliest proposed functional, also the most widely investigated one as of today, is the VRG functional [1416]. Even for it there are no conclusive results for B 0 in the literature. Here we examine this functional in detail for atoms in a B field, and show the problems inherent in it. The analysis indicates that the VRG functional is not cast in a suitable form, at least for magnetized atoms. The choice of vorticity 9 = V x JP as the basic variable to ensure gaugeinvariance, which is n the central result of references 1416, needs to be critically reexamined. The challenge to implementing CDFT is, somewhat paradoxically, that the current effect is presumably small. We do not expect that the current correction within CDFT will drastically change the DFT densities. Therefore the first question we encountered is which DFT functional should be used as a reference for the CDFT calculations. If the variation in outcomes that results from different choices of DFT functionals is much larger than the CDFT corrections, which seems to be the case in many situations, the predictability of the calculation is jeopardized. Of course, there is no easy answer to this question. Indeed, DFT functionals themselves are still undergoing improvement. Here we made the conventional DFT choice of using the LDA as the starting point. Even though not the most accurate one, the LDA is among the best understood DFT functionals. It is also easy to implement. Using selfconsistent KS orbitals obtained from LDA calculation for the helium atom ls2pl state in a field B = 1 au., I plotted various quantities that are important in CDFT along the z and p directions in Fig. 46. 104 102 100 o 4 102 0 10 4 10 106 108 1010 0 0.5 1 1.5 2 2.5 3 3.5 4 zor p (ao) FIG. 46 Various quantities (electron density n, paramagnetic current densityj, vorticity v, and the current correction to the exchangecorrelation energy density, gv2, in the VRG functional) for the helium atom ls2pl state in B = 1 au. All quantities are evaluated from the LDA KS orbitals and plotted along the z and p axes (cylindrical coordinates). Exponential decay of electron density was already seen in Fig. 42. Because the current density along the z axis is zero, it does not appear in Fig. 46. However, i is not zero on that axis. On the contrary, it diverges at the two poles of the atom. This divergence causes large values of g(n)v2, the energy density correction within the VRG functional (recall eqns. 1.39 and 2.9). If the prefactor g(n) does not decay fast enough to cancel this divergence, a convergence problem in the SCF solution of the CDFT KS equation will happen. Also notice that the electron density decays very rapidly along the z axis. At z = 3ao, the density is already smaller than 10 4a3. Remember the function g(n) was fitted to data points with rs < 10a thus in the lowdensity region it is not welldefined. Different fits to the same set of original data points vary considerably (refer to eqns. 1.41 through 1.46). Furthermore, even the accuracy of the original data set to be fitted is questionable at r~ 10a0. Even were these problems to be resolved, the underlying behavior shown in Fig. 46 would remain. The largest correction relative to ordinary DFT given by the VRG functional is at the places where the electron density and the current density are both almost zero, which is obviously peculiar if not outright unphysical. This peculiar (and difficult) behavior is rooted in the choice of v as the basic variable in the CDFT functional. To avoid the divergence problem, we introduced a rapidly decaying cutoff function. Details were given in chapter 2 (also recall eqns. 2.15 and 2.16). Using parameters ncutoff = 10 3a03, aCtff = 2ao' for the cutoff function, Table 43 lists some of the calculated results within the VRG approximation for the fully spinpolarized states of the helium, lithium, and beryllium atoms at several selected field strengths. An estimation of the current effect is to evaluate the VRG functional using the LDA KohnSham orbitals. Results for that estimation are listed in the third column of Table 43. This scheme can be thought as a nonselfconsistent postDFT calculation. Fully selfconsistent CDFT calculations were also accomplished when the B field is not too large, and they verified the LDAbased estimates. When the B field is larger than roughly 5 au., SCF convergence problems return because of the pathological behavior of VRG functional. Table 43 CDFT corrections to LDA results within VRG approximation (parameters c,,toff, = 10 a3, cutoff = 2ao0 are used for the cutoff function, AE in Hartree) Atom and State B (au) NonSCF AERG SCF AERG JHOMO He ls2p.1 0 0.0022 0.0021 0.0001 0.24 0.0031 0.0031 0.0013 0.5 0.0045 0.0047 0.0029 1 0.0077 0.0081 0.0071 5 0.036 10 0.074 100 0.81 Li 1s2p_13d2 0 0.0070 0.0071 0.0002 2 0.027 0.029 0.0077 5 0.065 10 0.129 Be ls2p_13d24f3 1 0.025 0.026 0.0017 5 0.085 10 0.166 Putting these concerns aside, consider Table 43. By design, the current correction given by the VRG functional is negative. It strongly depends on the particular atomic configuration. Within each configuration, the VRG contribution increases with increasing B field. Besides total atomic energies, the eigenvalues of the highest occupied KS orbitals are also slightly lowered by including the current term, but it helps little in bringing the HOMO energies closer to the ionization energies. This error, of course, is the wellknown selfinteraction problem. Because of the use of a cutoff function, these CDFT calculations can at best be thought of as semiquantitative. This is because the current corrections strongly depend on the chosen cutoff function. Use of different cutoff parameters gives quite different results (see Table 44), an outcome which is really undesirable. Of course, all of this is because the VRG functional does not provide a suitable form for either the low density or the highdensity regions, nor do we know its correct asymptotic expression. Table 44 Effect of cutoff function on CDFT corrections for the helium atom ls2p._ state in magnetic field B = 1 au. (energy in Hartree) ncuoff (ao a.CU (aO) NonSCF AERG 0.005 2.0 0.004 0.001 2.0 0.008 0.001 1.0 0.010 0.0001 2.0 0.025 0.00001 2.0 0.064 It is unsurprising that the VRG functional fails when applied to atomic systems in a strong magnetic field. It was developed from the study of the moderately dense to dense HEG in a weak magnetic field, for which Landau orbitals were used as approximations. This physical situation is quite different from a finite system. First, electron and paramagnetic current densities vary considerably within an atom, and the low density regions (r, > 10a0) are nonnegligible. Secondly, there is not a direct relationship between j (r) and the external B field as there is for the HEG. The question whether the electron gas remains homogeneous after imposing a substantial B field is even unclear. If the field induces some form of crystallization, the basic picture based on which the VRG functional proposed is completely lost. The analysis and numerical studies in this chapter suggest the picture of Landau orbitals used for the HEG may not be applicable at all for the atomiclike systems. Unlike the LDA, also based on the HEG, it seems that the VRG functional is too simple to encompass the essential physics of the atomic systems. A more fundamental question is whether i (r) is a suitable basic variable in gaugeinvariant CDFT as Vignale and Rasolt required [1416]. While it is appealing from a purely theoretical perspective, our numerical results on atomic systems in B fields suggest it is an inappropriate choice, or at least an awkward one, from the application perspective. Largely due to the choice of v (r) as the basic variable in the VRG functional, it gives unphysical results in our tests. Recently, Becke proposed a current dependent functional to resolve the discrepancy of atomic DFT energies of different multiplicity of openshell atoms [107]. Since this functional is based on the analysis of atomic systems, it may be more suitable for application to magnetized atoms than the VRG functional. There are significant technical barriers to its use. Nonetheless, we hope to investigate this functional in the future. Before attempting (sometimes in effect, guessing) better forms for the CDFT functional, we need to know some exact CDFT results to serve as touchstones for any possible proposed functional, This is the major task of the next chapter. Finally, one additional comment remains to be made about the results presented in this chapter. Relativistic effects and the effects due to finite nuclear mass are not considered. Those effects can be important for matter in superstrong fields (B > 104 au), in which regime the adiabatic approximation will be applicable. But for the field strengths involved in this chapter, both effects should be negligible. CHAPTER 5 HOOKE'S ATOM AS AN INSTRUCTIVE MODEL In DFT, the need for accurate approximations to the electronic exchange correlation energy Ex has motivated many studies of a model system often called Hooke's atom (HA) in the DFT literature [108121]. The basic HA is two electrons interacting by the Coulomb potential but confined by a harmonic potential rather than the nuclearelectron attraction. This system is significant for DFT because, for certain values of the confining constant, exact analytical solutions for various states of the HA are known [108, 111]. For other confining strengths, it can be solved numerically also with correlation effects fully included [113]. Since the DFT universal functional is independent of the external potential and the HA differs from atomic He (and isoelectronic ions) only by that potential, exact solutions of the HA allow construction of the exact Ex functional and comparative tests of approximate functionals for such two electron systems. Given that much less is known about the approximate functionals in CDFT than ordinary DFT, it would be of considerable value to the advancement of CDFT to have corresponding exact solutions for the HA in an external magnetic field. Hooke's Atom in Vanishing B Field There is a long history of investigating this problem. The system Hamiltonian reads to= lv1 2+v +(r+r2)+ (5.1) 2 2 r2 where ( (i = 1, 2) are the spatial coordinates of the electrons, and co the confinement frequency. Hartree atomic units are used throughout. By introducing center of mass (CM) and relative coordinates, 1 2 (5.2) r=r2 ' 1=2 (when I deal with the relative motion part, r always means r12), the Hamiltonian eqn. (5.1) becomes Hot = HCM + H (5.3) where H = V +co2R2 (5.4) 4 Hr + V+ 12r2 + (5.5) 4 r The solution to the threedimensional oscillator problem (5.4) can be found in any undergraduate QM textbook. It is the relative motion Schrodinger problem, defined by eqn. (5.5), that has been treated variously by different authors. Laufer and Krieger used the numerical solution to the relative motion problem to construct the exact DFT quantities, and found that, although most approximate functionals generate rather accurate total energies for this model system, the corresponding approximate XC potentials are significantly in error [113]. In 1989, Kais, Herschbach, and Levine found one analytical solution to the HA relative motion problem by dimensional scaling [108]. Samanta and Ghosh obtained solutions by adding an extra linear term to the Hamiltonian [110]. Later, Taut obtained a sequence of exact solutions for certain specific confinement frequencies [111], and used them in studies of DFT functionals [114116]. A basic observation about the HA follows from the Pauli principle, which requires the total wavefunction to be antisymmetric. Because the CM part is always symmetric under particle exchange i ++ r2 if the relative motion part is symmetric (e.g. s or dlike orbitals), the spin part must be antisymmetric, thus a spin singlet state; otherwise a spin triplet state. Thus we can concern ourselves with the spatial relative motion problem alone. Since the Hamiltonian (5.5) is spherically symmetric, the relative motion wave function can be written as the product of a spherical harmonic and a radial part. The radial part is in turn decomposed into a gaussian decaying part (ground state wavefunction of a harmonic oscillator) and a polynomial part. In some special conditions, the polynomial has only a finite number of terms, and thus the wavefunction is expressed explicitly in a closed form. Here I proceed slightly differently from the approach in reference 111. Insertion of the relative motion wavefunction lvr ()= Ylm (O,(p)e r2/4 akk (5.6) k=0 in Hir,/r (f) = EYf (F) and a little algebra give the recursion relation (k+2)(k+21+3)ak 2+ak++ k+L + j E ak= 0 (5.7) Suppose the polynomial part in eqn. (5.6) terminates at the nth term, e.g. a, 0 and ak>n = 0. The recursion relation (5.7) for k = n immediately gives E,, = 1+n+ (5.8) Repeatedly invoking eqn. (5.7) for k= n1, n2, ..., 0, 1, we get an expression for a.1 which by definition must be zero, in terms of co. Frequencies which make this expression be zero are the ones that correspond to analytical solutions with eigenvalues given by eqn. (5.8). Table 51 Confinement frequencies co for HA that have analytical solutions to eqn. (5.5) (see eqn. 5.8 for their eigenvalues) n /=0 1 0.50000000000000 2 0.10000000000000 3 0.03653726559926 0.38012940106740 4 0.01734620322217 0.08096840351940 5 0.00957842801556 0.03085793692937 0.31326733875878 6 0.00584170375528 0.01507863770249 0.06897467166559 7 0.00382334430066 0.00849974006449 0.02696238772621 0.26957696177107 8 0.00263809218012 0.00526419387919 0.01342801519820 0.06058986425144 9 0.00189655882218 0.00348659634110 0.00767969351968 0.02409197815100 0.23835310967398 10 0.00140897933719 0.00242861494144 0.00481042669358 0.01216213038015 0.05433349965263 = 1 0.25000000000000 0.05555555555556 0.02211227585113 0.20936920563036 0.01122668987403 0.04778618566245 0.00653448467629 0.01942406484507 0.18237478381198 0.00415579376716 0.01002629075547 0.04231138533718 0.00281378975218 0.00591291799966 0.01743843557070 0.16282427466688 0.00199650951781 0.00380045768734 0.00910586669888 0.03819659970201 0.00146924333165 0.00259554123244 0.00542189229787 0.01589809508448 0.14785508696009 0.00111335083551 0.00185491303393 0.00351289521069 0.00837269574743 0.03496458370680 /=2 0.16666666666667 0.03846153846154 0.01583274147996 0.14620429555708 0.00827862455572 0.03423838224700 0.00494304416061 0.01426990657388 0.13126853074700 0.00321380796521 0.00753956664388 0.03104720074689 0.00221804190308 0.00454144170886 0.01305357779085 0.11975836716865 0.00160027228709 0.00297472273003 0.00694983975395 0.02852875778009 0.00119499503998 0.00206608391617 0.00421416844040 0.01207319134740 0.11054467575052 0.00091728359398 0.00149877344939 0.00277638397339 0.00646564069324 0.02647749580751 /=3 0.12500000000000 0.02941176470588 0.01232668925503 0.11267331074497 0.00655187269690 0.02675808522737 0.00397054409092 0.01130595881607 0.10313090450042 0.00261635006133 0.00605214259197 0.02465640755471 0.00182765149628 0.00369051896040 0.01048109720372 0.09546288545482 0.00133306668779 0.00244506933776 0.00564110607268 0.02293923879074 0.00100531643239 0.00171616534853 0.00345660786555 0.00979689471802 0.08912851417778 0.00077860389150 0.00125702179525 0.00230004951343 0.00529550074157 0.02150263695791 An example may be helpful. Consider = 0, and n = 3. According to the previously 9 a a2 12a3 (1120)a3 prescribed procedure, we get E, = a a1 = 2 co 2 2co2 a,6a2 (124w)a3 a,2a, (130o +72C2)a3 a0 = 3 a I = 3a 603 4 a 2404 To ensure that the last expression vanishes requires that the confinement frequency be ) = [ or co = 0.3801294, 0.0365373. The solution corresponding to the smaller 24 frequency turns out to be a ground state, while the other one is an excited state. Confinement frequencies corresponding to analytical solutions for / = 0, 1, 2, 3 and n < 10 are compiled in Table 51. This tabulation includes more angular momentum quantum numbers and more significant figures than that presented by Taut [111]. For n > 3, there are several solutions. The smallest frequency corresponds to a ground state, the others are for excited states. Hooke's Atom in B Field, Analytical Solution When the HA is placed in an external magnetic field, its lateral confinement can exceed its vertical confinement. It is well known that the magnetic field can greatly complicate the motion of a columbic system. Even for the oneelectron system (H atom), substantial effort is required to get highly accurate results in aB field [7982]. Only recently have calculations on the He atom in a high field been pushed beyond the HF approximation [1113]. As far as I know, no exact solutions are reported in the literature for the 3D Hooke's atom in an external magnetic field. Taut only gave analytical solutions for a 2D HA in a perpendicular B field [117]. Here I present the exact analytical solutions to the magnetized HA [122]. When the nuclear attraction in the He atom is replaced by a harmonic potential, our exact analytical results can serve as a stringent check on the accuracy of the correlated calculations just mentioned. With an external magnetic field chosen along the z axis, the system Hamiltonian becomes H rot= (Vl +(A())2+ V2 + )2 (2 +r2) +1 + (5.9) 2L i 2 rl where Hs,,, = (s, + s )B is the spin part of the Hamiltonian, s, (i = 1,2) are the z components of the spin, and A(f) is the external vector potential. A similar separation of the CM and the relative motion parts is done as in the case of the B = 0 HA: Hto, = H + H, + H (5.10) 1 1 H cM=( ( + 2A(R)) + 2R2 = (V +40 2R2 +B2R2 sin20 + 2MB) (5.11) V. 11 21 1 2 1 1 1 H =( 1 +A(r))2 +0r += V 0 r +B2r2 sin2 0+mB+ i 2 4 r 4 16 2 r =V + p 2 2z2 + B+ (5.12) 2 r a0)2 B2 0) Here o = o) (5.13) 4 16 2 p = x + y2, m and M are magnetic quantum numbers for the relative and CM motion parts. In these expressions, the Coulomb gauge has been chosen. A(r) = Bx F (5.14) 2 The solution for the CM part is (unnormalized) Mc (R) =exp(mZ2 ) PM H ( 2oZ), F,(NP,, +1,LP2) exp(iMD) (5.15) with eigenvalue of EM=(N z+ )o+ B p+n +1 (5.16) 2 2 2 where P = R sin 0, Q = 42 + B2 H N is Nz th order Hermite polynomial., F, is the confluent hypergeometric function (or Kummer function) [123]. Nz,Np = 0, 1, 2,... are the quantum numbers. The relative motion eignvalue problem from eqn. (5.12) generally cannot be solved analytically in either spherical or cylindrical coordinates. The difficulty of solving the Schrodinger equation that corresponds to eqn. (5.12) lies in the different symmetries of the confining potential (cylindrical) and electronelectron interaction term (spherical). Since the effective potential in eqn. (5.12) V()= 22 + 2 expressed as a r combination of cylindrical coordinate variables (p, z) and the spherical coordinate variable r, it proves convenient also to express the relativemotion wavefunction in those combined, redundant variables y, (r) = v (p, z, r (p, z), p) In part motivated by the expected asymptotic behavior, we choose the form _Oz 2 p 2 V r()=e 2 2 pmzf u(r,Zz) e"m (5.17) where z = 0 for even z parity, 1 for odd z parity. Then (H E) V (r) = 0 yields { a2 82 2z 82 2r, j / \2 27 8 ( n 2 1 ~ o^  2z 0 1++l +(a o )z2 J +2 +2 z + u(r,z)=0 (5.18) &hr2 E2 r rr r2z I + (519r z where E =E mB(27, +l)c, 2(m +l)op. (5.19) 2 To avoid messy notation, no quantum numbers are appended to E. This differential equation is not easy to solve either. To proceed, we make a direct, double powerseries expansion u(r,z)= A,,rzr" z"= (5.20) n,,nz=0 to transform eqn. (5.18) into a recurrence relation, 2(n, + 2)(cO cz ) A+2,, 2(n, + 2)[n, + 3 + 2( m + Tz + n )AA,~ +,, + A4 +,n +[2(nrw ,+pnW ) Ac]4 (n +r+l1)(n, + +2),, 2= 0 (5.21) where A, = 0 for i < 0, or j < 0 or j = 2k +1. We seek values of E, o), a) for which the right side of eqn. (5.20) terminates at finite order. Assume the highest power of z that appears is N (A,,> = 0), where N, is an even number. For N, > 2, generally there is no solution to the set of equations that follow from eqn. (5.21). However, a judicious choice, op = 2ow, = 2Nw, (5.22) N allows us to set A, = 0 for 2i + j > N, since there are recurrencee relations of 2 eqn. (5.21) with 2n, + n, = N, that then are satisfied automatically. This is the special case, co, = 2o) = o which corresponds to imposition of an external field B upon a HA with magnitude B = 2,l (5.23) Now we find values of c, that correspond to analytical solutions. Repeated N application of eqn. (5.21) for each combination of 1< _< 2, 2 N 2(ln +1) nz > 2, allows us to express all the coefficients A O m i Nz N terms of tA ,0 homogenous linear equations involving {Ao,0 oN }. To have nontrivial solutions, the determinant of this set of equations must be zero, a requirement which reduces to finding the roots of a polynomial equation in co,. Energy eigenvalues can be easily found by substituting eqns. (5.22), (5.23) and (5.13) into (5.19), E, = (N, +,+2m)+5+3m C (5.24) Here I give an explicit example for m = rZ = 0, N, = 6. Its relative motion energy is E, = 6+ = 17o,, First we set all the coefficients A,,j =0 for 2i+ j > 6. The four equations derived from eqn. (5.21) are already satisfied for (n,nz) = (3,0), (2,2), (1,4), and (0,6). Repeatedly invoking the recursion relation (5.21) for n, = 1, n = 6,4,2, we find A1,4, A1,2, and A1,0, expressed in terms of Ao,2, A0,4, and A0,6. AO,6 A1, 4 2a 2A,4 1 A1,4 zwAO4 5AO,6 A1,2 2 22 20 202 ' Use eqn. (5.21) again twice for n, = 0,n = 4,2, A2 4 O4 A14 + 30A,6 16040 ,= A,6 A,4, 1442 A,2 +80A)z,2+12A,4 2 0z(4(0 +1)AO,4+(4200 17)Ao,6 A,0 = 16o0 2A',2 ; 4w 16w3 Employ eqn. (5.21) one more time for n, = 1, n = 2 to obtain S40A1, A2,2 +12A1,4 A4 1 280, A3, 0 + 80)2 S6w, 2w, 482 Now, all the coefficients are expressed in terms ofAo,2, A0,4, and A0,6. The next step is to apply eqn. (5.21) forn, = 0, n = 1,0,1,2. We have a set of four homogeneous equations, 4A, 2Al, = 0, Ao 12 Ao,o 6A2, 2A, = 0, 80, A, 12A3, + A2 2A,2 = 0, 4 A2,o + A3, 242, = 0. Substitute the expressions for A1,0, A2, 0, A 2,, A A2, and rearrange, [ O 0z,2 +30 ,4 15A,61=0, [ 9644A, +402(1+20w,)A4,60(3+40 )A44+(11112600w)A,6]=0, 1[96 ,34 A+20w (1+1400,)4 3(7+116 )4A,6=0, 16w t 1[384044A, 4802A,4 +(1+164 432002)A4,61=0. 480)3 Z Z z )AO,6 To have nontrivial solutions for Ao,o, Ao,2, A0,4, and Ao,6, the determinant for their coefficients must be zero. This requirement is equivalent to the polynomial equation 8m6 (206115840)4 1946560o3 +50256m2 420mo +) = 0 There is a standard procedure for solving the fourth power polynomial equation [123]. Here I give the nonzero solutions to the above equation 1 34506x,2,3 2 4345 m/)z = x ,, ' x + S8 2,3 1234762704 2,3 184032' 1/2 24759 1 5633902 34759 2( )r 358124231 where x { 425 cos cos +( 1) + 40257 3 1208188081 3 3704288112 Numerical evaluation gives ),Z = 0.0584428577856519844381713514636827195996701651, (third excitation) 0.0230491519033815661266886064985880747559374948, (second excitation) 0.0040457351480954861583832529737697350502295354, (ground state) 0.0089023525372406381159151962974840750107860006. (first excitation) The smallest frequency corresponds to a ground state, others correspond to excited states. Remember those states are not for the same confinement strength, hence not the same physical system. Table 52 lists all the frequencies that correspond to closedform analytical solutions for m = 0, 1, 2 and N. = 2, 4, 6, 8, 10, including both positive and negative z parties. For each frequency found in the previous step, the corresponding eigenvector {A0,0 N } determines the vector of all the coefficients A .N Table 53 gives explicitly some of the solutions to eqn. (5.12). Table 52 Confinement frequencies o, which have analytical solutions to eqn. (5.12) (o = co /2,B = 2,3o, see eqn. 5.24 for their eigenvalues) T N state m = 0 0 2 g 4 g e 6 g e e e 8 g e e e e e 10 g e e e e e e e e 2 g 4 g e 6 g e e e 8 g e e e e e 10 g e e e e e e e e 0.08333333333333(1) 0.01337996093554(5) 0.03958614075938 0.00404573514810 0.00890235253724 0.02304915190338 0.05844285778565 0.00169910717517 0.00326661504755 0.00563875253622 0.01040942255739 0.01653602158660 0.03180329263192 0.00086575722262 0.00147907803884 0.00241664809384 0.00334850037594 0.00397594182269 0.00736429574821 0.01303363861242 0.01973176390413 0.04608252623918 0.03571428571429 0.00707894326171 0.02581579358039 0.00254580870206 0.00563638108171 0.01883419465413 0.02725999959457 0.00119038526800 0.00217041875920 0.00436490208339 0.00573812529618 0.01453566471300 0.02126492440705 0.00064917903905 0.00105825552161 0.00178623911003 0.00215173389288 0.00344914023386 0.00477115358052 0.01179464597863 0.01632356905207 0.02247241353174 (a) g = ground state; e = excited state; (b) Numbers in parentheses are the listing number in Table 53. m=1 0.05000000000000 0.01000000000000 0.02500000000000(4) 0.00343014626071 0.00606707008623 0.01711506721549 0.03965895252427 0.00151575652301 0.00253598704899 0.00412145134063 0.00824400136608 0.01329159554766 0.02138083910939 0.00079244596296 0.00126081972751 0.00181507622152 0.00292405665076 0.00318249275883 0.00526294300755 0.01091689148853 0.01510846119750 0.03333142555505 0.02777777777778(2) 0.00591390023109 0.01964372058675 0.00223524705023 0.00451759211533 0.01489720978500 0.02265588831769 0.00107801657618 0.00183191437804 0.00352808023992 0.00498357626059 0.01198041559842 0.01697086960871 0.00059949900156 0.00093052390433 0.00148750195578 0.00193434188304 0.00286801805941 0.00395655886157 0.01002739489760 0.01328778382152 0.01936642452391 1 m=2 0.03571428571429 0.00778702514725 0.01938688789623 0.00291641684372 0.00468524599259 0.01417651105557 0.03017655861841 0.00135563312154 0.00203265499562 0.00339412325617 0.00670121318725 0.01137236016318 0.01706956908545 0.00072696957550 0.00107688670338 0.00147682860631 0.00255098333227 0.00273000792956 0.00417905341741 0.00954248742578 0.01275915629739 0.02636159496028 0.02272727272727(3) 0.00504676075803 0.01616733845868 0.00197793631977 0.00379720216639 0.01259718274291 0.01934687713121 0.00097964051831 0.00158115925765 0.00300865314483 0.00437993535596 0.01037221561899 0.01435335804110 0.00055475656842 0.00082578168452 0.00128288767589 0.00174595176146 0.00249502584282 0.00340206703840 0.00884019961249 0.01143529433090 0.01698248977049 83 104 x 3D HA in B, c =2co 10 3 12: 1=z1D/ 10  S3D QD in B, o) = 12 p z N P 10 p a 3D HA, B=0, co =co 100 2D  Av<: ]Av^ XCHdl>I: X 101 101 100 101 102 103 104 p FIG. 51 Confinement strengths subject to analytical solution to eqn. (5.12). For ID, C, = oc has been shifted to o, = 1. For 2D, o)Z = oc has been shifted to o) = 1. Hexagon, square, uptriangle, diamond, downtriangle, circle, lefttriangle, plus, righttriangle, and xmark stand for the highest order ofz in ID, p in 2D, r in the 3D spherical HA, t = ( r z) / 2 in QD, in the polynomial part of the relative motion wavefunctions being 1,2,...,10, respectively. For HA in a B field, they stand for N, + ~. Black, blue, and red symbols are for m = 0, 1, 2, respectively. For spherical HA, only r = 0 is included. Notice its odd parity (, = 1,m) and even parity states (; = 0,m +1) are degenerate. For another case of cap = o, /2, which can be thought as a twoelectron quantum dot (QD) in a suitable magnetic field, one can also find analytical solutions to eqn. (5.12) for some specific confinement frequencies. Together with two limiting cases of ID and 2D, they are summarized in reference 122. Figure 51 shows those frequencies subject to analytical solutions to the electron relative motion part. Table 53 Some # 6) 1 1/6 2 1/18 3 1/22 4 1/20 25 3 17 5 472 solutions to eqn. (5.12) for confinement potential o) = 2o3,, B = 2,f30 Relative Motion Wavefunction e (222)/24(l+r/2+2 /12) e (+2 ):72pz(1+r/6+z2/108)e' e (2 2P2)j88p2z(l+r/8+z2/176)e2z e (Z2 +2p2)8p(l+r/4/4z 2/40+p2/80rz2/160z4/3200)e 2 2p24 r 122C2 1+2c 2 118o2 11314cz4 e ++2 4 24 8 112 2 48 24 8 11328 ) Hooke's Atom in B Field, Numerical Solution For arbitrary o) and B values, eqn. (5.12) does not have an analytical solution. To have a clear picture for the dependence of the HA system behavior upon increasing cv or B field, more data points are essential. Expansion of the wavefunction in terms of spherical harmonics is satisfactory when the B field is not too large. For large B values, Landau orbitals are used for expansion. Consider the lowfield expansion first. Spherical expansion: 1 S(r) = (r)(0, ) 7r (5.25) Insertion of the foregoing expansion together with the Hamiltonian (5.12) in the relative motion Schrodinger equation gives a set of coupled differential equations, Erfmfi (r)r)Vf(()fir()) 0 dr2 (= 0, 2, 4,... or = 1, 3, 5,...) where the effective potential is (5.26) "ht ( 1(1 +21) 1 2 12 B2 Vj im(r)= ) m +12 B+ r r2sin2 'm( r2 r 2 4 16 SI (5.27) (l+1)m 2 2 B2 2 21+1 = s,1 ++B+r2 2 r + +(10,201 1'O)(m,20 I'm) r2 r 2 4 24 24 21'+ 1 and (lm, "m" l'm') is a ClebschGordon coefficient. This procedure is very similar to Ruder et al.'s method for treating the hydrogen atom in strong magnetic fields[78]. The numerical solver for eqn. (5.26) was obtained from reference 124, with appropriate modifications made to adapt it to this problem. Next turn to the strong field case, which requires a cylindrical expansion. The expansion used is: Cylindrical expansion. ,M (r) = I g, (z)O~,L" (p, A ) (5.28) Eg(z) 2 g, ()+ 'V,(z)gn(z)= 0 (n 0, 1,2, ...) (5.29) dz ' where the effective potential ;m (Z)=y m ) t 1 "(p,(0) + z2 + B+2 2 +1I 4 4 Calculating the effective potential eV"m (z) is not trivial. I followed Proschel et al.'s scheme [125]. Details are included in appendix D. By use of a similar argument as in reference 90, we can screen out the configurations pertaining to the HA global ground states in B field, which are (m = 0, r, = 0), (m = 1, rZ = 0), and (m = 3, zz = 0). Energies for the relative motion and spin parts are compiled in Table E1 and E.2 for two angular frequencies, co = 0.5, 0.1, respectively. States are labeled by their conserved quantities as v 2S+lm where (2S+1) is the spin multiplicity, and v is the degree of excitation within a given subspace. Their fieldfree notations are also included (e.g. Is, 2p, ...). The larger confinement frequency corresponds to the first analytical solution found by Kais, Herschbach, and Levine [108], and is also the most widely studied one. The smaller frequency has two analytical solutions, one for the Is state in B = 0 and another for the 2p1 state in B = 3/ 5 A sixteen spherical function expansion gives the relative motion energy of the latter state to be 0.4767949192445 Hartree, pleasingly accurate compared to the analytical result E,= L(N +7,z +2m)+5+3m + C = (2+0+2)+5 , *0.1= 0.47679491924311 For Tables E. 1 and E.2, numbers in parentheses denote the number of radial functions used in expansion (5.25); numbers in brackets are the number of Landau orbitals used in eqn. (5.28). It is easy to see that, in the low field regime (B < 1 au.), the spherical expansion outperforms the cylindrical expansion. However, its quality degrades as the B field increases. As Jones et al. have found and as is physically obvious, the high field regime is very demanding for a spherical basis [93]. Note in Table E2, for B = 10 au., the spherical expansion corresponds to /max = 48,49. Clearly, it cannot go much further on practical grounds. 1 The analytical solution for the singlet state of o = in vanishing B field is 10 1 I+ r(2 r2 2 (r lr5(r)= 1++ e 40 (5.31) 10 ;(240+61 ,) 2 20 