7. The Coupled Cluster Method#
7.1. 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.
Figure 1: Realistic computations of atomic nuclei with interactions from chiral EFT. The slow increase prior to 2015 is based on quantum Monte Carlo and the no-core shell model. These methods are exponentially expensive (in mass number \(A\)) and meet with exponentially increasing computer power (Moore’s law), thus leading to a progress that is linear in time. Methods such as coupled clusters and in-medium SRG carry a polynomial cost in mass number are transforming the field.
7.2. The normal-ordered Hamiltonian#
We start from the reference state
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 (1), 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
The reference state (1) 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
Here,
is the Fock matrix. We note that the Fock matrix is diagonal in the Hartree-Fock basis. The brackets \(\{\cdots\}\) in Eq. (3) denote normal ordering, i.e. all operators that annihilate the nontrivial vaccum (1) 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\).
7.3. Exercise 1: Practice in normal ordering#
Normal order the expression \(\sum\limits_{pq}\varepsilon_q^p a^\dagger_p a_q\).
Hint.
Answer. We have to move all operators that annihilate the reference state to the right of those that create on the reference state. Thus,
===== =====
We note that \(H = E_{HF} + H_N\), where
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.
7.4. Exercise 2: What does “good” mean?#
How do you know whether a Hartree-Fock state is a “good” reference? Which results of the Hartree-Fock computation will inform you?
Answer. Once the Hartree-Fock equations are solved, the Fock matrix (4) becomes diagonal, and its diagonal elements can be viewed as single-particle energies. Hopefully, there is a clear gap in the single-particle spectrum at the Fermi surface, i.e. after \(A\) orbitals are filled.
===== =====
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?
7.5. Exercise 3: How many nuclei are accessible with the coupled cluster method based on spherical mean fields?#
If one limits oneself to nuclei with mass number up to mass number \(A=60\), how many nuclei can potentially be described with the coupled-cluster method? Which of these nuclei are potentially interesting? Why?
Answer. Nuclear shell closures are at \(N,Z=2,8,20,28,50,82,126\), and subshell closures at \(N,Z=2,6,8,14,16,20,28,32,34,40,50,\ldots\).
In the physics of nuclei, the evolution of nuclear structure as neutrons are added (or removed) from an isotope is a key
interest. Examples are the rare isotopes of helium (He-8,10) oxygen (O-22,24,28), calcium (Ca-52,54,60), nickel (Ni-78) and tin (Sn-100,132). The coupled-cluster method has the potential to address questions regarding these nuclei, and in a several cases was used to make predictions before experimental data was available. In addition, the method can be used to compute neighbors of nuclei with closed subshells.
===== =====
7.6. 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
Here, \(T\) is an operator that induces correlations. We can now demand that the correlate state (11) 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
Here, \(E_c\) is the correlation energy, and the total energy is \(E=E_c+E_{HF}\). The similarity-transformed Hamiltonian is defined as
A more productive view on coupled-cluster theory thus emerges: This method seeks a similarity transformation such that the uncorrelated reference state (1) becomes an exact eigenstate of the similarity-transformed Hamiltonian (13).
7.7. Exercise 4: What \(T\) leads to Hermitian \(\overline{H_N}\) ?#
What are the conditions on \(T\) such that \(\overline{H_N}\) is Hermitian?
Answer. For a Hermitian \(\overline{H_N}\), we need a unitary \(e^T\), i.e. an anti-Hermitian \(T\) with \(T = -T^\dagger\)
===== =====
As we will see below, coupld-cluster theory employs a non-Hermitian Hamiltonian.
7.8. Exercise 5: Understanding (non-unitary) similarity transformations#
Show that \(\overline{H_N}\) has the same eigenvalues as \(H_N\) for arbitrary \(T\). What is the spectral decomposition of a non-Hermitian \(\overline{H_N}\) ?
Answer. Let \(H_N\vert E\rangle = E\vert E\rangle\). Thus
Thus, if \(\vert E\rangle\) is an eigenstate of \(H_N\) with eigenvalue \(E\), then \(e^{-T}\vert E\rangle\) is eigenstate of \(\overline{H_N}\) with the same eigenvalue.
A non-Hermitian \(\overline{H_N}\) has eigenvalues \(E_\alpha\) corresponding to left \(\langle L_\alpha\vert \) and right \(\vert R_\alpha \rangle\) eigenstates. Thus
with bi-orthonormal \(\langle L_\alpha\vert R_\beta\rangle = \delta_\alpha^\beta\).
===== =====
To make progress, we have to specify the cluster operator \(T\). In coupled cluster theory, this operator is
Thus, the operator (15) 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}\).
7.9. Exercise 6: How many unknowns?#
Show that the number of unknowns is as large as the FCI dimension of the problem, using the numbers \(A\) and \(n_u\).
Answer. We have to sum up all \(np-nh\) excitations, and there are \(\binom{n_u}{n}\) particle states and \(\binom{A}{A-n}\) hole states for each \(n\). Thus, we have for the total number
The right hand side are obviously all ways to distribute \(A\) fermions over \(n_0+A\) orbitals.
===== =====
Thus, the coupled-cluster method with the full cluster operator (15) 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
We need to determine the unknown cluster amplitudes that enter in CCSD. Let
be 1p-1h and 2p-2h excitations of the reference. Computing matrix elements of the Schroedinger Equation (12) yields
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. (20).
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.
7.10. Exercise 7: Why is CCD not exact?#
Above we argued that a similarity transformation preserves all eigenvalues. Nevertheless, the CCD correlation energy is not the exact correlation energy. Explain!
Answer. The CCD approximation does not make \(\vert\Phi_0\rangle\) an exact eigenstate of \(\overline{H_N}\); it is only an eigenstate when the similarity-transformed Hamiltonian is truncated to at most 2p-2h states. The full \(\overline{H_N}\), with \(T=T_2\), would involve six-body terms (do you understand this?), and this full Hamiltonian would reproduce the exact correlation energy. Thus CCD is a similarity transformation plus a truncation, which decouples the ground state only from 2p-2h states.
===== =====
7.11. Computing the similarity-transformed Hamiltonian#
The solution of the CCSD equations, i.e. the second and third line of Eq. (20), and the computation of the correlation energy requires us to compute matrix elements of the similarity-transformed Hamiltonian (13). This can be done with the Baker-Campbell-Hausdorff expansion
We now come to a key element of coupled-cluster theory: the cluster operator (15) 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 (23) can only be non-zero because each \(T\) must connect to \(H_N\) (but no \(T\) with another \(T\)). Thus, the expansion is finite.
7.12. Exercise 8: When does CCSD truncate?#
In CCSD and for two-body Hamiltonians, how many nested commutators yield nonzero results? Where does the Baker-Campbell-Hausdorff expansion terminate? What is the (many-body) rank of the resulting \(\overline{H_N}\)?
Answer. CCSD truncates for two-body operators at four-fold nested commutators, because each of the four annihilation and creation operators in \(\overline{H_N}\) can be knocked out with one term of \(T\).
===== =====
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
with
Thus, the CCSD Eqs. (20) for the amplitudes can be written as \(\overline{H}_i^a = 0\) and \(\overline{H}_{ij}^{ab}=0\).
7.13. Exercise 9: Compute the matrix element \(\overline{H}_{ab}^{ij}\equiv \langle ij\vert \overline{H_N}\vert ab\rangle\)#
Answer. This is a simple task. This matrix element is part of the operator \(\overline{H}_{ab}^{ij}a^\dagger_ia^\dagger_ja_ba_a\), i.e. particles are annihilated and holes are created. Thus, no contraction of the Hamiltonian \(H\) with any cluster operator \(T\) (remember that \(T\) annihilates holes and creates particles) can happen, and we simply have \(\overline{H}_{ab}^{ij} = \langle ij\vert V\vert ab\rangle\).
===== =====
We need to work out the similarity-transformed Hamiltonian of Eq. (23). 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.
7.14. Example: The contribution of \([F, T_2]\) to \(\overline{H_N}\)#
The commutator \([F, T_2]\) consists of two-body and one-body terms. Let us compute first the two-body term, as it results from a single contraction (i.e. a single application of \([a_p, a^\dagger_q] = \delta_p^q\)). We denote this as \([F, T_2]_{2b}\) and find
Here we exploited the antisymmetry \(t_{ij}^{ab} = -t_{ji}^{ab} = -t_{ij}^{ba} = t_{ji}^{ba}\) in the last step. Using \(a^\dagger_q a^\dagger_b a_j a_i = -a^\dagger_b a^\dagger_q a_j a_i \) and \(a^\dagger_a a^\dagger_b a_p a_i = a^\dagger_a a^\dagger_b a_i a_p\), we can make the expression manifest antisymmetric, i.e.
Thus, the contribution of \([F, T_2]_{2b}\) to the matrix element \(\overline{H}_{ij}^{ab}\) is
Here we used an arrow to indicate that this is just one contribution to this matrix element. We see that the derivation straight forward, but somewhat tedious. As no one likes to commute too much (neither in this example nor when going to and from work), and so we need a better approach. This is where diagramms come in handy.
===== =====
7.14.1. Diagrams#
The pictures in this Subsection are taken from Crawford and Schaefer.
By convention, hole lines (labels \(i, j, k,\ldots\)) are pointing down.
Figure 2: This is a hole line.
By convention, particle lines (labels \(a, b, c,\ldots\)) are pointing up.
Figure 3: This is a particle line.
Let us look at the one-body operator of the normal-ordered Hamiltonian, i.e. Fock matrix. Its diagrams are as follows.
Figure 4: The diagrams corresponding to \(f_a^b\). The dashed line with the ‘X’ denotes the interaction \(F\) between the incoming and outgoing lines. The labels \(a\) and \(b\) are not denoted, but you should label the outgoing and incoming lines accordingly.
Figure 5: The diagrams corresponding to \(f_i^j\). The dashed line with the ‘X’ denotes the interaction \(F\) between the incoming and outgoing lines.
Figure 6: The diagrams corresponding to \(f_a^i\). The dashed line with the ‘X’ denotes the interaction \(F\) between the incoming and outgoing lines.
Figure 7: The diagrams corresponding to \(f_i^a\). The dashed line with the ‘X’ denotes the interaction \(F\) between the incoming and outgoing lines.
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.
Figure 8: The diagrams corresponding to \(\langle ai\vert V\vert jb \rangle = - \langle ai\vert V\vert bj \rangle = -\langle ia\vert V\vert jb \rangle = \langle ia\vert V\vert bj\rangle\).
7.15. Exercise 10: Assign the correct matrix element \(\langle pq\vert V\vert rs\rangle\) to each of the following diagrams of the interaction#
Remember: \(\langle\rm{left-out, right-out}\vert V\vert \rm{left-in, right-in}\rangle\).
a)
Figure 9:
Answer. \(\langle ab\vert V\vert cd\rangle + \langle ij\vert V\vert kl\rangle + \langle ia\vert V\vert bj\rangle\)
b)
Figure 9:
Answer. \(\langle ai\vert V\vert bc\rangle + \langle ij\vert V\vert ka\rangle + \langle ab\vert V\vert ci\rangle\)
c)
Figure 9:
Answer. \(\langle ia\vert V\vert jk\rangle + \langle ab\vert V\vert ij\rangle + \langle ij\vert V\vert ab\rangle\)
===== =====
Finally, we have the following diagrams for the \(T_1\) and \(T_2\) amplitudes.
Figure 9: The horizontal full line is the cluster amplitude with incoming hole lines and outgoing particle lines as indicated.
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.
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. fig-symmetry are topologically identical, because they consist of incoming particle and hole lines and of outgoing particle and hole lines.
Write down the matrix elements that enter the diagram, and sum over all internal lines.
The overall sign is \((-1)\) to the power of [(number of hole lines) - (number of loops)].
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.
Antisymmetrize the outgoing and incoming lines as necessary.
Please note that this really works. You could derive these rules for yourself from the commutations and factors that enter the Baker-Campbell-Hausdorff expansion. The sign comes obviously from the arrangement of creation and annilhilation operators, while the symmetry factor stems from all the different ways, one can contract the cluster operator with the normal-ordered Hamiltonian.
7.16. Example: CCSD correlation energy#
The CCSD correlation energy, $E_c= \langle \Phi_0\vert \overline{H_N}\vert \Phi_0\rangle$, is the first of the CCSD equations ([20](#ccsd)). It is a vacuum expectation value and thus consists of all diagrams with no external legs. There are three such diagrams:
Figure 10: Three diagrams enter for the CCSD correlation energy, i.e. all diagrams that leave no external legs.
The correponding algebraic expression is \(E_c=\sum_{ia}f^i_a t_i^a +{1\over 4}\sum_{ijab} \langle ij\vert V\vert ab\rangle t_{ij}^{ab} + {1\over 2} \sum_{ijab} \langle ij\vert V\vert ab\rangle t_i^a t_j^b\).
The first algebraic expression is clear. We have one hole line and one loop, giving it a positive sign. There are no equivalent lines or vertices, giving it no symmetry factor. The second diagram has two loops and two hole lines, again leading to a positive sign. We have a pair of equivalent hole lines and a pair of equivalent particle lines, each giving a symmetry factor of \(1/2\). The third diagram has two loops and two hole lines, again leading to a positive sign. We have two indentical vertices (each connecting to a \(T_1\) in the same way) and thus a symmetry factor \(1/2\).
===== =====
7.17. 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).
7.18. Exercise 11: Derive the CCD equations!#
Let us consider the matrix element \(\overline{H}_{ij}^{ab}\). Clearly, it consists of all diagrams (i.e. all combinations of \(T_2\), and a single \(F\) or \(V\) that have two incoming hole lines and two outgoing particle lines. Write down all these diagrams.
Hint. Start systematically! Consider all combinations of \(F\) and \(V\) diagrams with 0, 1, and 2 cluster amplitudes \(T_2\).
Answer.
Figure 11: The diagrams for the \(T_2\) equation, i.e. the matrix elements of \(\overline{H}_{ij}^{ab}\). Taken from Baardsen et al (2013).
The corresponding algebraic expression is
===== =====
Let us now turn to the computational cost of a CCD computation.
7.19. Exercise 12: Computational scaling of CCD#
For each of the diagrams in Fig. fig-ccd write down the computational cost in terms of the number of occupied \(A\) and the number of unoccupied \(n_u\) orbitals.
Answer. The cost is \(A^2 n_u^2\), \(A^2 n_u^3\), \(A^3 n_u^2\), \(A^2 n_u^4\), \(A^4 n_u^2\), \(A^3 n_u^3\), \(A^4 n_u^4\), \(A^4 n_u^4\), \(A^4 n_u^4\), and \(A^4 n_u^4\) for the respective diagrams.
===== =====
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. fig-ccd. We break up the computation into two steps, computing first the intermediate
at a cost of \(A^4n_u^2\), and then
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.”
7.20. Exercise 13: Factorize the remaining diagrams of the CCD equation#
Diagrams 7, 8, and 9 of Fig. fig-ccd also need to be factorized.
Answer. For diagram number 7, we compute
at a cost of \(A^3 n_u^3\) and then compute
at the cost of \(A^3 n_u^3\).
For diagram number 8, we compute
at a cost of \(A^3 n_u^2\), and then compute
at the cost of \(A^3 n_u^2\).
For diagram number 9, we compute
at a cost of \(A^2 n_u^3\) and then compute
at the cost of \(A^3 n_u^3\).
===== =====
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}\).
7.21. Project 14: (Optional) Derive the CCSD equations!#
a) Let us consider the matrix element \(\overline{H}_i^a\) first. Clearly, it consists of all diagrams (i.e. all combinations of \(T_1\), \(T_2\), and a single \(F\) or \(V\) that have an incoming hole line and an outgoing particle line. Write down all these diagrams.
Answer.
Figure 12: The diagrams for the \(T_1\) equation, i.e. the matrix elements of \(\overline{H}_i^a\). Taken from Crawford and Schaefer. Here \(\langle pq\vert \vert rs\rangle \equiv \langle pq\vert V\vert rs\rangle\) and \(f_{pq}\equiv f^p_q\).
b) Let us now consider the matrix element \(\overline{H}_{ij}^{ab}\). Clearly, it consists of all diagrams (i.e. all combinations of \(T_1\), \(T_2\), and a single \(F\) or \(V\) that have two incoming hole lines and two outgoing particle lines. Write down all these diagrams and corresponding algebraic expressions.
Answer.
Figure 13: The diagrams for the \(T_2\) equation, i.e. the matrix elements of \(\overline{H}_{ij}^{ab}\). Taken from Crawford and Schaefer. Here \(\langle pq\vert \vert rs\rangle \equiv \langle pq\vert V\vert rs\rangle\), \(f_{pq}\equiv f^p_q\), and \(P(ab) = 1 - (a\leftrightarrow b)\) antisymmetrizes.
===== =====
We can now turn to the solution of the coupled-cluster equations.
7.22. Solving the CCD equations#
The CCD equations, depicted in Fig. fig-ccd, 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
Dividing by \((f_i^i +f_j^j -f_a^a -f_b^b)\) yields
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
as a starting point, compute \(\overline{H}_{ij}^{ab}\), and find a new \(t_{ij}^{ab}\) from the right-hand side of Eq. (36). We repeat this process until the amplitudes (or the CCD energy) converge.
7.23. 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
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\) ?
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?
Initialize the cluster amplitudes according to Eq. (37), and solve Eq. (36) by iteration. The cluster amplitudes \(T_1\) and \(T_2\) are two- and four-indexed arrays, respectively.
Please note that the contraction of tensors (i.e. the summation over
common indices in products of tensors) is very user friendly and
elegant in python when numpy.einsum is used.
7.24. Project 15: Solve the CCD equations for the pairing problem#
The Hamiltonian is
Check your results and reproduce Fig 8.5 and Table 8.12 from Lecture Notes in Physics 936.
Answer. Click for IPython notebook for FCI and CCD solutions
## 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)
parameters
delta = 1.0 , g = 0.5
MBPT2 energy = -0.2640986494237767 , compared to exact: -0.062393162393162394
iter= 0 erg= -0.36354426495308473 myeps= 1.0
iter= 1 erg= -0.37575624444298944 myeps= 0.03249973798308371
iter= 2 erg= -0.38585587217916006 myeps= 0.02617461198434523
iter= 3 erg= -0.3942075030185592 myeps= 0.021185874889362435
iter= 4 erg= -0.4011132425677644 myeps= 0.0172164337058469
iter= 5 erg= -0.40682329513318677 myeps= 0.014035706985640954
iter= 6 erg= -0.41154472066379555 myeps= 0.011472448299162774
iter= 7 erg= -0.415448787715655 myeps= 0.00939722817179461
iter= 8 erg= -0.41867711472018154 myeps= 0.007710779717883887
iter= 9 erg= -0.4213467805437514 myeps= 0.006336029956428388
iter= 10 erg= -0.4235545667262582 myeps= 0.005212518895903349
iter= 11 erg= -0.4253804729664931 myeps= 0.004292407283064816
iter= 12 erg= -0.42689062703898156 myeps= 0.003537566713430196
iter= 13 erg= -0.4281396916442187 myeps= 0.002917423050500827
iter= 14 erg= -0.4291728542011779 myeps= 0.0024073343568810036
iter= 15 erg= -0.4300274713630842 myeps= 0.0019873548059555165
iter= 16 erg= -0.4307344279442794 myeps= 0.0016412818092326651
iter= 17 erg= -0.4313192597672807 myeps= 0.0013559140004942823
iter= 18 erg= -0.43180308142773266 myeps= 0.001120468290435161
iter= 19 erg= -0.4322033528891356 myeps= 0.00092611836240335
iter= 20 erg= -0.43253451293787193 myeps= 0.0007656268779271082
iter= 21 erg= -0.43280850265722925 myeps= 0.0006330506856384872
iter= 22 erg= -0.4330351980489746 myeps= 0.0005235033844055304
iter= 23 erg= -0.4332227676000893 myeps= 0.00043296328157862334
iter= 24 erg= -0.4333779678407642 myeps= 0.0003581175144831452
iter= 25 erg= -0.4335063876676202 myeps= 0.00029623514326258533
iter= 26 erg= -0.4336126503302039 myeps= 0.0002450635665329708
iter= 27 erg= -0.43370058043009574 myeps= 0.00020274379113034994
iter= 28 erg= -0.4337733420034005 myeps= 0.0001677409980261048
iter= 29 erg= -0.43383355270178803 myeps= 0.00013878755576313752
iter= 30 erg= -0.4338833782156424 myeps= 0.00011483618952922357
iter= 31 erg= -0.43392461036313534 myeps= 9.50214542071232e-05
iter= 32 erg= -0.4339587316746199 myeps= 7.862800997893603e-05
iter= 33 erg= -0.43398696881076043 myeps= 6.506447928127211e-05
iter= 34 erg= -0.43401033674723327 myeps= 5.384188922313506e-05
iter= 35 erg= -0.4340296753237455 myeps= 4.455588548820297e-05
iter= 36 erg= -0.43404567947824346 myeps= 3.6872051156509844e-05
iter= 37 erg= -0.4340589242583793 myeps= 3.0513783718287337e-05
iter= 38 erg= -0.434069885513208 myeps= 2.525228124449486e-05
iter= 39 erg= -0.43407895701178667 myeps= 2.08982684650241e-05
iter= 40 erg= -0.43408646460614353 myeps= 1.729515884276886e-05
iter= 41 erg= -0.4340926779492775 myeps= 1.4313402297628368e-05
iter= 42 erg= -0.4340978201905307 myeps= 1.1845812197343312e-05
iter= 43 erg= -0.43410207599766415 myeps= 9.803701407424543e-06
iter= 44 erg= -0.4341055981945805 myeps= 8.113686925482528e-06
iter= 45 erg= -0.43410851325371014 myeps= 6.715047138304284e-06
iter= 46 erg= -0.43411092584077954 myeps= 5.557535933304656e-06
iter= 47 erg= -0.4341129225755269 myeps= 4.5995745427982914e-06
iter= 48 erg= -0.43411457514368107 myeps= 3.8067557478495853e-06
iter= 49 erg= -0.43411594287215505 myeps= 3.150606413873383e-06
iter= 50 erg= -0.4341170748600758 myeps= 2.607563687969726e-06
iter= 51 erg= -0.4341180117422817 myeps= 2.158127929698868e-06
iter= 52 erg= -0.4341187871486943 myeps= 1.7861618421804902e-06
iter= 53 erg= -0.4341194289120248 myeps= 1.478310547265082e-06
iter= 54 erg= -0.4341199600672291 myeps= 1.2235217292978426e-06
iter= 55 erg= -0.4341203996786241 myeps= 1.012648554028702e-06
iter= 56 erg= -0.43412076352439366 myeps= 8.381210947268386e-07
iter= 57 erg= -0.4341210646630742 myeps= 6.936744264669387e-07
iter= 58 erg= -0.43412131390237374 myeps= 5.741236183312693e-07
iter= 59 erg= -0.43412152018716577 myeps= 4.75177530791271e-07
iter= 60 erg= -0.43412169092059416 myeps= 3.9328472169422255e-07
iter= 61 erg= -0.43412183222982187 myeps= 3.2550592303140516e-07
iter= 62 erg= -0.43412194918596947 myeps= 2.6940851025684355e-07
iter= 63 erg= -0.43412204598613835 myeps= 2.2297915937601343e-07
iter= 64 erg= -0.43412212610405826 myeps= 1.8455156991990885e-07
iter= 65 erg= -0.43412219241476757 myeps= 1.5274664707545314e-07
iter= 66 erg= -0.4341222472978008 myeps= 1.2642299161312937e-07
iter= 67 erg= -0.4341222927225912 myeps= 1.0463593124585968e-07
iter= 68 erg= -0.43412233031915237 myeps= 8.660361048457761e-08
iter= 69 erg= -0.4341223614365778 myeps= 7.16789278285526e-08
iter= 70 erg= -0.43412238719145646 myeps= 5.9326308511289724e-08
iter= 71 erg= -0.4341224085079444 myeps= 4.910248249319024e-08
iter= 72 erg= -0.4341224261509296 myeps= 4.064057537950126e-08
iter= 73 erg= -0.43412244075348105 myeps= 3.3636942196453457e-08
iter= 74 erg= -0.4341224528395664 myeps= 2.7840267822341936e-08
iter= 75 erg= -0.43412246284285466 myeps= 2.304254930229631e-08
iter= 76 erg= -0.43412247112227875 myeps= 1.9071632181026577e-08
iter= 77 erg= -0.4341224779749146 myeps= 1.578502886931811e-08
iter= 78 erg= -0.43412248364664213 myeps= 1.306480946214148e-08
iter= 79 erg= -0.434122488340968 myeps= 1.0813367191382213e-08
Energy = -0.4341224922263279
8. Nucleonic Matter#
We want to compute nucleonic matter using coupled cluster or IMSRG methods, and start with considering the relevant symmetries.
8.1. Exercise 16: Which symmetries are relevant for nuclear matter?#
a) Enumerate continuous and discrete symmetries of nuclear matter.
Answer. The symmetries are the same as for nuclei. Continuous symmetries: translational and rotational invariance. Discrete symmetries: Parity and time reversal invariance.
b) What basis should we use to implement these symmetries? Why do we have to make a choice between the two continuous symmetries? Which basis is most convenient and why?
Answer. Angular momentum and momentum do not commute. Thus, there is no basis that respects both symmetries simultaneously. If we choose the spherical basis, we are computing a spherical blob of nuclear matter and have to contend with surface effects, i.e. with finite size effects. We also need a partial-wave decomposition of the nuclear interaction. This approach has been done followed in [“Coupled-cluster studies of infinite nuclear matter, “ G. Baardsen, A. Ekström, G. Hagen, M. Hjorth-Jensen, arXiv:1306.5681, Phys. Rev. C 88, 054312 (2013)]. If we choose a basis of discrete momentum states, translational invariance can be respected. This also facilitates the implementetation of modern nuclear interactions (which are often formulated in momentum space in effective field theories). However, we have to think about the finite size effects imposed by periodic boundary conditions (or generalized Bloch waves). This approch was implemented in [“Coupled-cluster calculations of nucleonic matter,” G. Hagen, T. Papenbrock, A. Ekström, K. A. Wendt, G. Baardsen, S. Gandolfi, M. Hjorth-Jensen, C. J. Horowitz, arXiv:1311.2925, Phys. Rev. C 89, 014319 (2014)].
===== =====
8.2. 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)\).
8.3. Exercise 17: Determine the basis states#
What are the discrete values of momenta admissable in \((k_x, k_y, k_z)\)?
Answer. In 1D position space \(\psi_k(x) \propto e^{i k x}\) with \(k = {2\pi n\over L}\) and \(n=0, \pm 1, \pm 2, \ldots\) fulfill \(\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.
8.4. Exercise 18: How large should the basis be?#
What values should be chosen for the box size \(L\) and for the maximum number \(n_{\rm max}\) , i.e. for the discrete momenta \(k = {2\pi n\over L}\) and \(n=0, \pm 1,\ \pm 2, \ldots, \pm n_{\rm max}\)?
Answer. Usually \(n_{\rm max}\) is fixed by computational cost, because we have \((2n_{\rm max}+1)^3\) basis states. We have used \(n_{\rm max}=4\) or up to \(n_{\rm max}=6\) in actual calculations to get converged results.
The maximum momentum must fulfill \(k_{n_{\rm max}}> \Lambda\), where \(\Lambda\) is the momemtum cutoff of the interaction. This then fixes \(L\) for a given \(n_{\rm max}\).
===== =====
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?
8.5. Exercise 19: Determine the lowest few magic numbers for a cubic lattice.#
Answer. As the spin-degeneracy is \(g_s=2\), we have the magic numbers \(g_s N\) with \(N=1, 7, 19, 27, 33, 57, 66, \ldots\) for neutrons.
===== =====
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.
8.6. Finite size effects#
We could also have considered the case of a more generalized boundary condition, i.e. \(\psi_k(x+L) =e^{i\theta}\psi_k(x)\). Admissable momenta that fulfill such a boundary condition are \(k_n(\theta) = {2\pi n +\theta \over L}\). Avering over the “twist” angle \(\theta\) removes finite size effects, because the discrete momenta are really drawn form a continuum. In three dimensions, there are three possible twist angles, and averaging over twist angles implies summing over many results corresponding to different angles. Thus, the removal of finite-size effects significantly increases the numerical expense. An example is shown in Figure 14, where we compute the kinetic energy per particle
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.
Figure 14: Relative finite-size corrections for the kinetic energy in pure neutron matter at the Fermi momentum \(k_F = 1.6795 {\rm fm}^{-1}\) versus the neutron number A. TABC10 are twist-averaged boundary conditions with 10 Gauss-Legendre points in each spatial direction.
8.7. 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
with \(\varepsilon_{\vec{k}, \sigma}^{\vec{k}, \sigma} = {k^2\over 2m}\). In Eq. (40) 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
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
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 (40) and (45) 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
Here,
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,
projects onto two-particle spin-isospin states as indicated, while
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
The radial functions are
and the parameters are as follows
| $\alpha$ | $V_\alpha$ | $\kappa_\alpha$ |
|---|---|---|
| $R$ | +200 MeV | 1.487 fm $^{-2}$ |
| $S$ | -91.85 MeV | 0.465 fm $^{-2}$ |
| $T$ | -178 MeV | 0.639 fm $^{-2}$ |
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
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\)
and it is understood that there is no contribution from \(V_T\). Please note that the matrix elements (57) 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.
8.8. Example: Channel structure and its usage#
We have single-particle states with momentum and spin, namely
Naively, two-body states are then
but using the center-of-mass transformation (41) we can rewrite
where \(P_{rs} = k_r +k_s\) is the total momentum and \(k_{rs}=(k_r-k_s)/2\) is the relative momentum. This representation of two-body states is well adapated to our problem, because the interaction and the \(T_2\) amplitudes preserve the total momentum. Thus, we store the cluster amplitudes \(t_{ij}^{ab}\) as matrices
and the conservation of total momentum is explicit.
Likewise, the pppp, pphh, and hhhh parts of the interaction can be written in this form, namely
and we also have
Using these objects, diagrams (1), (4), and (5) of Figure 11 can be done for each block of momentum \(P_{ij}\) as a copy and matrix-matrix multiplications, respectively
Similarly, the CCD correlation energy results from
These efficiences cannot be used for the sixth diagram of Figure 11. One could either change to a ph formulation, noting that \(k_a-k_i = k_k -k_c\) is also a preserved quantity in \(t_{ik}^{ac}\) and that \(k_k-k_c = k_j-k_b\) is preserved in \(V_{cj}^{kb}\). Thus \(\sum_{kc} t_{ik}^{ac}V_{cj}^{kb}\) has a conserved quantity \(k_k-k_c\) in the loop, and we can again use matrix-matrix multiplications for this diagram. This requires us to store the \(T_2\) amplitude in a phhp and in the usual pphh formulation. Alternatively, we could simply code this diagram with the loops over single-particle states. If this seems too tedious, one can also limit CCD to the first five diagrams in Figure 11 (these are the pp and hh ladders), which gives a good description for neutron matter, see the comparison between this and full CCD in [“Coupled-cluster calculations of nucleonic matter,” G. Hagen, T. Papenbrock, A. Ekström, K. A. Wendt, G. Baardsen, S. Gandolfi, M. Hjorth-Jensen, C. J. Horowitz, arXiv:1311.2925, Phys. Rev. C 89, 014319 (2014)].
===== =====
The steps towards the solution of the CCD equations for neutron matter are as follows
For a given density and UV cutoff, set up the lattice, i.e. determine the single-particle basis.
Determine the channels allowed by the (Minnesota) interaction, i.e. sets of two-body states that are connected by the interaction.
Exploit this channel structure when computing the diagrams.
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.
8.9. Exercise 20: Write a CCD code for neutron matter, focusing first on ladder approximation, i.e. including the first five diagrams in Figure 11.#
Answer. Click for IPython notebook
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)
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[2], line 227
217 return erg
220 ########################################################################################
221 # Functions for CCD of neutron matter
222 # Implementation uses only pp and hh ladders
223 #
224 ########################################################################################
--> 227 from numba import jit
228 # compile a few functions to gain speed; should probably done in Fortran or C++,
229 # and called from Python
231 @jit(nopython=True)
232 def minnesota_nn(p_out,s1_out,s2_out,p_in,s1_in,s2_in,Lbox):
File ~/miniforge3/envs/myenv/lib/python3.9/site-packages/numba/__init__.py:55
50 msg = ("Numba requires SciPy version 1.0 or greater. Got SciPy "
51 f"{scipy.__version__}.")
52 raise ImportError(msg)
---> 55 _ensure_critical_deps()
56 # END DO NOT MOVE
57 # ---------------------- WARNING WARNING WARNING ----------------------------
60 from ._version import get_versions
File ~/miniforge3/envs/myenv/lib/python3.9/site-packages/numba/__init__.py:42, in _ensure_critical_deps()
40 raise ImportError(msg)
41 elif numpy_version > (1, 24):
---> 42 raise ImportError("Numba needs NumPy 1.24 or less")
43 try:
44 import scipy
ImportError: Numba needs NumPy 1.24 or less
===== =====
8.10. 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.
9. From Structure to Reactions#
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.
9.1. 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
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
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
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
with \(\vert\tilde{\psi}_R(z)\rangle\) and \(\langle\tilde{\psi}_L(z^*)\vert\) fulfilling
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
for the right state, and
for the left state, make the Eq. (74) 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
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)].
9.2. Computing optical potentials from microscopic input#
The single-particle Green’s function
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
where \(\Sigma^*\) is the self energy and \(G^{(0)}\) is the Hartree-Fock Green’s function
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)]
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 (81). Computing the Green’s function (79) in coupled clusters, and inverting the Dyson equation (80) using
thus allows us to compute the optical potential. We note that the Green’s function (79) resembles in its structure the LIT (72), and we will indeed use a similar approach to compute this object with the coupled-cluster method. Indeed, we compute
by making the ansatz (76) and (77) 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
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.