Introduction

 
$$ \mid\Psi\rangle = \exp\left(-\hat{T}\right)\mid\Phi\rangle\ . $$

 
We will discuss how to define the operators later in this work. This simple ansatz carries enormous power. It leads to a non-perturbative many-body theory that includes summation of ladder diagrams , ring diagrams, and an infinite-order generalization of many-body perturbation theory.

Introduction

A non-practical way of solving the eigenvalue problem

Before we proceed with the derivation of the Coupled cluster equations, let us repeat some of the arguments we presented during our FCI lectures. In our FCI discussions, we rewrote the solution of the Schroedinger equation as a set of coupled equationsin the unknown coefficients \( C \). Let us repeat some of these arguments. To obtain the eigenstates and eigenvalues in terms of non-linear equations is not a very practical approach. However, it serves the scope of linking FCI theory with approximative solutions to the many-body problem like Coupled cluster (CC) theory

A non-practical way of solving the eigenvalue problem

 
$$ \langle \Phi_0 | \hat{H} -E| \Phi_0\rangle + \sum_{ai}\langle \Phi_0 | \hat{H} -E|\Phi_{i}^{a} \rangle C_{i}^{a}+ \sum_{abij}\langle \Phi_0 | \hat{H} -E|\Phi_{ij}^{ab} \rangle C_{ij}^{ab}=0, $$

 
or

 
$$ E-E_0 =\Delta E=\sum_{ai}\langle \Phi_0 | \hat{H}|\Phi_{i}^{a} \rangle C_{i}^{a}+ \sum_{abij}\langle \Phi_0 | \hat{H}|\Phi_{ij}^{ab} \rangle C_{ij}^{ab}, $$

 
where the energy \( E_0 \) is the reference energy and \( \Delta E \) defines the so-called correlation energy. The single-particle basis functions could be the results of a Hartree-Fock calculation or just the eigenstates of the non-interacting part of the Hamiltonian.

A non-practical way of solving the eigenvalue problem

 
$$ E-E_0 =\Delta E^{HF}=\sum_{abij}\langle \Phi_0 | \hat{H}|\Phi_{ij}^{ab} \rangle C_{ij}^{ab}. $$

A non-practical way of solving the eigenvalue problem

 
$$ \Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle C_{i}^{a}+ \sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}. $$

 
This equation determines the correlation energy but not the coefficients \( C \). We need more equations. Our next step is to set up

 
$$ \langle \Phi_i^a | \hat{H} -E| \Phi_0\rangle + \sum_{bj}\langle \Phi_i^a | \hat{H} -E|\Phi_{j}^{b} \rangle C_{j}^{b}+ \sum_{bcjk}\langle \Phi_i^a | \hat{H} -E|\Phi_{jk}^{bc} \rangle C_{jk}^{bc}+ \sum_{bcdjkl}\langle \Phi_i^a | \hat{H} -E|\Phi_{jkl}^{bcd} \rangle C_{jkl}^{bcd}=0, $$

 
as this equation will allow us to find an expression for the coefficents \( C_i^a \) since we can rewrite this equation as

 
$$ \langle i | \hat{f}| a\rangle +\langle \Phi_i^a | \hat{H}|\Phi_{i}^{a} \rangle C_{i}^{a}+ \sum_{bj\ne ai}\langle \Phi_i^a | \hat{H}|\Phi_{j}^{b} \rangle C_{j}^{b}+ \sum_{bcjk}\langle \Phi_i^a | \hat{H}|\Phi_{jk}^{bc} \rangle C_{jk}^{bc}+ \sum_{bcdjkl}\langle \Phi_i^a | \hat{H}|\Phi_{jkl}^{bcd} \rangle C_{jkl}^{bcd}=EC_i^a. $$

A non-practical way of solving the eigenvalue problem

 
$$ C_{i}^{a}=\frac{\langle i | \hat{f}| a\rangle}{\epsilon_i-\epsilon_a}. $$

The observant reader will however see that we need an equation for \( C_{jk}^{bc} \) and \( C_{jkl}^{bcd} \) as well. To find equations for these coefficients we need then to continue our multiplications from the left with the various \( \Phi_{H}^P \) terms.

A non-practical way of solving the eigenvalue problem

 
$$ \langle \Phi_{ij}^{ab} | \hat{H} -E| \Phi_0\rangle + \sum_{kc}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{k}^{c} \rangle C_{k}^{c}+ $$

 
$$ \sum_{cdkl}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{kl}^{cd} \rangle C_{kl}^{cd}+\sum_{cdeklm}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{klm}^{cde} \rangle C_{klm}^{cde}+\sum_{cdefklmn}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{klmn}^{cdef} \rangle C_{klmn}^{cdef}=0, $$

 
and we can isolate the coefficients \( C_{kl}^{cd} \) in a similar way as we did for the coefficients \( C_{i}^{a} \). A standard choice for the first iteration is to set

 
$$ C_{ij}^{ab} =\frac{\langle ij \vert \hat{v} \vert ab \rangle}{\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b}. $$

A non-practical way of solving the eigenvalue problem

Summarizing FCI and bringing in approximative methods

• It gives all eigenvalues, ground state and excited states
• The eigenvectors are obtained directly from the coefficients \( C_H^P \) which result from the diagonalization
• We can compute easily expectation values of other operators, as well as transition probabilities
• Correlations are easy to understand in terms of contributions to a given operator beyond the Hartree-Fock contribution. This is the standard approach in many-body theory.

Summarizing FCI and bringing in approximative methods

 
$$ \Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle C_{i}^{a}+ \sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}. $$

The coefficients \( C \) result from the solution of the eigenvalue problem. The energy of say the ground state is then

 
$$ E=E_{ref}+\Delta E, $$

 
where the so-called reference energy is the energy we obtain from a Hartree-Fock calculation, that is

 
$$ E_{ref}=\langle \Phi_0 \vert \hat{H} \vert \Phi_0 \rangle. $$

Summarizing FCI and bringing in approximative methods

Popular methods are

• Many-body perturbation theory (in essence a Taylor expansion)
• Coupled cluster theory (coupled non-linear equations)
• Green's function approaches (matrix inversion)
• Similarity group transformation methods (coupled ordinary differential equations)

All these methods start normally with a Hartree-Fock basis as the calculational basis.

A quick tour of Coupled Cluster theory

The ansatz for the wavefunction (ground state) is given by

 
$$ \begin{equation*} \vert \Psi\rangle = \vert \Psi_{CC}\rangle = e^{\hat{T}} \vert \Phi_0\rangle = \left( \sum_{n=1}^{A} \frac{1}{n!} \hat{T}^n \right) \vert \Phi_0\rangle, \end{equation*} $$

 
where \( A \) represents the maximum number of particle-hole excitations and \( \hat{T} \) is the cluster operator defined as

 
$$ \begin{align*} \hat{T} &= \hat{T}_1 + \hat{T}_2 + \ldots + \hat{T}_A \\ \hat{T}_n &= \left(\frac{1}{n!}\right)^2 \sum_{\substack{ i_1,i_2,\ldots i_n \\ a_1,a_2,\ldots a_n}} t_{i_1i_2\ldots i_n}^{a_1a_2\ldots a_n} a_{a_1}^\dagger a_{a_2}^\dagger \ldots a_{a_n}^\dagger a_{i_n} \ldots a_{i_2} a_{i_1}. \end{align*} $$

A quick tour of Coupled Cluster theory

 
$$ \begin{equation*} E_{\mathrm{CC}} = \langle\Phi_0\vert \overline{H}\vert \Phi_0\rangle, \end{equation*} $$

 
where \( \overline{H} \) is a similarity transformed Hamiltonian

 
$$ \begin{align*} \overline{H}&= e^{-\hat{T}} \hat{H}_N e^{\hat{T}} \\ \hat{H}_N &= \hat{H} - \langle\Phi_0\vert \hat{H} \vert \Phi_0\rangle. \end{align*} $$

A quick tour of Coupled Cluster theory

 
$$ \begin{equation*} 0 = \langle\Phi_{i_1 \ldots i_n}^{a_1 \ldots a_n}\vert \overline{H}\vert \Phi_0\rangle. \end{equation*} $$

 
The similarity transformed Hamiltonian \( \overline{H} \) is expanded using the Baker-Campbell-Hausdorff expression,

 
$$ \begin{align*} \overline{H}&= \hat{H}_N + \left[ \hat{H}_N, \hat{T} \right] + \frac{1}{2} \left[\left[ \hat{H}_N, \hat{T} \right], \hat{T}\right] + \ldots \\ & \quad \frac{1}{n!} \left[ \ldots \left[ \hat{H}_N, \hat{T} \right], \ldots \hat{T} \right] +\dots \end{align*} $$

 
and simplified using the connected cluster theorem

 
$$ \begin{equation*} \overline{H}= \hat{H}_N + \left( \hat{H}_N \hat{T}\right)_c + \frac{1}{2} \left( \hat{H}_N \hat{T}^2\right)_c + \dots + \frac{1}{n!} \left( \hat{H}_N \hat{T}^n\right)_c +\dots \end{equation*} $$

A quick tour of Coupled Cluster theory

The coupled cluster wavefunction is now given by

 
$$ \begin{equation*} \vert \Psi_{CC}\rangle = e^{\hat{T}_1 + \hat{T}_2} \vert \Phi_0\rangle \end{equation*} $$

 
where

 
$$ \begin{align*} \hat{T}_1 &= \sum_{ia} t_{i}^{a} a_{a}^\dagger a_i \\ \hat{T}_2 &= \frac{1}{4} \sum_{ijab} t_{ij}^{ab} a_{a}^\dagger a_{b}^\dagger a_{j} a_{i}. \end{align*} $$

A quick tour of Coupled Cluster theory

If we truncate our equations at the CCSD level, it corresponds to performing a transformation of the Hamiltonian matrix of the following type for a six particle problem (with a two-body Hamiltonian):

A quick tour of Coupled Cluster theory

 
$$ \Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle C_{i}^{a}+ \sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}. $$

In Coupled cluster theory it becomes (irrespective of level of truncation of \( T \))

 
$$ \Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle t_{i}^{a}+ \sum_{abij}\langle ij | \hat{v}| ab \rangle t_{ij}^{ab}. $$

A quick tour of Coupled Cluster theory

There are several interesting features:

• With a truncation like CCSD or CCSDT, we can include to infinite order correlations like \( 2p-2h \).
• We can include a large basis of single-particle states, not possible in standard FCI calculations

However, Coupled Cluster theory is

• non-variational
• if we want to find properties of excited states, additional calculations via for example equation of motion methods are needed
• if correlations are strong, a single-reference ansatz may not be the best starting point
• we cannot quantify properly the error we make when truncations are made in the cluster operator

The CCD approximation

We will now approximate the cluster operator \( \hat{T} \) to include only \( 2p-2h \) correlations. This leads to the so-called CCD approximation, that is

 
$$ \hat{T}\approx \hat{T}_2=\frac{1}{4}\sum_{abij}t_{ij}^{ab}a^{\dagger}_aa^{\dagger}_ba_ja_i, $$

 
meaning that we have

 
$$ \vert \Psi_0 \rangle \approx \vert \Psi_{CCD} \rangle = \exp{\left(\hat{T}_2\right)}\vert \Phi_0\rangle. $$

The CCD approximation

 
$$ \hat{H}=\hat{H}_N+E_{\mathrm{ref}}, $$

 
with

 
$$ \hat{H}_N=\sum_{pq}\langle p \vert \hat{f} \vert q \rangle a^{\dagger}_pa_q + \frac{1}{4}\sum_{pqrs}\langle pq \vert \hat{v} \vert rs \rangle a^{\dagger}_pa^{\dagger}_qa_sa_r, $$

 
we obtain that the energy can be written as

 
$$ \langle \Phi_0 \vert \exp{-\left(\hat{T}_2\right)}\hat{H}_N\exp{\left(\hat{T}_2\right)}\vert \Phi_0\rangle = \langle \Phi_0 \vert \hat{H}_N(1+\hat{T}_2)\vert \Phi_0\rangle = E_{CCD}. $$

The CCD approximation

 
$$ E_{CCD}=E_{\mathrm{ref}}+\frac{1}{4}\sum_{abij}\langle ij \vert \hat{v} \vert ab \rangle t_{ij}^{ab}, $$

 
where the latter is the correlation energy from this level of approximation of CC theory. Similarly, the expression for the amplitudes reads

 
$$ \langle \Phi_{ij}^{ab} \vert \exp{-\left(\hat{T}_2\right)}\hat{H}_N\exp{\left(\hat{T}_2\right)}\vert \Phi_0\rangle = 0. $$

The CCD approximation

 
$$ \begin{align} 0 = \langle ab \vert \hat{v} \vert ij \rangle + \left(\epsilon_a+\epsilon_b-\epsilon_i-\epsilon_j\right)t_{ij}^{ab} & \nonumber \\ +\frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle t_{ij}^{cd}+\frac{1}{2}\sum_{kl} \langle kl \vert \hat{v} \vert ij \rangle t_{kl}^{ab}+\hat{P}(ij\vert ab)\sum_{kc} \langle kb \vert \hat{v} \vert cj \rangle t_{ik}^{ac} & \nonumber \\ +\frac{1}{4}\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ij}^{cd}t_{kl}^{ab}+\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{ac}t_{jl}^{bd}& \nonumber \\ -\frac{1}{2}\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{dc}t_{lj}^{ab}-\frac{1}{2}\hat{P}(ab)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{lk}^{ac}t_{ij}^{db},& \tag{1} \end{align} $$

 
where we have defined

 
$$ \hat{P}\left(ab\right)= 1-\hat{P}_{ab}, $$

 
where \( \hat{P}_{ab} \) interchanges two particles occupying the quantum numbers \( a \) and \( b \).

The CCD approximation

 
$$ \hat{P}(ij\vert ab) = (1-\hat{P}_{ij})(1-\hat{P}_{ab}). $$

 
Recall also that the unknown amplitudes \( t_{ij}^{ab} \) represent anti-symmetrized matrix elements, meaning that they obey the same symmetry relations as the two-body interaction, that is

 
$$ t_{ij}^{ab}=-t_{ji}^{ab}=-t_{ij}^{ba}=t_{ji}^{ba}. $$

 
The two-body matrix elements are also anti-symmetrized, meaning that

 
$$ \langle ab \vert \hat{v} \vert ij \rangle = -\langle ab \vert \hat{v} \vert ji \rangle= -\langle ba \vert \hat{v} \vert ij \rangle=\langle ba \vert \hat{v} \vert ji \rangle. $$

 
The non-linear equations for the unknown amplitudes \( t_{ij}^{ab} \) are solved iteratively. We discuss the implementation of these equations below.

Approximations to the full CCD equations

It is useful to make approximations to the equations for the amplitudes. The standard method for solving these equations is to set up an iterative scheme where method's like Newton's method or similar root searching methods are used to find the amplitudes. Itreative solvers need a guess for the amplitudes. A good starting point is to use the correlated wave operator from perturbation theory to first order in the interaction. This means that we define the zeroth approximation to the amplitudes as

 
$$ t^{(0)}=\frac{\langle ab \vert \hat{v} \vert ij \rangle}{\left(\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b\right)}, $$

 
leading to our first approximation for the correlation energy at the CCD level to be equal to second-order perturbation theory without \( 1p-1h \) excitations, namely

 
$$ \Delta E_{\mathrm{CCD}}^{(0)}=\frac{1}{4}\sum_{abij} \frac{\langle ij \vert \hat{v} \vert ab \rangle \langle ab \vert \hat{v} \vert ij \rangle}{\left(\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b\right)}. $$

Approximations to the full CCD equations

 
$$ \begin{equation} 0 = \langle ab \vert \hat{v} \vert ij \rangle + \left(\epsilon_a+\epsilon_b-\epsilon_i-\epsilon_j\right)t_{ij}^{ab}+\frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle t_{ij}^{cd}. \tag{2} \end{equation} $$

Approximations to the full CCD equations

 
$$ \tau_{ij}^{ab}= \left(\omega-\epsilon_a-\epsilon_b\right)t_{ij}^{ab}, $$

 
and inserting

 
$$ 1=\frac{\left(\omega-\epsilon_c-\epsilon_d\right)}{\left(\omega-\epsilon_c-\epsilon_d\right)}, $$

 
in the intermediate sums over \( cd \) in Eq. (2), we can rewrite the latter equation as

 
$$ \tau_{ij}^{ab}(\omega)= \langle ab \vert \hat{v} \vert ij \rangle + \frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle \frac{1}{\omega-\epsilon_c-\epsilon_d}\tau_{ij}^{cd}(\omega), $$

 
where we have indicated an explicit energy dependence. This equation, transforming a two-particle configuration into a single index, can be transformed into a matrix inversion problem. Solving the equations for a fixed energy \( \omega \) allows us to compare directly with results from Green's function theory when only two-particle intermediate states are included.

Approximations to the full CCD equations

1. Start with a guess for the amplitudes and compute the zeroth approximation to the correlation energy
2. Use the ansatz for the amplitudes to solve Eq. (2) via for example your root-finding method of choice (Newton's method or modifications thereof can be used) and continue these iterations till the correlation energy does not change more than a prefixed quantity \( \lambda \); \( \Delta E_{\mathrm{CCD}}^{(i)}-\Delta E_{\mathrm{CCD}}^{(i-1)} \le \lambda \).
3. It is common during the iterations to scale the amplitudes with a parameter \( \alpha \), with \( \alpha \in (0,1] \) as \( t^{(i)}=\alpha t^{(i)}+(1-\alpha)t^{(i-1)} \).

Approximations to the full CCD equations

 
$$ \begin{equation} 0 = \langle ab \vert \hat{v} \vert ij \rangle + \left(\epsilon_a+\epsilon_b-\epsilon_i-\epsilon_j\right)t_{ij}^{ab}+\frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle t_{ij}^{cd}+\frac{1}{2}\sum_{kl} \langle kl \vert \hat{v} \vert ij \rangle t_{kl}^{ab}. \tag{3} \end{equation} $$

 
This equation is solved the same way as we would do for Eq. (2). The final step is then to include all terms in Eq. (1).

Introduction to studies of infinite matter

Studies of infinite nuclear matter play an important role in nuclear physics. The aim of this part of the lectures is to provide the necessary ingredients for perfoming studies of neutron star matter (or matter in \( \beta \)-equilibrium) and symmetric nuclear matter.

Here we will study infinite neutron matter

• at the Hartree-Fock with realistic nuclear forces and
• using many-body methods like coupled-cluster theory or many-body perturbation theory

Infinite nuclear matter and neutron star matter

Studies of dense baryonic matter are of central importance to our basic understanding of the stability of nuclear matter, spanning from matter at high densities and temperatures to matter as found within dense astronomical objects like neutron stars.

Neutron star matter at densities of 0.1 fm$^{-3}$ and greater, is often assumed to be made of mainly neutrons, protons, electrons and muons in beta equilibrium. However, other baryons like various hyperons may exist, as well as possible mesonic condensates and transitions to quark degrees of freedom at higher densities. Here we focus on specific definitions of various phases and focus on distinct phases of matter such as pure baryonic matter and/or quark matter. The composition of matter is then determined by the requirements of chemical and electrical equilibrium. Furthermore, we will also consider matter at temperatures much lower than the typical Fermi energies.

Properties of infinite nuclear matter

The equilibrium conditions are governed by the weak processes (normally referred to as the processes for \( \beta \)-equilibrium)

 
$$ \begin{equation} b_1 \rightarrow b_2 + l +\bar{\nu}_l \hspace{1cm} b_2 +l \rightarrow b_1 +\nu_l, \tag{4} \end{equation} $$

 
where \( b_1 \) and \( b_2 \) refer to for example the baryons being a neutron and a proton, respectively, \( l \) is either an electron or a muon and \( \bar{\nu}_l \) and \( \nu_l \) their respective anti-neutrinos and neutrinos. Muons typically appear at a density close to nuclear matter saturation density, the latter being

 
$$ n_0 \approx 0.16 \pm 0.02 \hspace{1cm} \mathrm{fm}^{-3}, $$

 
with a corresponding binding energy \( E_0 \) for symmetric nuclear matter (SNM) at saturation density of

 
$$ E_0 = B/A=-15.6\pm 0.2 \hspace{1cm} \mathrm{MeV}. $$

The infinite neutron gas as a homogenous system

This is a homogeneous system and the one-particle wave functions are given by plane wave functions normalized to a volume \( \Omega \) for a box with length \( L \) (the limit \( L\rightarrow \infty \) is to be taken after we have computed various expectation values)

 
$$ \psi_{\mathbf{k}\sigma}(\mathbf{r})= \frac{1}{\sqrt{\Omega}}\exp{(i\mathbf{kr})}\xi_{\sigma} $$

 
where \( \mathbf{k} \) is the wave number and \( \xi_{\sigma} \) is a spin function for either spin up or down

 
$$ \xi_{\sigma=+1/2}=\left(\begin{array}{c} 1 \\ 0 \end{array}\right) \hspace{0.5cm} \xi_{\sigma=-1/2}=\left(\begin{array}{c} 0 \\ 1 \end{array}\right). $$

Periodic boundary conditions and single-particle states

When using periodic boundary conditions, the discrete-momentum single-particle basis functions

 
$$ \phi_{\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot \mathbf{r}}/L^{d/2} $$

 
are associated with the single-particle energy

 
$$ \begin{align} \varepsilon_{n_{x}, n_{y}} = \frac{\hbar^{2}}{2m} \left( \frac{2\pi }{L}\right)^{2}\left( n_{x}^{2} + n_{y}^{2}\right) \tag{5} \end{align} $$

 
for two-dimensional sytems and

 
$$ \begin{align} \varepsilon_{n_{x}, n_{y}, n_{z}} = \frac{\hbar^{2}}{2m} \left( \frac{2\pi }{L}\right)^{2} \left( n_{x}^{2} + n_{y}^{2} + n_{z}^{2}\right) \tag{6} \end{align} $$

 
for three-dimensional systems.

More on periodic boundary conditions and single-particle states

The table on the next slide illustrates how single-particle energies fill energy shells in a two-dimensional neutron box. Here \( n_{x} \) and \( n_{y} \) are the momentum quantum numbers, \( n_{x}^{2} + n_{y}^{2} \) determines the single-particle energy level, \( N_{\uparrow \downarrow } \) represents the cumulated number of spin-orbitals in an unpolarized spin phase, and \( N_{\uparrow \uparrow } \) stands for the cumulated number of spin-orbitals in a spin-polarized system.

Magic numbers for the two-dimensional neutron (or electron) gas

Three-dimensional neutron gas

 
$$ \frac{\hbar^2}{2m}\left(k_{n_x}^2+k_{n_y}^2k_{n_z}^2\right), $$

 
and

 
$$ k_{n_i}=\frac{2\pi n_i}{L} \hspace{0.1cm} n_i = 0, \pm 1, \pm 2, \dots, $$

 
we can set up a similar table and obtain (assuming identical particles one and including spin up and spin down solutions) for energies less than or equal to \( n_{x}^{2}+n_{y}^{2}+n_{z}^{2}\le 3 \)

Single-particle states for the three-dimensional neutron gas

Continuing in this way we get for \( n_{x}^{2}+n_{y}^{2}+n_{z}^{2}=4 \) a total of 12 additional states, resulting in \( ? \) as a new magic number. We can continue like this by adding more shells.

When performing calculations based on many-body perturbation theory, Coupled cluster theory or other many-body methods, we need then to add states above the Fermi level in order to sum over single-particle states which are not occupied.

Input parameters

 
$$ \rho = g \frac{k_F^3}{6\pi^2}, $$

 
where you need to define \( k_F \) and the degeneracy \( g \), which is two for one type of spin-\( 1/2 \) particles and four for symmetric nuclear matter.

With the density we can define the length \( L \) of the box used with periodic boundary contributions, that is use the relation

 
$$ V= L^3= \frac{A}{\rho}. $$

 
Finally we can use \( L \) to define the spacing to set up the spacing between varipus \( k \)-values, that is

 
$$ \Delta k = \frac{2\pi}{L}. $$

 
Here, \( A \) can be the number of nucleons.

Potential model employed in code development

The interaction we will use for these calculations is a semirealistic nucleon-nucleon potential known as the Minnesota potential

 
$$ V_{\alpha}\left(r\right)=V_{\alpha}\exp{(-\alpha r^{2})}. $$

 
The spin and isospin dependence of the Minnesota potential is given by

 
$$ \begin{equation} V\left( r\right)=\frac{1}{2}\left( V_{R}+\frac{1}{2}\left(1+P_{12}^{\sigma}\right) V_{T}+\frac{1}{2}\left(1-P_{12}^{\sigma}\right) V_{S}\right)\left(1-P_{12}^{\sigma}P_{12}^{\tau}\right), \tag{7} \end{equation} $$

 
where

 
$$ P_{12}^{\sigma}=\frac{1}{2}\left(1+\sigma_{1}\cdot\sigma_{2}\right), $$

 
and

 
$$ P_{12}^{\tau}=\frac{1}{2}\left( 1+\tau_{1}\cdot\tau_{2}\right) $$

 
are the spin and isospin exchange operators, respectively.

Fourier transform

 
$$ V_{\alpha}\left(r\right) $$

 
is rather simple since the radial depends only on the magnitude of the relative distance and thereby the relative momentum

 
$$ \mathbf{q}=\frac{1}{2}\left(\mathbf{k}_{p}-\mathbf{k}_{q}-\mathbf{k}_{r}+\mathbf{k}_{s}\right) $$

 
Omitting spin and isospin dependencies, the momentum space version of the interaction reads

 
$$ \begin{equation} \langle \mathbf{k}_p \mathbf{k}_q \vert V_{\alpha}\vert\mathbf{k}_r\mathbf{k}_s\rangle= \frac{V_{\alpha}}{L^{3}}\left(\frac{\pi}{\alpha}\right)^{3/2}\exp{(\frac{-q^{2}}{4\alpha})}\delta_{\mathbf{k}_{p}+\mathbf{k}_{q},\mathbf{k}_{r}+\mathbf{k}_{s}} \tag{8} \end{equation} $$

Developing a program for infinite matter

• Structure a code in terms of functions.
• Modularize your codes.
• Be able to read input data flexibly.
• Write unit tests (test functions) and let your code undergo heavy testing.
• Refactor code in terms of classes (instead of functions only).
• Conduct and automate large-scale numerical experiments.
• New code is added in a modular fashion to a library (modules).
• Programs are run through convenient user interfaces.
• Use scripts in order to automatize tedious manual work.
• Make sure your scientific investigations are reproducible and document properly your results.
• Use version control software like for example git

Codes and reading material

• The codes are all available from the coming Lectures Notes in Physics volume An advanced course in computational nuclear physics: Bridging the scales from quarks to neutron stars, M. Hjorth-Jensen, M. P. Lombardo, U. van Kolck, Editors
• For more reading on CCD see chapter 8 of the coming Lectures Notes in Physics volume An advanced course in computational nuclear physics: Bridging the scales from quarks to neutron stars, M. Hjorth-Jensen, M. P. Lombardo, U. van Kolck, Editors

The CCD equation

 
$$ \begin{align} \left(\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b\right)t_{ij}^{ab} = \langle ab \vert \hat{v} \vert ij \rangle & \nonumber \\ +\frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle t_{ij}^{cd}+\frac{1}{2}\sum_{kl} \langle kl \vert \hat{v} \vert ij \rangle t_{kl}^{ab}+\hat{P}(ij\vert ab)\sum_{kc} \langle kb \vert \hat{v} \vert cj \rangle t_{ik}^{ac} & \nonumber \\ +\frac{1}{4}\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ij}^{cd}t_{kl}^{ab}+\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{ac}t_{jl}^{bd}& \nonumber \\ -\frac{1}{2}\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{dc}t_{lj}^{ab}-\frac{1}{2}\hat{P}(ab)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{lk}^{ac}t_{ij}^{db},& \tag{9} \end{align} $$

 
for all \( i < j \) and all \( a < b \), using the standard notation that \( a,b,... \) are particle states and \( i,j,... \) are hole states. With the CCD correlation energy given by

 
$$ \begin{equation} \Delta E_{CCD} = \frac{1}{4} \sum_{ijab} \langle ij \vert\hat{v}\vert ab\rangle t^{ab}_{ij}. \tag{10} \end{equation} $$

Solving the CCD equations

 
$$ \begin{align} t_{ij}^{ab}{}^{(n+1)} = \frac{1}{\epsilon^{ab}_{ij}} \bigg(\langle ab \vert \hat{v} \vert ij \rangle & \nonumber \\ +\frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle t_{ij}^{cd}{}^{(n)}+\frac{1}{2}\sum_{kl} \langle kl \vert \hat{v} \vert ij \rangle t_{kl}^{ab}{}^{(n)}+\hat{P}(ij\vert ab)\sum_{kc} \langle kb \vert \hat{v} \vert cj \rangle t_{ik}^{ac}{}^{(n)} & \nonumber \\ +\frac{1}{4}\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ij}^{cd}{}^{(n)}t_{kl}^{ab}{}^{(n)}+\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{ac}{}^{(n)}t_{jl}^{bd}{}^{(n)}& \nonumber \\ -\frac{1}{2}\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{dc}{}^{(n)}t_{lj}^{ab}{}^{(n)}-\frac{1}{2}\hat{P}(ab)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{lk}^{ac}{}^{(n)}t_{ij}^{db}{}^{(n)} \bigg),& \tag{11} \end{align} $$

 
for all \( i < j \) and all \( a < b \), where \( \epsilon^{ab}_{ij} =\left(\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b\right) \), and \( t_{ij}^{ab}{}^{(n)} \) is the \( t \) amplitude for the nth iteration of the series. This way, given some starting guess \( t_{ij}^{ab}{}^{(0)} \), we can generate subsequent \( t \) amplitudes that converges to some value.

Memory considerations

Care should thus be placed into how we store these objects. These are objects with four indices and a sensible first implementation of the CCD equations would be to create two four-dimensional arrays to store the objects. However, it is often more convenient to work with simple one-dimensional arrays instead.

The goal of our code is to calculate the correlation energy, \( \Delta E_{CCD} \), meaning that after each iteration of our equations, we use our newest \( t \) amplitudes to update the correlation energy

 
$$ \begin{equation} \Delta E_{CCD}^{(n)} = \frac{1}{4} \sum_{ijab} \langle ij\vert\hat{v}\vert ab\rangle t^{ab}_{ij}{}^{(n)}. \tag{12} \end{equation} $$

 
We check that our result is converged by testing whether the most recent iteration has changed the correlation energy by less than some tolerance threshold \( \eta \),

 
$$ \begin{equation} \eta > | \Delta E_{CCD}^{(n+1)} - \Delta E_{CCD}^{(n)} |. \tag{13} \end{equation} $$

More on memory

One limitation that will be ran into while trying to do realistic CCD calculations is that of memory. The four-indexed two-body matrix elements (TBMEs) and \( t \)-amplitudes have to store a lot of elements, and the size of these arrays can quickly exceed the available memory on a machine. If a calculation wants to use 500 single-particle basis states, then a structure like \( \langle pq\vert v\vert rs\rangle \) will need a length of 500 for each of its four indices, which means it will have \( 500^4 = 625\times 10^8 \) elements. To get a handle on how much memory is used, consider the elements as double-precision floating point type. One double takes up 8 bytes of memory. This specific array would take up \( 8\times 625\times 10^8 \) bytes = \( 5000 \times 10^8 \) bytes = \( 500 \) Gbytes of memory.

Using symmetries

Most personal computers in 2016 have 4-8 Gbytes of RAM, meaning that this calculation would be way out of reach. There are supercomputers that can handle applications using 500 Gbytes of memory, but we can quickly reduce the total memory required by applying some physical arguments. In addition to vanishing elements with repeated indices, mentioned above, elements that do not obey certain symmetries are also zero. Almost all realistic two-body forces preserve some quantities due to symmetries in the interaction. For example, an interaction with rotational symmetry will conserve angular momentum. This means that a two-body ket state \( \vert rs\rangle \), which has some set of quantum numbers, will retain quantum numbers corresponding to the interaction symmetries after being acted on by \( \hat{v} \). This state is then projected onto \( \vert pq\rangle \) with its own set of quantum numbers. Thus \( \langle pq|v|rs\rangle \) is only non-zero if \( \vert pq\rangle \) and \( \vert rs\rangle \) share the same quantum numbers that are preserved by \( \hat{v} \). In addition, because the cluster operators represent excitations due to the interaction, \( t_{ij}^{ab} \) is only non-zero if \( \vert ij\rangle \) has the same relevant quantum numbers as \( \vert ab\rangle \).

Using symmetries II

To take advantage of this, these two-body ket states can be organized into "channels" of shared quantum numbers. In the case of the pairing model, the interaction preserves the total spin projection of a two-body state, \( S_{z}=s_{z1}+s_{z2} \). The single particle states can have spin of +1/2 or -1/2, so there can be three two-body channels with \( S_{z}=-1,0,+1 \). These channels can then be indexed with a unique label in a similar way to the single particle index scheme. In more complicated systems, there will be many more channels involving multiple symmetries, so it is useful to create a data structure that stores the relevant two-body quantum numbers to keep track of the labeling scheme.

Using symmetries III

It is more efficient to use two-dimensional array data structures, where the first index refers to the channel number and the second refers to the element within that channel. So to access matrix elements or \( t \) amplitudes, you can loop over the channels first, then the indices within that channel. To get an idea of the savings using this block diagonal structure, let's look at a case with a plane wave basis, with three momentum and one spin quantum numbers, with an interaction that conserves linear momentum in all three dimensions, as well as the total spin projection. Using 502 basis states, the TBME's require about 0.23 Gb of memory in block diagonal form, which is an enormous saving from the 500 Gb mentioned earlier in the na\"ive storage scheme.

Using intermediates

Since the calculation of all zeros can now be avoided, improvements in speed and memory will now follow. To get a handle on how these CCD calculations are implemented we need only to look at the most expensive sum in equation (3). This corresponds to the sum over \( klcd \). Since this sum is repeated for all \( i < j \) and \( a < b \), it means that these equations will scale as \( \mathcal{O}(n_{p}^{4} n_{h}^{4}) \). However, they can be rewritten using intermediates as

 
$$ \begin{align} 0 = \langle ab|\hat{v}|ij \rangle + \hat{P}(ab) \sum_{c} \langle b| \chi |c\rangle \langle ac| t |ij\rangle - \hat{P}(ij) \sum_{k} \langle k| \chi |j\rangle \langle ab| t |ik\rangle & \nonumber \\ + \frac{1}{2}\sum_{cd} \langle ab| \chi |cd\rangle \langle cd| t |ij\rangle + \frac{1}{2} \sum_{kl} \langle ab| t |kl\rangle \langle kl| \chi |ij\rangle \tag{14}\\ + \hat{P}(ij)\hat{P}(ab) \sum_{kc} \langle ac| t |ik\rangle\langle kb| \chi |cj\rangle & \nonumber \end{align} $$

 
for all \( i,j,a,b \).

Defining intermediates

The intermediates \( \chi \) are defined as

 
$$ \begin{equation} \langle b| \chi |c \rangle =\langle b|f|c\rangle - \frac{1}{2} \sum_{kld} \langle bd|t|kl\rangle \langle kl|v|cd\rangle \tag{15} \end{equation} $$

 
$$ \begin{equation} \langle k| \chi |j\rangle = \langle k|f|j\rangle + \frac{1}{2} \sum_{cdl} \langle kl|v|cd\rangle \langle cd|t|jl\rangle \tag{16} \end{equation} $$

 
$$ \begin{equation} \langle kl| \chi |ij\rangle = \langle kl|v|ij\rangle + \frac{1}{2} \sum_{cd} \langle kl|v|cd\rangle \langle cd|t|ij\rangle \tag{17} \end{equation} $$

 
$$ \begin{equation} \langle kb| \chi |cj\rangle = \langle kb|v|cj\rangle + \frac{1}{2} \sum_{dl} \langle kl|v|cd\rangle \langle db|t|lj\rangle \tag{18} \end{equation} $$

 
$$ \begin{equation} \langle ab| \chi |cd\rangle = \langle ab|v|cd\rangle \tag{19} \end{equation} $$

With the introduction of the above intermediates, the CCD equations scale now as \( \mathcal{O}(n_{h}^{2}n_{p}^{4}) \).

Speed up

To further speed up these computations, we see that these sums can be written in terms of matrix-matrix multiplications. It is not obvious how to write all of these sums in such a way, but it is useful to first recall that the expression for the multiplication of two matrices \( \hat{C} =\hat{A}\times \hat{B} \) can be written as

 
$$ \begin{equation} C_{ij} = \sum_{k} A_{ik} \times B_{kj}. \tag{20} \end{equation} $$

 
We observe then that equation (17) can be written as \[ \langle K| \chi |I\rangle = \langle K|v|I\rangle + \frac{1}{2} \sum_{C} \langle K|v|C\rangle \langle C|t|I\rangle \] by mapping the two index pairs \( kl \to K, ij \to I, cd \to C \). The sum looks now like a matrix-matrix multiplication. This is useful because there are packages like BLAS (Basic Linear Algebra Subprograms) which have extremely fast implementations of matrix-matrix multiplication.

Testing brute force against block structures for MBPT(2)

Convergence properties for pure neutron matter I

Convergence properties for pure neutron matter II

Comparing CCD with Monte Carlo