The project is divided in three main parts. The first part deals with a simple pairing model and the development of a shell-model program related to this model. This program can then be developed into a more general shell-model program that allows you to study general nuclear structure problems. That is the second part of the project. In parallel, we will also use NushellX in order to perform more advanced shell-model studies and compare the results obtained with your own shell-model code to those of NushellX. We expect you to form working groups consisting of typically three (or more) participants. Every group should establish its own Github or Gitlab repository for the project.
In the first part of the project we will thus work with a simplified Hamiltonian consisting of a one-body operator and a so-called pairing interaction term. It is a model which to a large extent mimicks some central features of atomic nuclei, certain atoms and systems which exhibit superfluiditity or superconductivity. Pairing plays a central role in nuclear physics, in particular, for identical particles it makes up large fractions of the correlations among particles. The partial wave \( ^{1}S_0 \) of the nucleon-nucleon force plays a central role in setting up pairing correlations in nuclei. Without this particular partial wave, the \( J=0 \) ground state spin assignment for many nuclei with even numbers of particles would not be possible.
We define first the Hamiltonian, with a definition of the model space and the single-particle basis. Thereafter, we present the various steps which are needed to develop a shell-model program for studying the pairing problem.
The Hamiltonian acting in the complete Hilbert space (usually infinite dimensional) consists of an unperturbed one-body part, \( \hat{H}_0 \), and a perturbation \( \hat{H}_I \).
We limit ourselves to at most two-body interactions, our Hamiltonian is then represented by the following operators
$$
\begin{equation}
\hat{H} = \hat{H}_0 +\hat{H}_I=\sum_{pq}\langle p |h_0|q\rangle a_{p}^{\dagger}a_{q} +\frac{1}{4}\sum_{pqrs}\langle pq| V|rs\rangle a_{p}^{\dagger}a_{q}^{\dagger}a_{s}a_{r},
\tag{1}
\end{equation}
$$
where \( a_{p}^{\dagger} \) and \( a_{q} \) etc are standard fermion creation and annihilation operators, respectively,
and \( pqrs \) represent all possible single-particle quantum numbers.
The full single-particle space is defined by the completeness relation
\( \hat{1} = \sum_{p =1}^{\infty}|p \rangle \langle p| \).
In our calculations we will let the single-particle states \( |p\rangle \)
be eigenfunctions of the one-particle operator \( \hat{h}_0 \).
The above Hamiltonian acts in turn on various many-body Slater determinants constructed from the single-basis defined by the one-body operator \( \hat{h}_0 \).
Our specific model consists of only \( 2 \) doubly-degenerate and equally spaced single-particle levels labeled by \( p=1,2,\dots \) and spin \( \sigma=\pm 1 \). In Eq. (1) the labels \( pqrs \) could also include spin \( \sigma \). From now and for the rest of this project, labels like \( pqrs \) represent the states without spin. The spin quantum numbers need to be accounted for explicitely.
We write the Hamiltonian as
$$
\begin{equation*} \hat{H} = \hat{H}_0 +\hat{H}_I=\hat{H}_0 + \hat{V} , \end{equation*}
$$
where
$$
\begin{equation*}
\hat{H}_0=\xi\sum_{p\sigma}(p-1)a_{p\sigma}^{\dagger}a_{p\sigma}.
\end{equation*}
$$
Here, \( H_0 \) is the unperturbed Hamiltonian with a spacing between
successive single-particle states given by \( \xi \), which we will set to
a constant value \( \xi=1 \) without loss of generality.
The two-body operator \( \hat{V} \) has one term only. It represents the pairing contribution and carries a constant strength \( g \) and is given by
$$
\begin{equation*}
\langle q+q-| V|s+s-\rangle = -g
\end{equation*}
$$
where \( g \) is a constant. The above labeling means that for a general matrix elements
\( \langle pq| V|rs\rangle \) we require that the states \( p \) and \( q \) (and \( r \) and \( s \)) have the same number
quantum number \( q \) but opposite spins. The two spins values are
\( \sigma = \pm 1 \).
When setting up the Hamiltonian matrix you need to figure out how to make the two-body interaction antisymmetric.
The variables \( \sigma=\pm \) represent the two possible spin values. The
interaction can only couple pairs and excites therefore only two
particles at the time.
In our model we have kept both the interaction strength and the single-particle level as constants. In a realistic system like the atomic nucleus this is not the case.
The unperturbed Hamiltonian \( \hat{H}_0 \) and \( \hat{V} \) commute with the spin projection \( \hat{S}_z \) and the total spin \( \hat{S}^2 \). This is an important feature of our system that allows us to block-diagonalize the full Hamiltonian. In this project we will focus only on total spin \( S=0 \), this case is normally called the no-broken pair case.
$$
\begin{equation*}
\hat{S}_z := \frac{1}{2}\sum_{p\sigma} \sigma a^{\dagger}_{p\sigma}a_{p\sigma}
\end{equation*}
$$
and
$$
\begin{equation*}
\hat{S}^2 := \hat{S}_z^2 + \frac{1}{2}(\hat{S}_+\hat{S}_- +
\hat{S}_-\hat{S}_+),
\end{equation*}
$$
where
$$
\begin{equation*}
\hat{S}_\pm := \sum_{p} a^{\dagger}_{p\pm} a_{p\mp}.
\end{equation*}
$$
This is an important feature of our system that allows us to block-diagonalize the full Hamiltonian. We will focus on total spin \( S=0 \). In this case, it is convenient to define the so-called pair creation and pair annihilation operators
$$
\begin{equation*}
\hat{P}^{+}_p = a^{\dagger}_{p+}a^{\dagger}_{p-},
\end{equation*}
$$
and
$$
\begin{equation*}
\hat{P}^{-}_p = a_{p-}a_{p+},
\end{equation*}
$$
respectively.
The Hamiltonian (with \( \xi=1 \)) we will use can be written as
$$
\begin{equation*}
\hat{H}=\sum_{p\sigma}(p-1)a_{p\sigma}^{\dagger}a_{p\sigma}
-g\sum_{pq}\hat{P}^{+}_p\hat{P}^{-}_q.
\end{equation*}
$$
Show that Hamiltonian commutes with the product of the pair creation and annihilation operators.
This model corresponds to a system with no broken pairs. This means that the Hamiltonian can only link two-particle states in so-called spin-reversed states.
Find the eigenvalues by diagonalizing the Hamiltonian matrix. Vary your results for selected values of \( g\in [-1,1] \) and comment your results.
Codes to diagonalize in C++ or Fortran can be provided. For Python, numpy contains eigenvalue solvers based on for example Householder's and Givens' algorithms. These are topics which can we discuss separately. The lecture slides contain a rather detailed recipe on how to construct a Slater determinant basis and how to set up the Hamiltonian matrix to diagonalize.
$$
\begin{equation*}
E_0= -\frac{g}{4}n\left(\Omega-n+2\right).
\end{equation*}
$$
Enlarge now your system to six and eight fermions and to \( p=6 \) and \( p=8 \) single-particle states, respectively. Run your program for a degenerate single-particle state with degeneracy \( \Omega \) and test against the exact result for the ground state. Introduce thereafter a finite single-particle spacing and study the results as you vary \( g \), as done in b) and c). Comment your results.
The way we will set up the Slater determinants here follows a simple odometric recipe. The way it is done in more professional codes, is to use bitwise manipulations. The latter is a possible extension/challenge for those interested.
Part two of our project consists of at developing your own shell-model code that can perform shell-model studies of the oxygen isotopes using standard effective interactions (provided by us) using as example the \( 1s0d \) shell as model space. You may also need to consider a bit representation and manipulation of Slater determinants and to implement the Lanczos algorithm. These details will be discussed during our lectures.
The aim of this project is to study the structure of selected low-lying states of the oxygen and fluorine isotopes towards their respective dripline. These chains of isotopes have been studied extensively during the last years, with many efforts toward the understanding of their dripline properties, involving studies of low-lying excited states and electromagnetic transitions. For the oxygen isotopes, \( {}^{24}\mbox{O} \) is the last particle-stable nucleus, and for the fluorine isotopes \( {}^{31}\mbox{F} \) is assumed to the last stable one. This part can also be used to benchmark your shell-model program from the second part.
The task here is to study these isotopic chains, extract excitation energies and selected observables and compare with available data. To achieve this you will need to use an effective interaction designed for the \( 1s0d \) shell first and then, for nuclei beyond \( A=24 \) you may need to consider degrees of freedom from the \( 1p0f \) shell. Since a full calculation in these two major shells becomes quickly time-consuming for the fluorine isotopes, you will need to truncate the number of particles which can leave/occupy selected single-particle states. In the file which contains the single-particle data, you can reduce the size of the total space of Slater determinants by limiting the number of particles which can populate the \( 1p0f \) shell. Here you could limit yourselves to consider only the single-particle states \( 0f_{7/2} \) and \( 1p_{3/2} \).
In the build-up of a shell model 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
$$
\begin{equation*}
\vert \alpha_1\dots \alpha_n\rangle = a_{\alpha_1}^\dagger a_{\alpha_2}^\dagger \dots a_{\alpha_n}^\dagger \vert 0\rangle.
\end{equation*}
$$
The time consuming part stems from the action of the Hamiltonian
on the above determinant,
$$
\begin{equation*}
\left(\sum_{\alpha\beta} \langle \alpha\vert \hat{t}+\hat{u}\vert \beta\rangle a_\alpha^\dagger a_\beta + \frac{1}{4} \sum_{\alpha\beta\gamma\delta}
\langle\alpha \beta\vert \hat{V}\vert \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 \vert 0\rangle.
\end{equation*}
$$
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
$$
\begin{equation*}
\Phi_{\Lambda} = a_{\alpha_3}^\dagger a_{\alpha_6}^\dagger a_{\alpha_{10}}^\dagger a_{\alpha_{13}}^\dagger \vert 0 \rangle .
\end{equation*}
$$
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{equation*}
\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}
\end{equation*}
$$
which translates into the decimal number
$$
\begin{equation*}
2^3+2^6+2^{10}+2^{13}=9288.
\end{equation*}
$$
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
$$
\begin{equation*}
\mathrm{dim}(\mathcal{H}) = \left(\begin{array}{c} n \\N\end{array}\right).
\end{equation*}
$$
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
$$
\begin{equation*}
\Phi_{\Lambda} = a_{\alpha_3}^\dagger a_{\alpha_6}^\dagger a_{\alpha_{10}}^\dagger a_{\alpha_{13}}^\dagger \vert 0 \rangle ,
\end{equation*}
$$
in a more compact way as
$$
\begin{equation*}
\Phi_{3,6,10,13} = |0001001000100100\rangle.
\end{equation*}
$$
The action of a creation operator is thus
$$
\begin{equation*}
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 \vert 0 \rangle ,
\end{equation*}
$$
which becomes
$$
\begin{equation*}
-a_{\alpha_3}^\dagger a^\dagger_{\alpha_4} a_{\alpha_6}^\dagger a_{\alpha_{10}}^\dagger a_{\alpha_{13}}^\dagger \vert 0 \rangle =-|0001101000100100\rangle.
\end{equation*}
$$
Similarly
$$
\begin{equation*}
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 \vert 0 \rangle ,
\end{equation*}
$$
which becomes
$$
\begin{equation*}
-a^\dagger_{\alpha_4} (a_{\alpha_6}^\dagger)^ 2 a_{\alpha_{10}}^\dagger a_{\alpha_{13}}^\dagger \vert 0 \rangle =0!
\end{equation*}
$$
This gives a simple recipe:
Consider the action of \( a^\dagger_{\alpha_2} \) on various slater determinants:
$$
\begin{equation*}
\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}
\end{equation*}
$$
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
$$
\begin{equation*}
\Phi_{\Lambda} = a_{\alpha_0}^\dagger a_{\alpha_3}^\dagger a_{\alpha_6}^\dagger a_{\alpha_{10}}^\dagger a_{\alpha_{13}}^\dagger \vert 0 \rangle ,
\end{equation*}
$$
in a more compact way as
$$
\begin{equation*}
\Phi_{0,3,6,10,13} = |1001001000100100\rangle.
\end{equation*}
$$
The action of
$$
\begin{equation*}
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 \vert 0 \rangle ,
\end{equation*}
$$
which becomes
$$
\begin{equation*}
-a_{\alpha_3}^\dagger a^\dagger_{\alpha_4} a_{\alpha_6}^\dagger a_{\alpha_{10}}^\dagger a_{\alpha_{13}}^\dagger \vert 0 \rangle =-|0001101000100100\rangle.
\end{equation*}
$$
The action
$$
\begin{equation*}
a_{\alpha_0}\Phi_{0,3,6,10,13} = |0001001000100100\rangle,
\end{equation*}
$$
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
$$
\begin{equation*}
|1000000000000000\rangle,
\end{equation*}
$$
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
$$
\begin{equation*}
S=|0111000000000000\rangle
\end{equation*}
$$
$$
\begin{equation*}
\Phi_{0,3,6,10,13} = |1001001000100100\rangle,
\end{equation*}
$$
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.