Table of contents

The Coupled-Cluster Method
      Introduction
      The normal-ordered Hamiltonian
            Exercise 1: Practice in normal ordering
           
            Exercise 2: What does "good" mean?
           
            Exercise 3: How many nuclei are accessible with the coupled cluster method based on spherical mean fields?
           
      The similarity transformed Hamiltonian
            Exercise 4: What \( T \) leads to Hermitian \( \overline{H_N} \) ?
           
            Exercise 5: Understanding (non-unitary) similarity transformations
           
            Exercise 6: How many unknowns?
           
            Exercise 7: Why is CCD not exact?
           
      Computing the similarity-transformed Hamiltonian
            Exercise 8: When does CCSD truncate?
           
            Exercise 9: Compute the matrix element \( \overline{H}_{ab}^{ij}\equiv \langle ij\vert \overline{H_N}\vert ab\rangle \)
           
            Example: The contribution of \( [F, T_2] \) to \( \overline{H_N} \)
           
            Exercise 10: Assign the correct matrix element \( \langle pq\vert V\vert rs\rangle \) to each of the following diagrams of the interaction
           
            Example: CCSD correlation energy
           
      CCD Approximation
            Exercise 11: Derive the CCD equations!
           
            Exercise 12: Computational scaling of CCD
           
            Exercise 13: Factorize the remaining diagrams of the CCD equation
           
            Project 14: (Optional) Derive the CCSD equations!
           
            Solving the CCD equations
      CCD for the pairing Hamiltonian
            Project 15: Solve the CCD equations for the pairing problem
Nucleonic Matter
            Exercise 16: Which symmetries are relevant for nuclear matter?
           
      Basis states
            Exercise 17: Determine the basis states
           
            Exercise 18: How large should the basis be?
           
            Exercise 19: Determine the lowest few magic numbers for a cubic lattice.
           
            Finite size effects
      Channel structure of Hamiltonian and cluster amplitudes
            Example: Channel structure and its usage
           
            Exercise 20: Write a CCD code for neutron matter, focusing first on ladder approximation, i.e. including the first five diagrams in Figure 11.
           
      Benchmarks with the Minnesota potential
From Structure to Reactions
      Electroweak reactions
      Computing optical potentials from microscopic input

Introduction

The coupled-cluster method is an efficient tool to compute atomic nuclei with an effort that grows polynomial with system size. While this might still be expensive, it is now possible to compute nuclei with mass numbers about \( A\approx 100 \) with this method. Recall that full configuration interaction (FCI) such as the no-core shell model exhibits an exponential cost and is therefore limited to light nuclei.

The normal-ordered Hamiltonian

We start from the reference state $$ \begin{equation} \label{HFref} \vert\Phi_0\rangle = \prod_{i=1}^A a^\dagger_i \vert 0\rangle \end{equation} $$

for the description of a nucleus with mass number \( A \). Usually, this reference is the Hartree-Fock state, but that is not necessary. In the shell-model picture, it could also be a product state where the lowest \( A \) harmonic oscillator states are occupied. Here and in what follows, the indices \( i,j,k,\ldots \) run over hole states, i.e. orbitals occupied in the reference state \eqref{HFref}, while \( a,b,c,\ldots \) run over particle states, i.e. unoccupied orbitals. Indices \( p,q,r,s \) can identify any orbital. Let \( n_u \) be the number of unoccupied states, and \( A \) is of course the number of occupied states. We consider the Hamiltonian $$ \begin{equation} \label{Ham} H = \sum_{pq} \varepsilon^p_q a^\dagger_p a_q + \frac{1}{4}\sum_{pqrs}\langle pq\vert V\vert rs\rangle a^\dagger_pa^\dagger_q a_sa_r \end{equation} $$

The reference state \eqref{HFref} is a non-trivial vacuum of our theory. We normal order this Hamiltonian with respect to the nontrivial vacuum state given by the Hartree-Fock reference and obtain the normal-ordered Hamiltonian $$ \begin{equation} \label{HN} H_N = \sum_{pq} f_{pq} \left\{a^\dagger_p a_q\right\} + \frac{1}{4}\sum_{pqrs}\langle pq\vert V\vert rs\rangle \left\{a^\dagger_pa^\dagger_q a_sa_r\right\}. \end{equation} $$

Here, $$ \begin{equation} \label{Fock} f^p_q = \varepsilon^p_q + \sum_i \langle pi\vert V\vert qi\rangle \end{equation} $$

is the Fock matrix. We note that the Fock matrix is diagonal in the Hartree-Fock basis. The brackets \( \{\cdots\} \) in Eq. \eqref{HN} denote normal ordering, i.e. all operators that annihilate the nontrivial vaccum \eqref{HFref} are to the right of those operators that create with respect to that vaccum. Normal ordering implies that \( \langle \Phi_0\vert H_N\vert \Phi_0\rangle = 0 \).

We note that \( H = E_{HF} + H_N \), where $$ \begin{align} E_{HF} &\equiv \langle\Phi_0\vert H\vert \Phi_0\rangle = \sum_{i} \varepsilon^i_i +\frac{1}{2}\sum_{ij}\langle ij\vert V\vert ij\rangle \label{_auto6} \end{align} $$

is the Hartree-Fock energy. The coupled-cluster method is a very efficient tool to compute nuclei when a "good" reference state is available. Let us assume that the reference state results from a Hartree-Fock calculation.

If symmetry-restricted Hartree-Fock is used, one is limited to compute nuclei with closed subshells for neutrons and for protons. On a first view, this might seem as a severe limitation. But is it?

The similarity transformed Hamiltonian

There are several ways to view and understand the coupled-cluster method. A first simple view of coupled-cluster theory is that the method induces correlations into the reference state by expressing a correlated state as $$ \begin{equation} \label{psi} \vert\Psi\rangle = e^T \vert\Phi_0\rangle , \end{equation} $$

Here, \( T \) is an operator that induces correlations. We can now demand that the correlate state \eqref{psi} becomes and eigenstate of the Hamiltonian \( H_N \), i.e. \( H_N\vert \Psi\rangle = E\vert \Psi\rangle \). This view, while correct, is not the most productive one. Instead, we left-multiply the Schroedinger equation with \( e^{-T} \) and find $$ \begin{equation} \label{Schroedinger} \overline{H_N}\vert \Phi_0\rangle = E_c \vert \Phi_0\rangle . \end{equation} $$

Here, \( E_c \) is the correlation energy, and the total energy is \( E=E_c+E_{HF} \). The similarity-transformed Hamiltonian is defined as $$ \begin{equation} \label{Hsim} \overline{H_N} \equiv e^{-T} H_N e^T . \end{equation} $$

A more productive view on coupled-cluster theory thus emerges: This method seeks a similarity transformation such that the uncorrelated reference state \eqref{HFref} becomes an exact eigenstate of the similarity-transformed Hamiltonian \eqref{Hsim}.

As we will see below, coupld-cluster theory employs a non-Hermitian Hamiltonian.

To make progress, we have to specify the cluster operator \( T \). In coupled cluster theory, this operator is $$ \begin{equation} \label{Top} T \equiv \sum_{ia} t_i^a a^\dagger_a a_i + \frac{1}{4}\sum_{ijab}t_{ij}^{ab} a^\dagger_aa^\dagger_ba_ja_i + \cdots + \frac{1}{(A!)^2}\sum_{i_1\ldots i_A a_1 \ldots a_A} t_{i_1\ldots i_A}^{a_1\ldots a_A} a^\dagger_{a_1}\cdots a^\dagger_{a_A} a_{i_A}\cdots a_{i_1} . \end{equation} $$

Thus, the operator \eqref{Top} induces particle-hole (p-h) excitations with respect to the reference. In general, \( T \) generates up to \( Ap-Ah \) excitations, and the unknown parameters are the cluster amplitides \( t_i^a \), \( t_{ij}^{ab} \), ..., \( t_{i_1,\ldots,i_A}^{a_1,\ldots,a_A} \).

Thus, the coupled-cluster method with the full cluster operator \eqref{Top} is exponentially expensive, just as FCI. To make progress, we need to make an approximation by truncating the operator. Here, we will use the CCSD (coupled clusters singles doubles) approximation, where $$ \begin{equation} \label{Tccsd} T \equiv \sum_{ia} t_i^a a^\dagger_a a_i + \frac{1}{4}\sum_{ijab}t_{ij}^{ab} a^\dagger_aa^\dagger_ba_ja_i . \end{equation} $$

We need to determine the unknown cluster amplitudes that enter in CCSD. Let $$ \begin{align} \vert\Phi_i^a\rangle &= a^\dagger_a a_i \vert \Phi_0\rangle , \label{_auto9}\\ \vert\Phi_{ij}^{ab}\rangle &= a^\dagger_a a^\dagger_b a_j a_i \vert \Phi_0\rangle \label{_auto10} \end{align} $$

be 1p-1h and 2p-2h excitations of the reference. Computing matrix elements of the Schroedinger Equation \eqref{Schroedinger} yields $$ \begin{align} \label{ccsd} \langle \Phi_0\vert \overline{H_N}\vert \Phi_0\rangle &= E_c , \\ \langle \Phi_i^a\vert \overline{H_N}\vert \Phi_0\rangle &= 0 , \label{_auto11}\\ \langle \Phi_{ij}^{ab}\vert \overline{H_N}\vert \Phi_0\rangle &= 0 . \label{_auto12} \end{align} $$

The first equation states that the coupled-cluster correlation energy is an expectation value of the similarity-transformed Hamiltonian. The second and third equations state that the similarity-transformed Hamiltonian exhibits no 1p-1h and no 2p-2h excitations. These equations have to be solved to find the unknown amplitudes \( t_i^a \) and \( t_{ij}^{ab} \). Then one can use these amplitudes and compute the correlation energy from the first line of Eq. \eqref{ccsd}.

We note that in the CCSD approximation the reference state is not an exact eigenstates. Rather, it is decoupled from simple states but \( \overline{H} \) still connects this state to 3p-3h, and 4p-4h states etc.

At this point, it is important to recall that we assumed starting from a "good" reference state. In such a case, we might reasonably expect that the inclusion of 1p-1h and 2p-2h excitations could result in an accurate approximation. Indeed, empirically one finds that CCSD accounts for about 90% of the corelation energy, i.e. of the difference between the exact energy and the Hartree-Fock energy. The inclusion of triples (3p-3h excitations) typically yields 99% of the correlation energy.

We see that the coupled-cluster method in its CCSD approximation yields a similarity-transformed Hamiltonian that is of a two-body structure with respect to a non-trivial vacuum. When viewed in this light, the coupled-cluster method "transforms" an \( A \)-body problem (in CCSD) into a two-body problem, albeit with respect to a nontrivial vacuum.

Computing the similarity-transformed Hamiltonian

The solution of the CCSD equations, i.e. the second and third line of Eq. \eqref{ccsd}, and the computation of the correlation energy requires us to compute matrix elements of the similarity-transformed Hamiltonian \eqref{Hsim}. This can be done with the Baker-Campbell-Hausdorff expansion $$ \begin{align} \label{BCH} \overline{H_N} &= e^{-T} H_N e^T \\ &=H_N + \left[ H_N, T\right]+ \frac{1}{2!}\left[ \left[ H_N, T\right], T\right] + \frac{1}{3!}\left[\left[ \left[ H_N, T\right], T\right], T\right] +\ldots . \label{_auto13} \end{align} $$

We now come to a key element of coupled-cluster theory: the cluster operator \eqref{Top} consists of sums of terms that consist of particle creation and hole annihilation operators (but no particle annihilation or hole creation operators). Thus, all terms that enter \( T \) commute with each other. This means that the commutators in the Baker-Campbell-Hausdorff expansion \eqref{BCH} can only be non-zero because each \( T \) must connect to \( H_N \) (but no \( T \) with another \( T \)). Thus, the expansion is finite.

We see that the (disadvantage of having a) non-Hermitian Hamiltonian \( \overline{H_N} \) leads to the advantage that the Baker-Campbell-Hausdorff expansion is finite, thus leading to the possibility to compute \( \overline{H_N} \) exactly. In contrast, the IMSRG deals with a Hermitian Hamiltonian throughout, and the infinite Baker-Campbell-Hausdorff expansion is truncated at a high order when terms become very small.

We write the similarity-transformed Hamiltonian as $$ \begin{align} \overline{H_N}=\sum_{pq} \overline{H}^p_q a^\dagger_q a_p + {1\over 4} \sum_{pqrs} \overline{H}^{pq}_{rs} a^\dagger_p a^\dagger_q a_s a_r + \ldots \label{_auto14} \end{align} $$ with $$ \begin{align} \overline{H}^p_q &\equiv \langle p\vert \overline{H_N}\vert q\rangle , \label{_auto15}\\ \overline{H}^{pq}_{rs} &\equiv \langle pq\vert \overline{H_N}\vert rs\rangle . \label{_auto16} \end{align} $$

Thus, the CCSD Eqs. \eqref{ccsd} for the amplitudes can be written as \( \overline{H}_i^a = 0 \) and \( \overline{H}_{ij}^{ab}=0 \).

We need to work out the similarity-transformed Hamiltonian of Eq. \eqref{BCH}. To do this, we write \( T=T_1 +T_2 \) and \( H_N= F +V \), where \( T_1 \) and \( F \) are one-body operators, and \( T_2 \) and \( V \) are two-body operators.

Diagrams

By convention, hole lines (labels \( i, j, k,\ldots \)) are pointing down.

By convention, particle lines (labels \( a, b, c,\ldots \)) are pointing up.

Let us look at the one-body operator of the normal-ordered Hamiltonian, i.e. Fock matrix. Its diagrams are as follows.

We now turn to the two-body interaction. It is denoted as a horizontal dashed line with incoming and outgoing lines attached to it. We start by noting that the following diagrams of the interaction are all related by permutation symmetry.

Finally, we have the following diagrams for the \( T_1 \) and \( T_2 \) amplitudes.

We are now in the position to construct the diagrams of the similarity-transformed Hamiltonian, keeping in mind that these diagrams correspond to matrix elements of \( \overline{H_N} \). The rules are as follows.

1. Write down all topologically different diagrams corresponding to the desired matrix element. Topologically different diagrams differ in the number and type of lines (particle or hole) that connect the Fock matrix \( F \) or the interaction \( V \) to the cluster amplitudes \( T \), but not whether these connections are left or right (as those are related by antisymmetry). As an example, all diagrams in Fig. 8 are topologically identical, because they consist of incoming particle and hole lines and of outgoing particle and hole lines.
2. Write down the matrix elements that enter the diagram, and sum over all internal lines.
3. The overall sign is \( (-1) \) to the power of [(number of hole lines) – (number of loops)].
4. Symmetry factor: For each pair of equivalent lines (i.e. lines that connect the same two operators) multiply with a factor \( 1/2 \). For \( n \) identical vertices, multiply the algebraic expression by the symmery factor \( 1/n! \) to account properly for the number of ways the diagram can be constructed.
5. Antisymmetrize the outgoing and incoming lines as necessary.

CCD Approximation

In what follows, we will consider the coupled cluster doubles (CCD) approximation. This approximation is valid in cases where the system cannot exhibit any particle-hole excitations (such as nuclear matter when formulated on a momentum-space grid) or for the pairing model (as the pairing interactions only excites pairs of particles). In this case \( t_i^a=0 \) for all \( i, a \), and \( \overline{H}_i^a=0 \). The CCD approximation is also of some sort of leading order approximation in the Hartree-Fock basis (as the Hartree-Fock Hamiltonian exhibits no particle-hole excitations).

Let us now turn to the computational cost of a CCD computation.

Note that \( n_u\gg A \) in general. In textbooks, one reads that CCD (and CCSD) cost only \( A^2n_u^4 \). Our most expensive diagrams, however are \( A^4n_u^4 \). What is going on?

To understand this puzzle, let us consider the last diagram of Fig. 11. We break up the computation into two steps, computing first the intermediate $$ \begin{align} \chi_{ij}^{kl}\equiv {1\over 2} \sum_{cd} \langle kl\vert V\vert cd\rangle t_{ij}^{cd} \label{_auto17} \end{align} $$ at a cost of \( A^4n_u^2 \), and then $$ \begin{align} {1\over 2} \sum_{kl} \chi_{ij}^{kl} t_{kl}^{ab} \label{_auto18} \end{align} $$ at a cost of \( A^4n_u^2 \). This is affordable. The price to pay is the storage of the intermediate \( \chi_{ij}^{kl} \), i.e. we traded memory for computational cycles. This trick is known as "factorization."

We are now ready, to derive the full CCSD equations, i.e. the matrix elements of \( \overline{H}_i^a \) and \( \overline{H}_{ij}^{ab} \).

We can now turn to the solution of the coupled-cluster equations.

Solving the CCD equations

The CCD equations, depicted in Fig. 11, are nonlinear in the cluster amplitudes. How do we solve \( \overline{H}_{ij}^{ab}=0 \)? We subtract \( (f_a^a +f_b^b -f_i^i -f_j^j)t_{ij}^{ab} \) from both sides of \( \overline{H}_{ij}^{ab}=0 \) (because this term is contained in \( \overline{H}_{ij}^{ab} \)) and find $$ \begin{align*} (f_i^i +f_j^j -f_a^a -f_b^b)t_{ij}^{ab} &= (f_i^i +f_j^j -f_a^a -f_b^b)t_{ij}^{ab} +\overline{H}_{ij}^{ab} \end{align*} $$

Dividing by \( (f_i^i +f_j^j -f_a^a -f_b^b) \) yields $$ \begin{align} t_{ij}^{ab} &= t_{ij}^{ab} + \frac{\overline{H}_{ij}^{ab}}{f_i^i +f_j^j -f_a^a -f_b^b} \label{iter} \end{align} $$

This equation is of the type \( t=f(t) \), and we solve it by iteration, i.e. we start with a guess \( t_0 \) and iterate \( t_{n+1}=f(t_n) \), and hope that this will converge to a solution. We take the perturbative result $$ \begin{align} \label{pert} \left(t_{ij}^{ab}\right)_0 = \frac{\langle ab\vert V\vert ij\rangle}{f_i^i +f_j^j -f_a^a -f_b^b} \end{align} $$

as a starting point, compute \( \overline{H}_{ij}^{ab} \), and find a new \( t_{ij}^{ab} \) from the right-hand side of Eq. \eqref{iter}. We repeat this process until the amplitudes (or the CCD energy) converge.

CCD for the pairing Hamiltonian

You learned about the pairing Hamiltonian earlier in this school. Convince yourself that this Hamiltonian does not induce any 1p-1h excitations. Let us solve the CCD equations for this problem. This consists of the following steps

1. Write a function that compute the potential, i.e. it returns a four-indexed array (or tensor). We need \( \langle ab\vert V\vert cd\rangle \), \( \langle ij\vert V\vert kl\rangle \), and \( \langle ab\vert V\vert ij\rangle \). Why is there no \( \langle ab\vert V\vert id\rangle \) or \( \langle ai\vert V\vert jb\rangle \) ?
2. Write a function that computes the Fock matrix, i.e. a two-indexed array. We only need \( f_a^b \) and \( f_i^j \). Why?
3. Initialize the cluster amplitudes according to Eq. \eqref{pert}, and solve Eq. \eqref{iter} by iteration. The cluster amplitudes \( T_1 \) and \( T_2 \) are two- and four-indexed arrays, respectively.

Project 15: Solve the CCD equations for the pairing problem

The Hamiltonian is $$ \begin{align} H = \delta \sum_{p=1}^\Omega (p-1)\left(a^\dagger_{p+}a_{p+} + a^\dagger_{p-}a_{p-}\right) -{g \over 2} \sum_{p, q=1}^\Omega a^\dagger_{p+}a^\dagger_{p-} a_{q-} a_{q+} . \label{_auto25} \end{align} $$

Check your results and reproduce Fig 8.5 and Table 8.12 from Lecture Notes in Physics 936.

Answer.

## Coupled clusters in CCD approximation
## Implemented for the pairing model of Lecture Notes in Physics 936, Chapter 8.
## Thomas Papenbrock, June 2018

import numpy as np


def init_pairing_v(g,pnum,hnum):
    """
    returns potential matrices of the pairing model in three relevant channels
    
    param g: strength of the pairing interaction, as in Eq. (8.42)
    param pnum: number of particle states
    param hnum: number of hole states
    
    return v_pppp, v_pphh, v_hhhh: np.array(pnum,pnum,pnum,pnum), 
                                   np.array(pnum,pnum,hnum,hnum), 
                                   np.array(hnum,hnum,hnum,hnum), 
                                   The interaction as a 4-indexed tensor in three channels.
    """
    v_pppp=np.zeros((pnum,pnum,pnum,pnum))
    v_pphh=np.zeros((pnum,pnum,hnum,hnum))
    v_hhhh=np.zeros((hnum,hnum,hnum,hnum))
    
    gval=-0.5*g
    for a in range(0,pnum,2):
        for b in range(0,pnum,2):
            v_pppp[a,a+1,b,b+1]=gval
            v_pppp[a+1,a,b,b+1]=-gval
            v_pppp[a,a+1,b+1,b]=-gval
            v_pppp[a+1,a,b+1,b]=gval
            
    for a in range(0,pnum,2):
        for i in range(0,hnum,2):
            v_pphh[a,a+1,i,i+1]=gval
            v_pphh[a+1,a,i,i+1]=-gval
            v_pphh[a,a+1,i+1,i]=-gval
            v_pphh[a+1,a,i+1,i]=gval
    
    for j in range(0,hnum,2):
        for i in range(0,hnum,2):
            v_hhhh[j,j+1,i,i+1]=gval
            v_hhhh[j+1,j,i,i+1]=-gval
            v_hhhh[j,j+1,i+1,i]=-gval
            v_hhhh[j+1,j,i+1,i]=gval
        
    return v_pppp, v_pphh, v_hhhh
    
    
def init_pairing_fock(delta,g,pnum,hnum):
    """
    initializes the Fock matrix of the pairing model
    
    param delta: Single-particle spacing, as in Eq. (8.41)
    param g: pairing strength, as in eq. (8.42)
    param pnum: number of particle states
    param hnum: number of hole states
    
    return f_pp, f_hh: The Fock matrix in two channels as numpy arrays np.array(pnum,pnum), np.array(hnum,hnum). 
    """
# the Fock matrix for the pairing model. No f_ph needed, because we are in Hartree-Fock basis 
    deltaval=0.5*delta
    gval=-0.5*g
    f_pp = np.zeros((pnum,pnum))
    f_hh = np.zeros((hnum,hnum))

    for i in range(0,hnum,2):
        f_hh[i  ,i  ] = deltaval*i+gval
        f_hh[i+1,i+1] = deltaval*i+gval
        
    for a in range(0,pnum,2):
        f_pp[a  ,a  ] = deltaval*(hnum+a)
        f_pp[a+1,a+1] = deltaval*(hnum+a)
    
    return f_pp, f_hh


def init_t2(v_pphh,f_pp,f_hh):
    """
    Initializes t2 amlitudes as in MBPT2, see first equation on page 345
    
    param v_pphh: pairing tensor in pphh channel
    param f_pp:   Fock matrix in pp channel
    param f_hh:   Fock matrix in hh channel
    
    return t2: numpy array in pphh format, 4-indices tensor
    """
    pnum = len(f_pp)
    hnum = len(f_hh)
    t2_new = np.zeros((pnum,pnum,hnum,hnum))
    for i in range(hnum):
        for j in range(hnum):
            for a in range(pnum):
                for b in range(pnum):
                    t2_new[a,b,i,j] = v_pphh[a,b,i,j] / (f_hh[i,i]+f_hh[j,j]-f_pp[a,a]-f_pp[b,b])
    return t2_new


# CCD equations. Note that the "->abij" assignment is redundant, because indices are ordered alphabetically.
# Nevertheless, we retain it for transparency.
def ccd_iter(v_pppp,v_pphh,v_hhhh,f_pp,f_hh,t2):
    """
    Performs one iteration of the CCD equations (8.34), using also intermediates for the nonliniar terms
    
    param v_pppp: pppp-channel pairing tensor, numpy array
    param v_pphh: pphh-channel pairing tensor, numpy array
    param v_hhhh: hhhh-channel pairing tensor, numpy array
    param f_pp: Fock matrix in pp channel
    param f_hh: Fock matrix in hh channel
    param t2: Initial t2 amplitude, tensor in form of pphh channel
    
    return t2_new: new t2 amplitude, tensor in form of pphh channel
    """
    pnum = len(f_pp)
    hnum = len(f_hh)
    Hbar_pphh = (  v_pphh 
                 + np.einsum('bc,acij->abij',f_pp,t2) 
                 - np.einsum('ac,bcij->abij',f_pp,t2) 
                 - np.einsum('abik,kj->abij',t2,f_hh)
                 + np.einsum('abjk,ki->abij',t2,f_hh)
                 + 0.5*np.einsum('abcd,cdij->abij',v_pppp,t2) 
                 + 0.5*np.einsum('abkl,klij->abij',t2,v_hhhh)
                )

    # hh intermediate, see (8.47)
    chi_hh = 0.5* np.einsum('cdkl,cdjl->kj',v_pphh,t2)

    Hbar_pphh = Hbar_pphh - (  np.einsum('abik,kj->abij',t2,chi_hh) 
                             - np.einsum('abik,kj->abji',t2,chi_hh) )

    # pp intermediate, see (8.46)
    chi_pp = -0.5* np.einsum('cdkl,bdkl->cb',v_pphh,t2)

    Hbar_pphh = Hbar_pphh + (  np.einsum('acij,cb->abij',t2,chi_pp) 
                             - np.einsum('acij,cb->baij',t2,chi_pp) )

    # hhhh intermediate, see (8.48)
    chi_hhhh = 0.5 * np.einsum('cdkl,cdij->klij',v_pphh,t2)

    Hbar_pphh = Hbar_pphh + 0.5 * np.einsum('abkl,klij->abij',t2,chi_hhhh)

    # phph intermediate, see (8.49)
    chi_phph= + 0.5 * np.einsum('cdkl,dblj->bkcj',v_pphh,t2)


    Hbar_pphh = Hbar_pphh + (  np.einsum('bkcj,acik->abij',chi_phph,t2)
                             - np.einsum('bkcj,acik->baij',chi_phph,t2)
                             - np.einsum('bkcj,acik->abji',chi_phph,t2)
                             + np.einsum('bkcj,acik->baji',chi_phph,t2) )
                 
    t2_new=np.zeros((pnum,pnum,hnum,hnum))
    for i in range(hnum):
        for j in range(hnum):
            for a in range(pnum):
                for b in range(pnum):
                    t2_new[a,b,i,j] = (  t2[a,b,i,j] 
                                       + Hbar_pphh[a,b,i,j] / (f_hh[i,i]+f_hh[j,j]-f_pp[a,a]-f_pp[b,b]) )

    return t2_new


def ccd_energy(v_pphh,t2):
    """
    Computes CCD energy. Call as 
    energy = ccd_energy(v_pphh,t2)
    
    param v_pphh: pphh-channel pairing tensor, numpy array
    param t2: t2 amplitude, tensor in form of pphh channel
    
    return energy: CCD correlation energy
    """
    erg = 0.25*np.einsum('abij,abij',v_pphh,t2)
    return erg

###############################
######## Main Program

# set parameters as for model
pnum = 20 # number of particle states
hnum = 10 # number of hole states
delta = 1.0

g = 0.5

print("parameters")
print("delta =", delta, ", g =", g)


# Initialize pairing matrix elements and Fock matrix
v_pppp, v_pphh, v_hhhh = init_pairing_v(g,pnum,hnum)
f_pp, f_hh = init_pairing_fock(delta,g,pnum,hnum)

# Initialize T2 amplitudes from MBPT2
t2 = init_t2(v_pphh,f_pp,f_hh)
erg = ccd_energy(v_pphh,t2)

# Exact MBPT2 for comparison, see last equation on page 365 
exact_mbpt2 = -0.25*g**2*(1.0/(2.0+g) + 2.0/(4.0+g) + 1.0/(6.0+g))
print("MBPT2 energy =", erg, ", compared to exact:", exact_mbpt2)
    
    
# iterate CCD equations niter times
niter=200
mix=0.3
erg_old=0.0
eps=1.e-8
for iter in range(niter):
    t2_new = ccd_iter(v_pppp,v_pphh,v_hhhh,f_pp,f_hh,t2)
    erg = ccd_energy(v_pphh,t2_new)
    myeps = abs(erg-erg_old)/abs(erg)
    if myeps < eps: break
    erg_old=erg
    print("iter=", iter, "erg=", erg, "myeps=", myeps)
    t2 = mix*t2_new + (1.0-mix)*t2 
    
print("Energy = ", erg)

We want to compute nucleonic matter using coupled cluster or IMSRG methods, and start with considering the relevant symmetries.

Basis states

In what follows, we employ a basis made from discrete momentum states, i.e. those states \( \vert k_x, k_y, k_z\rangle \) in a cubic box of size \( L \) that exhibit periodic boundary conditions, i.e. \( \psi_k(x+L) =\psi_k(x) \).

Thus, we use a cubic lattice in momentum space. Note that the momentum states \( e^{i k x} \) are not invariant under time reversal (i.e. under \( k\to -k \)), and also do not exhibit good parity (\( x\to -x \)). The former implies that the Hamiltonian matrix will in general be complex Hermitian and that the cluster amplitudes will in general be complex.

Coupled cluster and IMSRG start from a Hartree-Fock reference state, and we need to think about this next. What are the magic numbers of a cubic lattice for neutron matter?

Given \( n_{\rm max} \) and \( L \) for the basis parameters, we can choose a magic neutron number \( N \). Clearly, the density of the system is then \( \rho=N/L^3 \). This summarizes the requirements for the basis. We choose \( n_{\rm max} \) as large as possible, i.e. as large as computationally feasible. Then \( L \) and \( N \) are constrained by the UV cutoff and density of the system.

Finite size effects

and compare with the infinite result \( T_{\rm inf} = {3\over 10} {\hbar^2 k_F^2\over m} N \) valid for the free Fermi gas. We clearly see strong shell effects (blue dashed line) and that averaging over the twist angles (red full line) very much reduces the shell oscillations. We also note that the neutron number 66 is quite attractive as it exhibits smaller finite size effects than other of the accessible magic numbers.

Channel structure of Hamiltonian and cluster amplitudes

Good quantum numbers for the nuclear interaction (i.e. operators that commute with the Hamiltonian and with each other) are total momentum, and the number of neutrons and protons, and – for simple interactions – also the spin (this is really spin, not orbital anular momentum or total angular momentum, as the latter two do not commute with momentum). Thus the Hamiltonian (and the cluster amplitudes) will consist of blocks, one for each set of quantum numbers. We call the set of quantum numbers that label each such block as a "channel." As the interaction is block diagonal, a numerically efficient implementation of nuclear matter has to take advantage of this channel structure. In fact, neutron matter cannot be computed in a numerically converged way (i.e. for large enough \( n_{\rm max} \)) if one does not exploit the channel structure.

The Hamiltonian is of the structure $$ \begin{align} \label{HQ} H = \sum_{\vec{k}, \sigma} \varepsilon_{\vec{k}, \sigma}^{\vec{k}, \sigma} a^\dagger_{\vec{k}, \sigma}a_{\vec{k}, \sigma} + \sum_{\vec{Q},\vec{p},\vec{k},\sigma_s} V_{\sigma_1\sigma_2}^{\sigma_3\sigma_4}(\vec{p},\vec{k}) a^\dagger_{\vec{Q/2}+\vec{p}, \sigma_3}a^\dagger_{\vec{Q/2}-\vec{p}, \sigma_4} a_{\vec{Q/2}-\vec{k}, \sigma_2}a_{\vec{Q/2}+\vec{k}, \sigma_1} \end{align} $$

with \( \varepsilon_{\vec{k}, \sigma}^{\vec{k}, \sigma} = {k^2\over 2m} \). In Eq. \eqref{HQ} we expressed the single-particle momenta in terms of center-of-mass momentum \( \vec{Q} \) and relative momenta \( (\vec{k},\vec{p}) \), i.e. the incoming momenta \( (\vec{k}_1, \vec{k}_2) \) and outgoing momenta \( (\vec{k}_3, \vec{k}_3) \) are $$ \begin{align} \label{CoM} \vec{k}_1 &= \vec{Q}/2 +\vec{k} ,\\ \vec{k}_2 &= \vec{Q}/2 -\vec{k} , \label{_auto27}\\ \vec{k}_3 &= \vec{Q}/2 +\vec{p} , \label{_auto28}\\ \vec{k}_4 &= \vec{Q}/2 -\vec{p} . \label{_auto29} \end{align} $$

The conservation of momentum is obvious in the two-body interaction as both the annihilation operators and the creation operators depend on the same center-of-mass momentum \( \vec{Q} \). We note that the two-body interaction \( V \) depends only on the relative momenta \( (\vec{k},\vec{p}) \) but not on the center-of-mass momentum. We also note that a local interaction (i.e. an interaction that is multiplicative in position space) depends only on the momentum transfer \( \vec{k}-\vec{p} \). The spin projections \( \pm 1/2 \) are denoted as \( \sigma \).

Thus, the \( T_2 \) operator is $$ \begin{align} \label{t2Q} T_2 = {1\over 4} \sum_{\vec{Q},\vec{p},\vec{k},\sigma_s} t_{\sigma_1\sigma_2}^{\sigma_3\sigma_4}(Q; \vec{p},\vec{k}) a^\dagger_{\vec{Q/2}+\vec{p}, \sigma_3}a^\dagger_{\vec{Q/2}-\vec{p}, \sigma_4} a_{\vec{Q/2}-\vec{k}, \sigma_2}a_{\vec{Q/2}+\vec{k}, \sigma_1} . \end{align} $$

We note that the amplitude \( t_{\sigma_1\sigma_2}^{\sigma_3\sigma_4}(Q; \vec{p},\vec{k}) \) depends on the center-of-mass momentum \( \vec{Q} \), in contrast to the potential matrix element \( V_{\sigma_1\sigma_2}^{\sigma_3\sigma_4}(\vec{p},\vec{k}) \).

In the expressions \eqref{HQ} and \eqref{t2Q} we supressed that \( \sigma_1+\sigma_2 = \sigma_3+\sigma_4 \). So, a channel is defined by \( \vec{Q} \) and total spin projection \( \sigma_1+\sigma_2 \).

Because of this channel structure, the simple solution we implemented for the pairing problem cannot be really re-used when computing neutron matter. Let us take a look at the Minnesota potential $$ \begin{align} V(r) = \left( V_R(r) + {1\over 2}(1+P_{12}^\sigma)V_T(r) + {1\over 2}(1-P_{12}^\sigma)V_S(r)\right) {1\over 2}(1-P_{12}^\sigma P_{12}^\tau). \label{_auto30} \end{align} $$

Here, $$ \begin{align} P^\sigma_{12}&= {1\over 2}(1+\vec{\sigma}_1\cdot\vec{\sigma}_2) , \nonumber\\ P^\tau_{12} &= {1\over 2}(1+\vec{\tau}_1\cdot\vec{\tau}_2) \label{_auto31} \end{align} $$

are spin and isospin exchange operators, respectively, and \( \vec{\sigma} \) and \( \vec{\tau} \) are vectors of Pauli matrices in spin and isospin space, respectively. Thus, $$ \begin{align} {1\over 2}(1-P_{12}^\sigma P_{12}^\tau) & = \vert S_{12}=0, T_{12}=1\rangle\langle S_{12}=0,T_{12}=1\vert + \vert S_{12}=1, T_{12}=0\rangle\langle S_{12}=1,T_{12}=0\vert \label{_auto32} \end{align} $$ projects onto two-particle spin-isospin states as indicated, while $$ \begin{align} {1\over 2}(1-P_{12}^\sigma) & = \vert S_{12}=0\rangle\langle S_{12}=0\vert , \label{_auto33}\\ {1\over 2}(1+P_{12}^\sigma) & = \vert S_{12}=1\rangle\langle S_{12}=1\vert \label{_auto34} \end{align} $$ project onto spin singlet and spin triplet combinations, respectively. For neutron matter, two-neutron states have isospin \( T_{12}=1 \), and the Minnesota potential has no triplet term \( V_T \). For the spin-exchange operator (and spins \( s_1, s_2=\pm 1/2 \)), we have \( P_{12}^\sigma\vert s_1s_2\rangle= \vert s_2s_1\rangle \). For neutron matter, \( P_{12}^\tau=1 \), because the states are symmetric under exchange of isospin. Thus, the Minnesota potential simplifies significantly for neutron matter as \( V_T \) does not contribute.

We note that the spin operator has matrix elements $$ \begin{align} \langle s_1' s_2'\vert {1\over 2}(1-P_{12}^{\sigma})\vert s_1 s_2\rangle = {1\over 2} \left( \delta_{s_1}^{s_1'}\delta_{s_2}^{s_2'} -\delta_{s_1}^{s_2'}\delta_{s_2}^{s_1'}\right) . \label{_auto35} \end{align} $$

The radial functions are $$ \begin{align} V_R(r) &=V_R e^{-\kappa_R r^2} , \label{_auto36}\\ V_S(r) &=V_S e^{-\kappa_S r^2} , \label{_auto37}\\ V_T(r) &=V_T e^{-\kappa_T r^2} , \label{_auto38}\\ \label{_auto39} \end{align} $$ and the parameters are as follows

Note that \( \kappa_\alpha^{1/2} \) sets the momentum scale of the Minnesota potential. We see that we deal with a short-ranged repulsive core (the \( V_R \) term) and longer ranged attractive terms in the singlet (the term \( V_S \)) and triplet (the term \( V_T \)) channels.

A Fourier transform (in the finite cube of length \( L \)) yields the momentum-space form of the potential $$ \begin{align} \langle k_p k_q \vert V_\alpha\vert k_r k_s \rangle = {V_\alpha\over L^3} \left({\pi\over\kappa_\alpha}\right)^{3/2} e^{- {q^2 \over 4\kappa_\alpha}} \delta_{k_p+k_q}^{k_r+k_s} . \label{_auto40} \end{align} $$ Here, \( q\equiv {1\over 2}(k_p-k_q-k_r+k_s) \) is the momentum transfer, and the momentum conservation \( k_p+k_q=k_r+k_s \) is explicit.

As we are dealing only with neutrons, the potential matrix elements (including spin) are for \( \alpha = R, S \) $$ \begin{align} \label{ME_of_V} \langle k_p s_p k_q s_q\vert V_\alpha\vert k_r s_r k_s s_s\rangle &= \langle k_p k_q \vert V_\alpha\vert k_r k_s \rangle {1\over 2}\left(\delta_{s_p}^{s_r}\delta_{s_q}^{s_s} - \delta_{s_p}^{s_s}\delta_{s_q}^{s_r}\right) , \end{align} $$ and it is understood that there is no contribution from \( V_T \). Please note that the matrix elements \eqref{ME_of_V} are not yet antisymmetric under exchange, but \( \langle k_p s_p k_q s_q\vert V_\alpha\vert k_r s_r k_s s_s\rangle - \langle k_p s_p k_q s_q\vert V_\alpha\vert k_s s_s k_r s_r\rangle \) is.

The steps towards the solution of the CCD equations for neutron matter are as follows

1. For a given density and UV cutoff, set up the lattice, i.e. determine the single-particle basis.
2. Determine the channels allowed by the (Minnesota) interaction, i.e. sets of two-body states that are connected by the interaction.
3. Exploit this channel structure when computing the diagrams.
4. Solve the coupled-cluster equations. Here we start first with the pp and hh ladders, i.e. using only the first five diagrams of Figure 11.

import numpy as np

##############################################################
# CCD Program for neutron matter with the Minnesota potential.
#
# Thomas Papenbrock, July/August 2018
#
# License: Free Software following Version 3 of GNU General Public License, 
# see https://www.gnu.org/licenses/gpl.html
#######


##########################
# Class for neutron matter 
#######

class MomSpaceBasis:
    """
    momentum-space basis class
    The constructor has the form
    
    MomSpaceBasis(Nmax,kmax)
    
    param Nmax: Number of lattice points in positive kx-direction
    param kmax: Highest lattice momentum (in 1/fm)
    
    return: MomSpaceBasis as a single-partcle basis. 
    attributes of MomSpaceBasis are
    
    dk : lattice spacing in 1/fm
    Lbox : linear dimension (in fm) of cubic box
    nvec : lattice vectors (integers)
    kvec : lattice momentum vectors (floats, in 1/fm)
    ngrid : total number of lattice points 
    """
    def __init__(self,Nmax,kmax,ordered=True):
        """
        the constructor
        
        Generates a cubic lattice in momentum space
        param Nmax: Number of lattice points in positive kx-direction
        param kmax: Highest lattice momentum (in 1/fm)
        param ordered: Optional parameter, True by default, will order lattice points by kinetic energy
        
        return MomSpaceBasis
        """
        self.Nmax = Nmax
        self.dim = 0
        self.ngrid = 0
        self._kvec=[]
        self._nvec=[]

        dk = kmax / Nmax
        self.dk = dk
        self.Lbox = 2.0*np.pi/dk
        
        nx=[]
        nvec=[]
        for i in range(-Nmax,Nmax+1):
            self.dim=self.dim+1
            nx.append(i)
        
        #print('nx=',nx)
        
        for i in nx:
            for j in nx:
                for k in nx:
                    nvec.append(np.array([i,j,k], dtype=int))
                    
        #print('nvec=',nvec)
        self.ngrid=len(nvec)
        
        if ordered:
            #print("ordered")
            norm=np.zeros(self.ngrid,dtype=int)
            for i, vec in enumerate(nvec):
                npvec=np.array(vec,dtype=int)
                norm[i]=np.dot(npvec,npvec)
               # print(i, vec, norm[i])
        
            index=np.argsort(norm)
            #print(index)
            self._nvec=[]
            for i, ind in enumerate(index):
                #print(i, ind, nvec[ind])
                self._nvec.append(nvec[ind])
                
        else: 
            self._nvec=nvec  # a list
            
        self._kvec = np.array(self._nvec)*dk # a numpy array

    
    def kvec(self,indx=-1):
        """
        MomSpaceBasis.kvec(i) returns ith momentum vector
        MomSpaceBasis.kvec() returns all momentum vectors
        
        param indx: index of k-vector to be returned, optional
        return 3-vector (if index non-negative), or all vectors if no index specified
        """
        if indx == -1:
            return self._kvec
        else:
            return self._kvec[indx]
        
    def nvec(self,indx=-1):
        """
        MomSpaceBasis.nvec(i) returns ith lattice vector
        MomSpaceBasis.nvec() returns all lattice vectors
        
        param indx: index of lattice vector to be returned, optional
        return 3-vector (if index non-negative), or all lattice vectors if no index specified
        """
        if indx == -1:
            return self._nvec
        else:
            return self._nvec[indx]
        
    def dens(self,num):
        """
        returns density of system if num particles are present
        param num: int, number of particles
        return dens: float
        """
        return num/(self.Lbox)**3
    
    def update(self,dk):
        """
        Uses dk as new lattice spacing and rescales existing lattice
        param dk: in 1/fm lattice spacing in momentum space
        """
        self.Lbox=2.0*np.pi/dk
        self._kvec = np.array(self._nvec)*dk
    
    def __len__(self):
        """
        overloading of the 'len' function
        """
        return self.ngrid
    
    
############
# useful functions

def magic_numbers(basis):
    """
    param basis: MomSpaceBasis object
    return magic: array of magic numbers
    """
    nvecs = basis.nvec()
    vec=np.array(nvecs[0],dtype=int)
    norm = np.dot(vec,vec)
    magic=[]
    for i in range(1,len(nvecs)):
        vec=np.array(nvecs[i],dtype=int)
        norm2 = np.dot(vec,vec)
        if norm2 > norm: 
            magic.append(2*i)
            norm=norm2
    return magic


def get_dk(rho,Num):
    """
    param rho: desired density
    param Num: magic number of particles
    return dk: grid spacing in momentum space (in 1/fm)
    """
    Lbox = (Num/rho)**(1.0/3.0)
    dk = 2.0*np.pi/Lbox
    return dk

def spbasis_from_MomSpaceBasis(lattice_vecs,st_degen):
    """
    converts a lattice to a single particle basis for spin-isospin degeneracy st_degen
    param lattice_vecs: list of lattice vectors for 1st particle
    param st_degen: spin-isospin degeneracy
    return: basis as a list of momenta
    """
    if st_degen != 2: # for now only neutron matter
        print("Unexpected parameter st_degen")
        return lattice_vecs
    
    basis=[]
    for vec in lattice_vecs:
        for st in range(st_degen):
            basis.append(np.array(vec,dtype=int))
            
    return basis



#########################################################
# Functions for comparisons with infinite free Fermi gas
                         
def kF_from_density(rho,st_degen=2):
    """
    Computes Fermi momentum for given density and spin/isospin degeneracy.
    
    param rho: density in inverse fm cubed
    param st_degen: spin-isospin degeneracy; default 2
    return: Fermi momentum in inverse fm
    """
    res = (6.0*(np.pi)**2*rho/st_degen)**(1.0/3.0)
    return res

def EnergyDensity_FermiGas(kF,st_degen=2):
    """
    Computes energy density of free Fermi gas at Fermi momentum and spin/isospin degeneracy
    param kF: Fermi momentum in inverse fm
    param st_degen: spin-isospin degeneracy; default 2
    return: Energy density in MeV/fm**3
    """
    pvec = np.array([kF,0.0,0.0])
    erg = (st_degen*kF**3/(10.0*np.pi**2)) * Tkin(pvec)
    return erg


########################################################################################
# Functions for CCD of neutron matter
# Implementation uses only pp and hh ladders
# 
########################################################################################


from numba import jit  
# compile a few functions to gain speed; should probably done in Fortran or C++, 
# and called from Python

@jit(nopython=True)
def minnesota_nn(p_out,s1_out,s2_out,p_in,s1_in,s2_in,Lbox):
    """
    The Minnesota potential between two neutrons, not yet anti-symmetrized 
    param p_out: relative out momentum
    param p_in : relative in momentum
    param s1_out, s2_out: spin projections of out particles 1 and 2
    param s1_in, s2_in  : spin projections of in particles 1 and 2
    Lbox : size of momentum box
    return: value of potential in MeV; not anti-symmetrized!
    """
    # parameters. VT is not active between two neutrons (no triplet)
    VR = 200.0
    VS = -91.85  # sign typo in Lecture Notes Physics 936, Chap. 8
    kappaR = 1.487
    kappaS = 0.465
    
    qvec=p_out-p_in
    q2=np.dot(qvec,qvec)
    
    s1_i =spin2spinor(s1_in)
    s2_i =spin2spinor(s2_in)
    s1_o =spin2spinor(s1_out)
    s2_o =spin2spinor(s2_out)
    
    spin_part = 0.5 * ( np.dot(s1_i,s1_o)*np.dot(s2_i,s2_o)
                       -np.dot(s1_i,s2_o)*np.dot(s2_i,s1_o) )
    
    
    pot = spin_part * (  VR*np.exp(-0.25*q2/kappaR) / (Lbox*np.sqrt(kappaR))**3 
                       + VS*np.exp(-0.25*q2/kappaS) / (Lbox*np.sqrt(kappaS))**3 )
    
    pot = pot*(np.sqrt(np.pi))**3 

    
    return pot

@jit
def spin_of_index(i):
    """
    Even indices of the lattive have spin up, odds have spin down
    param i: index of sp_basis
    return: spin as +/- 1
    """
    spin = 1-2*np.remainder(i,2)
    return spin

@jit
def spin2spinor(s):
    """
    Makes a two-component spinor of an integer s
    param s: spin = +/- 1
    return: two-component numpy array [1,0] for up and [0,1] for down
    """
    up  =np.array([1.0,0.0])
    down=np.array([0.0,1.0])
    if s == 1:
        return up
    else:
        return down

@jit
def Tkin(pvec):
    """
    Kinetic energy for a momentum vector
    param pvec: 3-component numpy array in inverse fm
    return: kinetic energy of that momentum in MeV
    """
    nucleon_mass = 938.92
    hbarc = 197.33
# More precise numbers for neutron mass and hbar follow.
# For N=14, this yields E_HF = 10.3337 MeV per nucleon in HF. Benchmarked with Ragnar Stroberg.
#     nucleon_mass = 939.56563
#     hbarc = 197.3269718
    p2 = np.dot(pvec,pvec)
    res = 0.5*hbarc**2*p2/nucleon_mass
    return res
    
@jit
def compute_total_Tkin(Nocc,sp_basis,dk):
    """
    Computes total kinetic energy of reference state
    param Nocc, sp_basis, dk: particle number, integer s.p. lattice, delta k 
    return: total kinetic energy
    """
    erg=0.0
    for i in range(Nocc):
        mom_vec = sp_basis[i]
        vec=np.array(mom_vec)*dk
        erg=erg+Tkin(vec)
        
    return erg



@jit
def Fock(pvec,s,sp_basis,Nocc,dk,Lbox):
    """
    Fock matrix of momentum pvec in hh space
    param pvec: 3-component numpy array in inverse fm
    param s: spin as +/- 1 of state
    param_sp_basis, Nocc, dk, Lbox : parameters of s.p. basis and system
    
    return: Fock matrix element = kinetic energy of that momentum in MeV
    """
    res = Tkin(pvec)
    
    dum=0.0
    for i in range(Nocc):
        vec=sp_basis[i]*dk
        si=spin_of_index(i)
        p_in = 0.5*(vec-pvec)
        p_out= p_in
        dum = dum + ( minnesota_nn(p_out,s,si, p_in, s,si,Lbox) 
                     -minnesota_nn(p_out,s,si,-p_in,si, s,Lbox) )  #antisymmetrized Minnesota
        
    res = res+dum
    return res

def compute_E_HF_simple(Nocc,sp_basis,dk):
    """
    Computes HF energy of reference state
    param Nocc, sp_basis, dk: particle number, integer s.p. lattice, delta k 
    return: total HF energy
    """
    erg=compute_total_Tkin(Nocc,sp_basis,dk)

    pot=0.0
    for i in range(Nocc):
        momi=sp_basis[i]*dk
        si = spin_of_index(i)
        for j in range(Nocc):
            momj=sp_basis[j]*dk
            sj = spin_of_index(j)
            p_rel = 0.5*(momi-momj)
            pot = pot + 0.5* (  minnesota_nn(p_rel,si,sj, p_rel,si,sj,Lbox)
                                - minnesota_nn(p_rel,si,sj,-p_rel,sj,si,Lbox) )
            
    erg = erg+pot
    return erg


def get_channels(sp_basis,start1,end1,start2,end2,identical,other_channels=None):
    """
    Returns channels for coupled cluster based on Minnesota potential
    param sp_Basis: A single-particle basis
    param start1: index to start for particle 1
    param end1: index to end for particle 1
    param start2: index to start for particle 2
    param end2: index to end for particle 2
    param identical: True for hh or pp, False for hp
    param other_channels: list of other channels to compare with
    return: channels, p_rel, t2amp. channels is a list of p12, where p12 is a momentum vector; 
            p_rel is a nested list with relative momenta and spins for each channel
    """
    channel=[]
    p_rel=[]
    for i, mom_vecs1 in enumerate(sp_basis[start1:end1]):
        #vec1=np.array(mom_vecs1,dtype=int)
        vec1=mom_vecs1
        spin1=spin_of_index(i)
            
        for j, mom_vecs2 in enumerate(sp_basis[start2:end2]):
            if identical and i==j: continue  #Fortran cycle
            #vec2=np.array(mom_vecs2,dtype=int)
            vec2=mom_vecs2
            spin2=spin_of_index(j)
                
            p12 = vec1+vec2
            prel= vec1-vec2
            spins=np.array([spin1,spin2],dtype=int)
            ps=[prel,spins]

            new=True
            needed=True
            if other_channels is not None: #check whether we need this channel
                needed=False
                for chan_o in other_channels:
                    if (chan_o==p12).all(): 
                        needed=True
                        break
            if needed: #check whether this channel exists already
                for ipos, chan in enumerate(channel):
                    if (chan==p12).all(): 
                        new=False
                        break
                    
            if needed and new: 
                channel.append(p12)
                p_rel.append([ps])
                
            if needed and not new:
                p_rel[ipos].append(ps)
                
    return channel, p_rel 
    

    
def setup_T2_amplitudes(sp_basis,NN,st_degen):
    """
    returns the t2 amplitudes and t2 channels
    param sp_basis: a sp_basis
    param NN: neutron number
    param st_degen: 2 for the moment, spin-isospin degeneracy
    return: hh_channels, pp_channels, p_rel_hh, p_rel_pp, t2amp
            these are the hh and pp channels of T2, lists of the relative momenta, 
            and t2amps as a list of numpy arrays set to zero 
    """
    num_states = len(sp_basis)
    
    hh_channels, p_rel_hh = get_channels(sp_basis,0,NN,0,NN,True)
    print('hh channels=', len(hh_channels))

    pp_channels, p_rel_pp = get_channels(sp_basis,NN,num_states,NN,num_states,True,hh_channels)
    print('pp channels=', len(pp_channels))
    
    if len(pp_channels) != len(hh_channels): print('pp and hh channels do not match')
     

    ordered_pp_channel=[]
    ordered_p_rel_pp=[]
    for i, chanhh in enumerate(hh_channels):
        for j, chanpp in enumerate(pp_channels):
            if (chanpp==chanhh).all():
                ordered_pp_channel.append(chanpp)
                ordered_p_rel_pp.append(p_rel_pp[j])    
                break
                
    pp_channels = ordered_pp_channel
    p_rel_pp = ordered_p_rel_pp
            
    # set t2 amplitudes to zero in each channel
    t2amp = fill_pot(Lbox, dk, pp_channels, hh_channels, p_rel_pp, p_rel_hh, True)
        
    return hh_channels, pp_channels, p_rel_hh, p_rel_pp, t2amp

def fill_pot(Lbox, dk, channels_out, channels_in, p_rel_out, p_rel_in, T2amp=False):
    """
    Fills lists of matrices such as Vhhhh, Vpphh, Vpppp, t2_pphh
    param Lbox: Lbox
    param dk: dk
    param channels_out, channels_in: the channels we have
    param p_rel_out, p_rel_in: the list of [prel, [s1,s2]]
    param T2amp=False: Set to True if t2_pphh needs to be computed
    return: The object of desire as a list of numpy matrices. 
            Contain matrix elements for potentials, zeros if T2amp=True is requested. 
    """
    Vpot=[]
    for i, chan_in in enumerate(channels_in):
        dim_in = len(p_rel_in[i])
        dim_out= len(p_rel_out[i])
        Vpot_chan=np.zeros((dim_out,dim_in))
        if not T2amp: 
            for ii, ps_i in enumerate(p_rel_in[i]):
                [pii, [s1, s2]] = ps_i
                pii = pii*dk*0.5
                for jj, ps_j in enumerate(p_rel_out[i]):
                    if dim_in == dim_out and jj > ii: continue
                    [pjj, [ss1, ss2]] = ps_j
                    pjj = pjj*dk*0.5
                    Vpot_chan[jj,ii] = ( minnesota_nn( pjj,ss1,ss2, pii,s1,s2,Lbox)
                                        -minnesota_nn(-pjj,ss2,ss1, pii,s1,s2,Lbox) )
                    if dim_in == dim_out : Vpot_chan[ii,jj] = Vpot_chan[jj,ii]
        
        Vpot.append(Vpot_chan)
    return Vpot


def init_V(Lbox, dk, hhchannels, ppchannels, p_relhh, p_relpp,zeros=False):
    """
    Sets up Vhhhh, Vpphh, and Vpppp. 
    
    return: Vhhhh, Vpphh, Vpppp as a lists of numpy arrays
    """
    Vhhhh = fill_pot(Lbox, dk, hhchannels, hhchannels, p_relhh, p_relhh,zeros)
    Vpphh = fill_pot(Lbox, dk, ppchannels, hhchannels, p_relpp, p_relhh,zeros)
    Vpppp = fill_pot(Lbox, dk, ppchannels, ppchannels, p_relpp, p_relpp,zeros)
    
    return Vhhhh, Vpphh, Vpppp  
            
@jit
def make_diagram(obj1,obj2,fac):
    """
    Makes diagrams for pp or hh ladders as matrix-matrix multiplications
    """
    hbar_pphh=[]
    dim1=len(obj1)
    for chan in range(dim1):
        mat1 = obj1[chan]
        mat2 = obj2[chan]
        hbar_pphh.append( fac*np.matmul(mat1,mat2) )
        
    return hbar_pphh


def make_diagrams2_3(t2_pphh,fabij):
    hbar_pphh=[]
    for i, t2_mat in enumerate(t2_pphh):
        f_mat = fabij[i]
        hbar_mat = t2_mat*f_mat
        hbar_pphh.append(hbar_mat)
    return hbar_pphh


def compute_hbar(v_pppp,v_pphh,v_hhhh,t2_pphh,fabij):
    diagram1 = v_pphh.copy()
    diagram23 = make_diagrams2_3(t2_pphh, fabij)
    diagram4 = make_diagram(v_pppp,t2_pphh,0.5)
    diagram5 = make_diagram(t2_pphh,v_hhhh,0.5)
    
    hbar_pphh=[]
    for i in range(len(t2_pphh)):
        mat = (  diagram1[i]
               + diagram23[i]
               + diagram4[i]
               + diagram5[i]  )
        hbar_pphh.append(mat)
    
    return hbar_pphh
    
def get_energy_denominator(hh_channels,p_rel_pp,p_rel_hh,sp_basis,Nocc,dk,Lbox):
    res=[]
    fabij=[]
    for i, Ptot in enumerate(hh_channels):
        dimhh=len(p_rel_hh[i])
        dimpp=len(p_rel_pp[i])
        res_mat = np.zeros((dimpp,dimhh))
        f_mat   = np.zeros((dimpp,dimhh))
        for ii, psh_rel in enumerate(p_rel_hh[i]):
            [pij, [si, sj]] = psh_rel  
            p_i = (Ptot+pij)//2
            p_i = p_i + np.array([Nmax,Nmax,Nmax],dtype=int)
            p_j = (Ptot-pij)//2
            p_j = p_j + np.array([Nmax,Nmax,Nmax],dtype=int)
            ssi = (1-si)//2
            ssj = (1-sj)//2
            fii = fock_mtx4[p_i[0],p_i[1],p_i[2],ssi]
            fjj = fock_mtx4[p_j[0],p_j[1],p_j[2],ssj]
            for jj, psp_rel in enumerate(p_rel_pp[i]):
                [pab, [sa, sb]] = psp_rel
                p_a = (Ptot+pab)//2
                p_a = p_a + np.array([Nmax,Nmax,Nmax],dtype=int)
                p_b = (Ptot-pab)//2
                p_b = p_b + np.array([Nmax,Nmax,Nmax],dtype=int)
                ssa = (1-sa)//2
                ssb = (1-sb)//2
                faa = fock_mtx4[p_a[0],p_a[1],p_a[2],ssa]
                fbb = fock_mtx4[p_b[0],p_b[1],p_b[2],ssb]
                
                res_mat[jj,ii] = 1.0 / (fii + fjj - faa - fbb) 
                f_mat[jj,ii] = faa + fbb - fii - fjj
        res.append(res_mat)
        fabij.append(f_mat)
    return res, fabij


def get_t2_from_mbpt(Vpphh,denom):
    """
    param Vpphh: Vpphh
    param denom: energy denominator in pphh format
    return t2: quotient of both, element for element 
    """
    res = []
    for i, vv in enumerate(Vpphh):
        dd = denom[i]
        res_mat = vv*dd  #how simple in python; element by element multiply
        res.append(res_mat)
    return res


def compute_E_CCD(Vpphh,T2pphh):
    erg=0.0
#     erg2=0.0
    for i, t2mat in enumerate(T2pphh):
        vmat  = Vpphh[i]
        erg = erg + 0.25*np.sum(vmat*t2mat)
    return erg


def compute_Fock_4(sp_basis,Nocc,dk,Lbox,Nmax):
    fock_mtx4=np.zeros(shape=(2*Nmax+1, 2*Nmax+1, 2*Nmax+1, 2))
    for i, vec in enumerate(sp_basis):
        pvec=vec*dk
        spin=spin_of_index(i)
        si = (1 - spin)//2
        px=vec[0]+Nmax
        py=vec[1]+Nmax
        pz=vec[2]+Nmax
        fock_mtx4[px,py,pz,si] = Fock(pvec,spin,sp_basis,Nocc,dk,Lbox)
    return fock_mtx4

#####################################
########### Main Program starts here    

from timeit import default_timer as timer
# for timing purposes

progstart=timer()

Nmax=1
kmax=1.0
mbase = MomSpaceBasis(Nmax,kmax)
lattice=mbase.nvec()

## set particle number
NN=14
st_degen=2  # spin up and down
print("chosen N =", NN)
print("magic numbers", magic_numbers(mbase))

## set density
rho=0.08

dk = get_dk(rho,NN)

mbase.update(dk)
Lbox = mbase.Lbox


## get single particle basis

sp_basis = spbasis_from_MomSpaceBasis(lattice,st_degen)
num_states = len(sp_basis)

print('number of s.p. states:', num_states)

# print out a few facts of the reference state
total_Tkin = compute_total_Tkin(NN,sp_basis,dk)
print('total Tkin per particle =', total_Tkin/NN )

k_fermi = kF_from_density(rho)

print("Fermi momentum =", k_fermi)

E_gas = EnergyDensity_FermiGas(k_fermi)

print("Energy per neutron of infinite free Fermi gas", E_gas/rho)

E_HF = compute_E_HF_simple(NN,sp_basis,dk)
E_HF = E_HF/NN
print("HF energy per neutron =", E_HF)

## now we start our business ...
## get all channels and two-body states within those channels; set T2 to zero
hh_channels, pp_channels, p_rel_hh, p_rel_pp, t2_pphh = setup_T2_amplitudes(sp_basis,NN,st_degen)

# get some insight in how big this all is
count=0
for i, channel in enumerate(p_rel_hh):
    dim=len(p_rel_hh[i])
    count=count+dim

print('hh number of total combinations', count)

count=0
for i, channel in enumerate(p_rel_pp):
    dim=len(p_rel_pp[i])
    count=count+dim
    
print('pp number of total combinations', count)


print("get v_hhhh, v_pphh, v_pppp")
start = timer()
v_hhhh, v_pphh, v_pppp = init_V(Lbox, dk, hh_channels, pp_channels, p_rel_hh, p_rel_pp)
end = timer()
print("what a hog!", end-start, 'seconds')


print("compute energy denominator")
start = timer()
fock_mtx4 = compute_Fock_4(sp_basis,NN,dk,Lbox,Nmax)
denom_pphh, f_abij = get_energy_denominator(pp_channels,p_rel_pp,p_rel_hh,sp_basis,NN,dk,Lbox)
end = timer()
print("that's faster", end-start, 'seconds')

print("Initialize T2 from MBPT2")
t2_pphh = get_t2_from_mbpt(v_pphh,denom_pphh)

erg = compute_E_CCD(v_pphh,t2_pphh)
print('MBPT2 correlation energy per neutron =', erg/NN)


print("start CCD iterations ...")

niter=200
mix=0.99
erg_old=0.0
eps=1.e-8
for iter in range(niter):
    
    start = timer()
    hbar_pphh = compute_hbar(v_pppp,v_pphh,v_hhhh,t2_pphh,f_abij)
    end = timer()
    print("time of making Hbar:", end-start, 'seconds')
    
    t2_new = get_t2_from_mbpt(hbar_pphh,denom_pphh)
    
    for i in range(len(t2_new)):
        t2_new[i] = t2_pphh[i] + t2_new[i]
    
    erg = compute_E_CCD(v_pphh,t2_new)
    
    myeps = abs(erg-erg_old)/abs(erg)
    if myeps < eps: break
    erg_old=erg

    print("iter=", iter, "Correlation energy per neutron=", erg/NN, ", epsilon=", myeps)
    
    for i in range(len(t2_pphh)):
        t2_pphh[i] = mix*t2_new[i] + (1.0-mix)*t2_pphh[i]
    
print("Correlation energy per neutron= ", erg/NN)

progend=timer()
print('total time in seconds', progend-progstart)

Benchmarks with the Minnesota potential

For the benchmarks, let us use a nucleon mass \( m=938.92 \) MeV, \( \hbar = 197.33 \) MeV fm, and \( c=1 \). For \( N=14 \) neutrons at a density \( \rho=0.08 {\rm fm}^{-3} \) one finds \( T_{\rm kin}(\vert\phi_0\rangle)/N = 22.427553 \) MeV, and \( E_{HF}/N = 10.3498 \) MeV. In model spaces with \( N_{\rm max}=1 \) and \( N_{\rm max}=2 \), and using only the first five diagrams of Figure 11 for the CCD calculation, yields the correlation energies per particle of \( E_{c}/N =-0.2118 \) MeV and \( E_{c}/N =-0.6923 \) MeV, respectively.

Nuclear coupled-cluster theory has also been used to describe aspects of nuclear reactions, namely photo reactions and computations of optical potentials. In what follows, we want to discuss these approaches.

Electroweak reactions

Let us assume we probe a nucleus with an electroweak probe (e.g. via photon or \( Z \)-boson exchange). The corresponding operator \( \hat{\Theta} \) introduces transitions in the nucleus. For photo reactions, the relevant operator is the dipole operator $$ \begin{align} \hat{\Theta} = \sum_{1=1}^A q_i (\vec{r}_i - \vec{R}_{CoM}) . \label{_auto53} \end{align} $$

Here \( q_i \) is the charge of nucleon \( i \), and \( \vec{R}_{CoM} \) is the position of the center of mass. The structure function or response function describing the reaction is $$ \begin{align} S(\omega) \equiv \sum_f \langle\psi_0|\hat{\Theta}^\dagger|\psi_f\rangle\langle \psi_f |\hat{\Theta}|\psi_0\rangle \delta(E_f-E_0-\omega). \label{_auto54} \end{align} $$

Here, the sum is over all final states. The structure function is difficult to compute because the sum is over (infinitely many) continuum states, and we seek a simpler formulation. The key idea is that the Lorentz integral transform (LIT) of the structure function $$ \begin{align} \label{LIT} L(\omega_0,\Gamma) &\equiv {\Gamma\over\pi} \int d\omega \frac{S(\omega)}{(\omega-\omega_0)^2+\Gamma^2} \\ &= {\Gamma\over\pi}\langle\psi_0|\hat{\Theta}^\dagger {1\over H-E_0-\omega_0+i\Gamma} {1\over H-E_0-\omega_0-i\Gamma}\hat{\Theta}|\psi_0\rangle \nonumber. \end{align} $$

We note that the LIT \( L(\omega_0,\Gamma) \) of the structure function is a ground-state expectation value and thus much easier to compute than the structure function itself. We also note that the LIT is not invertible (mathematically speaking), but making some assumptions about the structure function, and imposing a finite resolution \( \Gamma \) alleviates this problem in practical computation.

We next rewrite the LIT for coupled cluster using the shorthand \( z\equiv E_0+\omega_0+i\Gamma \) as $$ \begin{align} L(\omega_0,\Gamma) &\equiv {\Gamma\over\pi} \langle\tilde{\psi}_L(z^*) \vert \tilde{\psi}_R(z)\rangle , \label{_auto55} \end{align} $$

with \( \vert\tilde{\psi}_R(z)\rangle \) and \( \langle\tilde{\psi}_L(z^*)\vert \) fulfilling $$ \begin{align} \label{keyLR} \left(\overline{H}-z\right) \vert\tilde{\psi}_R(z)\rangle &= \overline{\Theta} \vert\phi_0\rangle , \\ \langle\tilde{\psi}_L(z^*)\vert \left(\overline{H}-z^*\right) &= \langle \phi_L\vert . \label{_auto56} \end{align} $$

Here, \( \langle \phi_L\vert \) is the left ground state, i.e. the left eigenstate of the similarity-transformed Hamiltonian. Note that in the coupled-cluster formulation we have distinguished between left and right states (as these are not adjoints of each other), and replaced all operators by their similarity transformations. Making the ansatz $$ \begin{align} \label{ansatzR} \vert\tilde{\psi}_R(z)\rangle &= \hat{R}\vert\phi_0\rangle \\ &= \left( r_0(z) + \sum_{ia} r_i^a(z) a_a^\dagger a_i + {1\over 4} \sum_{ijab} r_{ij}^{ab}(z) a_a^\dagger a^\dagger_b a_j a_i +\cdots\right)\vert\phi_0\rangle \nonumber \end{align} $$

for the right state, and $$ \begin{align} \label{ansatzL} \langle\tilde{\psi}_L(z^*)\vert&= \langle\phi_L\vert \hat{L} \\ &= \langle \phi_L\vert \left( l_0(z) + \sum_{ia} l^i_a(z) a_i^\dagger a_a + {1\over 4} \sum_{ijab} l^{ij}_{ab}(z) a_i^\dagger a^\dagger_j a_b a_a +\cdots\right) \nonumber \end{align} $$

for the left state, make the Eq. \eqref{keyLR} linear systems of equations that can be solved for the parameters \( r_0, r_i^a, r_{ij}^{ab} \) and \( l_0, l_a^i, l_{ab}^{ij} \) for each value of \( z \). The LIT then becomes $$ \begin{align} L(z) = l_0(z) r_0(z) + \sum_{ia} l_a^i(z) r_i^a(z) +{1\over 4} l_{ab}^{ij}(z) r_{ij}^{ab}(z) \label{_auto57} \end{align} $$ in the CCSD approximation. For details, please see ["Giant and pigmy dipole resonances in 4He, 16,22O, and 40Ca from chiral nucleon-nucleon interactions," S. Bacca, N. Barnea, G. Hagen, M. Miorelli, G. Orlandini, and T. Papenbrock, Phys. Rev. C 90, 064619 (2014)].

Computing optical potentials from microscopic input

The single-particle Green's function $$ \begin{align} \label{Green} G(\alpha, \beta, E) &\equiv \langle\psi_0\vert a_\alpha {1\over E-(H-E_0)+i\eta}a^\dagger_\beta\vert\psi_0\rangle + \langle\psi_0\vert a_\beta^\dagger {1\over E-(H-E_0)-i\eta}a_\alpha\vert\psi_0\rangle \end{align} $$

describes the propagation of particles and holes in the nucleus with a ground state \( \vert\psi_0\rangle \) and energy \( E_0 \). The Green's function fulfills the Dyson equation $$ \begin{align} \label{Dyson} G = G^{(0)} +G^{(0)} \Sigma^* G , \end{align} $$

where \( \Sigma^* \) is the self energy and \( G^{(0)} \) is the Hartree-Fock Green's function $$ \begin{align} \label{HFGreen} G^{(0)}(\alpha, \beta, E) &\equiv \langle\phi_0\vert a_\alpha {1\over E-(H_{HF}-E_{HF})+i\eta}a^\dagger_\beta\vert\phi_0\rangle + \langle\phi_0\vert a_\beta^\dagger {1\over E-(H_{HF}-E_{HF})-i\eta}a_\alpha\vert\phi_0\rangle \\ &= \delta_\alpha^\beta \left({\Theta(\alpha-F) \over E-\varepsilon_\alpha +i\eta} + {\Theta(F-\alpha) \over E-\varepsilon_\alpha -i\eta} \right) \nonumber. \end{align} $$

Here, \( \Theta() \) denotes the unit step function and \( F \) labels the index of the Fermi surface. The key point is now that the optical potential \( \Sigma' \) (which describes the reaction of a single nucleon with the nucleus) is related to the self energy and the Hartree-Fock potential \( U_{HF} \) by [F. Capuzzi and C. Mahaux, "Projection operator approach to the self-energy," Ann. Phys. (NY) 245, 147 (1996)] $$ \begin{align} \Sigma' \equiv \Sigma^* +U_{HF} . \label{_auto58} \end{align} $$

The idea is thus as follows. Starting from a Hartree-Fock calculations enables us to compute the Hartree-Fock potential \( U_{HF} \) and the Hartree-Fock Green's function \eqref{HFGreen}. Computing the Green's function \eqref{Green} in coupled clusters, and inverting the Dyson equation \eqref{Dyson} using $$ \begin{align} \Sigma^* = \left(G^{(0)}\right)^{-1} - G^{-1} \label{_auto59} \end{align} $$

thus allows us to compute the optical potential. We note that the Green's function \eqref{Green} resembles in its structure the LIT \eqref{LIT}, and we will indeed use a similar approach to compute this object with the coupled-cluster method. Indeed, we compute $$ \begin{align} \left(E-(\overline{H}-E_0)+ i\eta\right) \vert\tilde{\psi}_R\rangle &= \overline{a^\dagger_\beta} \vert\phi_0\rangle , \label{_auto60}\\ \langle\tilde{\psi}_L\vert \left(E-(E_0-\overline{H})-i\eta\right) &= \langle \phi_L\vert \overline{a_\beta^\dagger} , \label{_auto61} \end{align} $$

by making the ansatz \eqref{ansatzR} and \eqref{ansatzL} for the right and left states, respectively, solve the resulting linear systems, and then compute \( G(\alpha,\beta,E) = \langle\phi_L\vert \overline{a_\alpha} \vert\tilde{\psi}_R\rangle + \langle \tilde{\psi}_L\vert \overline{a_\alpha}\vert\phi_0\rangle \).

We note that a high quality optical potential can only be obtained if the structure of the nucleus is computed with a good accuracy, i.e. in good agreement to data on binding and separation energies, and charge and matter radii. For details, please see ["Optical potential from first principles," J. Rotureau, P. Danielewicz, G. Hagen, F. Nunes, and T. Papenbrock, Phys. Rev. C 95, 024315 (2017); arXiv:1611.04554]. We also note that this procedure must be repeated for any nucleus of interest, because the optical potential depends on the nucleus under consideration. Once the optical potential \( \Sigma' \) is computed, the single-particle Schroedinger equation $$ \begin{align} \left( -{\hbar^2\Delta\over 2m} +\Sigma'\right)\xi = E\xi \label{_auto62} \end{align} $$

describes the interaction of a single nucleon (and wave function \( \xi \)) with the nucleus.

Finally, we note that this approach does not depend on solving the nucleus \( A \) with the coupled-cluster method. Alternatives, such as the IMSRG or Self-Consistent Green's Function methods could also be used.