Coupled Cluster theory and infinite nuclear matter

Morten Hjorth-Jensen [1, 2]

[1] National Superconducting Cyclotron Laboratory and Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA
[2] Department of Physics, University of Oslo, N-0316 Oslo, Norway

28th Indian-Summer School on Ab Initio Methods in Nuclear Physics












Introduction

Coester and Kummel first developed the ideas that led to coupled-cluster theory in the late 1950s. The basic idea is that the correlated wave function of a many-body system \( \mid\Psi\rangle \) can be formulated as an exponential of correlation operators \( T \) acting on a reference state \( \mid\Phi\rangle \) $$ \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

Developments and applications of coupled-cluster theory took different routes in chemistry and nuclear physics. In quantum chemistry, coupled-cluster developments and applications have proven to be extremely useful, see for example the review by Barrett and Musial as well as the recent textbook by Shavitt and Barrett. Many previous applications to nuclear physics struggled with the repulsive character of the nuclear forces and limited basis sets used in the computations. Most of these problems have been overcome during the last decade and coupled-cluster theory is one of the computational methods of preference for doing nuclear physics, with applications ranging from light nuclei to medium-heavy nuclei, see for example the recent review by Hagen, Papenbrock, Hjorth-Jensen and Dean.

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

If we assume that we have a two-body operator at most, the Slater-Condon rule gives then an equation for the correlation energy in terms of \( C_i^a \) and \( C_{ij}^{ab} \) only. We get then $$ \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

In our notes on Hartree-Fock calculations, we have already computed the matrix \( \langle \Phi_0 | \hat{H}|\Phi_{i}^{a}\rangle \) and \( \langle \Phi_0 | \hat{H}|\Phi_{ij}^{ab}\rangle \). If we are using a Hartree-Fock basis, then the matrix elements \( \langle \Phi_0 | \hat{H}|\Phi_{i}^{a}\rangle=0 \) and we are left with a correlation energy given by $$ 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

Inserting the various matrix elements we can rewrite the previous equation as $$ \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

We see that on the right-hand side we have the energy \( E \). This leads to a non-linear equation in the unknown coefficients. These equations are normally solved iteratively ( that is we can start with a guess for the coefficients \( C_i^a \)). A common choice is to use perturbation theory for the first guess, setting thereby $$ 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

For \( C_{jk}^{bc} \) we need then $$ \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

At the end we can rewrite our solution of the Schroedinger equation in terms of \( n \) coupled equations for the coefficients \( C_H^P \). This is a very cumbersome way of solving the equation. However, by using this iterative scheme we can illustrate how we can compute the various terms in the wave operator or correlation operator \( \hat{C} \). We will later identify the calculation of the various terms \( C_H^P \) as parts of different many-body approximations to full CI. In particular, we can relate this non-linear scheme with Coupled Cluster theory and many-body perturbation theory.











Summarizing FCI and bringing in approximative methods

If we can diagonalize large matrices, FCI is the method of choice since:









Summarizing FCI and bringing in approximative methods

The correlation energy is defined as, with a two-body Hamiltonian, $$ \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

However, as we have seen, even for a small case like the four first major shells and a nucleus like oxygen-16, the dimensionality becomes quickly intractable. If we wish to include single-particle states that reflect weakly bound systems, we need a much larger single-particle basis. We need thus approximative methods that sum specific correlations to infinite order.

Popular methods are

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

The energy is given by $$ \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

The coupled cluster energy is a function of the unknown cluster amplitudes \( t_{i_1i_2\ldots i_n}^{a_1a_2\ldots a_n} \), given by the solutions to the amplitude equations $$ \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

A much used approximation is to truncate the cluster operator \( \hat{T} \) at the \( n=2 \) level. This defines the so-called singes and doubles approximation to the Coupled Cluster wavefunction, normally shortened to CCSD..

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

The amplutudes \( t \) play a role similar to the coefficients \( C \) in the shell-model calculations. They are obtained by solving a set of non-linear equations similar to those discussed above in connection withe FCI discussion.

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):

\( 0p-0h \) \( 1p-1h \) \( 2p-2h \) \( 3p-3h \) \( 4p-4h \) \( 5p-5h \) \( 6p-6h \)
\( 0p-0h \) \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) 0 0 0 0
\( 1p-1h \) 0 \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) 0 0 0
\( 2p-2h \) 0 \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) 0 0
\( 3p-3h \) 0 \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) 0
\( 4p-4h \) 0 0 \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \)
\( 5p-5h \) 0 0 0 \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \)
\( 6p-6h \) 0 0 0 0 \( \tilde{x} \) \( \tilde{x} \) \( \tilde{x} \)











A quick tour of Coupled Cluster theory

In our FCI discussion the correlation energy is defined as, with a two-body Hamiltonian, $$ \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

Coupled cluster theory has several interesting computational features and is the method of choice in quantum chemistry. The method was originally proposed by Coester and Kummel, two nuclear physicists (way back in the fifties). It came back in full strength in nuclear physics during the last decade.

There are several interesting features:

However, Coupled Cluster theory is









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

Inserting these equations in the expression for the computation of the energy we have, with a Hamiltonian defined with respect to a general vacuum (see the exercises in the second quantization part) $$ \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

This quantity becomes $$ 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

These equations can be reduced to (after several applications of Wick's theorem) to, for all \( i > j \) and all \( a > b \), $$ \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},& \label{eq:ccd} \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

The operator \( \hat{P}(ij\vert ab) \) is defined as $$ \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

With this starting point, we are now ready to solve Eq. \eqref{eq:ccd} iteratively. Before we attack the full equations, it is however instructive to study a truncated version of the equations. We will first study the following approximation where we take away all terms except the linear terms that involve the single-particle energies and the the two-particle intermediate excitations, that is $$ \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}. \label{eq:ccd1} \end{equation} $$











Approximations to the full CCD equations

Setting the single-particle energies for the hole states equal to an energy variable \( \omega = \epsilon_i+\epsilon_j \), Eq. \eqref{eq:ccd1} reduces to the well-known equations for the so-called \( G \)-matrix, widely used in infinite matter and finite nuclei studies. The equation can then be reordered and solved by matrix inversion. To see this let us define the following quantity $$ \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. \eqref{eq:ccd1}, 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

To solve Eq. \eqref{eq:ccd1}, we would thus start with a guess for the unknown amplitudes, typically using the wave operator defined by first order in perturbation theory, leading to a zeroth approximation to the energy given by second-order perturbation theory for the correlation energy. A simple approach to the solution of Eq. \eqref{eq:ccd1}, is to thus to
  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. \eqref{eq:ccd1} 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

The next approximation is to include the two-hole term in Eq. \eqref{eq:ccd}, a term which allow us to make a link with Green's function theory with two-particle and two-hole correlations. This means that we solve $$ \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}. \label{eq:ccd2} \end{equation} $$ This equation is solved the same way as we would do for Eq. \eqref{eq:ccd1}. The final step is then to include all terms in Eq. \eqref{eq:ccd}.











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











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, \label{eq:betadecay} \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) \label{_auto1} \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) \label{_auto2} \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

\( n_{x}^{2}+n_{y}^{2} \) \( n_{x} \) \( n_{y} \) \( N_{\uparrow \downarrow } \) \( N_{\uparrow \uparrow } \)
0 0 0 2 1
1 -1 0
1 0
0 -1
0 1 10 5
2 -1 -1
-1 1
1 -1
1 1 18 9
4 -2 0
2 0
0 -2
0 2 26 13
5 -2 -1
2 -1
-2 1
2 1
-1 -2
-1 2
1 -2
1 2 42 21











Three-dimensional neutron gas

Using the same approach as made with the two-dimensional electron gas with the single-particle kinetic energy defined as $$ \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

\( n_{x}^{2}+n_{y}^{2}+n_{z}^{2} \) \( n_{x} \) \( n_{y} \) \( n_{z} \) \( N_{\uparrow \downarrow } \)
0 0 0 0 2
1 -1 0 0
1 1 0 0
1 0 -1 0
1 0 1 0
1 0 0 -1
1 0 0 1 14
2 -1 -1 0
2 -1 1 0
2 1 -1 0
2 1 1 0
2 -1 0 -1
2 -1 0 1
2 1 0 -1
2 1 0 1
2 0 -1 -1
2 0 -1 1
2 0 1 -1
2 0 1 1 38
3 -1 -1 -1
3 -1 -1 1
3 -1 1 -1
3 -1 1 1
3 1 -1 -1
3 1 -1 1
3 1 1 -1
3 1 1 1 54

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

Every number of particles for filled shells defines also the number of particles to be used in a given calculation. We use the number of particles to define the density of the system $$ \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), \label{_auto3} \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

A Fourier transform to momentum space of the radial part $$ 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}} \label{_auto4} \end{equation} $$











Developing a program for infinite matter











Codes and reading material











The CCD equation

The CCD equations can be written as $$ \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},& \label{_auto5} \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}. \label{eq:ccdcorr} \end{equation} $$

Solving the CCD equations

One way to solve these equations, is to write equation \eqref{eq:ccd2} as a series of iterative nonlinear algebraic 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),& \label{eq:ccd3} \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)}. \label{_auto6} \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)} |. \label{_auto7} \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 \eqref{eq:ccd2}. 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 \label{_auto8}\\ + \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 \label{_auto9} \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 \label{_auto10} \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 \label{eq:mtxEx} \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 \label{_auto11} \end{equation} $$ $$ \begin{equation} \langle ab| \chi |cd\rangle = \langle ab|v|cd\rangle \label{_auto12} \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}. \label{_auto13} \end{equation} $$ We observe then that equation \eqref{eq:mtxEx} 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)


Figure 1: MBPT2 contribution to the correlation for pure neutron matter with \( N=14 \) neutrons and periodic boundary conditions. Up to approximately 1600 single-particle states have been included in the sums over intermediate states.

Convergence properties for pure neutron matter I


Figure 2: Energy per particle of pure neutron matter computed in the CCD approximation with the Minnesota potential for different numbers of particles with \( \mathrm{N_{max}=20} \).

Convergence properties for pure neutron matter II


Figure 3: Energy per particle of pure neutron matter computed in the CCD approximation with the Minnesota potential for different model space sizes with \( \mathrm{A=114} \).

Comparing CCD with Monte Carlo


Figure 4: CCD, Reference energy and Diffusion Monte Carlo results for pure neutron matter with \( 66 \) neutrons and \( \mathrm{N_{max}=36} \).

© 2013-2016, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license