Symmetry projection to coupled-cluster singles and doubles wave function through the Monte Carlo method

We begin to consider the coupled-cluster (CC) method FC58 ; CK60 ; KLZ78 ; HN78 . Based on a HF wave function $|\psi_{0})$, we define the CC wave function $|\psi_{cc})$ with the many-body $Z$ operator as $|\psi_{cc})=e^{Z}|\psi_{0}).$ Z is defined as $\sum\chi_{p_{1},p_{2},\cdots h_{1},h_{2},\cdots}a_{p_{1}}^{\dagger}a_{p_{2}}^{\dagger}\cdots a_{h_{1}}a_{h_{2}}\cdots$, where $p$’s ($h$’s) stand for unoccupied (occupied) orbits. The summation is arranged as $Z=Z_{1}+Z_{2}+\cdots$, where the $Z_{k}$ are defined as $Z_{k}=\sum_{i}x_{i}^{(k)}Z_{i}^{(k)}$ and the $Z_{i}^{(k)}$ are the different $k$-body operators and $x_{i}^{(k)}$ are their parameters. We limit $Z$ by up to two-body terms ($k=2)$ due to the increasing numerical complexity in the applications for higher values of $k$, which we call the CCSD approximation.

The parameters $x$ are determined in the following way. If $|\psi_{cc})$ is an exact wave function, we have

where $E$ is the exact energy and $H$ is the Hamiltonian expressed with the deformed operators. The equation corresponds to a non-Hermitian-type eigenvalue problem as follows

where the transformed Hamiltonian $\bar{H}=e^{-Z}He^{Z}$ is non-Hermitian. The energy is obtained by

where we refer to it as the CC energy hereafter. Multiplying Eq.(2) from the left with $(\psi^{*}|$ , the coefficients $x$ are determined by

where $(\psi^{*}|$ is any state orthogonal to $(\psi_{0}|$. Within the CCSD approximation, all operators $Z_{1}^{(k)}$, $Z_{2}^{(k)}$ and their parameters $x_{1}^{(k)}$, $x_{2}^{(k)}$ are renamed as $z_{i}(i=1,\cdots n_{0})$ and $x_{i}(i=1,\cdots n_{0})$ where $n_{0}$ is the number of parameters. The mean-field wave function $|\psi_{0})$ is simply shown by $|0)$ as $0$-th basis state. We evaluate the following matrix elements for an excited state $(k|$ as $a_{k}=(k|H|0)$, $b_{ki}=(k|[H,z_{i}]|0)$ and $c_{kij}=(k|[[H,z_{i}],z_{j}]|0)$ where $i,j,k=1,\cdots n_{0}$. Equation (4) is then rewritten as

The parameters $x$ can be obtained iteratively and we, thus, can evaluate the CC energy with the CCSD wave function.

We add the following iterative procedure to eliminate the one-body operator $Z_{1}$ by changing the mean-field for later convenience. The exponential of the $Z_{1}$ operator changes the mean-field determinant $|0)$ into another mean-field one $|0^{\prime})=e^{Z_{1}}|0)$, that is, the Slater determinant $D$ is modified by $e^{Z_{1}}$ and is changed into $D^{\prime}$. Then we recalculate all the matrix elements with this new $D^{\prime}$ and solve the coupled-cluster equations. This procedure is repeated iteratively. After convergence, we obtain a new mean-field $D$ and the corresponding CCSD wave function as $e^{Z_{2}}|0)$ where $|0)$ stands for the optimal mean-field state. We refer to this form as the CCD wave function, which is the starting point of this study.

In the coupled-cluster method, the energy is given by Eq.(3). This method, however, does not satisfy the Raleigh-Ritz variational principle. Therefore, we evaluate the energy directly as,

where $|\psi)$ stands for the CCSD or CCD wave function. This energy gives an upper limit to the exact energy, unlike the coupled-cluster method. Hereafter we refer to it as the RR energy to distinguish it from the CC energy.

Next, we introduce the symmetry projection operator $P^{L}$, which is a projection onto a state with a good quantum number $L$. For example, in nuclear structure physics, angular momentum projection has been often utilized. For a continuous symmetry, the symmetry projection operator is generally given by

where $\mathcal{N}$ is a normalization, $\mathcal{W}$ is a weight function and $R(\mu)$ is given by $e^{i\hat{O}\cdot\mu}$ where $\hat{O}$ is generally defined in terms of spherical operators, not by deformed ones. The integration is carried out over the parametrization $\mu$ of the continuous group. We will show an example below in section III. If $|\psi)$ is a superposition of several states with different $L$, the projection operator $P^{L}$ extracts $|\psi,L)$ with a definite quantum number $L$ as $|\psi,L)=P^{L}|\psi)=\frac{1}{\mathcal{N}}\int d\mu\mathcal{W}(L,\mu)R(\mu)|\psi)$.

The RR energy $E_{L}$ with the projected coupled-cluster wave function is given by

where we refer to it as projected energy. If we follow the terminology of projection method e.g., variation-after-projection (VAP) or projection-after-variation (PAV), we call this procedure projection-after-coupled-cluster (PACC). Note that this evaluation is usually carried out through spherical operators although the above equation is expressed in the deformed basis. To evaluate the projected energy, we need the matrix elements with the projected wave function as

where the rotated state $|\psi(\mu)\rangle$ is defined by

This is the standard way for projection calculations. In realistic applications with huge Hilbert spaces, this projection is generally not feasible because $|\psi)$ includes $e^{Z_{2}}$ in the CCD wave function, and evaluation of deformed operators through spherical operators additionally increases computational efforts.

To carry out the symmetry projection, we introduce the Monte Carlo procedure. Hereafter we choose the label $s$ to specify each basis of the spherical representation. Similarly, we choose the label $d$ to specify each basis in the deformed mean-field representation. The coupled-cluster wave function is given as $|\psi)$ with the deformed operators.

Inserting the unity operator $\sum_{s}|s\rangle\langle s|=1$, the projected energy in Eq.(8) can be rewritten as

For the spherical basis $s$, we define the projected local energy ${\cal E}_{L}(s)$ as

where $h_{s,s^{\prime}}$ is the Hamiltonian matrix element in the spherical representation and is generally very sparse. We can also define the projected density $\rho_{L}(s)$ as

where $\rho_{L}(s)\geq 0$ and $\sum_{s}\rho_{L}(s)=1$. This property allows us to stochastically generate the distribution of $s$ according to Eq.(13) by the Markov chain Monte Carlo (MCMC) method. Note that this technique has been presented in cases with pair-condensates BM06 ; TI08 ; VMC12 ; VMC18 . Therefore, by applying the Monte Carlo sampling, the symmetry projected energy can be estimated as

where $N_{0}$ is a number of Monte Carlo samples and we refer to it as the MC energy hereafter.

Next, we investigate the projected overlap $\langle s|P^{L}|\phi)$ in detail. By inserting $\sum_{d}|d)(d|=1$, the projected overlap is rewritten as

The correlation part of the overlap, $(d|\psi)=(d|e^{Z_{2}}|0)$ is easily computed because $Z_{2}$ is also represented with the deformed operators and $P^{L}$ does not operate on $|\psi)$. The symmetry projection is easily carried out in the form $\langle s|P^{L}|d)$, where $s$ and $d$ span, in principle, the whole Hilbert space, while on the other hand, for $s$, we use the Monte Carlo sampling so as to avoid the full use of the Hilbert space. For $d$, we also restrict the whole Hilbert space by introducing a truncation scheme due to the overlap $(d|\psi)$. For a certain $\epsilon$, we can truncate the coupled-cluster wave function as

where $\psi(d)=(d|\psi)$ is an amplitude of $|\psi)$ in the basis $d$. As the coupled-cluster wave function $|\psi)$ is an expansion around the optimized deformed mean-field wave function $|0)$, we can naturally expect that this truncation scheme works well as we will see in Sec. III C. Thus, we can also avoid the difficulty with the handling of the full Hilbert space for $d$.

First, we define the three-level Lipkin model with the levels $i=0,1,2$, each having a degeneracy of $N$. We consider $N$ spinless Fermions. The single particle energies, $\varepsilon_{i}$, are $\varepsilon_{0}=0$ and $\varepsilon_{1}=\varepsilon_{2}=1$. The creation and annihilation operators of the $k$-th particle on the $i$-th level are $c_{ik}^{\dagger}$ and $c_{ik}$, which we call, as before, spherical operators. The Hamiltonian is defined on this basis as

where $\varepsilon=1$ and $v$ is the strength of the two-body interaction. The $J$ operators are defined by $J_{pq}=\sum_{i}c_{p,i}^{\dagger}c_{q,i}$. The system is specified by the dimensionless parameter $\chi=v(N-1)/\varepsilon$. The spectra show two phases: One has a vibrational pattern around $\chi\sim 0$ and the other a rotational pattern for $\chi>2$. These two phases continuously change from one into the other when considering a finite number $N$.

The HF approximation is given by the mean-field state $|\phi)=a^{{\dagger}}_{0,1}\cdots a^{{\dagger}}_{0,N}|-\rangle$ with the deformed operators, $a^{{\dagger}}_{i,m}$ ($i=0\sim 2$, $m=1\sim N$), which are canonically related to the spherical ones with coefficients $D_{ij}$’s as $a^{{\dagger}}_{i,m}=\sum_{j}D_{ij}c^{{\dagger}}_{j,m}.$ The HF energy can be shown to be $E_{HF}=0$ for $\chi<1$ and $E_{HF}=\frac{N\varepsilon}{4}(2-\chi-\frac{1}{\chi})$ for $\chi>1$ HB00 .

For the degenerate two upper single-particle energies, there is a symmetry concerning the exchange of 1 and 2 levels. This can be embodied by the symmetry operator (identifiable with a rotational operator) $L=i(J_{21}-J_{12})$ which commutes with the Hamiltonian, $[H,L]=0$. Thereby, all wave functions have the quantum number $L_{0}=0,\pm 1,\pm 2,\cdots$. Except for $L_{0}=0$, the eigenstates with $L_{0}$ are doubly degenerate. With this symmetry operator, we can construct the projection operator as

which projects out the $L_{0}$ component from any wave function $|\psi)$ that is, $|\psi,L_{0})=P^{L_{0}}|\psi)$ where $L_{0}$ and $|\psi,L_{0})$ are eigenvalue and eigenstate of $L$, respectively, that is, $L|\psi,L_{0})=L_{0}|\psi,L_{0})$. Note that $L$ is defined by the $c$ and $c^{\dagger}$ operators and the projection should be carried out using the spherical operators.

The CCSD wave function $|\psi)$ is given by $|\psi)=e^{Z}|0)$ where $|0)$ is the deformed HF state, and $Z$ is organized into a sum of one-body and two-body operators

where the $K$ operators are defined as $K_{pq}=\sum_{i}a_{p,i}^{\dagger}a_{q,i}$. In applying the CCSD procedure described in Sec. II A, we take the components $|1,0)$, $|0,1)$, $|1,1)$, $|2,0)$ and $|0,2)$ in the deformed basis, defined by $|k_{1},k_{2})=\mathcal{N}_{k_{1},k_{2}}(K_{20})^{k_{2}}(K_{10})^{k_{1}}|0)$, with normalization factor $\mathcal{N}_{k_{1},k_{2}}$. With this, the CCSD equation can be solved easily. We can also carry out the CCD calculations by changing the mean-field.

Fig. 1(a) shows the CC energy differences $\Delta E=E-E_{exact}$ where this is relative to the exact ground state energy $E_{exact}$ as a function of $\chi$ for $N=20$. In Fig. 1(b), as a measure of improvement, we plot the ratio between $\Delta E$ and $\Delta E_{HF}=E_{HF}-E_{exact}$ which is the correlation energy. It shows that the coupled-cluster wave function contains a considerable amount of more correlations. Moreover, the CC energy of CCD is lower than that of CCSD. As the coupled-cluster method does not, however, satisfy the Raleigh-Ritz variational principle, we calculate the RR energies in Eq.(6) with the same CCSD and CCD wave functions. CCSD and CCD give almost the same RR energies. The CCD wave function may seem not to be so good. It shows, however, a distinct aspect if we apply the symmetry projection. In Fig. 1, projected ground state energies with CCSD and CCD relative to the exact ones $\Delta E$ and $\Delta E/\Delta E_{HF}$ are plotted. The projected energy of the ground state is vastly improved over the CC energy, especially for the CCD, where about 95% of correlation energy is taken into account. Note that VAP calculations can also be straightforwardly carried out in this simple model, and the calculated $\Delta E$ and $\Delta E/\Delta E_{HF}$ are also plotted as a reference result. The VAP energy with $\sim 10^{-4}-10^{-5}$ percentage error is almost perfect. We will, however, not further dwell on the VAP procedure, since it seems to be numerically intractable for realistic cases. On the other hand very recently the analog to the VAP calculation has been carried out for the pairing Hamiltonian in BD21 . We also will come back to the VAP approach in a follow-up paper in the near future.

In Fig. 2, the exact excitation energies for $L_{0}=1\sim 3$ are plotted as a function of $\chi$ for $N=20$. For $\chi\sim 0$, spectra show a vibrational feature, while $\chi$ around 5, the spectra show a deformed band. For $1<\chi<2$, a crossover between these two phases occurs. We plot the projected excitation energies of Eq.(8) with the CCSD and CCD wave functions, using the symmetry projection in Eq.(18) to show what results in the PACC procedure. The results show reasonably good excitation energies except for the crossover region.

For $0<\chi<1$, CCD and CCSD wave functions are the same by construction and coupled-cluster method gives a pure spherical mean-field. Thereby, in the projected calculations with CCD, we added a slight deformation to the pure spherical mean-field and succeeded in reproducing the excitation energies. Thus, the symmetry projection for the coupled-cluster wave function is quite promising. The remainder of the problem is how to calculate such a symmetry projection in realistic cases. Next, we introduce the Monte Carlo method and mixed representation for the projected overlap as in Eq.(15).

The symmetry projection in Eq.(9) requires heavy computations in realistic applications. So we introduce the Monte Carlo method as described in Sec. II C. The spherical basis $s$ is sampled by the MCMC, which stochastically generates the distribution obeyed to $\rho_{L}(s)$ in Eq.(13). To perform the MCMC, we take the basis states specified by $s$ as a random walker, and we move this $s$-basis to another nearby $s^{\prime}$-basis under the control of the Metropolis-Hasting (MH) algorithm, keeping the detailed balance as in VMC12 ; VMC18 .

As a benchmark test of the MC calculation, we take the projected energy in Eq.(8) with the CCD wave function, whose value for $L_{0}=0$ is -17.789. We investigate the same projected energy by the Monte Carlo method. As the MC calculations have statistical errors, we examine the convergence by taking several sampling numbers $N_{0}=$ 30, 100, 300, 1000, 3000, 10000, 30000. For each $N_{0}$, we carry out 100 sets with different random numbers. In Fig.3, the MC energy in Eq.(14) for each calculation is displayed. The average energy and its average variance over 100 sets are also plotted. Thus, the statistical errors of the Monte Carlo procedure are well-controlled. The sampling number is, in general, relatively constant for various quantum systems with a larger Hilbert space.

Finally, we discuss the remaining problem: the computation of the projected overlaps $\langle s|P^{L_{0}}|\psi)$. For the MC sampling in Eq.(12), the spherical basis is better than the deformed one due to the advantage of the sparsity in $h_{s,s^{\prime}}$. However, the projected overlap includes the exponential-type operators, of which exact projections become intractable for larger systems in general. Therefore we decompose the computation of the projected overlap as in Eq.(15), which allows us the following natural truncation of the deformed complete set. As the $|d)$ states are the ground state $|0)$ and all excited states thereof, the overlaps $\langle s|P^{L_{0}}|\psi)$ for highly excited states $|d)$ are expected to be negligible, which then can be truncated by the value of $(d|\psi)$. As the CCD wave function $|\psi)$ has no dependence of $\theta$ in Eq.(18) and is expressed by the deformed operators, the numerical evaluation of the values of $(d|\psi)$ is simple. Moreover, contributions of higher excited states $|d)$ are, in general, expected to be smaller. Therefore, the truncation of $d$ by $(d|\psi)$ is feasible.

In Fig.4, we present a benchmark test of such a truncation scheme for the states with $L_{0}=0\sim 3$. By setting a given threshold for the value $(d|\psi)$, we can truncate the summation for $d$ in Eq.(15). The MC energies in Eq.(14) with these projected states are plotted as a function of the ratio of the truncated basis number to the one of the whole Hilbert space dimension. The parameters for the three-level Lipkin model are $\chi=5$ and $N=20$. The calculation at the ratio of 0.02 shows that the truncation scheme works quite well. This ratio is expected to be increasingly smaller for larger systems. Its investigation for various quantum systems will, however, be our task for the future.