Nuclear Many-body Physics; a Computational Approach

May 4, 2015

Preface

Nearly all of physics is many-body physics at the most
microscopic level of understanding. A theoretical understanding of the behavior of quantum mechanical systems with many interacting particles, accompanied with experiments and simulations,
is a great challenge and provides fundamental insights into systems governed by quantum mechanical principles.

Many-body physics applications range from condensed matter physics (for example antiferromagnets, quantum liquids and solids, plasmas, metals, superconductors), nuclear and high energy physics (for example nuclear matter, finite nuclei, quark-gluon plasmas), dense matter astrophysics (for example neutron stars, white dwarfs), atoms and molecules (for example electron correlations, reactions, quantum chemistry), and elementary particles (for example quantum field theory, lattice gauge field models).

Microscopic many-body methodologies have attained a high degree of sophistication in the last decades, from both theoretical and applicative points of view. These developments have allowed researchers to adapt these technologies to a much wider range of physical systems, and to obtain quantitative descriptions of many observables.

Several of these methods have been developed independently of each other. In some cases the same methods have been applied and studied in different fields of research with little overlap and exchange of knowledge. A classic example is the development of coupled-cluster approaches, which originated from the nuclear many-body problem in the late fifties. These methods were adapted and refined to extremely high precision and predicability by quantum chemists over the next four decades. The developments comprised, and still do,
both methodological innovations and the usage of advances in high-performance computing.
As time elapsed,
the overlap with the original nuclear many-body problem vanished gradually, dwindling almost to a non-vanishing exchange of knowledge and methodologies. These methods were however revived in nuclear physics towards the turn of the last century and offer now a viable approach to ab initio We would like to come with a warning here, since it is our feeling that the labelling ab initio attached to many titles of scientific articles nowadays, has reached an inflationary stage. With ab initio we will mean solving Schroedinger's or Dirac's many-particle equations starting with the presumably best inter-particle Hamiltonian, making no approximations. The particles that make up the many-body problem will be approximated as the constituent ones. In the nuclear many-body problem these are various baryons such as protons, neutrons, isobars and hyperons and low-mass mesons. All many-body approaches discussed in this text involve some approximations, and only few of them, if any, can reach the glorious stage of being truely ab initio. Truely ab initio in nuclear physics means to solve the many-body problem starting with quarks and gluons. nuclear structure studies of nuclei, spanning the whole range of stable closed-shell and several open-shell nuclei.

Many of the developments are also not always readily accessible to a larger community of researchers and especially to the young and uninitiated ones. A notorious exception is again the quantum chemistry community, where highly effecient codes and methods are made available.
Numerical packages like Gamess, Dalton and Gaussian provide a reliable computational platform for many chemical applications and studies. In other fields this is not always the case.

Nuclear physics, where we have our basic research experience, is such a field. This book aims therefore at giving you an overview of widely used quantum mechanical many-body methods. We want also to provide you with computational tools in order to enable you to perform studies of realistic systems in nuclear physics. The methods we discuss are however not limited to applications in nuclear physics, rather, with minor changes one can apply them to other fields, from quantum chemistry, via studies of atoms and molecules to materials science and solid state devices such as quantum dots.

Typical methods we will discuss are configuration interaction and large-scale diagonalization.methods (shell-model in nuclear physics), perturbative many-body methods, coupled cluster theory, Green's function approaches, various Monte Carlo methods, extensions to weakly bound systems and links from ab initio methods to density functional and mean field theories. The methods we have chosen are those which enjoy a great degree of popularity in low and intermediate energy nuclear physics, but due to both our limited knowledge and due to space considerations, some of these will be treated in rather cavalier way.

Focussing on methods and tools only is however not necessarily a recipe for success, as we all know a a bad workman quarrels with his tools. Our aim is to provide you with a didactical description of some of these powerful methods and tools with the purpose of allowing you to gain further insights and a possible working experience in modern many-body methods. There are several textbooks on the market which cover specific many-body methods and/or applications, but almost none of them presents an updated and comparative description of various many-body methods. The reason for choosing such an approach is that we want you, the reader, to gain the experience and form your own opinion on the suitability and applicability of specific methods. Too often you may have met many-body practitioners claiming 'my method is better than your method because...'. We wish to avoid such a pitfal, although ours, as any other approach, is tainted by our biases, preferences and obviously ignorance about all possible details.

We also want to shed light, in a hopefully unbiased way, on possible interplays and differences between the various many-body methods with a focus on the physics content. We have therefore selected several physical systems, from simple models to realistic cases, whose properties can be described within the theoretical frameworks presented here. In this text we will expose you to problems which tour many nuclear physics applications, from the basic forces which bind nucleons together to compact objects such as neutron stars.
Our hope is that the material we have provided allows an eventual reader to assess by himself/herself the pro and cons of the methods discussed.

Moreover, the text can be used and read as a large project. Each chapter ends with a specific project. The code you end up developing for that particular project, can in turn be reused in the next chapter and its pertaining project. To give you an example, in chapter xx you will end up constructing your own nucleon-nucleon interaction and solve its Schr\"odinger equation using the Lippman-Schwinger equation. This solution results in the so-called $T$-matrix which in turn can be related to the experimental phase shifts. The Lippman-Schwinger equation can easily be modified to account for a given nuclear medium, resulting in the so-called $G$-matrix. In chapter xx you end up computing such an effective interaction, using similarity transformations as well. With these codes, you are in turn in the position where you can define effective interactions for say coupled-cluster calculations (the topic of chapter xx) or shell-model effective interactions (topics for chapter xx and xx).

Our hope is that this will enable you to assess at a much deeper level the pros and cons of the various methods.
We believe firmly that knowledge of the strengths and weaknesses of a given method allows you to realize where improvements can be made. The text has also a strong focus on computational issues, from how to build up large many-body codes to parallelization and high-performance computing topics.

This texts spans some four hundred pages, and with the promised wide scope we aim at, there have to be topics which will not be dealt with properly. In particular we will not cover reaction theories.
For reaction theories we provide all the inputs needed to compute onebody densities, spectroscopic factors and optical potentials. For density functional theories we discuss how one can constrain the exchange correlations based on ab initio methods. Furthermore, when it comes to the underlying nuclear forces, we will assume that various baryons and mesons are the effective degrees of freedom, limiting ourselves to low and intermediate energy nuclear physics. A discussion of recent advances in lattice quantumchromodynamics (QCD) and its link to effective field theories is also outside the scope of this book, although we will refer to advances in these fields as well and link the construction of nuclear forces to undergoing research in effective field theory and lattice QCD.
The material on infinite matter can be extended to finite temperatures, but we do not adventure into the realm of heavy-ion collisions and the search of the holy grail of quark-gluon plasma. We will throughout pay loyalty to a physical world governed by the degrees of freedom of baryons and mesons. We apologize for these shortcomings, which are mainly due to our lack of detailed research knowledge in the above fields. To be a jack of all trades leads normally to mastering none.

A final warning however. This text is not a text on many-body theories, rather it is a text on applications and implementations of many-body theories. This applies also to many of the algoritms we discuss. We do not go into details about for example Gaussian quadrature or the mathematical foundations of random walks and the Metropolis algorithm. We will often simply state results, but add enough references to tutorial texts so that you can look up the missing wonders of numerical mathematics yourself. Similarly, handling of angular momenta and their recouplings via $3j$, $6j$ or $9j$ symbols should not come as a bolt out from the blue. But again, don't despair, we'll guide you safely to the appropriate literature.

We have therefore taken the liberty to have certain expectations about you, the potential reader.
We assume that you are at least a graduate student who is embarking on studies in nuclear physics and that you have some familiarity with basic nuclear physics, angular momentum theory and many-body theory, typically at the level of texts like those of Talmi, Fetter and Walecka, Blaizot and Ripka, Dickhoff and Van Neck or similar monographs (add refs later). We will obviously state and repeat the basic rules and theorems, but will not go into derivations.
Some of the derivations will however be left to you via various exercises interspersed througout the text.

Friends/colleagues bla bla and sponsors bla bla. also add that errors, misunderstandings etc are mainly due to ourselves!!

Experimental data and definitions of the many-body problem

Introduction

To understand why matter is stable, and thereby shed light on the limits of nuclear stability, is one of the overarching aims and intellectual challenges of basic research in nuclear physics. To relate the stability of matter to the underlying fundamental forces and particles of nature as manifested in nuclear matter, is central to present and planned rare isotope facilities. Important properties of nuclear systems which can reveal information about these topics are for example masses, and thereby binding energies, and density distributions of nuclei.
These are quantities which convey important information on the shell structure of nuclei, with their pertinent magic numbers and shell closures or the eventual disappearence of the latter away from the valley of stability.

Neutron-rich nuclei are particularly interesting for the above endeavour. As a particular chain of isotopes becomes more and more neutron rich, one reaches finally the limit of stability, the so-called dripline, where one additional neutron makes the next isotopes unstable with respect to the previous ones. The appearence or not of magic numbers and shell structures, the formation of neutron skins and halos can thence be probed via investigations of quantities like the binding energy or the charge radii and neutron rms radii of neutron-rich nuclei. These quantities have in turn important consequences for theoretical models of nuclear structure and their application in astrophysics. For example, the neutron radius of $\,{}^{208}\mbox{Pb}$, recently extracted from the PREX experiment at Jefferson Laboratory can be used to constrain the equation of state of neutron matter. A related quantity to the neutron rms radius $r_n^{\mathrm{rms}}=\langle r^2\rangle_n^{1/2}$ is the neutron skin $r_{\mathrm{skin}}=r_n^{\mathrm{rms}}-r_p^{\mathrm{rms}}$, where $r_p^{\mathrm{rms}}$ is the corresponding proton rms radius.
There are several properties which relate the thickness of the neutron skin to quantities in nuclei and nuclear matter, such as the symmetry energy at the saturation point for nuclear matter, the slope of the equation of state for neutron matter or the low-energy electric dipole strength due to the pigmy dipole resonance. See Ref. http://iopscience.iop.org/1402-4896/2013/T152 for several interesting discussions of these topics.

Having access to precise measurements of masses, radii, and electromagnetic moments for a wide range of nuclei allows to study trends with varying neutron excess. A quantitative description of various experimental data with quantified uncertainty still remains a major challenge for nuclear structure theory. Global theoretical studies of isotopic chains, such as the Ca chain shown in the figure below here, make it possible to test systematic properties of effective interactions between nucleons. Such calculations also provide critical tests of limitations of many-body methods. As one approaches the particle emission thresholds, it becomes increasingly important to describe correctly the coupling to the continuum of decays and scattering channels. While the full treatment of antisymmetrization and short-range correlations has become routine in first principle approaches (to be defined later) to nuclear bound states, the many-body problem becomes more difficult when long-range correlations and continuum effects are considered.

Expected experimental information on the calcium isotopes that can be obtained at FRIB. The limits for detailed spectroscopic information are around $A\sim 60$.

The aim of this first section is to present some of the experimental data which can be used to extract information about correlations in nuclear systems. In particular, we will start with a theoretical analysis of a quantity called the separation energy for neutrons or protons. This quantity, to be discussed below, is defined as the difference between two binding energies (masses) of neighboring nuclei. As we will see from various figures below and exercises as well, the separation energies display a varying behavior as function of the number of neutrons or protons. These variations from one nucleus to another one, laid the foundation for the introduction of so-called magic numbers and a mean-field picture in order to describe nuclei theoretically.

With a mean- or average-field picture we mean that a given nucleon (either a proton or a neutron) moves in an average potential field which is set up by all other nucleons in the system. Consider for example a nucleus like $\,{}^{17}\mbox{O}$ with nine neutrons and eight protons. Many properties of this nucleus can be interpreted in terms of a picture where we can view it as one neutron on top of $\,{}^{16}\mbox{O}$. We infer from data and our theoretical interpretations that this additional neutron behaves almost as an individual neutron which sees an average interaction set up by the remaining 16 nucleons in $\,{}^{16}\mbox{O}$. A nucleus like $\,{}^{16}\mbox{O}$ is an example of what we in this course will denote as a good closed-shell nucleus. We will come back to what this means later.

Since we want to develop a theory capable of interpreting data in terms of our laws of motion and the pertinent forces, we can think of this neutron as a particle which moves in a potential field. We can hence attempt at solving our equations of motion (Schroedinger's equation in our case) for this system along the same lines as we did in atomic physics when we solved Schroedinger's equation for the hydrogen atom. We just need to define a model for our effective single-particle potential.

A simple potential model which enjoys quite some popularity in nuclear physics, is the three-dimensional harmonic oscillator. This potential model captures some of the physics of deeply bound single-particle states but fails in reproducing the less bound single-particle states. A parametrized, and more realistic, potential model which is widely used in nuclear physics, is the so-called Woods-Saxon potential. Both the harmonic oscillator and the Woods-Saxon potential models define computational problems that can easily be solved (see below), resulting (with the appropriate parameters) in a rather good reproduction of experiment for nuclei which can be approximated as one nucleon on top (or one nucleon removed) of a so-called closed-shell system.

To be able to interpret a nucleus in such a way requires at least that we are capable of parametrizing the abovementioned interactions in order to reproduce say the excitation spectrum of a nucleus like $\,{}^{17}\mbox{O}$.

With such a parametrized interaction we are able to solve Schroedinger's equation for the motion of one nucleon in a given field. A nucleus is however a true and complicated many-nucleon system, with extremely many degrees of freedom and complicated correlations, rendering the ideal solution of the many-nucleon Schroedinger equation an impossible enterprise. It is much easier to solve a single-particle problem with say a Woods-Saxon potential. Using such a potential hides however many of the complicated correlations and interactions which we see in nuclei. Such an effective single-nucleon potential is for example not capable of describing properties like the binding energy or the rms radius of a given nucleus.

An improvement to these simpler single-nucleon potentials is given by the Hartree-Fock method, where the variational principle is used to define a mean-field which the nucleons move in. There are many different classes of mean-field methods. An important difference between these methods and the simpler parametrized mean-field potentials like the harmonic oscillator and the Woods-Saxon potentials, is that the resulting equations contain information about the nuclear forces present in our models for solving Schroedinger's equation. Hartree-Fock and other mean-field methods like density functional theory form core topics in later lectures.

The aim of this section is to present some of the experimental data we will confront theory with. In particular, we will focus on separation and shell-gap energies and use these to build a picture of nuclei in terms of (from a philosophical stand we would call this a reductionist approach) a single-particle picture. The harmonic oscillator will serve as an excellent starting point in building nuclei from the bottom and up. Here we will neglect nuclear forces, these are introduced in the next section when we discuss the Hartree-Fock method.

The aim of this course is to develop our physics intuition of nuclear systems using a theoretical approach where we describe data in terms of the motion of individual nucleons and their mutual interactions.

How our theoretical pictures and models can be used to interpret data is in essence what this course is about. Our narrative will lead us along a path where we start with single-particle models and end with the theory of the nuclear shell-model. The latter will be used to understand and analyze excitation spectra and decay patterns of nuclei, linking our theoretical understanding with interpretations of experiment. The way we build up our theoretical descriptions and interpretations follows what we may call a standard reductionistic approach, that is we start with what we believe are our effective degrees of freedom (nucleons in our case) and interactions amongst these and solve thereafter the underlying equations of motions. This defines the nuclear many-body problem, and mean-field approaches like Hartree-Fock theory and the nuclear shell-model represent different approaches to our solutions of Schroedinger's equation.

We start our tour of experimental data and our interpretations by considering the chain of oxygen isotopes. In the exercises below you will be asked to perform similar analyses for other chains of isotopes.

The oxygen isotopes are the heaviest isotopes for which the drip line is well established. The drip line is defined as the point where adding one more nucleon leads to an unbound nucleus. Below we will see that we can define the dripline by studying the separation energy. Where the neutron (proton) separation energy changes sign as a function of the number of neutrons (protons) defines the neutron (proton) drip line.

The oxygen isotopes are simple enough to be described by some few selected single-particle degrees of freedom.

• Two out of four stable even-even isotopes exhibit a doubly magic nature, namely $\,{}^{22}\mbox{O}$ ($Z=8$, $N=14$) and $\,{}^{24}\mbox{O}$ ($Z=8$, $N=16$).

• The structure of $\,{}^{22}\mbox{O}$ and $\,{}^{24}\mbox{O}$ is assumed to be governed by the evolution of the $1s_{1/2}$ and $0d_{5/2}$ one-quasiparticle states.

• The isotopes $\,{}^{25}\mbox{O}$, $\,{}^{26}\mbox{O}$, $\,{}^{27}\mbox{O}$ and $\,{}^{28}\mbox{O}$ are outside the drip line, since the $0d_{3/2}$ orbit is not bound.

  Many experiments worldwide! These isotopes have been studied in series of recent experiments. Some of these experiments and theoretical interpretations are discussed in the following articles:

• $\,{}^{24}\mbox{O}$ and lighter: C. R. Hoffman et al., Phys. Lett. B 672, 17 (2009); R. Kanungo et al., Phys. Rev. Lett.~102, 152501 (2009); C. R. Hoffman et al., Phys. Rev. C 83, 031303(R) (2011); Stanoiu et al., Phys. Rev. C 69, 034312 (2004)

• $\,{}^{25}\mbox{O}$: C. R. Hoffman et al., Phys. Rev. Lett. 102,152501 (2009).

• $\,{}^{26}\mbox{O}$: E. Lunderberg et al., Phys. Rev. Lett. 108, 142503 (2012).

• $\,{}^{26}\mbox{O}$: Z. Kohley et al., Study of two-neutron radioactivity in the decay of 26O, Phys. Rev. Lett., 110, 152501 (2013).

• Theory: Oxygen isotopes with three-body forces, Otsuka et al., Phys. Rev. Lett. 105, 032501 (2010). Hagen et al., Phys. Rev. Lett., 108, 242501 (2012).

Masses and Binding energies

Our first approach in analyzing data theoretically, is to see if we can use experimental information to

• Extract information about a so-called single-particle behavior

• And interpret such a behavior in terms of the underlying forces and microscopic physics

The next step is to see if we could use these interpretations to say something about shell closures and magic numbers. Since we focus on single-particle properties, a quantity we can extract from experiment is the separation energy for protons and neutrons. Before we proceed, we need to define quantities like masses and binding energies. Two excellent reviews on recent trends in the determination of nuclear masses can be found in http://journals.aps.org/rmp/abstract/10.1103/RevModPhys.75.1021 and in http://iopscience.iop.org/1402-4896/2013/T152/014017/

A basic quantity which can be measured for the ground states of nuclei is the atomic mass $M(N, Z)$ of the neutral atom with atomic mass number $A$ and charge $Z$. The number of neutrons is $N$.

Atomic masses are usually tabulated in terms of the mass excess defined by

$$ \Delta M(N, Z) = M(N, Z) - uA, $$

where $u$ is the Atomic Mass Unit

$$ u = M(^{12}\mathrm{C})/12 = 931.49386 \hspace{0.1cm} \mathrm{MeV}/c^2. $$

In this course we will mainly use data from the 2003 compilation of Audi, Wapstra and Thibault, see http://www.sciencedirect.com/science/journal/03759474/729/1.

The nucleon masses are

$$ m_p = 938.27203(8)\hspace{0.1cm} \mathrm{MeV}/c^2 = 1.00727646688(13)u, $$

and

$$ m_n = 939.56536(8)\hspace{0.1cm} \mathrm{MeV}/c^2 = 1.0086649156(6)u. $$

In the 2003 mass evaluation there are 2127 nuclei measured with an accuracy of 0.2 MeV or better, and 101 nuclei measured with an accuracy of greater than 0.2 MeV. For heavy nuclei one observes several chains of nuclei with a constant $N-Z$ value whose masses are obtained from the energy released in $\alpha$-decay.

The nuclear binding energy is defined as the energy required to break up a given nucleus into its constituent parts of $N$ neutrons and $Z$ protons. In terms of the atomic masses $M(N, Z)$ the binding energy is defined by

$$ BE(N, Z) = ZM_H c^2 + Nm_n c^2 - M(N, Z)c^2 , $$

where $M_H$ is the mass of the hydrogen atom and $m_n$ is the mass of the neutron. In terms of the mass excess the binding energy is given by

$$ BE(N, Z) = Z\Delta_H c^2 + N\Delta_n c^2 -\Delta(N, Z)c^2 , $$

where $\Delta_H c^2 = 7.2890$ MeV and $\Delta_n c^2 = 8.0713$ MeV.

The following python program reads in the experimental data on binding energies and, stored in the file bindingenergies.dat, plots them as function of the mass number $A$. One notices clearly a saturation of the binding energy per nucleon at $A\approx 56$.

A popular and physically intuitive model which can be used to parametrize the experimental binding energies as function of $A$, is the so-called the liquid drop model. The ansatz is based on the following expression

$$ BE(N,Z) = a_1A-a_2A^{2/3}-a_3\frac{Z^2}{A^{1/3}}-a_4\frac{(N-Z)^2}{A}, $$

where $A$ stands for the number of nucleons and the $a_i$s are parameters which are determined by a fit to the experimental data.

To arrive at the above expression we have assumed that we can make the following assumptions:

• There is a volume term $a_1A$ proportional with the number of nucleons (the energy is also an extensive quantity). When an assembly of nucleons of the same size is packed together into the smallest volume, each interior nucleon has a certain number of other nucleons in contact with it. This contribution is proportional to the volume.

• There is a surface energy term $a_2A^{2/3}$. The assumption here is that a nucleon at the surface of a nucleus interacts with fewer other nucleons than one in the interior of the nucleus and hence its binding energy is less. This surface energy term takes that into account and is therefore negative and is proportional to the surface area.

• There is a Coulomb energy term $a_3\frac{Z^2}{A^{1/3}}$. The electric repulsion between each pair of protons in a nucleus yields less binding.

• There is an asymmetry term $a_4\frac{(N-Z)^2}{A}$. This term is associated with the Pauli exclusion principle and reflectd the fact that the proton-neutron interaction is more attractive on the average than the neutron-neutron and proton-proton interactions.

We could also add a so-called pairing term, which is a correction term that arises from the tendency of proton pairs and neutron pairs to occur. An even number of particles is more stable than an odd number. Performing a least-square fit to data, we obtain the following numerical values for the various constants

• $a_1=15.49$ MeV

• $a_2=17.23$ MeV

• $a_3=0.697$ MeV

• $a_4=22.6$ MeV

The python below here allows you to perform a fit of teh above parameters using nonlinear least squares curvefitting.

The following python program reads now in the experimental data on binding energies as well as the results from the above liquid drop model and plots these energies as function of the mass number $A$. One sees that for larger values of $A$, there is a better agreement with data.

This python program reads now in the experimental data on binding energies and performs a nonlinear least square fitting of the data. In the example here we use only the parameters $a_1$ and $a_2$, leaving it as an exercise to the reader to perform the fit for all four paramters. The results are plotted and compared with the experimental values. To read more about non-linear least square methods, see for example the text of M.J. Box, D. Davies and W.H. Swann, Non-Linear optimisation Techniques, Oliver & Boyd, 1969.

We are now interested in interpreting experimental binding energies in terms of a single-particle picture. In order to do so, we consider first energy conservation for nuclear transformations that include, for example, the fusion of two nuclei $a$ and $b$ into the combined system $c$

$$ {^{N_a+Z_a}}a+ {^{N_b+Z_b}}b\rightarrow {^{N_c+Z_c}}c $$

or the decay of nucleus $c$ into two other nuclei $a$ and $b$

$$ ^{N_c+Z_c}c \rightarrow ^{N_a+Z_a}a+ ^{N_b+Z_b}b $$

In general we have the reactions

$$ \sum_i {^{N_i+Z_i}}i \rightarrow \sum_f {^{N_f+Z_f}}f $$

We require also that the number of protons and neutrons (the total number of nucleons) is conserved in the initial stage and final stage, unless we have processes which violate baryon conservation,

$$ \sum_iN_i = \sum_f N_f \hspace{0.2cm}\mathrm{and} \hspace{0.2cm}\sum_iZ_i = \sum_f Z_f. $$

$Q$-values and separation energies

The above processes can be characterized by an energy difference called the $Q$ value, defined as

$$ Q=\sum_i M(N_i, Z_i)c^2-\sum_f M(N_f, Z_f)c^2=\sum_i BE(N_f, Z_f)-\sum_i BE(N_i, Z_i) $$

Spontaneous decay involves a single initial nuclear state and is allowed if $Q > 0$. In the decay, energy is released in the form of the kinetic energy of the final products. Reactions involving two initial nuclei are called endothermic (a net loss of energy) if $Q < 0$. The reactions are exothermic (a net release of energy) if $Q > 0$.

Let us study the Q values associated with the removal of one or two nucleons from a nucleus. These are conventionally defined in terms of the one-nucleon and two-nucleon separation energies. The neutron separation energy is defined as

$$ S_n= -Q_n= BE(N,Z)-BE(N-1,Z), $$

and the proton separation energy reads

$$ S_p= -Q_p= BE(N,Z)-BE(N,Z-1). $$

The two-neutron separation energy is defined as

$$ S_{2n}= -Q_{2n}= BE(N,Z)-BE(N-2,Z), $$

and the two-proton separation energy is given by

$$ S_{2p}= -Q_{2p}= BE(N,Z)-BE(N,Z-2), $$

Using say the neutron separation energies (alternatively the proton separation energies)

$$ S_n= -Q_n= BE(N,Z)-BE(N-1,Z), $$

we can define the so-called energy gap for neutrons (or protons) as

$$ \Delta S_n= BE(N,Z)-BE(N-1,Z)-\left(BE(N+1,Z)-BE(N,Z)\right), $$

or

$$ \Delta S_n= 2BE(N,Z)-BE(N-1,Z)-BE(N+1,Z). $$

This quantity can in turn be used to determine which nuclei are magic or not. For protons we would have

$$ \Delta S_p= 2BE(N,Z)-BE(N,Z-1)-BE(N,Z+1). $$

We leave it as an exercise to the reader to define and interpret the two-neutron or two-proton gaps.

The following python programs can now be used to plot the separation energies and the energy gaps for the oxygen isotopes. The following python code reads the separation energies from file for all oxygen isotopes from $A=13$ to $A=25$, The data are taken from the file snox.dat. This files contains the separation energies and the shell gap energies.

Here we display the python program for plotting the corresponding results for shell gaps for the oyxgen isotopes.

Since we will focus in the beginning on single-particle degrees of freedom and mean-field approaches before we start with nuclear forces and many-body approaches like the nuclear shell-model, there are some features to be noted

• In the discussion of the liquid drop model and binding energies, we note that the total binding energy is not that different from the sum of the individual neutron and proton masses.

One may thus infer that intrinsic properties of nucleons in a nucleus are close to those of free nucleons.

• In the discussion of the neutron separation energies for the oxygen isotopes, we note a clear staggering effect between odd and even isotopes with the even ones being more bound (larger separation energies). We will later link this to strong pairing correlations in nuclei.

• The neutron separation energy becomes negative at $\,{}^{25}\mbox{O}$, making this nucleus unstable with respect to the emission of one neutron. A nucleus like $\,{}^{24}\mbox{O}$ is thus the last stable oxygen isotopes which has been observed. Oyxgen-26 has been found to be unbound with respect to $\,{}^{24}\mbox{O}$, see <ournals.aps.org/prl/abstract/10.1103/PhysRevLett.108.142503> .

• We note also that there are large shell-gaps for some nuclei, meaning that more energy is needed to remove one nucleon. These gaps are used to define so-called magic numbers. For the oxygen isotopes we see a clear gap for $\,{}^{16}\mbox{O}$. We will interpret this gap as one of several experimental properties that define so-called magic numbers. In our discussion below we will make a first interpretation using single-particle states from the harmonic oscillator and the Woods-Saxon potential.

In the exercises below you will be asked to perform a similar analysis for other chains of isotopes and interpret the results.

Radii

The root-mean-square (rms) charge radius has been measured for the ground states of many nuclei. For a spherical charge density, $\rho(\boldsymbol{r})$, the mean-square radius is defined by

$$ \langle r^2\rangle = \frac{ \int d \boldsymbol{r} \rho(\boldsymbol{r}) r^2}{ \int d \boldsymbol{r} \rho(\boldsymbol{r})}, $$

and the rms radius is the square root of this quantity denoted by

$$ R =\sqrt{ \langle r^2\rangle}. $$

Radii for most stable nuclei have been deduced from electron scattering form factors and/or from the x-ray transition energies of muonic atoms. The relative radii for a series of isotopes can be extracted from the isotope shifts of atomic x-ray transitions. The rms radius for the nuclear point-proton density, $R_p$ is obtained from the rms charge radius by:

$$ R_p = \sqrt{R^2_{\mathrm{ch}}- R^2_{\mathrm{corr}}}, $$

where

$$ R^2_{\mathrm{corr}}= R^2_{\mathrm{op}}+(N/Z)R^2_{\mathrm{on}}+R^2_{\mathrm{rel}}, $$

where

$$ R_{\mathrm{op}}= 0.875(7) \mathrm{fm}. $$

is the rms radius of the proton, $R^2_{\mathrm{on}} = 0.116(2)$ $\mbox{fm}^{2}$ is the mean-square radius of the neutron and $R^2_{\mathrm{rel}} = 0.033$ $\mbox{fm}^{2}$ is the relativistic Darwin-Foldy correction. There are also smaller nucleus-dependent relativistic spin-orbit and mesonic-exchange corrections that should be included.

Definitions

We will now introduce the potential models we have discussex above, namely the harmonic oscillator and the Woods-Saxon potentials. In order to proceed, we need some definitions.

We define an operator as $\hat{O}$ throughout. Unless otherwise specified the total number of nucleons is always $A$ and $d$ is the dimension of the system. In nuclear physics we normally define the total number of particles to be $A=N+Z$, where $N$ is total number of neutrons and $Z$ the total number of protons. In case of other baryons such as isobars $\Delta$ or various hyperons such as $\Lambda$ or $\Sigma$, one needs to add their definitions. When we refer to a single neutron we will use the label $n$ and when we refer to a single proton we will use the label $p$. Unless otherwise specified, we will simply call these particles for nucleons.

The quantum numbers of a single-particle state in coordinate space are defined by the variables

$$ x=(\boldsymbol{r},\sigma), $$

where

$$ \boldsymbol{r}\in {\mathbb{R}}^{d}, $$

with $d=1,2,3$ represents the spatial coordinates and $\sigma$ is the eigenspin of the particle. For fermions with eigenspin $1/2$ this means that

$$ x\in {\mathbb{R}}^{d}\oplus (\frac{1}{2}), $$

and the integral

$$ \int dx = \sum_{\sigma}\int d^dr = \sum_{\sigma}\int d\boldsymbol{r}, $$

and

$$ \int d^Ax= \int dx_1\int dx_2\dots\int dx_A. $$

Since we are dealing with protons and neutrons we need to add isospin as a new degree of freedom.

Including isospin $\tau$ we have

$$ x=(\boldsymbol{r},\sigma,\tau), $$

where

$$ \boldsymbol{r}\in {\mathbb{R}}^{3}, $$

For nucleons, which are fermions with eigenspin $1/2$ and isospin $1/2$ this means that

$$ x\in {\mathbb{R}}^{d}\oplus (\frac{1}{2})\oplus (\frac{1}{2}), $$

and the integral

$$ \int dx = \sum_{\sigma\tau}\int d\boldsymbol{r}, $$

and

$$ \int d^Ax= \int dx_1\int dx_2\dots\int dx_A. $$

We will use the standard nuclear physics definition of isospin, resulting in $\tau_z=-1/2$ for protons and $\tau_z=1/2$ for neutrons.

The quantum mechanical wave function of a given state with quantum numbers $\lambda$ (encompassing all quantum numbers needed to specify the system), ignoring time, is

$$ \Psi_{\lambda}=\Psi_{\lambda}(x_1,x_2,\dots,x_A), $$

with $x_i=(\boldsymbol{r}_i,\sigma_i,\tau_i)$ and the projections of $\sigma_i$ and $\tau_i$ take the values $\{-1/2,+1/2\}$. We will hereafter always refer to $\Psi_{\lambda}$ as the exact wave function, and if the ground state is not degenerate we label it as

$$ \Psi_0=\Psi_0(x_1,x_2,\dots,x_A). $$

Since the solution $\Psi_{\lambda}$ seldomly can be found in closed form, approximations are sought. In this text we define an approximative wave function or an ansatz to the exact wave function as

$$ \Phi_{\lambda}=\Phi_{\lambda}(x_1,x_2,\dots,x_A), $$

with

$$ \Phi_{0}=\Phi_{0}(x_{1},x_{2},\dots,x_{A}), $$

being the ansatz for the ground state.

The wave function $\Psi_{\lambda}$ is sought in the Hilbert space of either symmetric or anti-symmetric $N$-body functions, namely

$$ \Psi_{\lambda}\in {\cal H}_A:= {\cal H}_1\oplus{\cal H}_1\oplus\dots\oplus{\cal H}_1, $$

where the single-particle Hilbert space $\hat{H}_1$ is the space of square integrable functions over $\in {\mathbb{R}}^{d}\oplus (\sigma)\oplus (\tau)$ resulting in

$$ {\cal H}_1:= L^2(\mathbb{R}^{d}\oplus (\sigma)\oplus (\tau)). $$

Our Hamiltonian is invariant under the permutation (interchange) of two particles. Since we deal with fermions however, the total wave function is antisymmetric. Let $\hat{P}$ be an operator which interchanges two particles. Due to the symmetries we have ascribed to our Hamiltonian, this operator commutes with the total Hamiltonian,

$$ [\hat{H},\hat{P}] = 0, $$

meaning that $\Psi_{\lambda}(x_1, x_2, \dots , x_A)$ is an eigenfunction of $\hat{P}$ as well, that is

$$ \hat{P}_{ij}\Psi_{\lambda}(x_1, x_2, \dots,x_i,\dots,x_j,\dots,x_A)= \beta\Psi_{\lambda}(x_1, x_2, \dots,x_j,\dots,x_i,\dots,x_A), $$

where $\beta$ is the eigenvalue of $\hat{P}$. We have introduced the suffix $ij$ in order to indicate that we permute particles $i$ and $j$. The Pauli principle tells us that the total wave function for a system of fermions has to be antisymmetric, resulting in the eigenvalue $\beta = -1$.

The Schrodinger equation reads

$$ \begin{equation} \hat{H}(x_1, x_2, \dots , x_A) \Psi_{\lambda}(x_1, x_2, \dots , x_A) = E_\lambda \Psi_\lambda(x_1, x_2, \dots , x_A), \label{eq:basicSE1} \tag{1} \end{equation} $$

where the vector $x_i$ represents the coordinates (spatial, spin and isospin) of particle $i$, $\lambda$ stands for all the quantum numbers needed to classify a given $A$-particle state and $\Psi_{\lambda}$ is the pertaining eigenfunction. Throughout this course, $\Psi$ refers to the exact eigenfunction, unless otherwise stated.

We write the Hamilton operator, or Hamiltonian, in a generic way

$$ \hat{H} = \hat{T} + \hat{V} $$

where $\hat{T}$ represents the kinetic energy of the system

$$ \hat{T} = \sum_{i=1}^A \frac{\mathbf{p}_i^2}{2m_i} = \sum_{i=1}^A \left( -\frac{\hbar^2}{2m_i} \mathbf{\nabla_i}^2 \right) = \sum_{i=1}^A t(x_i) $$

while the operator $\hat{V}$ for the potential energy is given by

$$ \begin{equation} \hat{V} = \sum_{i=1}^A \hat{u}_{\mathrm{ext}}(x_i) + \sum_{ji=1}^A v(x_i,x_j)+\sum_{ijk=1}^Av(x_i,x_j,x_k)+\dots \label{eq:firstv} \tag{2} \end{equation} $$

Hereafter we use natural units, viz. $\hbar=c=e=1$, with $e$ the elementary charge and $c$ the speed of light. This means that momenta and masses have dimension energy.

The potential energy part includes also an external potential $\hat{u}_{\mathrm{ext}}(x_i)$.

In a non-relativistic approach to atomic physics, this external potential is given by the attraction an electron feels from the atomic nucleus. The latter being much heavier than the involved electrons, is often used to define a natural center of mass. In nuclear physics there is no such external potential. It is the nuclear force which results in binding in nuclear systems. In a non-relativistic framework, the nuclear force contains two-body, three-body and more complicated degrees of freedom. The potential energy reads then

$$ \hat{V} = \sum_{ij}^A v(x_i,x_j)+\sum_{ijk}^Av(x_i,x_j,x_k)+\dots $$

Three-body and more complicated forces arise since we are dealing with protons and neutrons as effective degrees of freedom. We will come back to this topic later. Furthermore, in large parts of these lectures we will assume that the potential energy can be approximated by a two-body interaction only. Our Hamiltonian reads then

$$ \begin{equation} \hat{H} = \sum_{i=1}^A \frac{\mathbf{p}_i^2}{2m_i}+\sum_{ij}^A v(x_i,x_j). \label{eq:firstH} \tag{3} \end{equation} $$

A modified Hamiltonian

It is however, from a computational point of view, convenient to introduce an external potential $\hat{u}_{\mathrm{ext}}(x_i)$ by adding and substracting it to the original Hamiltonian. This means that our Hamiltonian can be rewritten as

$$ \hat{H} = \hat{H}_0 + \hat{H}_I = \sum_{i=1}^A \hat{h}_0(x_i) + \sum_{i < j=1}^A \hat{v}(x_{ij})-\sum_{i=1}^A\hat{u}_{\mathrm{ext}}(x_i), $$

with

$$ \hat{H}_0=\sum_{i=1}^A \hat{h}_0(x_i) = \sum_{i=1}^A\left(\hat{t}(x_i) + \hat{u}_{\mathrm{ext}}(x_i)\right). $$

The interaction (or potential energy term) reads now

$$ \hat{H}_I= \sum_{i < j=1}^A \hat{v}(x_{ij})-\sum_{i=1}^A\hat{u}_{\mathrm{ext}}(x_i). $$

In nuclear physics the one-body part $u_{\mathrm{ext}}(x_i)$ is often approximated by a harmonic oscillator potential or a Woods-Saxon potential. However, this is not fully correct, because as we have discussed, nuclei are self-bound systems and there is no external confining potential. As we will see later, the $\hat{H}_0$ part of the hamiltonian cannot be used to compute the binding energy of a nucleus since it is not based on a model for the nuclear forces. That is, the binding energy is not the sum of the individual single-particle energies.

Why do we introduce the Hamiltonian in the form

$$ \hat{H} = \hat{H}_0 + \hat{H}_I? $$

There are many reasons for this. Let us look at some of them, using the harmonic oscillator in three dimensions as our starting point. For the harmonic oscillator we know that

$$ \hat{h}_0(x_i)\psi_{\alpha}(x_i)=\varepsilon_{\alpha}\psi_{\alpha}(x_i), $$

where the eigenvalues are $\varepsilon_{\alpha}$ and the eigenfunctions are $\psi_{\alpha}(x_i)$. The subscript $\alpha$ represents quantum numbers like the orbital angular momentum $l_{\alpha}$, its projection $m_{l_{\alpha}}$ and the
principal quantum number $n_{\alpha}=0,1,2,\dots$.

The eigenvalues are

$$ \varepsilon_{\alpha} = \hbar\omega \left(2n_{\alpha}+l_{\alpha}+\frac{3}{2}\right). $$

The following mathematical properties of the harmonic oscillator are handy.

• First of all we have a complete basis of orthogonal eigenvectors. These have well-know expressions and can be easily be encoded.

• With a complete basis $\psi_{\alpha}(x_i)$, we can construct a new basis $\phi_{\tau}(x_i)$ by expanding in terms of a harmonic oscillator basis, that is

$$ \phi_{\tau}(x_i)=\sum_{\alpha} C_{\tau\alpha}\psi_{\alpha}(x_i), $$

where $C_{\tau\alpha}$ represents the overlap between the two basis sets.

• As we will see later, the harmonic oscillator basis allows us to compute in an expedient way matrix elements of the interactions between two nucleons. Using the above expansion we can in turn represent nuclear forces in terms of new basis, for example the Woods-Saxon basis to be discussed later here.

The harmonic oscillator (a shifted one by a negative constant) provides also a very good approximation to most bound single-particle states. Furthermore, it serves as a starting point in building up our picture of nuclei, in particular how we define magic numbers and systems with one nucleon added to (or removed from) a closed-shell core nucleus. The figure here shows the various harmonic oscillator states, with those obtained with a Woods-Saxon potential as well, including a spin-orbit splitting (to be discussed below).

Single-particle spectrum and quantum numbers for a harmonic oscillator potential and a Woods-Saxon potential with and without a spin-orbit force.

In nuclear physics the one-body part $u_{\mathrm{ext}}(x_i)$ is often approximated by a harmonic oscillator potential. However, as we also noted with the Woods-Saxon potential there is no external confining potential in nuclei.

What many people do then, is to add and subtract a harmonic oscillator potential, with

$$ \hat{u}_{\mathrm{ext}}(x_i)=\hat{u}_{\mathrm{ho}}(x_i)= \frac{1}{2}m\omega^2 r_i^2, $$

where $\omega$ is the oscillator frequency. This leads to

$$ \hat{H} = \hat{H_0} + \hat{H_I} = \sum_{i=1}^A \hat{h}_0(x_i) + \sum_{i < j=1}^A \hat{v}(x_{ij})-\sum_{i=1}^A\hat{u}_{\mathrm{ho}}(x_i), $$

with

$$ H_0=\sum_{i=1}^A \hat{h}_0(x_i) = \sum_{i=1}^A\left(\hat{t}(x_i) + \hat{u}_{\mathrm{ho}}(x_i)\right). $$

Many practitioners use this as the standard Hamiltonian when doing nuclear structure calculations. This is ok if the number of nucleons is large, but still with this Hamiltonian, we do not obey translational invariance. How can we cure this?

In setting up a translationally invariant Hamiltonian
the following expressions are helpful. The center-of-mass (CoM) momentum is

$$ P=\sum_{i=1}^A\boldsymbol{p}_i, $$

and we have that

$$ \sum_{i=1}^A\boldsymbol{p}_i^2 = \frac{1}{A}\left[\boldsymbol{P}^2+\sum_{i < j}(\boldsymbol{p}_i-\boldsymbol{p}_j)^2\right] $$

meaning that

$$ \left[\sum_{i=1}^A\frac{\boldsymbol{p}_i^2}{2m} -\frac{\boldsymbol{P}^2}{2mA}\right] =\frac{1}{2mA}\sum_{i < j}(\boldsymbol{p}_i-\boldsymbol{p}_j)^2. $$

In a similar fashion we can define the CoM coordinate

$$ \boldsymbol{R}=\frac{1}{A}\sum_{i=1}^{A}\boldsymbol{r}_i, $$

which yields

$$ \sum_{i=1}^A\boldsymbol{r}_i^2 = \frac{1}{A}\left[A^2\boldsymbol{R}^2+\sum_{i < j}(\boldsymbol{r}_i-\boldsymbol{r}_j)^2\right]. $$

If we then introduce the harmonic oscillator one-body Hamiltonian

$$ H_0= \sum_{i=1}^A\left(\frac{\boldsymbol{p}_i^2}{2m}+ \frac{1}{2}m\omega^2\boldsymbol{r}_i^2\right), $$

with $\omega$ the oscillator frequency, we can rewrite the latter as

$$ H_{\mathrm{HO}}= \frac{\boldsymbol{P}^2}{2mA}+\frac{mA\omega^2\boldsymbol{R}^2}{2} +\frac{1}{2mA}\sum_{i < j}(\boldsymbol{p}_i-\boldsymbol{p}_j)^2 +\frac{m\omega^2}{2A}\sum_{i < j}(\boldsymbol{r}_i-\boldsymbol{r}_j)^2. \label{eq:obho} \tag{4} $$

Alternatively, we could write it as

$$ H_{\mathrm{HO}}= H_{\mathrm{CoM}}+\frac{1}{2mA}\sum_{i < j}(\boldsymbol{p}_i-\boldsymbol{p}_j)^2 +\frac{m\omega^2}{2A}\sum_{i < j}(\boldsymbol{r}_i-\boldsymbol{r}_j)^2, $$

The center-of-mass term is defined as

$$ H_{\mathrm{CoM}}= \frac{\boldsymbol{P}^2}{2mA}+\frac{mA\omega^2\boldsymbol{R}^2}{2}. $$

The translationally invariant one- and two-body Hamiltonian reads for an A-nucleon system,

$$ \label{eq:ham} \tag{5} \hat{H}=\left[\sum_{i=1}^A\frac{\boldsymbol{p}_i^2}{2m} -\frac{\boldsymbol{P}^2}{2mA}\right] +\sum_{i < j}^A V_{ij} \; , $$

where $V_{ij}$ is the nucleon-nucleon interaction. Adding zero as here

$$ \sum_{i=1}^A\frac{1}{2}m\omega^2\boldsymbol{r}_i^2- \frac{m\omega^2}{2A}\left[\boldsymbol{R}^2+\sum_{i < j}(\boldsymbol{r}_i-\boldsymbol{r}_j)^2\right]=0. $$

we can then rewrite the Hamiltonian as

$$ \hat{H}=\sum_{i=1}^A \left[ \frac{\boldsymbol{p}_i^2}{2m} +\frac{1}{2}m\omega^2 \boldsymbol{r}^2_i \right] + \sum_{i < j}^A \left[ V_{ij}-\frac{m\omega^2}{2A} (\boldsymbol{r}_i-\boldsymbol{r}_j)^2 \right]-H_{\mathrm{CoM}}. $$

The Woods-Saxon potential is a mean field potential for the nucleons (protons and neutrons) inside an atomic nucleus. It represent an average potential that a given nucleon feels from the forces applied on each nucleon. The parametrization is

$$ \hat{u}_{\mathrm{ext}}(r)=-\frac{V_0}{1+\exp{(r-R)/a}}, $$

with $V_0\approx 50$ MeV representing the potential well depth, $a\approx 0.5$ fm length representing the "surface thickness" of the nucleus and $R=r_0A^{1/3}$, with $r_0=1.25$ fm and $A$ the number of nucleons. The value for $r_0$ can be extracted from a fit to data, see for example M. Kirson, http://www.sciencedirect.com/science/article/pii/S037594740600769X.

The following python code produces a plot of the Woods-Saxon potential with the above parameters.

From the plot we notice that the potential

• rapidly approaches zero as $r$ goes to infinity, reflecting the short-distance nature of the strong nuclear force.

• For large $A$, it is approximately flat in the center.

• Nucleons near the surface of the nucleus experience a large force towards the center.

We have introduced a single-particle Hamiltonian

$$ H_0=\sum_{i=1}^A \hat{h}_0(x_i) = \sum_{i=1}^A\left(\hat{t}(x_i) + \hat{u}_{\mathrm{ext}}(x_i)\right), $$

with an external and central symmetric potential $u_{\mathrm{ext}}(x_i)$, which is often approximated by a harmonic oscillator potential or a Woods-Saxon potential. Being central symmetric leads to a degeneracy in energy which is not observed experimentally. We see this from for example our discussion of separation energies and magic numbers. There are, in addition to the assumed magic numbers from a harmonic oscillator basis of $2,8,20,40,70\dots$ magic numbers like $28$, $50$, $82$ and $126$.

To produce these additional numbers, we need to add a phenomenological spin-orbit force which lifts the degeneracy, that is

$$ \hat{h}(x_i) = \hat{t}(x_i) + \hat{u}_{\mathrm{ext}}(x_i) +\xi(\boldsymbol{r})\boldsymbol{ls}=\hat{h}_0(x_i)+\xi(\boldsymbol{r})\boldsymbol{ls}. $$

We have introduced a modified single-particle Hamiltonian

$$ \hat{h}(x_i) = \hat{t}(x_i) + \hat{u}_{\mathrm{ext}}(x_i) +\xi(\boldsymbol{r})\boldsymbol{ls}=\hat{h}_0(x_i)+\xi(\boldsymbol{r})\boldsymbol{ls}. $$

We can calculate the expectation value of the latter using the fact that

$$ \xi(\boldsymbol{r})\boldsymbol{ls}=\frac{1}{2}\xi(\boldsymbol{r})\left(\boldsymbol{j}^2-\boldsymbol{l}^2-\boldsymbol{s}^2\right). $$

For a single-particle state with quantum numbers $nlj$ (we suppress $s$ and $m_j$), with $s=1/2$, we obtain the single-particle energies

$$ \varepsilon_{nlj} = \varepsilon_{nlj}^{(0)}+\Delta\varepsilon_{nlj}, $$

with $\varepsilon_{nlj}^{(0)}$ being the single-particle energy obtained with $\hat{h}_0(x)$ and

$$ \Delta\varepsilon_{nlj}=\frac{C}{2}\left(j(j+1)-l(l+1)-\frac{3}{4}\right). $$

The spin-orbit force gives thus an additional contribution to the energy

$$ \Delta\varepsilon_{nlj}=\frac{C}{2}\left(j(j+1)-l(l+1)-\frac{3}{4}\right), $$

which lifts the degeneracy we have seen before in the harmonic oscillator or Woods-Saxon potentials. The value $C$ is the radial integral involving $\xi(\boldsymbol{r})$. Depending on the value of $j=l\pm 1/2$, we obtain

$$ \Delta\varepsilon_{nlj=l-1/2}=\frac{C}{2}l, $$

or

$$ \Delta\varepsilon_{nlj=l+1/2}=-\frac{C}{2}(l+1), $$

clearly lifting the degeneracy. Note well that till now we have simply postulated the spin-orbit force in ad hoc way. Later, we will see how this term arises from the two-nucleon force in a natural way.

With the spin-orbit force, we can modify our Woods-Saxon potential to

$$ \hat{u}_{\mathrm{ext}}(r)=-\frac{V_0}{1+\exp{(r-R)/a}}+V_{so}(r)\boldsymbol{ls}, $$

with

$$ V_{so}(r) = V_{so}\frac{1}{r}\frac{d f_{so}(r)}{dr}, $$

where we have

$$ f_{so}(r) = \frac{1}{1+\exp{(r-R_{so})/a_{so}}}. $$

We can also add, in case of proton, a Coulomb potential. The Woods-Saxon potential has been widely used in parametrizations of effective single-particle potentials. However, as was the case with the harmonic oscillator, none of these potentials are linked directly to the nuclear forces. Our next step is to build a mean field based on the nucleon-nucleon interaction. This will lead us to our first and simplest many-body theory, Hartree-Fock theory.

The Woods-Saxon potential does allow for closed-form or analytical solutions of the eigenvalue problem

$$ \hat{h}_0(x_i)\psi_{\alpha}(x_i)=\varepsilon_{\alpha}\psi_{\alpha}(x_i). $$

For the harmonic oscillator in three dimensions we have closed-form expressions for the energies and analytical solutions for the eigenstates, with the latter given by either Hermite polynomials (cartesian coordinates) or Laguerre polynomials (spherical coordinates).

To solve the above equation is however rather straightforward numerically.

Numerical solution of the single-particle Schroedinger equation

We will illustrate the numerical solution of Schroedinger's equation by solving it for the harmonic oscillator in three dimensions. It is straightforward to change the harmonic oscillator potential with a Woods-Saxon potential, or any other type of potentials.

We are interested in the solution of the radial part of Schroedinger's equation for one nucleon. The angular momentum part is given by the so-called Spherical harmonics.

The radial equation reads

$$ -\frac{\hbar^2}{2 m} \left ( \frac{1}{r^2} \frac{d}{dr} r^2 \frac{d}{dr} - \frac{l (l + 1)}{r^2} \right )R(r) + V(r) R(r) = E R(r). $$

In our case $V(r)$ is the harmonic oscillator potential $(1/2)kr^2$ with $k=m\omega^2$ and $E$ is the energy of the harmonic oscillator in three dimensions. The oscillator frequency is $\omega$ and the energies are

$$ E_{nl}= \hbar \omega \left(2n+l+\frac{3}{2}\right), $$

with $n=0,1,2,\dots$ and $l=0,1,2,\dots$.

Since we have made a transformation to spherical coordinates it means that $r\in [0,\infty)$.
The quantum number $l$ is the orbital momentum of the nucleon. Then we substitute $R(r) = (1/r) u(r)$ and obtain

$$ -\frac{\hbar^2}{2 m} \frac{d^2}{dr^2} u(r) + \left ( V(r) + \frac{l (l + 1)}{r^2}\frac{\hbar^2}{2 m} \right ) u(r) = E u(r) . $$

The boundary conditions are $u(0)=0$ and $u(\infty)=0$.

We introduce a dimensionless variable $\rho = (1/\alpha) r$ where $\alpha$ is a constant with dimension length and get

$$ -\frac{\hbar^2}{2 m \alpha^2} \frac{d^2}{d\rho^2} u(\rho) + \left ( V(\rho) + \frac{l (l + 1)}{\rho^2} \frac{\hbar^2}{2 m\alpha^2} \right ) u(\rho) = E u(\rho) . $$

Let us specialize to $l=0$. Inserting $V(\rho) = (1/2) k \alpha^2\rho^2$ we end up with

$$ -\frac{\hbar^2}{2 m \alpha^2} \frac{d^2}{d\rho^2} u(\rho) + \frac{k}{2} \alpha^2\rho^2u(\rho) = E u(\rho) . $$

We multiply thereafter with $2m\alpha^2/\hbar^2$ on both sides and obtain

$$ -\frac{d^2}{d\rho^2} u(\rho) + \frac{mk}{\hbar^2} \alpha^4\rho^2u(\rho) = \frac{2m\alpha^2}{\hbar^2}E u(\rho) . $$

We have thus

$$ -\frac{d^2}{d\rho^2} u(\rho) + \frac{mk}{\hbar^2} \alpha^4\rho^2u(\rho) = \frac{2m\alpha^2}{\hbar^2}E u(\rho) . $$

The constant $\alpha$ can now be fixed so that

$$ \frac{mk}{\hbar^2} \alpha^4 = 1, $$

or

$$ \alpha = \left(\frac{\hbar^2}{mk}\right)^{1/4}. $$

Defining

$$ \lambda = \frac{2m\alpha^2}{\hbar^2}E, $$

we can rewrite Schroedinger's equation as

$$ -\frac{d^2}{d\rho^2} u(\rho) + \rho^2u(\rho) = \lambda u(\rho) . $$

This is the first equation to solve numerically. In three dimensions the eigenvalues for $l=0$ are $\lambda_0=3,\lambda_1=7,\lambda_2=11,\dots .$

We use the standard expression for the second derivative of a function $u$

$$ \begin{equation} u''=\frac{u(\rho+h) -2u(\rho) +u(\rho-h)}{h^2} +O(h^2), \label{eq:diffoperation} \tag{6} \end{equation} $$

where $h$ is our step. Next we define minimum and maximum values for the variable $\rho$, $\rho_{\mathrm{min}}=0$ and $\rho_{\mathrm{max}}$, respectively. You need to check your results for the energies against different values $\rho_{\mathrm{max}}$, since we cannot set $\rho_{\mathrm{max}}=\infty$.

With a given number of steps, $n_{\mathrm{step}}$, we then define the step $h$ as

$$ h=\frac{\rho_{\mathrm{max}}-\rho_{\mathrm{min}} }{n_{\mathrm{step}}}. $$

Define an arbitrary value of $\rho$ as

$$ \rho_i= \rho_{\mathrm{min}} + ih \hspace{1cm} i=0,1,2,\dots , n_{\mathrm{step}} $$

we can rewrite the Schroedinger equation for $\rho_i$ as

$$ -\frac{u(\rho_i+h) -2u(\rho_i) +u(\rho_i-h)}{h^2}+\rho_i^2u(\rho_i) = \lambda u(\rho_i), $$

or in a more compact way

$$ -\frac{u_{i+1} -2u_i +u_{i-1}}{h^2}+\rho_i^2u_i=-\frac{u_{i+1} -2u_i +u_{i-1} }{h^2}+V_iu_i = \lambda u_i, $$

where $V_i=\rho_i^2$ is the harmonic oscillator potential.

Define first the diagonal matrix element

$$ d_i=\frac{2}{h^2}+V_i, $$

and the non-diagonal matrix element

$$ e_i=-\frac{1}{h^2}. $$

In this case the non-diagonal matrix elements are given by a mere constant. All non-diagonal matrix elements are equal.

With these definitions the Schroedinger equation takes the following form

$$ d_iu_i+e_{i-1}u_{i-1}+e_{i+1}u_{i+1} = \lambda u_i, $$

where $u_i$ is unknown. We can write the latter equation as a matrix eigenvalue problem

$$ \begin{equation} \left( \begin{array}{ccccccc} d_1 & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & d_2 & e_2 & 0 & \dots &0 &0 \\ 0 & e_2 & d_3 & e_3 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &d_{n_{\mathrm{step}}-2} & e_{n_{\mathrm{step}}-1}\\ 0 & \dots & \dots & \dots &\dots &e_{n_{\mathrm{step}}-1} & d_{n_{\mathrm{step}}-1} \end{array} \right) \left( \begin{array}{c} u_{1} \\ u_{2} \\ \dots\\ \dots\\ \dots\\ u_{n_{\mathrm{step}}-1} \end{array} \right)=\lambda \left( \begin{array}{c} u_{1} \\ u_{2} \\ \dots\\ \dots\\ \dots\\ u_{n_{\mathrm{step}}-1} \end{array} \right) \label{eq:sematrix} \tag{7} \end{equation} $$

or if we wish to be more detailed, we can write the tridiagonal matrix as

$$ \begin{equation} \left( \begin{array}{ccccccc} \frac{2}{h^2}+V_1 & -\frac{1}{h^2} & 0 & 0 & \dots &0 & 0 \\ -\frac{1}{h^2} & \frac{2}{h^2}+V_2 & -\frac{1}{h^2} & 0 & \dots &0 &0 \\ 0 & -\frac{1}{h^2} & \frac{2}{h^2}+V_3 & -\frac{1}{h^2} &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &\frac{2}{h^2}+V_{n_{\mathrm{step}}-2} & -\frac{1}{h^2}\\ 0 & \dots & \dots & \dots &\dots &-\frac{1}{h^2} & \frac{2}{h^2}+V_{n_{\mathrm{step}}-1} \end{array} \right) \label{eq:matrixse} \tag{8} \end{equation} $$

Recall that the solutions are known via the boundary conditions at $i=n_{\mathrm{step}}$ and at the other end point, that is for $\rho_0$. The solution is zero in both cases.

The following python program is an example of how one can obtain the eigenvalues for a single-nucleon moving in a harmonic oscillator potential. It is rather easy to change the onebody-potential with ones like a Woods-Saxon potential.

• The c++ and Fortran versions of this program can be found at https://github.com/NuclearStructure/PHY981/tree/master/doc/pub/spdata/programs.

• The c++ program uses the c++ library armadillo, see http://arma.sourceforge.net/.

• To install armadillo see the guidelines at http://www.uio.no/studier/emner/matnat/fys/FYS4411/v14/guides/installing-armadillo/.

• For mac users I recommend using brew, see http://brew.sh/.

• If you use ipython notebook, you can run c++ programs following the instructions at http://nbviewer.ipython.org/github/dragly/cppmagic/blob/master/example.ipynb

The code sets up the Hamiltonian matrix by defining the the minimun and maximum values of $r$ with a maximum value of integration points. These are set in the initialization function. It plots the eigenfunctions of the three lowest eigenstates.

Exercise 1: Masses and binding energies

The data on binding energies can be found in the file bedata.dat at the github address of the course, see https://github.com/NuclearStructure/PHY981/tree/master/doc/pub/spdata/programs

a) Write a small program which reads in the proton and neutron numbers and the binding energies and make a plot of all neutron separation energies for the chain of oxygen (O), calcium (Ca), nickel (Ni), tin (Sn) and lead (Pb) isotopes, that is you need to plot

$$ S_n= BE(N,Z)-BE(N-1,Z). $$

Comment your results.

b) In the same figures, you should also include the liquid drop model results of Eq. (2.17) of Alex Brown's text, namely

$$ BE(N,Z)= \alpha_1A-\alpha_2A^{2/3}-\alpha_3\frac{Z^2}{A^{1/3}}-\alpha_4\frac{(N-Z)^2}{A}, $$

with $\alpha_1=15.49$ MeV, $\alpha_2=17.23$ MeV, $\alpha_3=0.697$ MeV and $\alpha_4=22.6$ MeV. Again, comment your results.

c) Make also a plot of the binding energies as function of $A$ using the data in the file on binding energies and the above liquid drop model. Make a figure similar to figure 2.5 of Alex Brown where you set the various parameters $\alpha_i=0$. Comment your results.

d) Use the liquid drop model to find the neutron drip lines for Z values up to 120. Analyze then the fluorine isotopes and find, where available the corresponding experimental data, and compare the liquid drop model predicition with experiment. Comment your results. A program example in C++ and the input data file bedata.dat can be found found at the github repository for the course

Exercise 2: Eigenvalues and eigenvectors for various single-particle potentials

The program for finding the eigenvalues of the harmonic oscillator are in the github folder https://github.com/NuclearStructure/PHY981/tree/master/doc/pub/spdata/programs.

You can use this program to solve the exercises below, or write your own using your preferred programming language, be it python, fortran or c++ or other languages. Here I will mainly provide fortran, python and c++.

a) Compute the eigenvalues of the five lowest states with a given orbital momentum and oscillator frequency $\omega$. Study these results as functions of the the maximum value of $r$ and the number of integration points $n$, starting with $r_{\mathrm{max}}=10$. Compare the computed ones with the exact values and comment your results.

b) Plot thereafter the eigenfunctions as functions of $r$ for the lowest-lying state with a given orbital momentum $l$.

c) Replace thereafter the harmonic oscillator potential with a Woods-Saxon potential using the parameters discussed above. Compute the lowest five eigenvalues and plot the eigenfunction of the lowest-lying state. How does this compare with the harmonic oscillator? Comment your results and possible implications for nuclear physics studies.

Exercise 3: Operators and Slater determinants

Consider the Slater determinant

$$ \Phi_{\lambda}^{AS}(x_{1}x_{2}\dots x_{N};\alpha_{1}\alpha_{2}\dots\alpha_{N}) =\frac{1}{\sqrt{N!}}\sum_{p}(-)^{p}P\prod_{i=1}^{N}\psi_{\alpha_{i}}(x_{i}). $$

where $P$ is an operator which permutes the coordinates of two particles. We have assumed here that the number of particles is the same as the number of available single-particle states, represented by the greek letters $\alpha_{1}\alpha_{2}\dots\alpha_{N}$.

a) Write out $\Phi^{AS}$ for $N=3$.

b) Show that

$$ \int dx_{1}dx_{2}\dots dx_{N}\left\vert \Phi_{\lambda}^{AS}(x_{1}x_{2}\dots x_{N};\alpha_{1}\alpha_{2}\dots\alpha_{N}) \right\vert^{2} = 1. $$

c) Define a general onebody operator $\hat{F} = \sum_{i}^N\hat{f}(x_{i})$ and a general twobody operator $\hat{G}=\sum_{i>j}^N\hat{g}(x_{i},x_{j})$ with $g$ being invariant under the interchange of the coordinates of particles $i$ and $j$. Calculate the matrix elements for a two-particle Slater determinant

$$ \langle\Phi_{\alpha_{1}\alpha_{2}}^{AS}|\hat{F}|\Phi_{\alpha_{1}\alpha_{2}}^{AS}\rangle, $$

and

$$ \langle\Phi_{\alpha_{1}\alpha_{2}}^{AS}|\hat{G}|\Phi_{\alpha_{1}\alpha_{2}}^{AS}\rangle. $$

Explain the short-hand notation for the Slater determinant. Which properties do you expect these operators to have in addition to an eventual permutation symmetry?

Exercise 4: First simple shell-model calculation

We will now consider a simple three-level problem, depicted in the figure below. This is our first and very simple model of a possible many-nucleon (or just fermion) problem and the shell-model. The single-particle states are labelled by the quantum number $p$ and can accomodate up to two single particles, viz., every single-particle state is doubly degenerate (you could think of this as one state having spin up and the other spin down). We let the spacing between the doubly degenerate single-particle states be constant, with value $d$. The first state has energy $d$. There are only three available single-particle states, $p=1$, $p=2$ and $p=3$, as illustrated in the figure.

a) How many two-particle Slater determinants can we construct in this space?

We limit ourselves to a system with only the two lowest single-particle orbits and two particles, $p=1$ and $p=2$. We assume that we can write the Hamiltonian as

$$ \hat{H}=\hat{H}_0+\hat{H}_I, $$

and that the onebody part of the Hamiltonian with single-particle operator $\hat{h}_0$ has the property

$$ \hat{h}_0\psi_{p\sigma} = p\times d \psi_{p\sigma}, $$

where we have added a spin quantum number $\sigma$. We assume also that the only two-particle states that can exist are those where two particles are in the same state $p$, as shown by the two possibilities to the left in the figure. The two-particle matrix elements of $\hat{H}_I$ have all a constant value, $-g$.

b) Show then that the Hamiltonian matrix can be written as

$$ \left(\begin{array}{cc}2d-g &-g \\ -g &4d-g \end{array}\right), $$

c) Find the eigenvalues and eigenvectors. What is mixing of the state with two particles in $p=2$ to the wave function with two-particles in $p=1$? Discuss your results in terms of a linear combination of Slater determinants.

d) Add the possibility that the two particles can be in the state with $p=3$ as well and find the Hamiltonian matrix, the eigenvalues and the eigenvectors. We still insist that we only have two-particle states composed of two particles being in the same level $p$. You can diagonalize numerically your $3\times 3$ matrix. This simple model catches several birds with a stone. It demonstrates how we can build linear combinations of Slater determinants and interpret these as different admixtures to a given state. It represents also the way we are going to interpret these contributions. The two-particle states above $p=1$ will be interpreted as excitations from the ground state configuration, $p=1$ here. The reliability of this ansatz for the ground state, with two particles in $p=1$, depends on the strength of the interaction $g$ and the single-particle spacing $d$. Finally, this model is a simple schematic ansatz for studies of pairing correlations and thereby superfluidity/superconductivity
in fermionic systems.

Schematic plot of the possible single-particle levels with double degeneracy. The filled circles indicate occupied particle states. The spacing between each level $p$ is constant in this picture. We show some possible two-particle states.

Mean-field approaches

Derivation of Hartree-Fock equations in coordinate space

Hartree-Fock (HF) theory is an algorithm for finding an approximative expression for the ground state of a given Hamiltonian. The basic ingredients are

• Define a single-particle basis $\{\psi_{\alpha}\}$ so that

$$ \hat{h}^{\mathrm{HF}}\psi_{\alpha} = \varepsilon_{\alpha}\psi_{\alpha} $$

with the Hartree-Fock Hamiltonian defined as

$$ \hat{h}^{\mathrm{HF}}=\hat{t}+\hat{u}_{\mathrm{ext}}+\hat{u}^{\mathrm{HF}} $$

• The term $\hat{u}^{\mathrm{HF}}$ is a single-particle potential to be determined by the HF algorithm.

  • The HF algorithm means to choose $\hat{u}^{\mathrm{HF}}$ in order to have

$$ \langle \hat{H} \rangle = E^{\mathrm{HF}}= \langle \Phi_0 | \hat{H}|\Phi_0 \rangle $$

that is to find a local minimum with a Slater determinant $\Phi_0$ being the ansatz for the ground state.

• The variational principle ensures that $E^{\mathrm{HF}} \ge E_0$, with $E_0$ the exact ground state energy.

We will show that the Hartree-Fock Hamiltonian $\hat{h}^{\mathrm{HF}}$ equals our definition of the operator $\hat{f}$ discussed in connection with the new definition of the normal-ordered Hamiltonian (see later lectures), that is we have, for a specific matrix element

$$ \langle p |\hat{h}^{\mathrm{HF}}| q \rangle =\langle p |\hat{f}| q \rangle=\langle p|\hat{t}+\hat{u}_{\mathrm{ext}}|q \rangle +\sum_{i\le F} \langle pi | \hat{V} | qi\rangle_{AS}, $$

meaning that

$$ \langle p|\hat{u}^{\mathrm{HF}}|q\rangle = \sum_{i\le F} \langle pi | \hat{V} | qi\rangle_{AS}. $$

The so-called Hartree-Fock potential $\hat{u}^{\mathrm{HF}}$ brings an explicit medium dependence due to the summation over all single-particle states below the Fermi level $F$. It brings also in an explicit dependence on the two-body interaction (in nuclear physics we can also have complicated three- or higher-body forces). The two-body interaction, with its contribution from the other bystanding fermions, creates an effective mean field in which a given fermion moves, in addition to the external potential $\hat{u}_{\mathrm{ext}}$ which confines the motion of the fermion. For systems like nuclei, there is no external confining potential. Nuclei are examples of self-bound systems, where the binding arises due to the intrinsic nature of the strong force. For nuclear systems thus, there would be no external one-body potential in the Hartree-Fock Hamiltonian.

The calculus of variations involves problems where the quantity to be minimized or maximized is an integral.

In the general case we have an integral of the type

$$ E[\Phi]= \int_a^b f(\Phi(x),\frac{\partial \Phi}{\partial x},x)dx, $$

where $E$ is the quantity which is sought minimized or maximized. The problem is that although $f$ is a function of the variables $\Phi$, $\partial \Phi/\partial x$ and $x$, the exact dependence of $\Phi$ on $x$ is not known. This means again that even though the integral has fixed limits $a$ and $b$, the path of integration is not known. In our case the unknown quantities are the single-particle wave functions and we wish to choose an integration path which makes the functional $E[\Phi]$ stationary. This means that we want to find minima, or maxima or saddle points. In physics we search normally for minima. Our task is therefore to find the minimum of $E[\Phi]$ so that its variation $\delta E$ is zero subject to specific constraints. In our case the constraints appear as the integral which expresses the orthogonality of the single-particle wave functions. The constraints can be treated via the technique of Lagrangian multipliers

Let us specialize to the expectation value of the energy for one particle in three-dimensions. This expectation value reads

$$ E=\int dxdydz \psi^*(x,y,z) \hat{H} \psi(x,y,z), $$

with the constraint

$$ \int dxdydz \psi^*(x,y,z) \psi(x,y,z)=1, $$

and a Hamiltonian

$$ \hat{H}=-\frac{1}{2}\nabla^2+V(x,y,z). $$

We will, for the sake of notational convenience, skip the variables $x,y,z$ below, and write for example $V(x,y,z)=V$.

The integral involving the kinetic energy can be written as, with the function $\psi$ vanishing strongly for large values of $x,y,z$ (given here by the limits $a$ and $b$),

$$ \int_a^b dxdydz \psi^* \left(-\frac{1}{2}\nabla^2\right) \psi dxdydz = \psi^*\nabla\psi|_a^b+\int_a^b dxdydz\frac{1}{2}\nabla\psi^*\nabla\psi. $$

We will drop the limits $a$ and $b$ in the remaining discussion. Inserting this expression into the expectation value for the energy and taking the variational minimum we obtain

$$ \delta E = \delta \left\{\int dxdydz\left( \frac{1}{2}\nabla\psi^*\nabla\psi+V\psi^*\psi\right)\right\} = 0. $$

The constraint appears in integral form as

$$ \int dxdydz \psi^* \psi=\mathrm{constant}, $$

and multiplying with a Lagrangian multiplier $\lambda$ and taking the variational minimum we obtain the final variational equation

$$ \delta \left\{\int dxdydz\left( \frac{1}{2}\nabla\psi^*\nabla\psi+V\psi^*\psi-\lambda\psi^*\psi\right)\right\} = 0. $$

We introduce the function $f$

$$ f = \frac{1}{2}\nabla\psi^*\nabla\psi+V\psi^*\psi-\lambda\psi^*\psi= \frac{1}{2}(\psi^*_x\psi_x+\psi^*_y\psi_y+\psi^*_z\psi_z)+V\psi^*\psi-\lambda\psi^*\psi, $$

where we have skipped the dependence on $x,y,z$ and introduced the shorthand $\psi_x$, $\psi_y$ and $\psi_z$ for the various derivatives.

For $\psi^*$ the Euler-Lagrange equations yield

$$ \frac{\partial f}{\partial \psi^*}- \frac{\partial }{\partial x}\frac{\partial f}{\partial \psi^*_x}-\frac{\partial }{\partial y}\frac{\partial f}{\partial \psi^*_y}-\frac{\partial }{\partial z}\frac{\partial f}{\partial \psi^*_z}=0, $$

which results in

$$ -\frac{1}{2}(\psi_{xx}+\psi_{yy}+\psi_{zz})+V\psi=\lambda \psi. $$

We can then identify the Lagrangian multiplier as the energy of the system. The last equation is nothing but the standard Schroedinger equation and the variational approach discussed here provides a powerful method for obtaining approximate solutions of the wave function.

Definitions and notations

Before we proceed we need some definitions. We will assume that the interacting part of the Hamiltonian can be approximated by a two-body interaction. This means that our Hamiltonian is written as the sum of some onebody part and a twobody part

$$ \begin{equation} \hat{H} = \hat{H}_0 + \hat{H}_I = \sum_{i=1}^A \hat{h}_0(x_i) + \sum_{i < j}^A \hat{v}(r_{ij}), \label{Hnuclei} \tag{9} \end{equation} $$

with

$$ \begin{equation} H_0=\sum_{i=1}^A \hat{h}_0(x_i). \label{hinuclei} \tag{10} \end{equation} $$

The onebody part $u_{\mathrm{ext}}(x_i)$ is normally approximated by a harmonic oscillator potential or the Coulomb interaction an electron feels from the nucleus. However, other potentials are fully possible, such as one derived from the self-consistent solution of the Hartree-Fock equations to be discussed here.

Our Hamiltonian is invariant under the permutation (interchange) of two particles. Since we deal with fermions however, the total wave function is antisymmetric. Let $\hat{P}$ be an operator which interchanges two particles. Due to the symmetries we have ascribed to our Hamiltonian, this operator commutes with the total Hamiltonian,

$$ [\hat{H},\hat{P}] = 0, $$

meaning that $\Psi_{\lambda}(x_1, x_2, \dots , x_A)$ is an eigenfunction of $\hat{P}$ as well, that is

$$ \hat{P}_{ij}\Psi_{\lambda}(x_1, x_2, \dots,x_i,\dots,x_j,\dots,x_A)= \beta\Psi_{\lambda}(x_1, x_2, \dots,x_i,\dots,x_j,\dots,x_A), $$

where $\beta$ is the eigenvalue of $\hat{P}$. We have introduced the suffix $ij$ in order to indicate that we permute particles $i$ and $j$. The Pauli principle tells us that the total wave function for a system of fermions has to be antisymmetric, resulting in the eigenvalue $\beta = -1$.

In our case we assume that we can approximate the exact eigenfunction with a Slater determinant

$$ \begin{equation} \Phi(x_1, x_2,\dots ,x_A,\alpha,\beta,\dots, \sigma)=\frac{1}{\sqrt{A!}} \left| \begin{array}{ccccc} \psi_{\alpha}(x_1)& \psi_{\alpha}(x_2)& \dots & \dots & \psi_{\alpha}(x_A)\\ \psi_{\beta}(x_1)&\psi_{\beta}(x_2)& \dots & \dots & \psi_{\beta}(x_A)\\ \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots \\ \psi_{\sigma}(x_1)&\psi_{\sigma}(x_2)& \dots & \dots & \psi_{\sigma}(x_A)\end{array} \right|, \label{eq:HartreeFockDet} \tag{11} \end{equation} $$

where $x_i$ stand for the coordinates and spin values of a particle $i$ and $\alpha,\beta,\dots, \gamma$ are quantum numbers needed to describe remaining quantum numbers.

The single-particle function $\psi_{\alpha}(x_i)$ are eigenfunctions of the onebody Hamiltonian $h_i$, that is

$$ \hat{h}_0(x_i)=\hat{t}(x_i) + \hat{u}_{\mathrm{ext}}(x_i), $$

with eigenvalues

$$ \hat{h}_0(x_i) \psi_{\alpha}(x_i)=\left(\hat{t}(x_i) + \hat{u}_{\mathrm{ext}}(x_i)\right)\psi_{\alpha}(x_i)=\varepsilon_{\alpha}\psi_{\alpha}(x_i). $$

The energies $\varepsilon_{\alpha}$ are the so-called non-interacting single-particle energies, or unperturbed energies. The total energy is in this case the sum over all single-particle energies, if no two-body or more complicated many-body interactions are present.

Let us denote the ground state energy by $E_0$. According to the variational principle we have

$$ E_0 \le E[\Phi] = \int \Phi^*\hat{H}\Phi d\mathbf{\tau} $$

where $\Phi$ is a trial function which we assume to be normalized

$$ \int \Phi^*\Phi d\mathbf{\tau} = 1, $$

where we have used the shorthand $d\mathbf{\tau}=dx_1dr_2\dots dr_A$.

In the Hartree-Fock method the trial function is the Slater determinant of Eq. (11) which can be rewritten as

$$ \Phi(x_1,x_2,\dots,x_A,\alpha,\beta,\dots,\nu) = \frac{1}{\sqrt{A!}}\sum_{P} (-)^P\hat{P}\psi_{\alpha}(x_1) \psi_{\beta}(x_2)\dots\psi_{\nu}(x_A)=\sqrt{A!}\hat{A}\Phi_H, $$

where we have introduced the antisymmetrization operator $\hat{A}$ defined by the summation over all possible permutations of two particles.

It is defined as

$$ \begin{equation} \hat{A} = \frac{1}{A!}\sum_{p} (-)^p\hat{P}, \label{antiSymmetryOperator} \tag{12} \end{equation} $$

with $p$ standing for the number of permutations. We have introduced for later use the so-called Hartree-function, defined by the simple product of all possible single-particle functions

$$ \Phi_H(x_1,x_2,\dots,x_A,\alpha,\beta,\dots,\nu) = \psi_{\alpha}(x_1) \psi_{\beta}(x_2)\dots\psi_{\nu}(x_A). $$

Both $\hat{H}_0$ and $\hat{H}_I$ are invariant under all possible permutations of any two particles and hence commute with $\hat{A}$

$$ \begin{equation} [H_0,\hat{A}] = [H_I,\hat{A}] = 0. \label{commutionAntiSym} \tag{13} \end{equation} $$

Furthermore, $\hat{A}$ satisfies

$$ \begin{equation} \hat{A}^2 = \hat{A}, \label{AntiSymSquared} \tag{14} \end{equation} $$

since every permutation of the Slater determinant reproduces it.

The expectation value of $\hat{H}_0$

$$ \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau} = A! \int \Phi_H^*\hat{A}\hat{H}_0\hat{A}\Phi_H d\mathbf{\tau} $$

is readily reduced to

$$ \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau} = A! \int \Phi_H^*\hat{H}_0\hat{A}\Phi_H d\mathbf{\tau}, $$

where we have used Eqs. (13) and (14). The next step is to replace the antisymmetrization operator by its definition and to replace $\hat{H}_0$ with the sum of one-body operators

$$ \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau} = \sum_{i=1}^A \sum_{p} (-)^p\int \Phi_H^*\hat{h}_0\hat{P}\Phi_H d\mathbf{\tau}. $$

The integral vanishes if two or more particles are permuted in only one of the Hartree-functions $\Phi_H$ because the individual single-particle wave functions are orthogonal. We obtain then

$$ \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau}= \sum_{i=1}^A \int \Phi_H^*\hat{h}_0\Phi_H d\mathbf{\tau}. $$

Orthogonality of the single-particle functions allows us to further simplify the integral, and we arrive at the following expression for the expectation values of the sum of one-body Hamiltonians

$$ \begin{equation} \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau} = \sum_{\mu=1}^A \int \psi_{\mu}^*(x)\hat{h}_0\psi_{\mu}(x)dx d\mathbf{r}. \label{H1Expectation} \tag{15} \end{equation} $$

We introduce the following shorthand for the above integral

$$ \langle \mu | \hat{h}_0 | \mu \rangle = \int \psi_{\mu}^*(x)\hat{h}_0\psi_{\mu}(x)dx, $$

and rewrite Eq. (15) as

$$ \begin{equation} \int \Phi^*\hat{H}_0\Phi d\tau = \sum_{\mu=1}^A \langle \mu | \hat{h}_0 | \mu \rangle. \label{H1Expectation1} \tag{16} \end{equation} $$

The expectation value of the two-body part of the Hamiltonian is obtained in a similar manner. We have

$$ \int \Phi^*\hat{H}_I\Phi d\mathbf{\tau} = A! \int \Phi_H^*\hat{A}\hat{H}_I\hat{A}\Phi_H d\mathbf{\tau}, $$

which reduces to

$$ \int \Phi^*\hat{H}_I\Phi d\mathbf{\tau} = \sum_{i\le j=1}^A \sum_{p} (-)^p\int \Phi_H^*\hat{v}(r_{ij})\hat{P}\Phi_H d\mathbf{\tau}, $$

by following the same arguments as for the one-body Hamiltonian.

Because of the dependence on the inter-particle distance $r_{ij}$, permutations of any two particles no longer vanish, and we get

$$ \int \Phi^*\hat{H}_I\Phi d\mathbf{\tau} = \sum_{i < j=1}^A \int \Phi_H^*\hat{v}(r_{ij})(1-P_{ij})\Phi_H d\mathbf{\tau}. $$

where $P_{ij}$ is the permutation operator that interchanges particle $i$ and particle $j$. Again we use the assumption that the single-particle wave functions are orthogonal.

We obtain

$$ \begin{equation} \begin{split} \int \Phi^*\hat{H}_I\Phi d\mathbf{\tau} = \frac{1}{2}\sum_{\mu=1}^A\sum_{\nu=1}^A &\left[ \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\mu}(x_i)\psi_{\nu}(x_j) dx_idx_j \right.\\ &\left. - \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j) \hat{v}(r_{ij})\psi_{\nu}(x_i)\psi_{\mu}(x_j) dx_idx_j \right]. \label{H2Expectation} \tag{17} \end{split} \end{equation} $$

The first term is the so-called direct term. It is frequently also called the Hartree term, while the second is due to the Pauli principle and is called the exchange term or just the Fock term. The factor $1/2$ is introduced because we now run over all pairs twice.

The last equation allows us to introduce some further definitions.
The single-particle wave functions $\psi_{\mu}(x)$, defined by the quantum numbers $\mu$ and $x$ are defined as the overlap

$$ \psi_{\alpha}(x) = \langle x | \alpha \rangle . $$

We introduce the following shorthands for the above two integrals

$$ \langle \mu\nu|\hat{v}|\mu\nu\rangle = \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\mu}(x_i)\psi_{\nu}(x_j) dx_idx_j, $$

and

$$ \langle \mu\nu|\hat{v}|\nu\mu\rangle = \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j) \hat{v}(r_{ij})\psi_{\nu}(x_i)\psi_{\mu}(x_j) dx_idx_j. $$

Derivation of Hartree-Fock equations in coordinate space

Let us denote the ground state energy by $E_0$. According to the variational principle we have

$$ E_0 \le E[\Phi] = \int \Phi^*\hat{H}\Phi d\mathbf{\tau} $$

where $\Phi$ is a trial function which we assume to be normalized

$$ \int \Phi^*\Phi d\mathbf{\tau} = 1, $$

where we have used the shorthand $d\mathbf{\tau}=dx_1dx_2\dots dx_A$.

In the Hartree-Fock method the trial function is a Slater determinant which can be rewritten as

$$ \Psi(x_1,x_2,\dots,x_A,\alpha,\beta,\dots,\nu) = \frac{1}{\sqrt{A!}}\sum_{P} (-)^PP\psi_{\alpha}(x_1) \psi_{\beta}(x_2)\dots\psi_{\nu}(x_A)=\sqrt{A!}\hat{A}\Phi_H, $$

where we have introduced the anti-symmetrization operator $\hat{A}$ defined by the summation over all possible permutations p of two fermions. It is defined as

$$ \hat{A} = \frac{1}{A!}\sum_{p} (-)^p\hat{P}, $$

with the the Hartree-function given by the simple product of all possible single-particle function

$$ \Phi_H(x_1,x_2,\dots,x_A,\alpha,\beta,\dots,\nu) = \psi_{\alpha}(x_1) \psi_{\beta}(x_2)\dots\psi_{\nu}(x_A). $$

Our functional is written as

$$ E[\Phi] = \sum_{\mu=1}^A \int \psi_{\mu}^*(x_i)\hat{h}_0(x_i)\psi_{\mu}(x_i) dx_i + \frac{1}{2}\sum_{\mu=1}^A\sum_{\nu=1}^A \left[ \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\mu}(x_i)\psi_{\nu}(x_j)dx_idx_j- \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j) \hat{v}(r_{ij})\psi_{\nu}(x_i)\psi_{\mu}(x_j)dx_idx_j\right] $$

The more compact version reads

$$ E[\Phi] = \sum_{\mu}^A \langle \mu | \hat{h}_0 | \mu\rangle+ \frac{1}{2}\sum_{\mu\nu}^A\left[\langle \mu\nu |\hat{v}|\mu\nu\rangle-\langle \nu\mu |\hat{v}|\mu\nu\rangle\right]. $$

Since the interaction is invariant under the interchange of two particles it means for example that we have

$$ \langle \mu\nu|\hat{v}|\mu\nu\rangle = \langle \nu\mu|\hat{v}|\nu\mu\rangle, $$

or in the more general case

$$ \langle \mu\nu|\hat{v}|\sigma\tau\rangle = \langle \nu\mu|\hat{v}|\tau\sigma\rangle. $$

The direct and exchange matrix elements can be brought together if we define the antisymmetrized matrix element

$$ \langle \mu\nu|\hat{v}|\mu\nu\rangle_{AS}= \langle \mu\nu|\hat{v}|\mu\nu\rangle-\langle \mu\nu|\hat{v}|\nu\mu\rangle, $$

or for a general matrix element

$$ \langle \mu\nu|\hat{v}|\sigma\tau\rangle_{AS}= \langle \mu\nu|\hat{v}|\sigma\tau\rangle-\langle \mu\nu|\hat{v}|\tau\sigma\rangle. $$

It has the symmetry property

$$ \langle \mu\nu|\hat{v}|\sigma\tau\rangle_{AS}= -\langle \mu\nu|\hat{v}|\tau\sigma\rangle_{AS}=-\langle \nu\mu|\hat{v}|\sigma\tau\rangle_{AS}. $$

The antisymmetric matrix element is also hermitian, implying

$$ \langle \mu\nu|\hat{v}|\sigma\tau\rangle_{AS}= \langle \sigma\tau|\hat{v}|\mu\nu\rangle_{AS}. $$

With these notations we rewrite the Hartree-Fock functional as

$$ \begin{equation} \int \Phi^*\hat{H_I}\Phi d\mathbf{\tau} = \frac{1}{2}\sum_{\mu=1}^A\sum_{\nu=1}^A \langle \mu\nu|\hat{v}|\mu\nu\rangle_{AS}. \label{H2Expectation2} \tag{18} \end{equation} $$

Adding the contribution from the one-body operator $\hat{H}_0$ to (18) we obtain the energy functional

$$ \begin{equation} E[\Phi] = \sum_{\mu=1}^A \langle \mu | h | \mu \rangle + \frac{1}{2}\sum_{{\mu}=1}^A\sum_{{\nu}=1}^A \langle \mu\nu|\hat{v}|\mu\nu\rangle_{AS}. \label{FunctionalEPhi} \tag{19} \end{equation} $$

In our coordinate space derivations below we will spell out the Hartree-Fock equations in terms of their integrals.

If we generalize the Euler-Lagrange equations to more variables and introduce $N^2$ Lagrange multipliers which we denote by $\epsilon_{\mu\nu}$, we can write the variational equation for the functional of $E$

$$ \delta E - \sum_{\mu\nu}^A \epsilon_{\mu\nu} \delta \int \psi_{\mu}^* \psi_{\nu} = 0. $$

For the orthogonal wave functions $\psi_{i}$ this reduces to

$$ \delta E - \sum_{\mu=1}^A \epsilon_{\mu} \delta \int \psi_{\mu}^* \psi_{\mu} = 0. $$

Variation with respect to the single-particle wave functions $\psi_{\mu}$ yields then

1 7 6

< < < ! ! M A T H _ B L O C K

$$ \sum_{\mu=1}^A \int \psi_{\mu}^*\hat{h_0}(x_i)\delta\psi_{\mu} dx_i + \frac{1}{2}\sum_{{\mu}=1}^A\sum_{{\nu}=1}^A \left[ \int \psi_{\mu}^*\psi_{\nu}^*\hat{v}(r_{ij})\delta\psi_{\mu}\psi_{\nu} dx_idx_j- \int \psi_{\mu}^*\psi_{\nu}^*\hat{v}(r_{ij})\psi_{\nu}\delta\psi_{\mu} dx_idx_j \right]- \sum_{{\mu}=1}^A E_{\mu} \int \delta\psi_{\mu}^* \psi_{\mu}dx_i - \sum_{{\mu}=1}^A E_{\mu} \int \psi_{\mu}^* \delta\psi_{\mu}dx_i = 0. $$

Although the variations $\delta\psi$ and $\delta\psi^*$ are not independent, they may in fact be treated as such, so that the terms dependent on either $\delta\psi$ and $\delta\psi^*$ individually may be set equal to zero. To see this, simply replace the arbitrary variation $\delta\psi$ by $i\delta\psi$, so that $\delta\psi^*$ is replaced by $-i\delta\psi^*$, and combine the two equations. We thus arrive at the Hartree-Fock equations

$$ \begin{equation} \left[ -\frac{1}{2}\nabla_i^2+ \sum_{\nu=1}^A\int \psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\nu}(x_j)dx_j \right]\psi_{\mu}(x_i) - \left[ \sum_{{\nu}=1}^A \int\psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\mu}(x_j) dx_j\right] \psi_{\nu}(x_i) = \epsilon_{\mu} \psi_{\mu}(x_i). \label{eq:hartreefockcoordinatespace} \tag{20} \end{equation} $$

Notice that the integration $\int dx_j$ implies an integration over the spatial coordinates $\mathbf{r_j}$ and a summation over the spin-coordinate of fermion $j$. We note that the factor of $1/2$ in front of the sum involving the two-body interaction, has been removed. This is due to the fact that we need to vary both $\delta\psi_{\mu}^*$ and $\delta\psi_{\nu}^*$. Using the symmetry properties of the two-body interaction and interchanging $\mu$ and $\nu$ as summation indices, we obtain two identical terms.

The two first terms in the last equation are the one-body kinetic energy and the electron-nucleus potential. The third or direct term is the averaged electronic repulsion of the other electrons. As written, the term includes the self-interaction of electrons when $\mu=\nu$. The self-interaction is cancelled in the fourth term, or the exchange term. The exchange term results from our inclusion of the Pauli principle and the assumed determinantal form of the wave-function. Equation (20), in addition to the kinetic energy and the attraction from the atomic nucleus that confines the motion of a single electron, represents now the motion of a single-particle modified by the two-body interaction. The additional contribution to the Schroedinger equation due to the two-body interaction, represents a mean field set up by all the other bystanding electrons, the latter given by the sum over all single-particle states occupied by $N$ electrons.

The Hartree-Fock equation is an example of an integro-differential equation. These equations involve repeated calculations of integrals, in addition to the solution of a set of coupled differential equations. The Hartree-Fock equations can also be rewritten in terms of an eigenvalue problem. The solution of an eigenvalue problem represents often a more practical algorithm and the solution of coupled integro-differential equations. This alternative derivation of the Hartree-Fock equations is given below.

Analysis of Hartree-Fock equations in coordinate space

A theoretically convenient form of the Hartree-Fock equation is to regard the direct and exchange operator defined through

$$ V_{\mu}^{d}(x_i) = \int \psi_{\mu}^*(x_j) \hat{v}(r_{ij})\psi_{\mu}(x_j) dx_j $$

and

$$ V_{\mu}^{ex}(x_i) g(x_i) = \left(\int \psi_{\mu}^*(x_j) \hat{v}(r_{ij})g(x_j) dx_j \right)\psi_{\mu}(x_i), $$

respectively.

The function $g(x_i)$ is an arbitrary function, and by the substitution $g(x_i) = \psi_{\nu}(x_i)$ we get

$$ V_{\mu}^{ex}(x_i) \psi_{\nu}(x_i) = \left(\int \psi_{\mu}^*(x_j) \hat{v}(r_{ij})\psi_{\nu}(x_j) dx_j\right)\psi_{\mu}(x_i). $$

We may then rewrite the Hartree-Fock equations as

$$ \hat{h}^{HF}(x_i) \psi_{\nu}(x_i) = \epsilon_{\nu}\psi_{\nu}(x_i), $$

with

$$ \hat{h}^{HF}(x_i)= \hat{h}_0(x_i) + \sum_{\mu=1}^AV_{\mu}^{d}(x_i) - \sum_{\mu=1}^AV_{\mu}^{ex}(x_i), $$

and where $\hat{h}_0(i)$ is the one-body part. The latter is normally chosen as a part which yields solutions in closed form. The harmonic oscilltor is a classical problem thereof. We normally rewrite the last equation as

$$ \hat{h}^{HF}(x_i)= \hat{h}_0(x_i) + \hat{u}^{HF}(x_i). $$

Hartree-Fock by varying the coefficients of a wave function expansion

Another possibility is to expand the single-particle functions in a known basis and vary the coefficients, that is, the new single-particle wave function is written as a linear expansion in terms of a fixed chosen orthogonal basis (for example the well-known harmonic oscillator functions or the hydrogen-like functions etc). We define our new Hartree-Fock single-particle basis by performing a unitary transformation on our previous basis (labelled with greek indices) as

$$ \begin{equation} \psi_p^{HF} = \sum_{\lambda} C_{p\lambda}\phi_{\lambda}. \label{eq:newbasis} \tag{21} \end{equation} $$

In this case we vary the coefficients $C_{p\lambda}$. If the basis has infinitely many solutions, we need to truncate the above sum. We assume that the basis $\phi_{\lambda}$ is orthogonal. A unitary transformation keeps the orthogonality, as discussed in exercise 1 below.

It is normal to choose a single-particle basis defined as the eigenfunctions of parts of the full Hamiltonian. The typical situation consists of the solutions of the one-body part of the Hamiltonian, that is we have

$$ \hat{h}_0\phi_{\lambda}=\epsilon_{\lambda}\phi_{\lambda}. $$

The single-particle wave functions $\phi_{\lambda}({\bf r})$, defined by the quantum numbers $\lambda$ and ${\bf r}$ are defined as the overlap

$$ \phi_{\lambda}({\bf r}) = \langle {\bf r} | \lambda \rangle . $$

In our discussions hereafter we will use our definitions of single-particle states above and below the Fermi ($F$) level given by the labels $ijkl\dots \le F$ for so-called single-hole states and $abcd\dots > F$ for so-called particle states. For general single-particle states we employ the labels $pqrs\dots$.

In Eq. (19), restated here

$$ E[\Phi] = \sum_{\mu=1}^A \langle \mu | h | \mu \rangle + \frac{1}{2}\sum_{{\mu}=1}^A\sum_{{\nu}=1}^A \langle \mu\nu|\hat{v}|\mu\nu\rangle_{AS}, $$

we found the expression for the energy functional in terms of the basis function $\phi_{\lambda}({\bf r})$. We then varied the above energy functional with respect to the basis functions $|\mu \rangle$. Now we are interested in defining a new basis defined in terms of a chosen basis as defined in Eq. (21). We can then rewrite the energy functional as

$$ \begin{equation} E[\Phi^{HF}] = \sum_{i=1}^A \langle i | h | i \rangle + \frac{1}{2}\sum_{ij=1}^A\langle ij|\hat{v}|ij\rangle_{AS}, \label{FunctionalEPhi2} \tag{22} \end{equation} $$

where $\Phi^{HF}$ is the new Slater determinant defined by the new basis of Eq. (21).

Using Eq. (21) we can rewrite Eq. (22) as

$$ \begin{equation} E[\Psi] = \sum_{i=1}^A \sum_{\alpha\beta} C^*_{i\alpha}C_{i\beta}\langle \alpha | h | \beta \rangle + \frac{1}{2}\sum_{ij=1}^A\sum_{{\alpha\beta\gamma\delta}} C^*_{i\alpha}C^*_{j\beta}C_{i\gamma}C_{j\delta}\langle \alpha\beta|\hat{v}|\gamma\delta\rangle_{AS}. \label{FunctionalEPhi3} \tag{23} \end{equation} $$

We wish now to minimize the above functional. We introduce again a set of Lagrange multipliers, noting that since $\langle i | j \rangle = \delta_{i,j}$ and $\langle \alpha | \beta \rangle = \delta_{\alpha,\beta}$, the coefficients $C_{i\gamma}$ obey the relation

$$ \langle i | j \rangle=\delta_{i,j}=\sum_{\alpha\beta} C^*_{i\alpha}C_{i\beta}\langle \alpha | \beta \rangle= \sum_{\alpha} C^*_{i\alpha}C_{i\alpha}, $$

which allows us to define a functional to be minimized that reads

$$ \begin{equation} F[\Phi^{HF}]=E[\Phi^{HF}] - \sum_{i=1}^A\epsilon_i\sum_{\alpha} C^*_{i\alpha}C_{i\alpha}. \end{equation} $$

Minimizing with respect to $C^*_{i\alpha}$, remembering that the equations for $C^*_{i\alpha}$ and $C_{i\alpha}$ can be written as two independent equations, we obtain

$$ \frac{d}{dC^*_{i\alpha}}\left[ E[\Phi^{HF}] - \sum_{j}\epsilon_j\sum_{\alpha} C^*_{j\alpha}C_{j\alpha}\right]=0, $$

which yields for every single-particle state $i$ and index $\alpha$ (recalling that the coefficients $C_{i\alpha}$ are matrix elements of a unitary (or orthogonal for a real symmetric matrix) matrix) the following Hartree-Fock equations

$$ \sum_{\beta} C_{i\beta}\langle \alpha | h | \beta \rangle+ \sum_{j=1}^A\sum_{\beta\gamma\delta} C^*_{j\beta}C_{j\delta}C_{i\gamma}\langle \alpha\beta|\hat{v}|\gamma\delta\rangle_{AS}=\epsilon_i^{HF}C_{i\alpha}. $$

We can rewrite this equation as (changing dummy variables)

$$ \sum_{\beta} \left\{\langle \alpha | h | \beta \rangle+ \sum_{j}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}\right\}C_{i\beta}=\epsilon_i^{HF}C_{i\alpha}. $$

Note that the sums over greek indices run over the number of basis set functions (in principle an infinite number).

Defining

$$ h_{\alpha\beta}^{HF}=\langle \alpha | h | \beta \rangle+ \sum_{j=1}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}, $$

we can rewrite the new equations as

$$ \begin{equation} \sum_{\gamma}h_{\alpha\beta}^{HF}C_{i\beta}=\epsilon_i^{HF}C_{i\alpha}. \label{eq:newhf} \tag{24} \end{equation} $$

The latter is nothing but a standard eigenvalue problem. Compared with Eq. (20), we see that we do not need to compute any integrals in an iterative procedure for solving the equations. It suffices to tabulate the matrix elements $\langle \alpha | h | \beta \rangle$ and $\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}$ once and for all. Successive iterations require thus only a look-up in tables over one-body and two-body matrix elements. These details will be discussed below when we solve the Hartree-Fock equations numerical.

Our Hartree-Fock matrix is thus

$$ \hat{h}_{\alpha\beta}^{HF}=\langle \alpha | \hat{h}_0 | \beta \rangle+ \sum_{j=1}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}. $$

The Hartree-Fock equations are solved in an iterative waym starting with a guess for the coefficients $C_{j\gamma}=\delta_{j,\gamma}$ and solving the equations by diagonalization till the new single-particle energies $\epsilon_i^{\mathrm{HF}}$ do not change anymore by a prefixed quantity.

Normally we assume that the single-particle basis $|\beta\rangle$ forms an eigenbasis for the operator $\hat{h}_0$, meaning that the Hartree-Fock matrix becomes

$$ \hat{h}_{\alpha\beta}^{HF}=\epsilon_{\alpha}\delta_{\alpha,\beta}+ \sum_{j=1}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}. $$

The Hartree-Fock eigenvalue problem

$$ \sum_{\beta}\hat{h}_{\alpha\beta}^{HF}C_{i\beta}=\epsilon_i^{\mathrm{HF}}C_{i\alpha}, $$

can be written out in a more compact form as

$$ \hat{h}^{HF}\hat{C}=\epsilon^{\mathrm{HF}}\hat{C}. $$

Hartree-Fock algorithm

The Hartree-Fock equations are, in their simplest form, solved in an iterative way, starting with a guess for the coefficients $C_{i\alpha}$. We label the coefficients as $C_{i\alpha}^{(n)}$, where the subscript $n$ stands for iteration $n$. To set up the algorithm we can proceed as follows:

• We start with a guess $C_{i\alpha}^{(0)}=\delta_{i,\alpha}$. Alternatively, we could have used random starting values as long as the vectors are normalized. Another possibility is to give states below the Fermi level a larger weight.

• The Hartree-Fock matrix simplifies then to (assuming that the coefficients $C_{i\alpha} $ are real)

$$ \hat{h}_{\alpha\beta}^{HF}=\epsilon_{\alpha}\delta_{\alpha,\beta}+ \sum_{j = 1}^A\sum_{\gamma\delta} C_{j\gamma}^{(0)}C_{j\delta}^{(0)}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}. $$

Solving the Hartree-Fock eigenvalue problem yields then new eigenvectors $C_{i\alpha}^{(1)}$ and eigenvalues $\epsilon_i^{HF(1)}$.

• With the new eigenvalues we can set up a new Hartree-Fock potential

$$ \sum_{j = 1}^A\sum_{\gamma\delta} C_{j\gamma}^{(1)}C_{j\delta}^{(1)}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}. $$

The diagonalization with the new Hartree-Fock potential yields new eigenvectors and eigenvalues. This process is continued till for example

$$ \frac{\sum_{p} |\epsilon_i^{(n)}-\epsilon_i^{(n-1)}|}{m} \le \lambda, $$

where $\lambda$ is a user prefixed quantity ($\lambda \sim 10^{-8}$ or smaller) and $p$ runs over all calculated single-particle energies and $m$ is the number of single-particle states.

Analysis of Hartree-Fock equations and Koopman's theorem

We can rewrite the ground state energy by adding and subtracting $\hat{u}^{HF}(x_i)$

$$ E_0^{HF} =\langle \Phi_0 | \hat{H} | \Phi_0\rangle = \sum_{i\le F}^A \langle i | \hat{h}_0 +\hat{u}^{HF}| j\rangle+ \frac{1}{2}\sum_{i\le F}^A\sum_{j \le F}^A\left[\langle ij |\hat{v}|ij \rangle-\langle ij|\hat{v}|ji\rangle\right]-\sum_{i\le F}^A \langle i |\hat{u}^{HF}| i\rangle, $$

which results in

$$ E_0^{HF} = \sum_{i\le F}^A \varepsilon_i^{HF} + \frac{1}{2}\sum_{i\le F}^A\sum_{j \le F}^A\left[\langle ij |\hat{v}|ij \rangle-\langle ij|\hat{v}|ji\rangle\right]-\sum_{i\le F}^A \langle i |\hat{u}^{HF}| i\rangle. $$

Our single-particle states $ijk\dots$ are now single-particle states obtained from the solution of the Hartree-Fock equations.

Using our definition of the Hartree-Fock single-particle energies we obtain then the following expression for the total ground-state energy

$$ E_0^{HF} = \sum_{i\le F}^A \varepsilon_i - \frac{1}{2}\sum_{i\le F}^A\sum_{j \le F}^A\left[\langle ij |\hat{v}|ij \rangle-\langle ij|\hat{v}|ji\rangle\right]. $$

This form will be used in our discussion of Koopman's theorem.

We have

$$ E[\Phi^{\mathrm{HF}}(N)] = \sum_{i=1}^H \langle i | \hat{h}_0 | i \rangle + \frac{1}{2}\sum_{ij=1}^N\langle ij|\hat{v}|ij\rangle_{AS}, $$

where $\Phi^{\mathrm{HF}}(N)$ is the new Slater determinant defined by the new basis of Eq. (21) for $N$ electrons (same $Z$). If we assume that the single-particle wave functions in the new basis do not change when we remove one electron or add one electron, we can then define the corresponding energy for the $N-1$ systems as

$$ E[\Phi^{\mathrm{HF}}(N-1)] = \sum_{i=1; i\ne k}^N \langle i | \hat{h}_0 | i \rangle + \frac{1}{2}\sum_{ij=1;i,j\ne k}^N\langle ij|\hat{v}|ij\rangle_{AS}, $$

where we have removed a single-particle state $k\le F$, that is a state below the Fermi level.

Calculating the difference

$$ E[\Phi^{\mathrm{HF}}(N)]- E[\Phi^{\mathrm{HF}}(N-1)] = \langle k | \hat{h}_0 | k \rangle + \frac{1}{2}\sum_{i=1;i\ne k}^N\langle ik|\hat{v}|ik\rangle_{AS} \frac{1}{2}\sum_{j=1;j\ne k}^N\langle kj|\hat{v}|kj\rangle_{AS}, $$

we obtain

$$ E[\Phi^{\mathrm{HF}}(N)]- E[\Phi^{\mathrm{HF}}(N-1)] = \langle k | \hat{h}_0 | k \rangle + \frac{1}{2}\sum_{j=1}^N\langle kj|\hat{v}|kj\rangle_{AS} $$

which is just our definition of the Hartree-Fock single-particle energy

$$ E[\Phi^{\mathrm{HF}}(N)]- E[\Phi^{\mathrm{HF}}(N-1)] = \epsilon_k^{\mathrm{HF}} $$

Similarly, we can now compute the difference (we label the single-particle states above the Fermi level as $abcd > F$)

$$ E[\Phi^{\mathrm{HF}}(N+1)]- E[\Phi^{\mathrm{HF}}(N)]= \epsilon_a^{\mathrm{HF}}. $$

These two equations can thus be used to the electron affinity or ionization energies, respectively. Koopman's theorem states that for example the ionization energy of a closed-shell system is given by the energy of the highest occupied single-particle state. If we assume that changing the number of electrons from $N$ to $N+1$ does not change the Hartree-Fock single-particle energies and eigenfunctions, then Koopman's theorem simply states that the ionization energy of an atom is given by the single-particle energy of the last bound state. In a similar way, we can also define the electron affinities.

As an example, consider a simple model for atomic sodium, Na. Neutral sodium has eleven electrons, with the weakest bound one being confined the $3s$ single-particle quantum numbers. The energy needed to remove an electron from neutral sodium is rather small, 5.1391 eV, a feature which pertains to all alkali metals. Having performed a Hartree-Fock calculation for neutral sodium would then allows us to compute the ionization energy by using the single-particle energy for the $3s$ states, namely $\epsilon_{3s}^{\mathrm{HF}}$.

From these considerations, we see that Hartree-Fock theory allows us to make a connection between experimental observables (here ionization and affinity energies) and the underlying interactions between particles.
In this sense, we are now linking the dynamics and structure of a many-body system with the laws of motion which govern the system. Our approach is a reductionistic one, meaning that we want to understand the laws of motion in terms of the particles or degrees of freedom which we believe are the fundamental ones. Our Slater determinant, being constructed as the product of various single-particle functions, follows this philosophy.

With similar arguments as in atomic physics, we can now use Hartree-Fock theory to make a link between nuclear forces and separation energies. Changing to nuclear system, we define

$$ E[\Phi^{\mathrm{HF}}(A)] = \sum_{i=1}^A \langle i | \hat{h}_0 | i \rangle + \frac{1}{2}\sum_{ij=1}^A\langle ij|\hat{v}|ij\rangle_{AS}, $$

where $\Phi^{\mathrm{HF}}(A)$ is the new Slater determinant defined by the new basis of Eq. (21) for $A$ nucleons, where $A=N+Z$, with $N$ now being the number of neutrons and $Z$ th enumber of protons. If we assume again that the single-particle wave functions in the new basis do not change from a nucleus with $A$ nucleons to a nucleus with $A-1$ nucleons, we can then define the corresponding energy for the $A-1$ systems as

$$ E[\Phi^{\mathrm{HF}}(A-1)] = \sum_{i=1; i\ne k}^A \langle i | \hat{h}_0 | i \rangle + \frac{1}{2}\sum_{ij=1;i,j\ne k}^A\langle ij|\hat{v}|ij\rangle_{AS}, $$

where we have removed a single-particle state $k\le F$, that is a state below the Fermi level.

Calculating the difference

$$ E[\Phi^{\mathrm{HF}}(A)]- E[\Phi^{\mathrm{HF}}(A-1)] = \langle k | \hat{h}_0 | k \rangle + \frac{1}{2}\sum_{i=1;i\ne k}^A\langle ik|\hat{v}|ik\rangle_{AS} \frac{1}{2}\sum_{j=1;j\ne k}^A\langle kj|\hat{v}|kj\rangle_{AS}, $$

which becomes

$$ E[\Phi^{\mathrm{HF}}(A)]- E[\Phi^{\mathrm{HF}}(A-1)] = \langle k | \hat{h}_0 | k \rangle + \frac{1}{2}\sum_{j=1}^A\langle kj|\hat{v}|kj\rangle_{AS} $$

which is just our definition of the Hartree-Fock single-particle energy

$$ E[\Phi^{\mathrm{HF}}(A)]- E[\Phi^{\mathrm{HF}}(A-1)] = \epsilon_k^{\mathrm{HF}} $$

Similarly, we can now compute the difference (recall that the single-particle states $abcd > F$)

$$ E[\Phi^{\mathrm{HF}}(A+1)]- E[\Phi^{\mathrm{HF}}(A)]= \epsilon_a^{\mathrm{HF}}. $$

If we then recall that the binding energy differences

$$ BE(A)-BE(A-1) \hspace{0.5cm} \mathrm{and} \hspace{0.5cm} BE(A+1)-BE(A), $$

define the separation energies, we see that the Hartree-Fock single-particle energies can be used to define separation energies. We have thus our first link between nuclear forces (included in the potential energy term) and an observable quantity defined by differences in binding energies.

We have thus the following interpretations (if the single-particle field do not change)

$$ BE(A)-BE(A-1)\approx E[\Phi^{\mathrm{HF}}(A)]- E[\Phi^{\mathrm{HF}}(A-1)] = \epsilon_k^{\mathrm{HF}}, $$

and

$$ BE(A+1)-BE(A)\approx E[\Phi^{\mathrm{HF}}(A+1)]- E[\Phi^{\mathrm{HF}}(A)] = \epsilon_a^{\mathrm{HF}}. $$

If we use ${}^{16}\mbox{O}$ as our closed-shell nucleus, we could then interpret the separation energy

$$ BE(^{16}\mathrm{O})-BE(^{15}\mathrm{O})\approx \epsilon_{0p^{\nu}_{1/2}}^{\mathrm{HF}}, $$

and

$$ BE(^{16}\mathrm{O})-BE(^{15}\mathrm{N})\approx \epsilon_{0p^{\pi}_{1/2}}^{\mathrm{HF}}. $$

Similalry, we could interpret

$$ BE(^{17}\mathrm{O})-BE(^{16}\mathrm{O})\approx \epsilon_{0d^{\nu}_{5/2}}^{\mathrm{HF}}, $$

and

$$ BE(^{17}\mathrm{F})-BE(^{16}\mathrm{O})\approx\epsilon_{0d^{\pi}_{5/2}}^{\mathrm{HF}}. $$

We can continue like this for all $A\pm 1$ nuclei where $A$ is a good closed-shell (or subshell closure) nucleus. Examples are ${}^{22}\mbox{O}$, ${}^{24}\mbox{O}$, ${}^{40}\mbox{Ca}$, ${}^{48}\mbox{Ca}$, ${}^{52}\mbox{Ca}$, ${}^{54}\mbox{Ca}$, ${}^{56}\mbox{Ni}$, ${}^{68}\mbox{Ni}$, ${}^{78}\mbox{Ni}$, ${}^{90}\mbox{Zr}$, ${}^{88}\mbox{Sr}$, ${}^{100}\mbox{Sn}$, ${}^{132}\mbox{Sn}$ and ${}^{208}\mbox{Pb}$, to mention some possile cases.

We can thus make our first interpretation of the separation energies in terms of the simplest possible many-body theory. If we also recall that the so-called energy gap for neutrons (or protons) is defined as

$$ \Delta S_n= 2BE(N,Z)-BE(N-1,Z)-BE(N+1,Z), $$

for neutrons and the corresponding gap for protons

$$ \Delta S_p= 2BE(N,Z)-BE(N,Z-1)-BE(N,Z+1), $$

we can define the neutron and proton energy gaps for ${}^{16}\mbox{O}$ as

$$ \Delta S_{\nu}=\epsilon_{0d^{\nu}_{5/2}}^{\mathrm{HF}}-\epsilon_{0p^{\nu}_{1/2}}^{\mathrm{HF}}, $$

and

$$ \Delta S_{\pi}=\epsilon_{0d^{\pi}_{5/2}}^{\mathrm{HF}}-\epsilon_{0p^{\pi}_{1/2}}^{\mathrm{HF}}. $$

Exercise 5: Derivation of Hartree-Fock equations

Consider a Slater determinant built up of single-particle orbitals $\psi_{\lambda}$, with $\lambda = 1,2,\dots,N$.

The unitary transformation

$$ \psi_a = \sum_{\lambda} C_{a\lambda}\phi_{\lambda}, $$

brings us into the new basis.
The new basis has quantum numbers $a=1,2,\dots,N$.

a) Show that the new basis is orthonormal.

b) Show that the new Slater determinant constructed from the new single-particle wave functions can be written as the determinant based on the previous basis and the determinant of the matrix $C$.

c) Show that the old and the new Slater determinants are equal up to a complex constant with absolute value unity.

Hint. Hint, $C$ is a unitary matrix.

Exercise 6: Derivation of Hartree-Fock equations

Consider the Slater determinant

$$ \Phi_{0}=\frac{1}{\sqrt{n!}}\sum_{p}(-)^{p}P \prod_{i=1}^{n}\psi_{\alpha_{i}}(x_{i}). $$

A small variation in this function is given by

$$ \delta\Phi_{0}=\frac{1}{\sqrt{n!}}\sum_{p}(-)^{p}P \psi_{\alpha_{1}}(x_{1})\psi_{\alpha_{2}}(x_{2})\dots \psi_{\alpha_{i-1}}(x_{i-1})(\delta\psi_{\alpha_{i}}(x_{i})) \psi_{\alpha_{i+1}}(x_{i+1})\dots\psi_{\alpha_{n}}(x_{n}). $$

a) Show that

$$ \langle \delta\Phi_{0}|\sum_{i=1}^{n}\left\{t(x_{i})+u(x_{i}) \right\}+\frac{1}{2} \sum_{i\neq j=1}^{n}v(x_{i},x_{j})|\Phi_{0}\rangle =\sum_{i=1}^{n}\langle \delta\psi_{\alpha_{i}}|\hat{t}+\hat{u} |\phi_{\alpha_{i}}\rangle $$

$$ +\sum_{i\neq j=1}^{n}\left\{\langle\delta\psi_{\alpha_{i}} \psi_{\alpha_{j}}|\hat{v}|\psi_{\alpha_{i}}\psi_{\alpha_{j}}\rangle- \langle\delta\psi_{\alpha_{i}}\psi_{\alpha_{j}}|\hat{v} |\psi_{\alpha_{j}}\psi_{\alpha_{i}}\rangle\right\} $$

Second quantization

Second quantization

• Second quantization and operators, two-body operator with examples

• Wick's theorem (missing in html and ipython notebook but in pdf files)

• Thouless' theorem and analysis of the Hartree-Fock equations using second quantization

• Examples on how to use bit representations for Slater determinants

• TO DO: add material about Wick's theorem and diagrammatic representation

• TO DO: Link determinantal analysis with Thouless' theorem

We introduce the time-independent operators $a_\alpha^{\dagger}$ and $a_\alpha$ which create and annihilate, respectively, a particle in the single-particle state $\varphi_\alpha$. We define the fermion creation operator $a_\alpha^{\dagger}$

$$ \begin{equation} a_\alpha^{\dagger}|0\rangle \equiv |\alpha\rangle \label{eq:2-1a} \tag{25}, \end{equation} $$

and

$$ \begin{equation} a_\alpha^{\dagger}|\alpha_1\dots \alpha_n\rangle_{\mathrm{AS}} \equiv |\alpha\alpha_1\dots \alpha_n\rangle_{\mathrm{AS}} \label{eq:2-1b} \tag{26} \end{equation} $$

In Eq. (25) the operator $a_\alpha^{\dagger}$ acts on the vacuum state $|0\rangle$, which does not contain any particles. Alternatively, we could define a closed-shell nucleus or atom as our new vacuum, but then we need to introduce the particle-hole formalism, see the discussion to come.

In Eq. (26) $a_\alpha^{\dagger}$ acts on an antisymmetric $n$-particle state and creates an antisymmetric $(n+1)$-particle state, where the one-body state $\varphi_\alpha$ is occupied, under the condition that $\alpha \ne \alpha_1, \alpha_2, \dots, \alpha_n$. It follows that we can express an antisymmetric state as the product of the creation operators acting on the vacuum state.

$$ \begin{equation} |\alpha_1\dots \alpha_n\rangle_{\mathrm{AS}} = a_{\alpha_1}^{\dagger} a_{\alpha_2}^{\dagger} \dots a_{\alpha_n}^{\dagger} |0\rangle \label{eq:2-2} \tag{27} \end{equation} $$

It is easy to derive the commutation and anticommutation rules for the fermionic creation operators $a_\alpha^{\dagger}$. Using the antisymmetry of the states (27)

$$ \begin{equation} |\alpha_1\dots \alpha_i\dots \alpha_k\dots \alpha_n\rangle_{\mathrm{AS}} = - |\alpha_1\dots \alpha_k\dots \alpha_i\dots \alpha_n\rangle_{\mathrm{AS}} \label{eq:2-3a} \tag{28} \end{equation} $$

we obtain

$$ \begin{equation} a_{\alpha_i}^{\dagger} a_{\alpha_k}^{\dagger} = - a_{\alpha_k}^{\dagger} a_{\alpha_i}^{\dagger} \label{eq:2-3b} \tag{29} \end{equation} $$

Using the Pauli principle

$$ \begin{equation} |\alpha_1\dots \alpha_i\dots \alpha_i\dots \alpha_n\rangle_{\mathrm{AS}} = 0 \label{eq:2-4a} \tag{30} \end{equation} $$

it follows that

$$ \begin{equation} a_{\alpha_i}^{\dagger} a_{\alpha_i}^{\dagger} = 0. \label{eq:2-4b} \tag{31} \end{equation} $$

If we combine Eqs. (29) and (31), we obtain the well-known anti-commutation rule

$$ \begin{equation} a_{\alpha}^{\dagger} a_{\beta}^{\dagger} + a_{\beta}^{\dagger} a_{\alpha}^{\dagger} \equiv \{a_{\alpha}^{\dagger},a_{\beta}^{\dagger}\} = 0 \label{eq:2-5} \tag{32} \end{equation} $$

The hermitian conjugate of $a_\alpha^{\dagger}$ is

$$ \begin{equation} a_{\alpha} = ( a_{\alpha}^{\dagger} )^{\dagger} \label{eq:2-6} \tag{33} \end{equation} $$

If we take the hermitian conjugate of Eq. (32), we arrive at

$$ \begin{equation} \{a_{\alpha},a_{\beta}\} = 0 \label{eq:2-7} \tag{34} \end{equation} $$

What is the physical interpretation of the operator $a_\alpha$ and what is the effect of $a_\alpha$ on a given state $|\alpha_1\alpha_2\dots\alpha_n\rangle_{\mathrm{AS}}$? Consider the following matrix element

$$ \begin{equation} \langle\alpha_1\alpha_2 \dots \alpha_n|a_\alpha|\alpha_1'\alpha_2' \dots \alpha_m'\rangle \label{eq:2-8} \tag{35} \end{equation} $$

where both sides are antisymmetric. We distinguish between two cases. The first (1) is when $\alpha \in \{\alpha_i\}$. Using the Pauli principle of Eq. (30) it follows

$$ \begin{equation} \langle\alpha_1\alpha_2 \dots \alpha_n|a_\alpha = 0 \label{eq:2-9a} \tag{36} \end{equation} $$

The second (2) case is when $\alpha \notin \{\alpha_i\}$. It follows that an hermitian conjugation

$$ \begin{equation} \langle \alpha_1\alpha_2 \dots \alpha_n|a_\alpha = \langle\alpha\alpha_1\alpha_2 \dots \alpha_n| \label{eq:2-9b} \tag{37} \end{equation} $$

Eq. (37) holds for case (1) since the lefthand side is zero due to the Pauli principle. We write Eq. (35) as

$$ \begin{equation} \langle\alpha_1\alpha_2 \dots \alpha_n|a_\alpha|\alpha_1'\alpha_2' \dots \alpha_m'\rangle = \langle \alpha_1\alpha_2 \dots \alpha_n|\alpha\alpha_1'\alpha_2' \dots \alpha_m'\rangle \label{eq:2-10} \tag{38} \end{equation} $$

Here we must have $m = n+1$ if Eq. (38) has to be trivially different from zero.

For the last case, the minus and plus signs apply when the sequence $\alpha ,\alpha_1, \alpha_2, \dots, \alpha_n$ and $\alpha_1', \alpha_2', \dots, \alpha_{n+1}'$ are related to each other via even and odd permutations. If we assume that $\alpha \notin \{\alpha_i\}$ we obtain

$$ \begin{equation} \langle\alpha_1\alpha_2 \dots \alpha_n|a_\alpha|\alpha_1'\alpha_2' \dots \alpha_{n+1}'\rangle = 0 \label{eq:2-12} \tag{39} \end{equation} $$

when $\alpha \in \{\alpha_i'\}$. If $\alpha \notin \{\alpha_i'\}$, we obtain

$$ \begin{equation} a_\alpha\underbrace{|\alpha_1'\alpha_2' \dots \alpha_{n+1}'}\rangle_{\neq \alpha} = 0 \label{eq:2-13a} \tag{40} \end{equation} $$

and in particular

$$ \begin{equation} a_\alpha |0\rangle = 0 \label{eq:2-13b} \tag{41} \end{equation} $$

If $\{\alpha\alpha_i\} = \{\alpha_i'\}$, performing the right permutations, the sequence $\alpha ,\alpha_1, \alpha_2, \dots, \alpha_n$ is identical with the sequence $\alpha_1', \alpha_2', \dots, \alpha_{n+1}'$. This results in

$$ \begin{equation} \langle\alpha_1\alpha_2 \dots \alpha_n|a_\alpha|\alpha\alpha_1\alpha_2 \dots \alpha_{n}\rangle = 1 \label{eq:2-14} \tag{42} \end{equation} $$

and thus

$$ \begin{equation} a_\alpha |\alpha\alpha_1\alpha_2 \dots \alpha_{n}\rangle = |\alpha_1\alpha_2 \dots \alpha_{n}\rangle \label{eq:2-15} \tag{43} \end{equation} $$

The action of the operator $a_\alpha$ from the left on a state vector is to to remove one particle in the state $\alpha$. If the state vector does not contain the single-particle state $\alpha$, the outcome of the operation is zero. The operator $a_\alpha$ is normally called for a destruction or annihilation operator.

The next step is to establish the commutator algebra of $a_\alpha^{\dagger}$ and $a_\beta$.

The action of the anti-commutator $\{a_\alpha^{\dagger}$,$a_\alpha\}$ on a given $n$-particle state is

$$ \begin{equation} a_\alpha^{\dagger} a_\alpha \underbrace{|\alpha_1\alpha_2 \dots \alpha_{n}\rangle}_{\neq \alpha} = 0 \nonumber \end{equation} $$

$$ \begin{equation} a_\alpha a_\alpha^{\dagger} \underbrace{|\alpha_1\alpha_2 \dots \alpha_{n}\rangle}_{\neq \alpha} = a_\alpha \underbrace{|\alpha \alpha_1\alpha_2 \dots \alpha_{n}\rangle}_{\neq \alpha} = \underbrace{|\alpha_1\alpha_2 \dots \alpha_{n}\rangle}_{\neq \alpha} \label{eq:2-16a} \tag{44} \end{equation} $$

if the single-particle state $\alpha$ is not contained in the state.

If it is present we arrive at

$$ \begin{equation} a_\alpha^{\dagger} a_\alpha |\alpha_1\alpha_2 \dots \alpha_{k}\alpha \alpha_{k+1} \dots \alpha_{n-1}\rangle = a_\alpha^{\dagger} a_\alpha (-1)^k |\alpha \alpha_1\alpha_2 \dots \alpha_{n-1}\rangle \nonumber \end{equation} $$

$$ \begin{equation} = (-1)^k |\alpha \alpha_1\alpha_2 \dots \alpha_{n-1}\rangle = |\alpha_1\alpha_2 \dots \alpha_{k}\alpha \alpha_{k+1} \dots \alpha_{n-1}\rangle \nonumber \end{equation} $$

$$ \begin{equation} a_\alpha a_\alpha^{\dagger}|\alpha_1\alpha_2 \dots \alpha_{k}\alpha \alpha_{k+1} \dots \alpha_{n-1}\rangle = 0 \label{eq:2-16b} \tag{45} \end{equation} $$

From Eqs. (44) and (45) we arrive at

$$ \begin{equation} \{a_\alpha^{\dagger} , a_\alpha \} = a_\alpha^{\dagger} a_\alpha + a_\alpha a_\alpha^{\dagger} = 1 \label{eq:2-17} \tag{46} \end{equation} $$

The action of $\left\{a_\alpha^{\dagger}, a_\beta\right\}$, with $\alpha \ne \beta$ on a given state yields three possibilities. The first case is a state vector which contains both $\alpha$ and $\beta$, then either $\alpha$ or $\beta$ and finally none of them.

The first case results in

$$ \begin{eqnarray} a_\alpha^{\dagger} a_\beta |\alpha\beta\alpha_1\alpha_2 \dots \alpha_{n-2}\rangle = 0 \nonumber \\ a_\beta a_\alpha^{\dagger} |\alpha\beta\alpha_1\alpha_2 \dots \alpha_{n-2}\rangle = 0 \label{eq:2-18a} \tag{47} \end{eqnarray} $$

while the second case gives

$$ \begin{eqnarray} a_\alpha^{\dagger} a_\beta |\beta \underbrace{\alpha_1\alpha_2 \dots \alpha_{n-1}}_{\neq \alpha}\rangle &=& |\alpha \underbrace{\alpha_1\alpha_2 \dots \alpha_{n-1}}_{\neq \alpha}\rangle \nonumber \\ a_\beta a_\alpha^{\dagger} |\beta \underbrace{\alpha_1\alpha_2 \dots \alpha_{n-1}}_{\neq \alpha}\rangle &=& a_\beta |\alpha\beta\underbrace{\beta \alpha_1\alpha_2 \dots \alpha_{n-1}}_{\neq \alpha}\rangle \nonumber \\ &=& - |\alpha\underbrace{\alpha_1\alpha_2 \dots \alpha_{n-1}}_{\neq \alpha}\rangle \label{eq:2-18b} \tag{48} \end{eqnarray} $$

Finally if the state vector does not contain $\alpha$ and $\beta$

$$ \begin{eqnarray} a_\alpha^{\dagger} a_\beta |\underbrace{\alpha_1\alpha_2 \dots \alpha_{n}}_{\neq \alpha,\beta}\rangle &=& 0 \nonumber \\ a_\beta a_\alpha^{\dagger} |\underbrace{\alpha_1\alpha_2 \dots \alpha_{n}}_{\neq \alpha,\beta}\rangle &=& a_\beta |\alpha \underbrace{\alpha_1\alpha_2 \dots \alpha_{n}}_{\neq \alpha,\beta}\rangle = 0 \label{eq:2-18c} \tag{49} \end{eqnarray} $$

For all three cases we have

$$ \begin{equation} \{a_\alpha^{\dagger},a_\beta \} = a_\alpha^{\dagger} a_\beta + a_\beta a_\alpha^{\dagger} = 0, \quad \alpha \neq \beta \label{eq:2-19} \tag{50} \end{equation} $$

We can summarize our findings in Eqs. (46) and (50) as

$$ \begin{equation} \{a_\alpha^{\dagger},a_\beta \} = \delta_{\alpha\beta} \label{eq:2-20} \tag{51} \end{equation} $$

with $\delta_{\alpha\beta}$ is the Kroenecker $\delta$-symbol.

The properties of the creation and annihilation operators can be summarized as (for fermions)

$$ a_\alpha^{\dagger}|0\rangle \equiv |\alpha\rangle, $$

and

$$ a_\alpha^{\dagger}|\alpha_1\dots \alpha_n\rangle_{\mathrm{AS}} \equiv |\alpha\alpha_1\dots \alpha_n\rangle_{\mathrm{AS}}. $$

from which follows

$$ |\alpha_1\dots \alpha_n\rangle_{\mathrm{AS}} = a_{\alpha_1}^{\dagger} a_{\alpha_2}^{\dagger} \dots a_{\alpha_n}^{\dagger} |0\rangle. $$

The hermitian conjugate has the folowing properties

$$ a_{\alpha} = ( a_{\alpha}^{\dagger} )^{\dagger}. $$

Finally we found

$$ a_\alpha\underbrace{|\alpha_1'\alpha_2' \dots \alpha_{n+1}'}\rangle_{\neq \alpha} = 0, \quad \textrm{in particular } a_\alpha |0\rangle = 0, $$

and

$$ a_\alpha |\alpha\alpha_1\alpha_2 \dots \alpha_{n}\rangle = |\alpha_1\alpha_2 \dots \alpha_{n}\rangle, $$

and the corresponding commutator algebra

$$ \{a_{\alpha}^{\dagger},a_{\beta}^{\dagger}\} = \{a_{\alpha},a_{\beta}\} = 0 \hspace{0.5cm} \{a_\alpha^{\dagger},a_\beta \} = \delta_{\alpha\beta}. $$

One-body operators in second quantization

A very useful operator is the so-called number-operator. Most physics cases we will study in this text conserve the total number of particles. The number operator is therefore a useful quantity which allows us to test that our many-body formalism conserves the number of particles. In for example $(d,p)$ or $(p,d)$ reactions it is important to be able to describe quantum mechanical states where particles get added or removed. A creation operator $a_\alpha^{\dagger}$ adds one particle to the single-particle state $\alpha$ of a give many-body state vector, while an annihilation operator $a_\alpha$ removes a particle from a single-particle state $\alpha$.

Let us consider an operator proportional with $a_\alpha^{\dagger} a_\beta$ and $\alpha=\beta$. It acts on an $n$-particle state resulting in

$$ \begin{equation} a_\alpha^{\dagger} a_\alpha |\alpha_1\alpha_2 \dots \alpha_{n}\rangle = \begin{cases} 0 &\alpha \notin \{\alpha_i\} \\ \\ |\alpha_1\alpha_2 \dots \alpha_{n}\rangle & \alpha \in \{\alpha_i\} \end{cases} \end{equation} $$

Summing over all possible one-particle states we arrive at

$$ \begin{equation} \left( \sum_\alpha a_\alpha^{\dagger} a_\alpha \right) |\alpha_1\alpha_2 \dots \alpha_{n}\rangle = n |\alpha_1\alpha_2 \dots \alpha_{n}\rangle \label{eq:2-21} \tag{52} \end{equation} $$

The operator

$$ \begin{equation} \hat{N} = \sum_\alpha a_\alpha^{\dagger} a_\alpha \label{eq:2-22} \tag{53} \end{equation} $$

is called the number operator since it counts the number of particles in a give state vector when it acts on the different single-particle states. It acts on one single-particle state at the time and falls therefore under category one-body operators. Next we look at another important one-body operator, namely $\hat{H}_0$ and study its operator form in the occupation number representation.

We want to obtain an expression for a one-body operator which conserves the number of particles. Here we study the one-body operator for the kinetic energy plus an eventual external one-body potential. The action of this operator on a particular $n$-body state with its pertinent expectation value has already been studied in coordinate space. In coordinate space the operator reads

$$ \begin{equation} \hat{H}_0 = \sum_i \hat{h}_0(x_i) \label{eq:2-23} \tag{54} \end{equation} $$

and the anti-symmetric $n$-particle Slater determinant is defined as

$$ \Phi(x_1, x_2,\dots ,x_n,\alpha_1,\alpha_2,\dots, \alpha_n)= \frac{1}{\sqrt{n!}} \sum_p (-1)^p\hat{P}\psi_{\alpha_1}(x_1)\psi_{\alpha_2}(x_2) \dots \psi_{\alpha_n}(x_n). $$

Defining

$$ \begin{equation} \hat{h}_0(x_i) \psi_{\alpha_i}(x_i) = \sum_{\alpha_k'} \psi_{\alpha_k'}(x_i) \langle\alpha_k'|\hat{h}_0|\alpha_k\rangle \label{eq:2-25} \tag{55} \end{equation} $$

we can easily evaluate the action of $\hat{H}_0$ on each product of one-particle functions in Slater determinant. From Eq. (55) we obtain the following result without permuting any particle pair

$$ \begin{eqnarray} && \left( \sum_i \hat{h}_0(x_i) \right) \psi_{\alpha_1}(x_1)\psi_{\alpha_2}(x_2) \dots \psi_{\alpha_n}(x_n) \nonumber \\ & =&\sum_{\alpha_1'} \langle \alpha_1'|\hat{h}_0|\alpha_1\rangle \psi_{\alpha_1'}(x_1)\psi_{\alpha_2}(x_2) \dots \psi_{\alpha_n}(x_n) \nonumber \\ &+&\sum_{\alpha_2'} \langle \alpha_2'|\hat{h}_0|\alpha_2\rangle \psi_{\alpha_1}(x_1)\psi_{\alpha_2'}(x_2) \dots \psi_{\alpha_n}(x_n) \nonumber \\ &+& \dots \nonumber \\ &+&\sum_{\alpha_n'} \langle \alpha_n'|\hat{h}_0|\alpha_n\rangle \psi_{\alpha_1}(x_1)\psi_{\alpha_2}(x_2) \dots \psi_{\alpha_n'}(x_n) \label{eq:2-26} \tag{56} \end{eqnarray} $$

If we interchange particles $1$ and $2$ we obtain

$$ \begin{eqnarray} && \left( \sum_i \hat{h}_0(x_i) \right) \psi_{\alpha_1}(x_2)\psi_{\alpha_1}(x_2) \dots \psi_{\alpha_n}(x_n) \nonumber \\ & =&\sum_{\alpha_2'} \langle \alpha_2'|\hat{h}_0|\alpha_2\rangle \psi_{\alpha_1}(x_2)\psi_{\alpha_2'}(x_1) \dots \psi_{\alpha_n}(x_n) \nonumber \\ &+&\sum_{\alpha_1'} \langle \alpha_1'|\hat{h}_0|\alpha_1\rangle \psi_{\alpha_1'}(x_2)\psi_{\alpha_2}(x_1) \dots \psi_{\alpha_n}(x_n) \nonumber \\ &+& \dots \nonumber \\ &+&\sum_{\alpha_n'} \langle \alpha_n'|\hat{h}_0|\alpha_n\rangle \psi_{\alpha_1}(x_2)\psi_{\alpha_1}(x_2) \dots \psi_{\alpha_n'}(x_n) \label{eq:2-27} \tag{57} \end{eqnarray} $$

We can continue by computing all possible permutations. We rewrite also our Slater determinant in its second quantized form and skip the dependence on the quantum numbers $x_i.$ Summing up all contributions and taking care of all phases $(-1)^p$ we arrive at

$$ \begin{eqnarray} \hat{H}_0|\alpha_1,\alpha_2,\dots, \alpha_n\rangle &=& \sum_{\alpha_1'}\langle \alpha_1'|\hat{h}_0|\alpha_1\rangle |\alpha_1'\alpha_2 \dots \alpha_{n}\rangle \nonumber \\ &+& \sum_{\alpha_2'} \langle \alpha_2'|\hat{h}_0|\alpha_2\rangle |\alpha_1\alpha_2' \dots \alpha_{n}\rangle \nonumber \\ &+& \dots \nonumber \\ &+& \sum_{\alpha_n'} \langle \alpha_n'|\hat{h}_0|\alpha_n\rangle |\alpha_1\alpha_2 \dots \alpha_{n}'\rangle \label{eq:2-28} \tag{58} \end{eqnarray} $$

In Eq. (58) we have expressed the action of the one-body operator of Eq. (54) on the $n$-body state in its second quantized form. This equation can be further manipulated if we use the properties of the creation and annihilation operator on each primed quantum number, that is

$$ \begin{equation} |\alpha_1\alpha_2 \dots \alpha_k' \dots \alpha_{n}\rangle = a_{\alpha_k'}^{\dagger} a_{\alpha_k} |\alpha_1\alpha_2 \dots \alpha_k \dots \alpha_{n}\rangle \label{eq:2-29} \tag{59} \end{equation} $$

Inserting this in the right-hand side of Eq. (58) results in

$$ \begin{eqnarray} \hat{H}_0|\alpha_1\alpha_2 \dots \alpha_{n}\rangle &=& \sum_{\alpha_1'}\langle \alpha_1'|\hat{h}_0|\alpha_1\rangle a_{\alpha_1'}^{\dagger} a_{\alpha_1} |\alpha_1\alpha_2 \dots \alpha_{n}\rangle \nonumber \\ &+& \sum_{\alpha_2'} \langle \alpha_2'|\hat{h}_0|\alpha_2\rangle a_{\alpha_2'}^{\dagger} a_{\alpha_2} |\alpha_1\alpha_2 \dots \alpha_{n}\rangle \nonumber \\ &+& \dots \nonumber \\ &+& \sum_{\alpha_n'} \langle \alpha_n'|\hat{h}_0|\alpha_n\rangle a_{\alpha_n'}^{\dagger} a_{\alpha_n} |\alpha_1\alpha_2 \dots \alpha_{n}\rangle \nonumber \\ &=& \sum_{\alpha, \beta} \langle \alpha|\hat{h}_0|\beta\rangle a_\alpha^{\dagger} a_\beta |\alpha_1\alpha_2 \dots \alpha_{n}\rangle \label{eq:2-30a} \tag{60} \end{eqnarray} $$

In the number occupation representation or second quantization we get the following expression for a one-body operator which conserves the number of particles

$$ \begin{equation} \hat{H}_0 = \sum_{\alpha\beta} \langle \alpha|\hat{h}_0|\beta\rangle a_\alpha^{\dagger} a_\beta \label{eq:2-30b} \tag{61} \end{equation} $$

Obviously, $\hat{H}_0$ can be replaced by any other one-body operator which preserved the number of particles. The stucture of the operator is therefore not limited to say the kinetic or single-particle energy only.

The opearator $\hat{H}_0$ takes a particle from the single-particle state $\beta$ to the single-particle state $\alpha$ with a probability for the transition given by the expectation value $\langle \alpha|\hat{h}_0|\beta\rangle$.

It is instructive to verify Eq. (61) by computing the expectation value of $\hat{H}_0$ between two single-particle states

$$ \begin{equation} \langle \alpha_1|\hat{h}_0|\alpha_2\rangle = \sum_{\alpha\beta} \langle \alpha|\hat{h}_0|\beta\rangle \langle 0|a_{\alpha_1}a_\alpha^{\dagger} a_\beta a_{\alpha_2}^{\dagger}|0\rangle \label{eq:2-30c} \tag{62} \end{equation} $$

Using the commutation relations for the creation and annihilation operators we have

$$ \begin{equation} a_{\alpha_1}a_\alpha^{\dagger} a_\beta a_{\alpha_2}^{\dagger} = (\delta_{\alpha \alpha_1} - a_\alpha^{\dagger} a_{\alpha_1} )(\delta_{\beta \alpha_2} - a_{\alpha_2}^{\dagger} a_{\beta} ), \label{eq:2-30d} \tag{63} \end{equation} $$

which results in

$$ \begin{equation} \langle 0|a_{\alpha_1}a_\alpha^{\dagger} a_\beta a_{\alpha_2}^{\dagger}|0\rangle = \delta_{\alpha \alpha_1} \delta_{\beta \alpha_2} \label{eq:2-30e} \tag{64} \end{equation} $$

and

$$ \begin{equation} \langle \alpha_1|\hat{h}_0|\alpha_2\rangle = \sum_{\alpha\beta} \langle \alpha|\hat{h}_0|\beta\rangle\delta_{\alpha \alpha_1} \delta_{\beta \alpha_2} = \langle \alpha_1|\hat{h}_0|\alpha_2\rangle \label{eq:2-30f} \tag{65} \end{equation} $$

Two-body operators in second quantization

Let us now derive the expression for our two-body interaction part, which also conserves the number of particles. We can proceed in exactly the same way as for the one-body operator. In the coordinate representation our two-body interaction part takes the following expression

$$ \begin{equation} \hat{H}_I = \sum_{i < j} V(x_i,x_j) \label{eq:2-31} \tag{66} \end{equation} $$

where the summation runs over distinct pairs. The term $V$ can be an interaction model for the nucleon-nucleon interaction or the interaction between two electrons. It can also include additional two-body interaction terms.

The action of this operator on a product of two single-particle functions is defined as

$$ \begin{equation} V(x_i,x_j) \psi_{\alpha_k}(x_i) \psi_{\alpha_l}(x_j) = \sum_{\alpha_k'\alpha_l'} \psi_{\alpha_k}'(x_i)\psi_{\alpha_l}'(x_j) \langle \alpha_k'\alpha_l'|\hat{v}|\alpha_k\alpha_l\rangle \label{eq:2-32} \tag{67} \end{equation} $$

We can now let $\hat{H}_I$ act on all terms in the linear combination for $|\alpha_1\alpha_2\dots\alpha_n\rangle$. Without any permutations we have

$$ \begin{eqnarray} && \left( \sum_{i < j} V(x_i,x_j) \right) \psi_{\alpha_1}(x_1)\psi_{\alpha_2}(x_2)\dots \psi_{\alpha_n}(x_n) \nonumber \\ &=& \sum_{\alpha_1'\alpha_2'} \langle \alpha_1'\alpha_2'|\hat{v}|\alpha_1\alpha_2\rangle \psi_{\alpha_1}'(x_1)\psi_{\alpha_2}'(x_2)\dots \psi_{\alpha_n}(x_n) \nonumber \\ & +& \dots \nonumber \\ &+& \sum_{\alpha_1'\alpha_n'} \langle \alpha_1'\alpha_n'|\hat{v}|\alpha_1\alpha_n\rangle \psi_{\alpha_1}'(x_1)\psi_{\alpha_2}(x_2)\dots \psi_{\alpha_n}'(x_n) \nonumber \\ & +& \dots \nonumber \\ &+& \sum_{\alpha_2'\alpha_n'} \langle \alpha_2'\alpha_n'|\hat{v}|\alpha_2\alpha_n\rangle \psi_{\alpha_1}(x_1)\psi_{\alpha_2}'(x_2)\dots \psi_{\alpha_n}'(x_n) \nonumber \\ & +& \dots \label{eq:2-33} \tag{68} \end{eqnarray} $$

where on the rhs we have a term for each distinct pairs.

For the other terms on the rhs we obtain similar expressions and summing over all terms we obtain

$$ \begin{eqnarray} H_I |\alpha_1\alpha_2\dots\alpha_n\rangle &=& \sum_{\alpha_1', \alpha_2'} \langle \alpha_1'\alpha_2'|\hat{v}|\alpha_1\alpha_2\rangle |\alpha_1'\alpha_2'\dots\alpha_n\rangle \nonumber \\ &+& \dots \nonumber \\ &+& \sum_{\alpha_1', \alpha_n'} \langle \alpha_1'\alpha_n'|\hat{v}|\alpha_1\alpha_n\rangle |\alpha_1'\alpha_2\dots\alpha_n'\rangle \nonumber \\ &+& \dots \nonumber \\ &+& \sum_{\alpha_2', \alpha_n'} \langle \alpha_2'\alpha_n'|\hat{v}|\alpha_2\alpha_n\rangle |\alpha_1\alpha_2'\dots\alpha_n'\rangle \nonumber \\ &+& \dots \label{eq:2-34} \tag{69} \end{eqnarray} $$

We introduce second quantization via the relation

$$ \begin{eqnarray} && a_{\alpha_k'}^{\dagger} a_{\alpha_l'}^{\dagger} a_{\alpha_l} a_{\alpha_k} |\alpha_1\alpha_2\dots\alpha_k\dots\alpha_l\dots\alpha_n\rangle \nonumber \\ &=& (-1)^{k-1} (-1)^{l-2} a_{\alpha_k'}^{\dagger} a_{\alpha_l'}^{\dagger} a_{\alpha_l} a_{\alpha_k} |\alpha_k\alpha_l \underbrace{\alpha_1\alpha_2\dots\alpha_n}_{\neq \alpha_k,\alpha_l}\rangle \nonumber \\ &=& (-1)^{k-1} (-1)^{l-2} |\alpha_k'\alpha_l' \underbrace{\alpha_1\alpha_2\dots\alpha_n}_{\neq \alpha_k',\alpha_l'}\rangle \nonumber \\ &=& |\alpha_1\alpha_2\dots\alpha_k'\dots\alpha_l'\dots\alpha_n\rangle \label{eq:2-35} \tag{70} \end{eqnarray} $$

Inserting this in (69) gives

$$ \begin{eqnarray} H_I |\alpha_1\alpha_2\dots\alpha_n\rangle &=& \sum_{\alpha_1', \alpha_2'} \langle \alpha_1'\alpha_2'|\hat{v}|\alpha_1\alpha_2\rangle a_{\alpha_1'}^{\dagger} a_{\alpha_2'}^{\dagger} a_{\alpha_2} a_{\alpha_1} |\alpha_1\alpha_2\dots\alpha_n\rangle \nonumber \\ &+& \dots \nonumber \\ &=& \sum_{\alpha_1', \alpha_n'} \langle \alpha_1'\alpha_n'|\hat{v}|\alpha_1\alpha_n\rangle a_{\alpha_1'}^{\dagger} a_{\alpha_n'}^{\dagger} a_{\alpha_n} a_{\alpha_1} |\alpha_1\alpha_2\dots\alpha_n\rangle \nonumber \\ &+& \dots \nonumber \\ &=& \sum_{\alpha_2', \alpha_n'} \langle \alpha_2'\alpha_n'|\hat{v}|\alpha_2\alpha_n\rangle a_{\alpha_2'}^{\dagger} a_{\alpha_n'}^{\dagger} a_{\alpha_n} a_{\alpha_2} |\alpha_1\alpha_2\dots\alpha_n\rangle \nonumber \\ &+& \dots \nonumber \\ &=& \sum_{\alpha, \beta, \gamma, \delta} ' \langle \alpha\beta|\hat{v}|\gamma\delta\rangle a^{\dagger}_\alpha a^{\dagger}_\beta a_\delta a_\gamma |\alpha_1\alpha_2\dots\alpha_n\rangle \label{eq:2-36} \tag{71} \end{eqnarray} $$

Here we let $\sum'$ indicate that the sums running over $\alpha$ and $\beta$ run over all single-particle states, while the summations $\gamma$ and $\delta$ run over all pairs of single-particle states. We wish to remove this restriction and since

$$ \begin{equation} \langle \alpha\beta|\hat{v}|\gamma\delta\rangle = \langle \beta\alpha|\hat{v}|\delta\gamma\rangle \label{eq:2-37} \tag{72} \end{equation} $$

we get

$$ \begin{eqnarray} \sum_{\alpha\beta} \langle \alpha\beta|\hat{v}|\gamma\delta\rangle a^{\dagger}_\alpha a^{\dagger}_\beta a_\delta a_\gamma &=& \sum_{\alpha\beta} \langle \beta\alpha|\hat{v}|\delta\gamma\rangle a^{\dagger}_\alpha a^{\dagger}_\beta a_\delta a_\gamma \label{eq:2-38a} \tag{73} \\ &=& \sum_{\alpha\beta}\langle \beta\alpha|\hat{v}|\delta\gamma\rangle a^{\dagger}_\beta a^{\dagger}_\alpha a_\gamma a_\delta \label{eq:2-38b} \tag{74} \end{eqnarray} $$

where we have used the anti-commutation rules.

Changing the summation indices $\alpha$ and $\beta$ in (74) we obtain

$$ \begin{equation} \sum_{\alpha\beta} \langle \alpha\beta|\hat{v}|\gamma\delta\rangle a^{\dagger}_\alpha a^{\dagger}_\beta a_\delta a_\gamma = \sum_{\alpha\beta} \langle \alpha\beta|\hat{v}|\delta\gamma\rangle a^{\dagger}_\alpha a^{\dagger}_\beta a_\gamma a_\delta \label{eq:2-38c} \tag{75} \end{equation} $$

From this it follows that the restriction on the summation over $\gamma$ and $\delta$ can be removed if we multiply with a factor $\frac{1}{2}$, resulting in

$$ \begin{equation} \hat{H}_I = \frac{1}{2} \sum_{\alpha\beta\gamma\delta} \langle \alpha\beta|\hat{v}|\gamma\delta\rangle a^{\dagger}_\alpha a^{\dagger}_\beta a_\delta a_\gamma \label{eq:2-39} \tag{76} \end{equation} $$

where we sum freely over all single-particle states $\alpha$, $\beta$, $\gamma$ og $\delta$.

With this expression we can now verify that the second quantization form of $\hat{H}_I$ in Eq. (76) results in the same matrix between two anti-symmetrized two-particle states as its corresponding coordinate space representation. We have

$$ \begin{equation} \langle \alpha_1 \alpha_2|\hat{H}_I|\beta_1 \beta_2\rangle = \frac{1}{2} \sum_{\alpha\beta\gamma\delta} \langle \alpha\beta|\hat{v}|\gamma\delta\rangle \langle 0|a_{\alpha_2} a_{\alpha_1} a^{\dagger}_\alpha a^{\dagger}_\beta a_\delta a_\gamma a_{\beta_1}^{\dagger} a_{\beta_2}^{\dagger}|0\rangle. \label{eq:2-40} \tag{77} \end{equation} $$

Using the commutation relations we get

$$ \begin{eqnarray} && a_{\alpha_2} a_{\alpha_1}a^{\dagger}_\alpha a^{\dagger}_\beta a_\delta a_\gamma a_{\beta_1}^{\dagger} a_{\beta_2}^{\dagger} \nonumber \\ &=& a_{\alpha_2} a_{\alpha_1}a^{\dagger}_\alpha a^{\dagger}_\beta ( a_\delta \delta_{\gamma \beta_1} a_{\beta_2}^{\dagger} - a_\delta a_{\beta_1}^{\dagger} a_\gamma a_{\beta_2}^{\dagger} ) \nonumber \\ &=& a_{\alpha_2} a_{\alpha_1}a^{\dagger}_\alpha a^{\dagger}_\beta (\delta_{\gamma \beta_1} \delta_{\delta \beta_2} - \delta_{\gamma \beta_1} a_{\beta_2}^{\dagger} a_\delta - a_\delta a_{\beta_1}^{\dagger}\delta_{\gamma \beta_2} + a_\delta a_{\beta_1}^{\dagger} a_{\beta_2}^{\dagger} a_\gamma ) \nonumber \\ &=& a_{\alpha_2} a_{\alpha_1}a^{\dagger}_\alpha a^{\dagger}_\beta (\delta_{\gamma \beta_1} \delta_{\delta \beta_2} - \delta_{\gamma \beta_1} a_{\beta_2}^{\dagger} a_\delta \nonumber \\ && \qquad - \delta_{\delta \beta_1} \delta_{\gamma \beta_2} + \delta_{\gamma \beta_2} a_{\beta_1}^{\dagger} a_\delta + a_\delta a_{\beta_1}^{\dagger} a_{\beta_2}^{\dagger} a_\gamma ) \label{eq:2-41} \tag{78} \end{eqnarray} $$

The vacuum expectation value of this product of operators becomes

$$ \begin{eqnarray} && \langle 0|a_{\alpha_2} a_{\alpha_1} a^{\dagger}_\alpha a^{\dagger}_\beta a_\delta a_\gamma a_{\beta_1}^{\dagger} a_{\beta_2}^{\dagger}|0\rangle \nonumber \\ &=& (\delta_{\gamma \beta_1} \delta_{\delta \beta_2} - \delta_{\delta \beta_1} \delta_{\gamma \beta_2} ) \langle 0|a_{\alpha_2} a_{\alpha_1}a^{\dagger}_\alpha a^{\dagger}_\beta|0\rangle \nonumber \\ &=& (\delta_{\gamma \beta_1} \delta_{\delta \beta_2} -\delta_{\delta \beta_1} \delta_{\gamma \beta_2} ) (\delta_{\alpha \alpha_1} \delta_{\beta \alpha_2} -\delta_{\beta \alpha_1} \delta_{\alpha \alpha_2} ) \label{eq:2-42b} \tag{79} \end{eqnarray} $$

Insertion of Eq. (79) in Eq. (77) results in

$$ \begin{eqnarray} \langle \alpha_1\alpha_2|\hat{H}_I|\beta_1\beta_2\rangle &=& \frac{1}{2} \big[ \langle \alpha_1\alpha_2|\hat{v}|\beta_1\beta_2\rangle- \langle \alpha_1\alpha_2|\hat{v}|\beta_2\beta_1\rangle \nonumber \\ && - \langle \alpha_2\alpha_1|\hat{v}|\beta_1\beta_2\rangle + \langle \alpha_2\alpha_1|\hat{v}|\beta_2\beta_1\rangle \big] \nonumber \\ &=& \langle \alpha_1\alpha_2|\hat{v}|\beta_1\beta_2\rangle - \langle \alpha_1\alpha_2|\hat{v}|\beta_2\beta_1\rangle \nonumber \\ &=& \langle \alpha_1\alpha_2|\hat{v}|\beta_1\beta_2\rangle_{\mathrm{AS}}. \label{eq:2-43b} \tag{80} \end{eqnarray} $$

The two-body operator can also be expressed in terms of the anti-symmetrized matrix elements we discussed previously as

$$ \begin{eqnarray} \hat{H}_I &=& \frac{1}{2} \sum_{\alpha\beta\gamma\delta} \langle \alpha \beta|\hat{v}|\gamma \delta\rangle a_\alpha^{\dagger} a_\beta^{\dagger} a_\delta a_\gamma \nonumber \\ &=& \frac{1}{4} \sum_{\alpha\beta\gamma\delta} \left[ \langle \alpha \beta|\hat{v}|\gamma \delta\rangle - \langle \alpha \beta|\hat{v}|\delta\gamma \rangle \right] a_\alpha^{\dagger} a_\beta^{\dagger} a_\delta a_\gamma \nonumber \\ &=& \frac{1}{4} \sum_{\alpha\beta\gamma\delta} \langle \alpha \beta|\hat{v}|\gamma \delta\rangle_{\mathrm{AS}} a_\alpha^{\dagger} a_\beta^{\dagger} a_\delta a_\gamma \label{eq:2-45} \tag{81} \end{eqnarray} $$

The factors in front of the operator, either $\frac{1}{4}$ or $\frac{1}{2}$ tells whether we use antisymmetrized matrix elements or not.

We can now express the Hamiltonian operator for a many-fermion system in the occupation basis representation as

$$ \begin{equation} H = \sum_{\alpha, \beta} \langle \alpha|\hat{t}+\hat{u}_{\mathrm{ext}}|\beta\rangle a_\alpha^{\dagger} a_\beta + \frac{1}{4} \sum_{\alpha\beta\gamma\delta} \langle \alpha \beta|\hat{v}|\gamma \delta\rangle a_\alpha^{\dagger} a_\beta^{\dagger} a_\delta a_\gamma. \label{eq:2-46b} \tag{82} \end{equation} $$

This is the form we will use in the rest of these lectures, assuming that we work with anti-symmetrized two-body matrix elements.

Particle-hole formalism

Second quantization is a useful and elegant formalism for constructing many-body states and quantum mechanical operators. One can express and translate many physical processes into simple pictures such as Feynman diagrams. Expecation values of many-body states are also easily calculated. However, although the equations are seemingly easy to set up, from a practical point of view, that is the solution of Schroedinger's equation, there is no particular gain. The many-body equation is equally hard to solve, irrespective of representation. The cliche that there is no free lunch brings us down to earth again.
Note however that a transformation to a particular basis, for cases where the interaction obeys specific symmetries, can ease the solution of Schroedinger's equation.

But there is at least one important case where second quantization comes to our rescue. It is namely easy to introduce another reference state than the pure vacuum $|0\rangle $, where all single-particle states are active. With many particles present it is often useful to introduce another reference state than the vacuum state$|0\rangle $. We will label this state $|c\rangle$ ($c$ for core) and as we will see it can reduce considerably the complexity and thereby the dimensionality of the many-body problem. It allows us to sum up to infinite order specific many-body correlations. The particle-hole representation is one of these handy representations.

In the original particle representation these states are products of the creation operators $a_{\alpha_i}^\dagger$ acting on the true vacuum $|0\rangle $. Following Eq. (27) we have

$$ \begin{eqnarray} |\alpha_1\alpha_2\dots\alpha_{n-1}\alpha_n\rangle &=& a_{\alpha_1}^\dagger a_{\alpha_2}^\dagger \dots a_{\alpha_{n-1}}^\dagger a_{\alpha_n}^\dagger |0\rangle \label{eq:2-47a} \tag{83} \\ |\alpha_1\alpha_2\dots\alpha_{n-1}\alpha_n\alpha_{n+1}\rangle &=& a_{\alpha_1}^\dagger a_{\alpha_2}^\dagger \dots a_{\alpha_{n-1}}^\dagger a_{\alpha_n}^\dagger a_{\alpha_{n+1}}^\dagger |0\rangle \label{eq:2-47b} \tag{84} \\ |\alpha_1\alpha_2\dots\alpha_{n-1}\rangle &=& a_{\alpha_1}^\dagger a_{\alpha_2}^\dagger \dots a_{\alpha_{n-1}}^\dagger |0\rangle \label{eq:2-47c} \tag{85} \end{eqnarray} $$

If we use Eq. (83) as our new reference state, we can simplify considerably the representation of this state

$$ \begin{equation} |c\rangle \equiv |\alpha_1\alpha_2\dots\alpha_{n-1}\alpha_n\rangle = a_{\alpha_1}^\dagger a_{\alpha_2}^\dagger \dots a_{\alpha_{n-1}}^\dagger a_{\alpha_n}^\dagger |0\rangle \label{eq:2-48a} \tag{86} \end{equation} $$

The new reference states for the $n+1$ and $n-1$ states can then be written as

$$ \begin{eqnarray} |\alpha_1\alpha_2\dots\alpha_{n-1}\alpha_n\alpha_{n+1}\rangle &=& (-1)^n a_{\alpha_{n+1}}^\dagger |c\rangle \equiv (-1)^n |\alpha_{n+1}\rangle_c \label{eq:2-48b} \tag{87} \\ |\alpha_1\alpha_2\dots\alpha_{n-1}\rangle &=& (-1)^{n-1} a_{\alpha_n} |c\rangle \equiv (-1)^{n-1} |\alpha_{n-1}\rangle_c \label{eq:2-48c} \tag{88} \end{eqnarray} $$

The first state has one additional particle with respect to the new vacuum state $|c\rangle $ and is normally referred to as a one-particle state or one particle added to the many-body reference state. The second state has one particle less than the reference vacuum state $|c\rangle $ and is referred to as a one-hole state. When dealing with a new reference state it is often convenient to introduce new creation and annihilation operators since we have from Eq. (88)

$$ \begin{equation} a_\alpha |c\rangle \neq 0 \label{eq:2-49} \tag{89} \end{equation} $$

since $\alpha$ is contained in $|c\rangle $, while for the true vacuum we have $a_\alpha |0\rangle = 0$ for all $\alpha$.

The new reference state leads to the definition of new creation and annihilation operators which satisfy the following relations

$$ \begin{eqnarray} b_\alpha |c\rangle &=& 0 \label{eq:2-50a} \tag{90} \\ \{b_\alpha^\dagger , b_\beta^\dagger \} = \{b_\alpha , b_\beta \} &=& 0 \nonumber \\ \{b_\alpha^\dagger , b_\beta \} &=& \delta_{\alpha \beta} \label{eq:2-50c} \tag{91} \end{eqnarray} $$

We assume also that the new reference state is properly normalized

$$ \begin{equation} \langle c | c \rangle = 1 \label{eq:2-51} \tag{92} \end{equation} $$

The physical interpretation of these new operators is that of so-called quasiparticle states. This means that a state defined by the addition of one extra particle to a reference state $|c\rangle $ may not necesseraly be interpreted as one particle coupled to a core. We define now new creation operators that act on a state $\alpha$ creating a new quasiparticle state

$$ \begin{equation} b_\alpha^\dagger|c\rangle = \Bigg\{ \begin{array}{ll} a_\alpha^\dagger |c\rangle = |\alpha\rangle, & \alpha > F \\ \\ a_\alpha |c\rangle = |\alpha^{-1}\rangle, & \alpha \leq F \end{array} \label{eq:2-52} \tag{93} \end{equation} $$

where $F$ is the Fermi level representing the last occupied single-particle orbit of the new reference state $|c\rangle $.

The annihilation is the hermitian conjugate of the creation operator

$$ b_\alpha = (b_\alpha^\dagger)^\dagger, $$

resulting in

$$ \begin{equation} b_\alpha^\dagger = \Bigg\{ \begin{array}{ll} a_\alpha^\dagger & \alpha > F \\ \\ a_\alpha & \alpha \leq F \end{array} \qquad b_\alpha = \Bigg\{ \begin{array}{ll} a_\alpha & \alpha > F \\ \\ a_\alpha^\dagger & \alpha \leq F \end{array} \label{eq:2-54} \tag{94} \end{equation} $$

With the new creation and annihilation operator we can now construct many-body quasiparticle states, with one-particle-one-hole states, two-particle-two-hole states etc in the same fashion as we previously constructed many-particle states. We can write a general particle-hole state as

$$ \begin{equation} |\beta_1\beta_2\dots \beta_{n_p} \gamma_1^{-1} \gamma_2^{-1} \dots \gamma_{n_h}^{-1}\rangle \equiv \underbrace{b_{\beta_1}^\dagger b_{\beta_2}^\dagger \dots b_{\beta_{n_p}}^\dagger}_{>F} \underbrace{b_{\gamma_1}^\dagger b_{\gamma_2}^\dagger \dots b_{\gamma_{n_h}}^\dagger}_{\leq F} |c\rangle \label{eq:2-56} \tag{95} \end{equation} $$

We can now rewrite our one-body and two-body operators in terms of the new creation and annihilation operators. The number operator becomes

$$ \begin{equation} \hat{N} = \sum_\alpha a_\alpha^\dagger a_\alpha= \sum_{\alpha > F} b_\alpha^\dagger b_\alpha + n_c - \sum_{\alpha \leq F} b_\alpha^\dagger b_\alpha \label{eq:2-57b} \tag{96} \end{equation} $$

where $n_c$ is the number of particle in the new vacuum state $|c\rangle $.
The action of $\hat{N}$ on a many-body state results in

$$ \begin{equation} N |\beta_1\beta_2\dots \beta_{n_p} \gamma_1^{-1} \gamma_2^{-1} \dots \gamma_{n_h}^{-1}\rangle = (n_p + n_c - n_h) |\beta_1\beta_2\dots \beta_{n_p} \gamma_1^{-1} \gamma_2^{-1} \dots \gamma_{n_h}^{-1}\rangle \label{2-59} \tag{97} \end{equation} $$

Here $n=n_p +n_c - n_h$ is the total number of particles in the quasi-particle state of Eq. (95). Note that $\hat{N}$ counts the total number of particles present

$$ \begin{equation} N_{qp} = \sum_\alpha b_\alpha^\dagger b_\alpha, \label{eq:2-60} \tag{98} \end{equation} $$

gives us the number of quasi-particles as can be seen by computing

$$ \begin{equation} N_{qp}= |\beta_1\beta_2\dots \beta_{n_p} \gamma_1^{-1} \gamma_2^{-1} \dots \gamma_{n_h}^{-1}\rangle = (n_p + n_h)|\beta_1\beta_2\dots \beta_{n_p} \gamma_1^{-1} \gamma_2^{-1} \dots \gamma_{n_h}^{-1}\rangle \label{eq:2-61} \tag{99} \end{equation} $$

where $n_{qp} = n_p + n_h$ is the total number of quasi-particles.

We express the one-body operator $\hat{H}_0$ in terms of the quasi-particle creation and annihilation operators, resulting in

$$ \begin{eqnarray} \hat{H}_0 &=& \sum_{\alpha\beta > F} \langle \alpha|\hat{h}_0|\beta\rangle b_\alpha^\dagger b_\beta + \sum_{\begin{array}{c} \alpha > F \\ \beta \leq F \end{array}} \left[ \langle \alpha|\hat{h}_0|\beta\rangle b_\alpha^\dagger b_\beta^\dagger + \langle \beta|\hat{h}_0|\alpha\rangle b_\beta b_\alpha \right] \nonumber \\ &+& \sum_{\alpha \leq F} \langle \alpha|\hat{h}_0|\alpha\rangle - \sum_{\alpha\beta \leq F} \langle \beta|\hat{h}_0|\alpha\rangle b_\alpha^\dagger b_\beta \label{eq:2-63b} \tag{100} \end{eqnarray} $$

The first term gives contribution only for particle states, while the last one contributes only for holestates. The second term can create or destroy a set of quasi-particles and the third term is the contribution from the vacuum state $|c\rangle$.

Before we continue with the expressions for the two-body operator, we introduce a nomenclature we will use for the rest of this text. It is inspired by the notation used in quantum chemistry. We reserve the labels $i,j,k,\dots$ for hole states and $a,b,c,\dots$ for states above $F$, viz. particle states. This means also that we will skip the constraint $\leq F$ or $> F$ in the summation symbols. Our operator $\hat{H}_0$ reads now

$$ \begin{eqnarray} \hat{H}_0 &=& \sum_{ab} \langle a|\hat{h}|b\rangle b_a^\dagger b_b + \sum_{ai} \left[ \langle a|\hat{h}|i\rangle b_a^\dagger b_i^\dagger + \langle i|\hat{h}|a\rangle b_i b_a \right] \nonumber \\ &+& \sum_{i} \langle i|\hat{h}|i\rangle - \sum_{ij} \langle j|\hat{h}|i\rangle b_i^\dagger b_j \label{eq:2-63b} \tag{101} \end{eqnarray} $$

The two-particle operator in the particle-hole formalism is more complicated since we have to translate four indices $\alpha\beta\gamma\delta$ to the possible combinations of particle and hole states. When performing the commutator algebra we can regroup the operator in five different terms

$$ \begin{equation} \hat{H}_I = \hat{H}_I^{(a)} + \hat{H}_I^{(b)} + \hat{H}_I^{(c)} + \hat{H}_I^{(d)} + \hat{H}_I^{(e)} \label{eq:2-65} \tag{102} \end{equation} $$

Using anti-symmetrized matrix elements, bthe term $\hat{H}_I^{(a)}$ is

$$ \begin{equation} \hat{H}_I^{(a)} = \frac{1}{4} \sum_{abcd} \langle ab|\hat{V}|cd\rangle b_a^\dagger b_b^\dagger b_d b_c \label{eq:2-66} \tag{103} \end{equation} $$

The next term $\hat{H}_I^{(b)}$ reads

$$ \begin{equation} \hat{H}_I^{(b)} = \frac{1}{4} \sum_{abci}\left(\langle ab|\hat{V}|ci\rangle b_a^\dagger b_b^\dagger b_i^\dagger b_c +\langle ai|\hat{V}|cb\rangle b_a^\dagger b_i b_b b_c\right) \label{eq:2-67b} \tag{104} \end{equation} $$

This term conserves the number of quasiparticles but creates or removes a three-particle-one-hole state. For $\hat{H}_I^{(c)}$ we have

$$ \begin{eqnarray} \hat{H}_I^{(c)}& =& \frac{1}{4} \sum_{abij}\left(\langle ab|\hat{V}|ij\rangle b_a^\dagger b_b^\dagger b_j^\dagger b_i^\dagger + \langle ij|\hat{V}|ab\rangle b_a b_b b_j b_i \right)+ \nonumber \\ && \frac{1}{2}\sum_{abij}\langle ai|\hat{V}|bj\rangle b_a^\dagger b_j^\dagger b_b b_i + \frac{1}{2}\sum_{abi}\langle ai|\hat{V}|bi\rangle b_a^\dagger b_b. \label{eq:2-68c} \tag{105} \end{eqnarray} $$

The first line stands for the creation of a two-particle-two-hole state, while the second line represents the creation to two one-particle-one-hole pairs while the last term represents a contribution to the particle single-particle energy from the hole states, that is an interaction between the particle states and the hole states within the new vacuum state. The fourth term reads

$$ \begin{eqnarray} \hat{H}_I^{(d)}& = &\frac{1}{4} \sum_{aijk}\left(\langle ai|\hat{V}|jk\rangle b_a^\dagger b_k^\dagger b_j^\dagger b_i+ \langle ji|\hat{V}|ak\rangle b_k^\dagger b_j b_i b_a\right)+\nonumber \\ &&\frac{1}{4}\sum_{aij}\left(\langle ai|\hat{V}|ji\rangle b_a^\dagger b_j^\dagger+ \langle ji|\hat{V}|ai\rangle - \langle ji|\hat{V}|ia\rangle b_j b_a \right). \label{eq:2-69d} \tag{106} \end{eqnarray} $$

The terms in the first line stand for the creation of a particle-hole state interacting with hole states, we will label this as a two-hole-one-particle contribution. The remaining terms are a particle-hole state interacting with the holes in the vacuum state. Finally we have

$$ \begin{equation} \hat{H}_I^{(e)} = \frac{1}{4} \sum_{ijkl} \langle kl|\hat{V}|ij\rangle b_i^\dagger b_j^\dagger b_l b_k+ \frac{1}{2}\sum_{ijk}\langle ij|\hat{V}|kj\rangle b_k^\dagger b_i +\frac{1}{2}\sum_{ij}\langle ij|\hat{V}|ij\rangle \label{eq:2-70d} \tag{107} \end{equation} $$

The first terms represents the interaction between two holes while the second stands for the interaction between a hole and the remaining holes in the vacuum state. It represents a contribution to single-hole energy to first order. The last term collects all contributions to the energy of the ground state of a closed-shell system arising from hole-hole correlations.

Summarizing and defining a normal-ordered Hamiltonian

$$ \Phi_{AS}(\alpha_1, \dots, \alpha_A; x_1, \dots x_A)= \frac{1}{\sqrt{A}} \sum_{\hat{P}} (-1)^P \hat{P} \prod_{i=1}^A \psi_{\alpha_i}(x_i), $$

which is equivalent with $|\alpha_1 \dots \alpha_A\rangle= a_{\alpha_1}^{\dagger} \dots a_{\alpha_A}^{\dagger} |0\rangle$. We have also

3 2 8

< < < ! ! M A T H _ B L O C K

$$ \delta_{pq} = \left\{a_p, a_q^\dagger \right\}, $$

and

3 3 0

< < < ! ! M A T H _ B L O C K

$$ |\Phi_0\rangle = |\alpha_1 \dots \alpha_A\rangle, \quad \alpha_1, \dots, \alpha_A \leq \alpha_F $$

3 3 2

< < < ! ! M A T H _ B L O C K

$$ \left\{a_p, a_q^\dagger \right\} = \delta_{pq}, p, q > \alpha_F $$

with $i,j,\ldots \leq \alpha_F, \quad a,b,\ldots > \alpha_F, \quad p,q, \ldots - \textrm{any}$

$$ a_i|\Phi_0\rangle = |\Phi_i\rangle, \hspace{0.5cm} a_a^\dagger|\Phi_0\rangle = |\Phi^a\rangle $$

and

$$ a_i^\dagger|\Phi_0\rangle = 0 \hspace{0.5cm} a_a|\Phi_0\rangle = 0 $$

The one-body operator is defined as

$$ \hat{F} = \sum_{pq} \langle p|\hat{f}|q\rangle a_p^\dagger a_q $$

while the two-body opreator is defined as

$$ \hat{V} = \frac{1}{4} \sum_{pqrs} \langle pq|\hat{v}|rs\rangle_{AS} a_p^\dagger a_q^\dagger a_s a_r $$

where we have defined the antisymmetric matrix elements

$$ \langle pq|\hat{v}|rs\rangle_{AS} = \langle pq|\hat{v}|rs\rangle - \langle pq|\hat{v}|sr\rangle. $$

We can also define a three-body operator

$$ \hat{V}_3 = \frac{1}{36} \sum_{pqrstu} \langle pqr|\hat{v}_3|stu\rangle_{AS} a_p^\dagger a_q^\dagger a_r^\dagger a_u a_t a_s $$

with the antisymmetrized matrix element

$$ \begin{equation} \langle pqr|\hat{v}_3|stu\rangle_{AS} = \langle pqr|\hat{v}_3|stu\rangle + \langle pqr|\hat{v}_3|tus\rangle + \langle pqr|\hat{v}_3|ust\rangle- \langle pqr|\hat{v}_3|sut\rangle - \langle pqr|\hat{v}_3|tsu\rangle - \langle pqr|\hat{v}_3|uts\rangle. \end{equation} $$

Hartree-Fock in second quantization and stability of HF solution

We wish now to derive the Hartree-Fock equations using our second-quantized formalism and study the stability of the equations. Our ansatz for the ground state of the system is approximated as (this is our representation of a Slater determinant in second quantization)

$$ |\Phi_0\rangle = |c\rangle = a^{\dagger}_i a^{\dagger}_j \dots a^{\dagger}_l|0\rangle. $$

We wish to determine $\hat{u}^{HF}$ so that $E_0^{HF}= \langle c|\hat{H}| c\rangle$ becomes a local minimum.

In our analysis here we will need Thouless' theorem, which states that an arbitrary Slater determinant $|c'\rangle$ which is not orthogonal to a determinant $| c\rangle ={\displaystyle\prod_{i=1}^{n}} a_{\alpha_{i}}^{\dagger}|0\rangle$, can be written as

$$ |c'\rangle=exp\left\{\sum_{a>F}\sum_{i\le F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| c\rangle $$

Let us give a simple proof of Thouless' theorem. The theorem states that we can make a linear combination av particle-hole excitations with respect to a given reference state $\vert c\rangle$. With this linear combination, we can make a new Slater determinant $\vert c'\rangle $ which is not orthogonal to $\vert c\rangle$, that is

$$ \langle c|c'\rangle \ne 0. $$

To show this we need some intermediate steps. The exponential product of two operators $\exp{\hat{A}}\times\exp{\hat{B}}$ is equal to $\exp{(\hat{A}+\hat{B})}$ only if the two operators commute, that is

$$ [\hat{A},\hat{B}] = 0. $$

Thouless' theorem

If the operators do not commute, we need to resort to the Baker-Campbell-Hauersdorf. This relation states that

$$ \exp{\hat{C}}=\exp{\hat{A}}\exp{\hat{B}}, $$

with

$$ \hat{C}=\hat{A}+\hat{B}+\frac{1}{2}[\hat{A},\hat{B}]+\frac{1}{12}[[\hat{A},\hat{B}],\hat{B}]-\frac{1}{12}[[\hat{A},\hat{B}],\hat{A}]+\dots $$

From these relations, we note that in our expression for $|c'\rangle$ we have commutators of the type

$$ [a_{a}^{\dagger}a_{i},a_{b}^{\dagger}a_{j}], $$

and it is easy to convince oneself that these commutators, or higher powers thereof, are all zero. This means that we can write out our new representation of a Slater determinant as

$$ |c'\rangle=exp\left\{\sum_{a>F}\sum_{i\le F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| c\rangle=\prod_{i}\left\{1+\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}+\left(\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\right)^2+\dots\right\}| c\rangle $$

We note that

$$ \prod_{i}\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\sum_{b>F}C_{bi}a_{b}^{\dagger}a_{i}| c\rangle =0, $$

and all higher-order powers of these combinations of creation and annihilation operators disappear due to the fact that $(a_i)^n| c\rangle =0$ when $n > 1$. This allows us to rewrite the expression for $|c'\rangle $ as

$$ |c'\rangle=\prod_{i}\left\{1+\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| c\rangle, $$

which we can rewrite as

$$ |c'\rangle=\prod_{i}\left\{1+\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| a^{\dagger}_{i_1} a^{\dagger}_{i_2} \dots a^{\dagger}_{i_n}|0\rangle. $$

The last equation can be written as

$$ \begin{equation} |c'\rangle=\prod_{i}\left\{1+\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| a^{\dagger}_{i_1} a^{\dagger}_{i_2} \dots a^{\dagger}_{i_n}|0\rangle=\left(1+\sum_{a>F}C_{ai_1}a_{a}^{\dagger}a_{i_1}\right)a^{\dagger}_{i_1} \end{equation} $$

$$ \begin{equation} \times\left(1+\sum_{a>F}C_{ai_2}a_{a}^{\dagger}a_{i_2}\right)a^{\dagger}_{i_2} \dots |0\rangle=\prod_{i}\left(a^{\dagger}_{i}+\sum_{a>F}C_{ai}a_{a}^{\dagger}\right)|0\rangle. \end{equation} $$

New operators

If we define a new creation operator

$$ \begin{equation} b^{\dagger}_{i}=a^{\dagger}_{i}+\sum_{a>F}C_{ai}a_{a}^{\dagger}, \label{eq:newb} \tag{108} \end{equation} $$

we have

$$ |c'\rangle=\prod_{i}b^{\dagger}_{i}|0\rangle=\prod_{i}\left(a^{\dagger}_{i}+\sum_{a>F}C_{ai}a_{a}^{\dagger}\right)|0\rangle, $$

meaning that the new representation of the Slater determinant in second quantization, $|c'\rangle$, looks like our previous ones. However, this representation is not general enough since we have a restriction on the sum over single-particle states in Eq. (108). The single-particle states have all to be above the Fermi level. The question then is whether we can construct a general representation of a Slater determinant with a creation operator

$$ \tilde{b}^{\dagger}_{i}=\sum_{p}f_{ip}a_{p}^{\dagger}, $$

where $f_{ip}$ is a matrix element of a unitary matrix which transforms our creation and annihilation operators $a^{\dagger}$ and $a$ to $\tilde{b}^{\dagger}$ and $\tilde{b}$. These new operators define a new representation of a Slater determinant as

$$ |\tilde{c}\rangle=\prod_{i}\tilde{b}^{\dagger}_{i}|0\rangle. $$

Showing that $|\tilde{c}\rangle= |c'\rangle$

We need to show that $|\tilde{c}\rangle= |c'\rangle$. We need also to assume that the new state is not orthogonal to $|c\rangle$, that is $\langle c| \tilde{c}\rangle \ne 0$. From this it follows that

$$ \langle c| \tilde{c}\rangle=\langle 0| a_{i_n}\dots a_{i_1}\left(\sum_{p=i_1}^{i_n}f_{i_1p}a_{p}^{\dagger} \right)\left(\sum_{q=i_1}^{i_n}f_{i_2q}a_{q}^{\dagger} \right)\dots \left(\sum_{t=i_1}^{i_n}f_{i_nt}a_{t}^{\dagger} \right)|0\rangle, $$

which is nothing but the determinant $det(f_{ip})$ which we can, using the intermediate normalization condition, normalize to one, that is

$$ det(f_{ip})=1, $$

meaning that $f$ has an inverse defined as (since we are dealing with orthogonal, and in our case unitary as well, transformations)

$$ \sum_{k} f_{ik}f^{-1}_{kj} = \delta_{ij}, $$

and

$$ \sum_{j} f^{-1}_{ij}f_{jk} = \delta_{ik}. $$

Using these relations we can then define the linear combination of creation (and annihilation as well) operators as

$$ \sum_{i}f^{-1}_{ki}\tilde{b}^{\dagger}_{i}=\sum_{i}f^{-1}_{ki}\sum_{p=i_1}^{\infty}f_{ip}a_{p}^{\dagger}=a_{k}^{\dagger}+\sum_{i}\sum_{p=i_{n+1}}^{\infty}f^{-1}_{ki}f_{ip}a_{p}^{\dagger}. $$

Defining

$$ c_{kp}=\sum_{i \le F}f^{-1}_{ki}f_{ip}, $$

we can redefine

$$ a_{k}^{\dagger}+\sum_{i}\sum_{p=i_{n+1}}^{\infty}f^{-1}_{ki}f_{ip}a_{p}^{\dagger}=a_{k}^{\dagger}+\sum_{p=i_{n+1}}^{\infty}c_{kp}a_{p}^{\dagger}=b_k^{\dagger}, $$

our starting point. We have shown that our general representation of a Slater determinant

$$ |\tilde{c}\rangle=\prod_{i}\tilde{b}^{\dagger}_{i}|0\rangle=|c'\rangle=\prod_{i}b^{\dagger}_{i}|0\rangle, $$

with

$$ b_k^{\dagger}=a_{k}^{\dagger}+\sum_{p=i_{n+1}}^{\infty}c_{kp}a_{p}^{\dagger}. $$

This means that we can actually write an ansatz for the ground state of the system as a linear combination of terms which contain the ansatz itself $|c\rangle$ with an admixture from an infinity of one-particle-one-hole states. The latter has important consequences when we wish to interpret the Hartree-Fock equations and their stability. We can rewrite the new representation as

$$ |c'\rangle = |c\rangle+|\delta c\rangle, $$

where $|\delta c\rangle$ can now be interpreted as a small variation. If we approximate this term with contributions from one-particle-one-hole (1p-1h) states only, we arrive at

$$ |c'\rangle = \left(1+\sum_{ai}\delta C_{ai}a_{a}^{\dagger}a_i\right)|c\rangle. $$

In our derivation of the Hartree-Fock equations we have shown that

$$ \langle \delta c| \hat{H} | c\rangle =0, $$

which means that we have to satisfy

$$ \langle c|\sum_{ai}\delta C_{ai}\left\{a_{a}^{\dagger}a_i\right\} \hat{H} | c\rangle =0. $$

With this as a background, we are now ready to study the stability of the Hartree-Fock equations.

Hartree-Fock in second quantization and stability of HF solution

The variational condition for deriving the Hartree-Fock equations guarantees only that the expectation value $\langle c | \hat{H} | c \rangle$ has an extreme value, not necessarily a minimum. To figure out whether the extreme value we have found is a minimum, we can use second quantization to analyze our results and find a criterion for the above expectation value to a local minimum. We will use Thouless' theorem and show that

$$ \frac{\langle c' |\hat{H} | c'\rangle}{\langle c' |c'\rangle} \ge \langle c |\hat{H} | c\rangle= E_0, $$

with

$$ {|c'\rangle} = {|c\rangle + |\delta c\rangle}. $$

Using Thouless' theorem we can write out ${|c'\rangle}$ as

$$ \begin{equation} {|c'\rangle}=\exp\left\{\sum_{a > F}\sum_{i \le F}\delta C_{ai}a_{a}^{\dagger}a_{i}\right\}| c\rangle \end{equation} $$

$$ \begin{equation} =\left\{1+\sum_{a > F}\sum_{i \le F}\delta C_{ai}a_{a}^{\dagger} a_{i}+\frac{1}{2!}\sum_{ab > F}\sum_{ij \le F}\delta C_{ai}\delta C_{bj}a_{a}^{\dagger}a_{i}a_{b}^{\dagger}a_{j}+\dots\right\} \end{equation} $$

where the amplitudes $\delta C$ are small.

The norm of $|c'\rangle$ is given by (using the intermediate normalization condition $\langle c' |c\rangle=1$)

$$ \langle c' | c'\rangle = 1+\sum_{a>F} \sum_{i\le F}|\delta C_{ai}|^2+O(\delta C_{ai}^3). $$

The expectation value for the energy is now given by (using the Hartree-Fock condition)

3 7 6

< < < ! ! M A T H _ B L O C K

$$ \frac{1}{2!}\sum_{ab>F} \sum_{ij\le F}\delta C_{ai}\delta C_{bj}\langle c |\hat{H}a_{a}^{\dagger}a_{i}a_{b}^{\dagger}a_{j}|c\rangle+\frac{1}{2!}\sum_{ab>F} \sum_{ij\le F}\delta C_{ai}^*\delta C_{bj}^*\langle c|a_{j}^{\dagger}a_{b}a_{i}^{\dagger}a_{a}\hat{H}|c\rangle +\dots $$

We have already calculated the second term on the right-hand side of the previous equation

$$ \begin{equation} \langle c | \left(\{a^\dagger_i a_a\} \hat{H} \{a^\dagger_b a_j\} \right) | c\rangle=\sum_{pq} \sum_{ijab}\delta C_{ai}^*\delta C_{bj} \langle p|\hat{h}_0 |q\rangle \langle c | \left(\{a^{\dagger}_i a_a\}\{a^{\dagger}_pa_q\} \{a^{\dagger}_b a_j\} \right)| c\rangle \end{equation} $$

$$ \begin{equation} +\frac{1}{4} \sum_{pqrs} \sum_{ijab}\delta C_{ai}^*\delta C_{bj} \langle pq| \hat{v}|rs\rangle \langle c | \left(\{a^\dagger_i a_a\}\{a^{\dagger}_p a^{\dagger}_q a_s a_r\} \{a^{\dagger}_b a_j\} \right)| c\rangle , \end{equation} $$

resulting in

$$ E_0\sum_{ai}|\delta C_{ai}|^2+\sum_{ai}|\delta C_{ai}|^2(\varepsilon_a-\varepsilon_i)-\sum_{ijab} \langle aj|\hat{v}| bi\rangle \delta C_{ai}^*\delta C_{bj}. $$

$$ \frac{1}{2!}\langle c |\left(\{a^\dagger_j a_b\} \{a^\dagger_i a_a\} \hat{V}_N \right) | c\rangle = \frac{1}{2!}\langle c |\left( \hat{V}_N \{a^\dagger_a a_i\} \{a^\dagger_b a_j\} \right)^{\dagger} | c\rangle $$

which is nothing but

$$ \frac{1}{2!}\langle c | \left( \hat{V}_N \{a^\dagger_a a_i\} \{a^\dagger_b a_j\} \right) | c\rangle^* =\frac{1}{2} \sum_{ijab} (\langle ij|\hat{v}|ab\rangle)^*\delta C_{ai}^*\delta C_{bj}^* $$

or

$$ \frac{1}{2} \sum_{ijab} (\langle ab|\hat{v}|ij\rangle)\delta C_{ai}^*\delta C_{bj}^* $$

where we have used the relation

$$ \langle a |\hat{A} | b\rangle = (\langle b |\hat{A}^{\dagger} | a\rangle)^* $$

due to the hermiticity of $\hat{H}$ and $\hat{V}$.

We define two matrix elements

$$ A_{ai,bj}=-\langle aj|\hat{v} bi\rangle $$

and

$$ B_{ai,bj}=\langle ab|\hat{v}|ij\rangle $$

both being anti-symmetrized.

With these definitions we write out the energy as

$$ \begin{equation} \langle c'|H|c'\rangle = \left(1+\sum_{ai}|\delta C_{ai}|^2\right)\langle c |H|c\rangle+\sum_{ai}|\delta C_{ai}|^2(\varepsilon_a^{HF}-\varepsilon_i^{HF})+\sum_{ijab}A_{ai,bj}\delta C_{ai}^*\delta C_{bj}+ \end{equation} $$

$$ \begin{equation} \frac{1}{2} \sum_{ijab} B_{ai,bj}^*\delta C_{ai}\delta C_{bj}+\frac{1}{2} \sum_{ijab} B_{ai,bj}\delta C_{ai}^*\delta C_{bj}^* +O(\delta C_{ai}^3), \end{equation} $$

which can be rewritten as

$$ \langle c'|H|c'\rangle = \left(1+\sum_{ai}|\delta C_{ai}|^2\right)\langle c |H|c\rangle+\Delta E+O(\delta C_{ai}^3), $$

and skipping higher-order terms we arrived

$$ \frac{\langle c' |\hat{H} | c'\rangle}{\langle c' |c'\rangle} =E_0+\frac{\Delta E}{\left(1+\sum_{ai}|\delta C_{ai}|^2\right)}. $$

We have defined

$$ \Delta E = \frac{1}{2} \langle \chi | \hat{M}| \chi \rangle $$

with the vectors

$$ \chi = \left[ \delta C\hspace{0.2cm} \delta C^*\right]^T $$

and the matrix

$$ \hat{M}=\left(\begin{array}{cc} \Delta + A & B \\ B^* & \Delta + A^*\end{array}\right), $$

with $\Delta_{ai,bj} = (\varepsilon_a-\varepsilon_i)\delta_{ab}\delta_{ij}$.

The condition

$$ \Delta E = \frac{1}{2} \langle \chi | \hat{M}| \chi \rangle \ge 0 $$

for an arbitrary vector

$$ \chi = \left[ \delta C\hspace{0.2cm} \delta C^*\right]^T $$

means that all eigenvalues of the matrix have to be larger than or equal zero. A necessary (but no sufficient) condition is that the matrix elements (for all $ai$ )

$$ (\varepsilon_a-\varepsilon_i)\delta_{ab}\delta_{ij}+A_{ai,bj} \ge 0. $$

This equation can be used as a first test of the stability of the Hartree-Fock equation.

Operators in second quantization

In the build-up of a shell-model or FCI code that is meant to tackle large dimensionalities is the action of the Hamiltonian $\hat{H}$ on a Slater determinant represented in second quantization as

$$ |\alpha_1\dots \alpha_n\rangle = a_{\alpha_1}^{\dagger} a_{\alpha_2}^{\dagger} \dots a_{\alpha_n}^{\dagger} |0\rangle. $$

The time consuming part stems from the action of the Hamiltonian on the above determinant,

$$ \left(\sum_{\alpha\beta} \langle \alpha|t+u|\beta\rangle a_\alpha^{\dagger} a_\beta + \frac{1}{4} \sum_{\alpha\beta\gamma\delta} \langle \alpha \beta|\hat{v}|\gamma \delta\rangle a_\alpha^{\dagger} a_\beta^{\dagger} a_\delta a_\gamma\right)a_{\alpha_1}^{\dagger} a_{\alpha_2}^{\dagger} \dots a_{\alpha_n}^{\dagger} |0\rangle. $$

A practically useful way to implement this action is to encode a Slater determinant as a bit pattern.

Assume that we have at our disposal $n$ different single-particle orbits $\alpha_0,\alpha_2,\dots,\alpha_{n-1}$ and that we can distribute among these orbits $N\le n$ particles.

A Slater determinant can then be coded as an integer of $n$ bits. As an example, if we have $n=16$ single-particle states $\alpha_0,\alpha_1,\dots,\alpha_{15}$ and $N=4$ fermions occupying the states $\alpha_3$, $\alpha_6$, $\alpha_{10}$ and $\alpha_{13}$ we could write this Slater determinant as

$$ \Phi_{\Lambda} = a_{\alpha_3}^{\dagger} a_{\alpha_6}^{\dagger} a_{\alpha_{10}}^{\dagger} a_{\alpha_{13}}^{\dagger} |0\rangle. $$

The unoccupied single-particle states have bit value $0$ while the occupied ones are represented by bit state $1$. In the binary notation we would write this 16 bits long integer as

$$ \begin{array}{cccccccccccccccc} {\alpha_0}&{\alpha_1}&{\alpha_2}&{\alpha_3}&{\alpha_4}&{\alpha_5}&{\alpha_6}&{\alpha_7} & {\alpha_8} &{\alpha_9} & {\alpha_{10}} &{\alpha_{11}} &{\alpha_{12}} &{\alpha_{13}} &{\alpha_{14}} & {\alpha_{15}} \\ {0} & {0} &{0} &{1} &{0} &{0} &{1} &{0} &{0} &{0} &{1} &{0} &{0} &{1} &{0} & {0} \\ \end{array} $$

which translates into the decimal number

$$ 2^3+2^6+2^{10}+2^{13}=9288. $$

We can thus encode a Slater determinant as a bit pattern.

With $N$ particles that can be distributed over $n$ single-particle states, the total number of Slater determinats (and defining thereby the dimensionality of the system) is

$$ \mathrm{dim}(\mathcal{H}) = \left(\begin{array}{c} n \\N\end{array}\right). $$

The total number of bit patterns is $2^n$.

We assume again that we have at our disposal $n$ different single-particle orbits $\alpha_0,\alpha_2,\dots,\alpha_{n-1}$ and that we can distribute among these orbits $N\le n$ particles. The ordering among these states is important as it defines the order of the creation operators. We will write the determinant

$$ \Phi_{\Lambda} = a_{\alpha_3}^{\dagger} a_{\alpha_6}^{\dagger} a_{\alpha_{10}}^{\dagger} a_{\alpha_{13}}^{\dagger} |0\rangle, $$

in a more compact way as

$$ \Phi_{3,6,10,13} = |0001001000100100\rangle. $$

The action of a creation operator is thus

$$ a^{\dagger}_{\alpha_4}\Phi_{3,6,10,13} = a^{\dagger}_{\alpha_4}|0001001000100100\rangle=a^{\dagger}_{\alpha_4}a_{\alpha_3}^{\dagger} a_{\alpha_6}^{\dagger} a_{\alpha_{10}}^{\dagger} a_{\alpha_{13}}^{\dagger} |0\rangle, $$

which becomes

$$ -a_{\alpha_3}^{\dagger} a^{\dagger}_{\alpha_4} a_{\alpha_6}^{\dagger} a_{\alpha_{10}}^{\dagger} a_{\alpha_{13}}^{\dagger} |0\rangle=-|0001101000100100\rangle. $$

Similarly

$$ a^{\dagger}_{\alpha_6}\Phi_{3,6,10,13} = a^{\dagger}_{\alpha_6}|0001001000100100\rangle=a^{\dagger}_{\alpha_6}a_{\alpha_3}^{\dagger} a_{\alpha_6}^{\dagger} a_{\alpha_{10}}^{\dagger} a_{\alpha_{13}}^{\dagger} |0\rangle, $$

which becomes

$$ -a^{\dagger}_{\alpha_4} (a_{\alpha_6}^{\dagger})^ 2 a_{\alpha_{10}}^{\dagger} a_{\alpha_{13}}^{\dagger} |0\rangle=0! $$

This gives a simple recipe:

• If one of the bits $b_j$ is $1$ and we act with a creation operator on this bit, we return a null vector

• If $b_j=0$, we set it to $1$ and return a sign factor $(-1)^l$, where $l$ is the number of bits set before bit $j$.

Consider the action of $a^{\dagger}_{\alpha_2}$ on various slater determinants:

$$ \begin{array}{ccc} a^{\dagger}_{\alpha_2}\Phi_{00111}& = a^{\dagger}_{\alpha_2}|00111\rangle&=0\times |00111\rangle\\ a^{\dagger}_{\alpha_2}\Phi_{01011}& = a^{\dagger}_{\alpha_2}|01011\rangle&=(-1)\times |01111\rangle\\ a^{\dagger}_{\alpha_2}\Phi_{01101}& = a^{\dagger}_{\alpha_2}|01101\rangle&=0\times |01101\rangle\\ a^{\dagger}_{\alpha_2}\Phi_{01110}& = a^{\dagger}_{\alpha_2}|01110\rangle&=0\times |01110\rangle\\ a^{\dagger}_{\alpha_2}\Phi_{10011}& = a^{\dagger}_{\alpha_2}|10011\rangle&=(-1)\times |10111\rangle\\ a^{\dagger}_{\alpha_2}\Phi_{10101}& = a^{\dagger}_{\alpha_2}|10101\rangle&=0\times |10101\rangle\\ a^{\dagger}_{\alpha_2}\Phi_{10110}& = a^{\dagger}_{\alpha_2}|10110\rangle&=0\times |10110\rangle\\ a^{\dagger}_{\alpha_2}\Phi_{11001}& = a^{\dagger}_{\alpha_2}|11001\rangle&=(+1)\times |11101\rangle\\ a^{\dagger}_{\alpha_2}\Phi_{11010}& = a^{\dagger}_{\alpha_2}|11010\rangle&=(+1)\times |11110\rangle\\ \end{array} $$

What is the simplest way to obtain the phase when we act with one annihilation(creation) operator on the given Slater determinant representation?

We have an SD representation

$$ \Phi_{\Lambda} = a_{\alpha_0}^{\dagger} a_{\alpha_3}^{\dagger} a_{\alpha_6}^{\dagger} a_{\alpha_{10}}^{\dagger} a_{\alpha_{13}}^{\dagger} |0\rangle, $$

in a more compact way as

$$ \Phi_{0,3,6,10,13} = |1001001000100100\rangle. $$

The action of

$$ a^{\dagger}_{\alpha_4}a_{\alpha_0}\Phi_{0,3,6,10,13} = a^{\dagger}_{\alpha_4}|0001001000100100\rangle=a^{\dagger}_{\alpha_4}a_{\alpha_3}^{\dagger} a_{\alpha_6}^{\dagger} a_{\alpha_{10}}^{\dagger} a_{\alpha_{13}}^{\dagger} |0\rangle, $$

which becomes

$$ -a_{\alpha_3}^{\dagger} a^{\dagger}_{\alpha_4} a_{\alpha_6}^{\dagger} a_{\alpha_{10}}^{\dagger} a_{\alpha_{13}}^{\dagger} |0\rangle=-|0001101000100100\rangle. $$

The action

$$ a_{\alpha_0}\Phi_{0,3,6,10,13} = |0001001000100100\rangle, $$

can be obtained by subtracting the logical sum (AND operation) of $\Phi_{0,3,6,10,13}$ and a word which represents only $\alpha_0$, that is

$$ |1000000000000000\rangle, $$

from $\Phi_{0,3,6,10,13}= |1001001000100100\rangle$.

This operation gives $|0001001000100100\rangle$.

Similarly, we can form $a^{\dagger}_{\alpha_4}a_{\alpha_0}\Phi_{0,3,6,10,13}$, say, by adding $|0000100000000000\rangle$ to $a_{\alpha_0}\Phi_{0,3,6,10,13}$, first checking that their logical sum is zero in order to make sure that orbital $\alpha_4$ is not already occupied.

It is trickier however to get the phase $(-1)^l$. One possibility is as follows

• Let $S_1$ be a word that represents the $1-$bit to be removed and all others set to zero.

In the previous example $S_1=|1000000000000000\rangle$

• Define $S_2$ as the similar word that represents the bit to be added, that is in our case

$S_2=|0000100000000000\rangle$.

• Compute then $S=S_1-S_2$, which here becomes

$$ S=|0111000000000000\rangle $$

• Perform then the logical AND operation of $S$ with the word containing

$$ \Phi_{0,3,6,10,13} = |1001001000100100\rangle, $$

which results in $|0001000000000000\rangle$. Counting the number of $1-$bits gives the phase. Here you need however an algorithm for bitcounting. Several efficient ones available.

Exercise 7: Relation between basis functions

This exercise serves to convince you about the relation between two different single-particle bases, where one could be our new Hartree-Fock basis and the other a harmonic oscillator basis.

Consider a Slater determinant built up of single-particle orbitals $\psi_{\lambda}$, with $\lambda = 1,2,\dots,A$. The unitary transformation

$$ \psi_a = \sum_{\lambda} C_{a\lambda}\phi_{\lambda}, $$

brings us into the new basis.
The new basis has quantum numbers $a=1,2,\dots,A$. Show that the new basis is orthonormal. Show that the new Slater determinant constructed from the new single-particle wave functions can be written as the determinant based on the previous basis and the determinant of the matrix $C$. Show that the old and the new Slater determinants are equal up to a complex constant with absolute value unity. (Hint, $C$ is a unitary matrix).

Starting with the second quantization representation of the Slater determinant

$$ \Phi_{0}=\prod_{i=1}^{n}a_{\alpha_{i}}^{\dagger}|0\rangle, $$

use Wick's theorem to compute the normalization integral $\langle\Phi_{0}|\Phi_{0}\rangle$.

Exercise 8: Matrix elements

Calculate the matrix elements

$$ \langle \alpha_{1}\alpha_{2}|\hat{F}|\alpha_{1}\alpha_{2}\rangle $$

and

$$ \langle \alpha_{1}\alpha_{2}|\hat{G}|\alpha_{1}\alpha_{2}\rangle $$

with

4 2 2

< < < ! ! M A T H _ B L O C K

4 2 3

< < < ! ! M A T H _ B L O C K

4 2 4

< < < ! ! M A T H _ B L O C K

$$ \hat{G} = \frac{1}{2}\sum_{\alpha\beta\gamma\delta} \langle \alpha\beta |\hat{g}|\gamma\delta\rangle a_{\alpha}^{\dagger}a_{\beta}^{\dagger}a_{\delta}a_{\gamma} , $$

and

$$ \langle \alpha\beta |\hat{g}|\gamma\delta\rangle= \int\int \psi_{\alpha}^{*}(x_{1})\psi_{\beta}^{*}(x_{2})g(x_{1}, x_{2})\psi_{\gamma}(x_{1})\psi_{\delta}(x_{2})dx_{1}dx_{2} $$

Compare these results with those from exercise 3c).

Exercise 9: Normal-ordered one-body operator

Show that the onebody part of the Hamiltonian

$$ \hat{H}_0 = \sum_{pq} \langle p|\hat{h}_0|q\rangle a^{\dagger}_p a_q, $$

can be written, using standard annihilation and creation operators, in normal-ordered form as

$$ \hat{H}_0 = \sum_{pq} \langle p|\hat{h}_0|q\rangle \left\{a^\dagger_p a_q\right\} + \sum_i \langle i|\hat{h}_0|i\rangle. $$

Explain the meaning of the various symbols. Which reference vacuum has been used?

Exercise 10: Normal-ordered two-body operator

Show that the twobody part of the Hamiltonian

$$ \hat{H}_I = \frac{1}{4} \sum_{pqrs} \langle pq|\hat{v}|rs\rangle a^\dagger_p a^\dagger_q a_s a_r, $$

can be written, using standard annihilation and creation operators, in normal-ordered form as

$$ \hat{H}_I =\frac{1}{4} \sum_{pqrs} \langle pq|\hat{v}|rs\rangle \left\{a^\dagger_p a^\dagger_q a_s a_r\right\} + \sum_{pqi} \langle pi|\hat{v}|qi\rangle \left\{a^\dagger_p a_q\right\} + \frac{1}{2} \sum_{ij}\langle ij|\hat{v}|ij\rangle. $$

Explain again the meaning of the various symbols.

This exercise is optional: Derive the normal-ordered form of the threebody part of the Hamiltonian.

$$ \hat{H}_3 = \frac{1}{36} \sum_{\substack{pqr \\ stu}} \langle pqr|\hat{v}_3|stu\rangle a^\dagger_p a^\dagger_q a^\dagger_r a_u a_t a_s, $$

and specify the contributions to the twobody, onebody and the scalar part.

Exercise 11: Matrix elements using the Slater-Condon rule

The aim of this exercise is to set up specific matrix elements that will turn useful when we start our discussions of the nuclear shell model. In particular you will notice, depending on the character of the operator, that many matrix elements will actually be zero.

Consider three $N$-particle Slater determinants $|\Phi_0$, $|\Phi_i^a\rangle$ and $|\Phi_{ij}^{ab}\rangle$, where the notation means that Slater determinant $|\Phi_i^a\rangle$ differs from $|\Phi_0\rangle$ by one single-particle state, that is a single-particle state $\psi_i$ is replaced by a single-particle state $\psi_a$. It is often interpreted as a so-called one-particle-one-hole excitation. Similarly, the Slater determinant $|\Phi_{ij}^{ab}\rangle$ differs by two single-particle states from $|\Phi_0\rangle$ and is normally thought of as a two-particle-two-hole excitation. We assume also that $|\Phi_0\rangle$ represents our new vacuum reference state and the labels $ijk\dots$ represent single-particle states below the Fermi level and $abc\dots$ represent states above the Fermi level, so-called particle states. We define thereafter a general onebody normal-ordered (with respect to the new vacuum state) operator as

$$ \hat{F}_N=\sum_{pq}\langle p |f |\beta\rangle \left\{a_{p}^{\dagger}a_{q}\right\} , $$

with

$$ \langle p |f| q\rangle=\int \psi_{p}^{*}(x)f(x)\psi_{q}(x)dx , $$

and a general normal-ordered two-body operator

$$ \hat{G}_N = \frac{1}{4}\sum_{pqrs} \langle pq |g| rs\rangle_{AS} \left\{a_{p}^{\dagger}a_{q}^{\dagger}a_{s}a_{r}\right\} , $$

with for example the direct matrix element given as

$$ \langle pq |g| rs\rangle= \int\int \psi_{p}^{*}(x_{1})\psi_{q}^{*}(x_{2})g(x_{1}, x_{2})\psi_{r}(x_{1})\psi_{s}(x_{2})dx_{1}dx_{2} $$

with $g$ being invariant under the interchange of the coordinates of two particles. The single-particle states $\psi_i$ are not necessarily eigenstates of $\hat{f}$. The curly brackets mean that the operators are normal-ordered with respect to the new vacuum reference state.

How would you write the above Slater determinants in a second quantization formalism, utilizing the fact that we have defined a new reference state?

Use thereafter Wick's theorem to find the expectation values of

$$ \langle \Phi_0 \vert\hat{F}_N\vert\Phi_0\rangle, $$

and

$$ \langle \Phi_0\hat{G}_N|\Phi_0\rangle. $$

Find thereafter

$$ \langle \Phi_0 |\hat{F}_N|\Phi_i^a\rangle, $$

and

$$ \langle \Phi_0|\hat{G}_N|\Phi_i^a\rangle, $$

Finally, find

$$ \langle \Phi_0 |\hat{F}_N|\Phi_{ij}^{ab}\rangle, $$

and

$$ \langle \Phi_0|\hat{G}_N|\Phi_{ij}^{ab}\rangle. $$

What happens with the two-body operator if we have a transition probability of the type

$$ \langle \Phi_0|\hat{G}_N|\Phi_{ijk}^{abc}\rangle, $$

where the Slater determinant to the right of the operator differs by more than two single-particle states?

Exercise 12: Program to set up Slater determinants

Write a program which sets up all possible Slater determinants given $N=4$ eletrons which can occupy the atomic single-particle states defined by the $1s$, $2s2p$ and $3s3p3d$ shells. How many single-particle states $n$ are there in total? Include the spin degrees of freedom as well.

Exercise 13: Using sympy to compute matrix elements

Compute the matrix element

$$ \langle\alpha_{1}\alpha_{2}\alpha_{3}|\hat{G}|\alpha_{1}'\alpha_{2}'\alpha_{3}'\rangle, $$

using Wick's theorem and express the two-body operator $G$ in the occupation number (second quantization) representation.

Exercise 14: Using sympy to compute matrix elements

The last exercise can be solved using the symbolic Python package called SymPy. SymPy is a Python package for general purpose symbolic algebra. There is a physics module with several interesting submodules. Among these, the submodule called secondquant, contains several functionalities that allow us to test our algebraic manipulations using Wick's theorem and operators for second quantization.

The code defines single-particle states above and below the Fermi level, in addition to the genereal symbols $pq$ which can refer to any type of state below or above the Fermi level. Wick's theorem is implemented between the creation and annihilation operators Fd and F, respectively. Using the simplify option, one can lump together several Kronecker-$\delta$ functions.

Exercise 15: Using sympy to compute matrix elements

We can expand the above Python code by defining one-body and two-body operators using the following SymPy code

Here we have used the AntiSymmetricTensor functionality, together with normal-ordering defined by the NO function. Using the latex option, this program produces the following output

$$ f^{p}_{q} \left\{a^\dagger_{p} a_{q}\right\} - \frac{1}{4} v^{qp}_{sr} \left\{a^\dagger_{p} a^\dagger_{q} a_{r} a_{s}\right\} $$

Exercise 16: Using sympy to compute matrix elements

We can now use this code to compute the matrix elements between two two-body Slater determinants using Wick's theorem.

The result is as expected,

$$ \delta_{a c} f^{b}_{d} - \delta_{a d} f^{b}_{c} - \delta_{b c} f^{a}_{d} + \delta_{b d} f^{a}_{c} + v^{ab}_{cd}. $$

Exercise 17: Using sympy to compute matrix elements

We can continue along these lines and define a normal-ordered Hamiltonian with respect to a given reference state. In our first step we just define the Hamiltonian

which results in

$$ f^{q}_{p} a^\dagger_{q} a_{p} + \frac{1}{4} v^{sr}_{qp} a^\dagger_{s} a^\dagger_{r} a_{p} a_{q} $$

Exercise 18: Using sympy to compute matrix elements

In our next step we define the reference energy $E_0$ and redefine the Hamiltonian by subtracting the reference energy and collecting the coefficients for all normal-ordered products (given by the NO function).

which gives us

$$ - f^{i}_{i} + f^{q}_{p} a^\dagger_{q} a_{p} - \frac{1}{4} v^{ii}_{ii} - \frac{1}{4} v^{ii}_{ii} + \frac{1}{4} v^{sr}_{qp} a^\dagger_{r} a^\dagger_{s} a_{q} a_{p}, $$

again as expected, with the reference energy to be subtracted.

Exercise 19: Using sympy to compute matrix elements

We can now go back to exercise 7 and define the Hamiltonian and the second-quantized representation of a three-body Slater determinant.