## Abstract

Pairing properties of a trapped fermion gas consisting of atoms with magnetic dipole moment are studied using the Hartree–Fock–Bogoliubov theory (HFB) and the results are compared with those obtained from the exact diagonalization approach (EDA). It is shown that HFB underestimates two-body correlations and gives quite different atom–atom correlations than the results in EDA. Extended time-dependent HFB (ETDHFB), which is derived from the truncation of equations of motion for reduced density matrices, is also applied and it is found that the difficulties of HFB can be remedied by ETDHFB.

## Introduction

Degenerate Fermi gases of $$^{161}$$Dy [1] and $$^{167}$$Er [2] have recently been achieved following the realization of Bose–Einstein condensation of $$^{164}$$Dy [3]. These isotopes have large magnetic moments: $$^{161}$$Dy has 10$$\mu _B$$ and $$^{167}$$Er 7$$\mu _B$$ ($$\mu _B$$ being the Bohr magneton). The progress of these experiments provides an opportunity to study exotic many-body physics with magnetic dipolar moments. One interesting subject is to investigate how superfluid phases are realized in dipolar fermion gases [4]. The standard approach to studying pairing correlations in fermion systems is the Bardeen–Cooper–Schrieffer (BCS) theory, or its generalization, the Hartree–Fock–Bogoliubov theory (HFB) [5]. It has been shown [6, 7] that BCS predicts a $$p$$-wave spin triplet Cooper pairing in magnetic dipolar fermion gases. BCS and HFB have also been used to study the $$p$$-wave pairing in neutron stars and neutron matter [8–10]. The dipole–dipole interaction is a spatially anisotropic and long-range interaction, which is quite different from the attractive contact interaction to which BCS and HFB are usually applied. Therefore, it is desirable to know how accurately BCS or HFB can describe atom–atom correlations induced by the dipole–dipole interaction. In this paper we apply HFB to a trapped magnetic dipolar fermion gas and compare the results with those from the exact diagonalization approach (EDA). We are restricted to a system consisting of a rather small number of atoms. Such systems have often been used for theoretical investigations of dipolar Fermi gases [11] and may be realized in an array of microtraps or optical lattices, as discussed in Refs. [11–13]. It is shown that HFB underestimates atom–atom correlations and gives a pair correlation function which differs from the result in EDA. We also use an extended time-dependent Hartree–Fock–Bogoliubov theory (ETDHFB) which includes higher-order effects neglected in the time-dependent Hartree–Fock–Bogoliubov theory (TDHFB) and HFB. ETDHFB has been formulated [14] using a truncation scheme similar to that used in the time-dependent density-matrix theory (TDDM) in the normal-fluid regime [15, 16], where higher-order reduced density matrices are approximated by lower-order density matrices to truncate the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) hierarchy for reduced density matrices. TDDM has in the past demonstrated its effectiveness in various applications [16–18], and it can reasonably be assumed that its extension to the superfluid case will show equally good performance. The advantages of ETDHFB are that it has a direct connection to TDHFB and that various correction terms are expressed explicitly, contrary to EDA and Monte Carlo approaches. We show that the difficulties of HFB can be removed by using ETDHFB. The paper is organized as follows: the formulation is given in Sect. 2, the results are shown in Sect. 3, and Sect. 4 is devoted to a summary.

## Formulation

### Hamiltonian

We consider a magnetic dipolar gas of fermions with spin one half, which is trapped in a spherically symmetric harmonic potential with frequency $$\omega$$. The system is described by the Hamiltonian

(1)
$$H=\sum_\alpha\epsilon_\alpha a^\dagger_\alpha a_\alpha +\frac{1}{2}\sum_{\alpha\beta\alpha'\beta'}\left\langle\alpha\beta|v|\alpha'\beta'\right\rangle a^\dagger_{\alpha}a^\dagger_\beta a_{\beta'}a_{\alpha'},$$
where $$a^\dagger _\alpha$$ and $$a_\alpha$$ are the creation and annihilation operators of an atom at a harmonic oscillator state $$\alpha$$ corresponding to the trapping potential $$V(r)=m\omega ^2r^2/2$$ and $$\epsilon _\alpha =\omega (n+3/2)$$ with $$n=0,1,2, \ldots$$. We use units such that $$\hbar =1$$ and assume that $$\alpha$$ contains the spin quantum number. In Eq. (1), $$\langle \alpha \beta |v|\alpha '\beta '\rangle$$ is the matrix element of a pure magnetic dipole–dipole interaction [19]
(2)
$$v({\boldsymbol r})=-\frac{1}{r^3}\left(3\left({\boldsymbol d}_1\cdot\hat{{\boldsymbol r}}\right)\left({\boldsymbol d}_2\cdot\hat{{\boldsymbol r}}\right) -{\boldsymbol d}_1\cdot{\boldsymbol d}_2\right),$$
where $${\boldsymbol d}$$ is the magnetic dipole moment, $${\boldsymbol r}={\boldsymbol r}_1-{\boldsymbol r}_2$$, and $$\hat {{\boldsymbol r}}={\boldsymbol r}/r$$. The magnetic dipole moment for spin 1/2 is given by $${\boldsymbol d}=d\boldsymbol {\sigma }$$, where $$\boldsymbol {\sigma }$$ is a Pauli matrix. The dipole–dipole interaction contains a contact term $$-8\pi {\boldsymbol d}_1\cdot {\boldsymbol d}_2\delta ^3\left ({\boldsymbol r}\right )/3$$ [19], which is usually omitted in the study of dipolar gases. In this study we also neglect it, assuming that its effect could be canceled by the contact interaction $$g\delta ^3\left ({\boldsymbol r}\right )$$, which is usually additionally included in the study of cold atoms. The dipole–dipole interaction Eq. (2) is attractive in the $$^3P_1$$ channel [8] and can induce the $$p$$-wave pairing correlation. To compare the $$p$$-wave pairing with the $$s$$-wave pairing, we also consider an attractive contact interaction $$g\delta ^3\left ({\boldsymbol r}\right )$$.

### Time-dependent Hartree–Fock–Bogoliubov equations

We first consider the equation of motion for the density matrix $${n}_{\alpha \alpha '}$$, which is defined as $${n}_{\alpha \alpha '}=\left \langle \Phi (t)\left |a^\dagger _{\alpha '}a_{\alpha }\right |\Phi (t)\right \rangle$$. Here, $$\left |\Phi (t)\right \rangle$$ is the time-dependent total wavefunction $$\left |\Phi (t)\right \rangle =\exp \left (-i\left (H-\mu \hat {N}\right )\right )\left |\Phi (0)\right \rangle$$, where $$\hat {N}$$ is the number operator and $$\mu$$ is the chemical potential. In the equation of motion for the density matrix $$i \dot {n}_{\alpha \alpha '}=\left \langle \Phi (t)\left |\left [a^\dagger _{\alpha '}a_{\alpha },H-\mu \hat {N}\right ]\right |\Phi (t)\right \rangle$$, there appears a two-body density matrix $$\rho _{\alpha \beta \alpha '\beta '}=\left \langle \Phi (t)\left |a^\dagger _{\alpha '}a^\dagger _{\beta '}a_{\beta }a_{\alpha }\right |\Phi (t)\right \rangle$$. In the derivation of the TDHFB equations it is decomposed as

(3)
$$\rho_{\alpha\beta\alpha'\beta'}\approx n_{\alpha\alpha'}n_{\beta\beta'}-n_{\alpha\beta'}n_{\beta\alpha'} +\kappa_{\alpha\beta}\kappa^*_{\alpha'\beta'}.$$
Here, $$\kappa _{\alpha \beta }$$ is the pairing tensor given by $$\kappa _{\alpha \beta }=\left \langle \Phi (t)\left |a_{\beta }a_{\alpha }\right |\Phi (t)\right \rangle$$. Similarly, in the equation of motion for the pairing tensor $$i\dot {\kappa }_{\alpha \beta }=\left \langle \Phi (t)\left |\left [a_{\beta }a_{\alpha },H-\mu \hat {N}\right ]\right |\Phi (t)\right \rangle$$, there appears a matrix given by $$\left \langle \Phi (t)\left |a^\dagger _{\alpha '}a_{\gamma }a_{\beta }a_{\alpha }\right |\Phi (t)\right \rangle$$. In the derivation of the TDHFB equations it is decomposed as
(4)
$$\left\langle\Phi(t)\left|a^\dagger_{\alpha'}a_{\gamma}a_{\beta}a_{\alpha}\right|\Phi(t)\right\rangle \approx n_{\gamma\alpha'}\kappa_{\alpha\beta} -n_{\beta\alpha'}\kappa_{\alpha\gamma} +n_{\alpha\alpha'}\kappa_{\beta\gamma}.$$

The TDHFB equation for the density matrix $$n_{\alpha \alpha '}$$ is given by

(5)
$$i \dot{n}_{\alpha\alpha'}= \sum_{\lambda}\Big(\epsilon_{\alpha\lambda}{n}_{\lambda\alpha'}-{n}_{\alpha\lambda}\epsilon_{\lambda\alpha'}\Big) +\sum_{\lambda} \Big(\Delta_{\alpha\lambda}\kappa^*_{\alpha'\lambda}-\Delta^*_{\alpha'\lambda}\kappa_{\alpha\lambda}\Big),$$
where $$\epsilon _{\alpha \alpha '}$$ is defined as
(6)
$$\epsilon_{\alpha\alpha'}=\epsilon_{\alpha}\delta_{\alpha\alpha'} +\sum_{\lambda_1\lambda_2} \left\langle\alpha\lambda_1|v|\alpha'\lambda_2\right\rangle_A n_{\lambda_2\lambda_1},$$
and the pairing potential $$\Delta _{\alpha \beta }$$ by
(7)
$$\Delta_{\alpha\beta}=\frac{1}{2}\sum_{\lambda_1\lambda_2}\langle\alpha\beta|v|\lambda_1\lambda_2\rangle_A\kappa_{\lambda_1\lambda_2}.$$
Here, the subscript $$A$$ means that the corresponding matrix is antisymmetrized. The TDHFB equation for the pairing tensor is given by
(8)
$$i \dot{\kappa}_{\alpha\beta}=\sum_{\lambda}\left(\tilde{\epsilon}_{\alpha\lambda}{\kappa}_{\lambda\beta}+\tilde{\epsilon}_{\beta\lambda}{\kappa}_{\alpha\lambda}\right) +\Delta_{\alpha\beta} +\sum_{\lambda} \left(\Delta_{\beta\lambda}n_{\alpha\lambda}-\Delta_{\alpha\lambda}n_{\beta\lambda}\right),$$
where $$\tilde {\epsilon }_{\alpha \alpha '}$$ includes the chemical potential $$\mu$$ as $$\tilde {\epsilon }_{\alpha \alpha '}=\epsilon _{\alpha \alpha '}-\mu \delta _{\alpha \alpha '}$$.

The ground state in HFB is given as a stationary solution of the TDHFB equations [Eqs. (5) and (8)]:

(9)
\begin{align} 0&=\sum_{\lambda}\left(\epsilon_{\alpha\lambda}{n}_{\lambda\alpha'}-{n}_{\alpha\lambda}\epsilon_{\lambda\alpha'}\right) +\sum_{\lambda} \left(\Delta_{\alpha\lambda}\kappa^*_{\alpha'\lambda}-\Delta^*_{\alpha'\lambda}\kappa_{\alpha\lambda}\right), \end{align}

(10)
\begin{align} 0&=\sum_{\lambda}\left(\tilde{\epsilon}_{\alpha\lambda}{\kappa}_{\lambda\beta}+\tilde{\epsilon}_{\beta\lambda}{\kappa}_{\alpha\lambda}\right) +\Delta_{\alpha\beta} +\sum_{\lambda} \left(\Delta_{\beta\lambda}n_{\alpha\lambda}-\Delta_{\alpha\lambda}n_{\beta\lambda}\right). \end{align}
The solution of Eqs. (9) and (10) can be obtained using the usual iteration method.

### Extension of time-dependent Hartree–Fock–Bogoliubov equations

The equations in ETDHFB can be derived by extending the decompositions as

(11)
$$\rho_{\alpha\beta\alpha'\beta'}=n_{\alpha\alpha'}n_{\beta\beta'}-n_{\alpha\beta'}n_{\beta\alpha'} +\kappa_{\alpha\beta}\kappa^*_{\alpha'\beta'} +{\mathcal C}_{\alpha\beta\alpha'\beta'}$$
and
(12)
$$\left\langle\Phi(t)\left|a^\dagger_{\alpha'}a_{\gamma}a_{\beta}a_{\alpha}\right|\Phi(t)\right\rangle =n_{\gamma\alpha'}\kappa_{\alpha\beta} -n_{\beta\alpha'}\kappa_{\alpha\gamma} +n_{\alpha\alpha'}\kappa_{\beta\gamma}+K_{\alpha\beta\gamma:\alpha'}.$$
The matrices $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ and $$K_{\alpha \beta \gamma :\alpha '}$$ describe higher-order effects. The equation for the density matrix is now extended as
(13)
\begin{align} i \dot{n}_{\alpha\alpha'}&= \sum_{\lambda}\left(\epsilon_{\alpha\lambda}{n}_{\lambda\alpha'}-{n}_{\alpha\lambda}\epsilon_{\lambda\alpha'}\right) +\sum_{\lambda} \left(\Delta_{\alpha\lambda}\kappa^*_{\alpha'\lambda}-\Delta^*_{\alpha'\lambda}\kappa_{\alpha\lambda}\right),\nonumber \\ &\quad + \sum_{\lambda_1\lambda_2\lambda_3} \left[\left\langle\alpha\lambda_1|v|\lambda_2\lambda_3\right\rangle {\mathcal C}_{\lambda_2\lambda_3\alpha'\lambda_1} -{\mathcal C}_{\alpha\lambda_1\lambda_2\lambda_3}\left\langle\lambda_2\lambda_3|v|\alpha'\lambda_1\right\rangle\right]. \end{align}
The equation of motion for $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ is given by [14]
(14)
\begin{align} i\dot{\mathcal C}_{\alpha\beta\alpha'\beta'}&= \sum_{\lambda}\left(\epsilon_{\alpha\lambda}{\mathcal C}_{\lambda\beta\alpha'\beta'} +\epsilon_{\beta\lambda}{\mathcal C}_{\alpha\lambda\alpha'\beta'} -\epsilon_{\lambda\alpha'}{\mathcal C}_{\alpha\beta\lambda\beta'} -\epsilon_{\lambda\beta'}{\mathcal C}_{\alpha\beta\alpha'\lambda}\right)\nonumber \\ &\quad + B_{\alpha\beta\alpha'\beta'} +P_{\alpha\beta\alpha'\beta'}+H_{\alpha\beta\alpha'\beta'} +S_{\alpha\beta\alpha'\beta'}+T_{\alpha\beta\alpha'\beta'}. \end{align}
In order to close the coupled chain of equations of motion, we approximated the matrix $$\left \langle \Phi (t)\left |a^\dagger _{\alpha '}a^\dagger _{\beta '}a^\dagger _{\gamma '}a_{\gamma }a_{\beta }a_{\alpha }\right |\Phi (t)\right \rangle$$ by antisymmetrized product combinations of $$n_{\alpha \alpha '}$$, $$\kappa _{\alpha \beta }$$, $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$, and $$K_{\alpha \beta \gamma :\alpha '}$$ such as $$n_{\alpha \alpha '}n_{\beta \beta '}n_{\gamma \gamma '}$$, $$n_{\alpha \alpha '}C_{\beta \gamma \beta '\gamma '}$$, $$n_{\alpha \alpha '}\kappa _{\beta \gamma }\kappa ^*_{\beta '\gamma '}$$, $$\kappa _{\alpha \beta }K^*_{\alpha '\beta '\gamma ':\gamma }$$, and $$\kappa ^*_{\alpha '\beta '}K_{\alpha \beta \gamma :\gamma '}$$. In Eq. (14), $$B_{\alpha \beta \alpha '\beta '}$$ describes the two particle (2p)–two hole (2h) and 2h–2p excitations, $$P_{\alpha \beta \alpha '\beta '}$$ the p–p (and h–h) correlations which are not included in the pairing tensor, and $$H_{\alpha \beta \alpha '\beta '}$$ the p–h correlations. The terms in $$S_{\alpha \beta \alpha '\beta '}$$ and $$T_{\alpha \beta \alpha '\beta '}$$ express the coupling to $$\kappa _{\alpha \beta }$$ and $$K_{\alpha \beta \gamma :\alpha '}$$, respectively. The expressions for the matrices in Eq. (14) are given in Ref. [14], but shown in Appendix A for completeness. The equation for $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ without $$S_{\alpha \beta \alpha '\beta '}$$ and $$T_{\alpha \beta \alpha '\beta '}$$ is the same as that in TDDM [16]. Since the total wavefunction $$\left |\Phi (t)\right \rangle$$ is not an eigenstate of the number operator, the couplings to $$\kappa _{\alpha \beta }$$ and $$K_{\alpha \beta \gamma :\alpha '}$$ appear in Eq. (14).

The equation for the pairing tensor is also extended such that

(15)
\begin{align} i\dot{\kappa}_{\alpha\beta}&=\sum_{\lambda}\left(\tilde{\epsilon}_{\alpha\lambda}{\kappa}_{\lambda\beta}+\tilde{\epsilon}_{\beta\lambda}{\kappa}_{\alpha\lambda}\right) +\Delta_{\alpha\beta} +\sum_{\lambda} \left(\Delta_{\beta\lambda}n_{\alpha\lambda}-\Delta_{\alpha\lambda}n_{\beta\lambda}\right)\nonumber \\ &\quad -\sum_{\lambda_1\lambda_2\lambda_3}\left(\left\langle\alpha\lambda_1|v|\lambda_2\lambda_3\right\rangle K_{\beta\lambda_2\lambda_3:\lambda_1} +\left\langle\lambda_1\beta|v|\lambda_2\lambda_3\right\rangle K_{\alpha\lambda_2\lambda_3:\lambda_1}\right). \end{align}
The equation for $$K_{\alpha \beta \gamma :\alpha '}$$ is written as [14]
(16)
\begin{align} i\dot{K}_{\alpha\beta\gamma:\alpha'}&=\sum_{\lambda}\left(\tilde{\epsilon}_{\alpha\lambda}{K}_{\lambda\beta\gamma:\alpha'}+\tilde{\epsilon}_{\beta\lambda}{K}_{\alpha\lambda\gamma:\alpha'} +\tilde{\epsilon}_{\gamma\lambda}{K}_{\alpha\beta\lambda:\alpha'}-\tilde{\epsilon}_{\lambda\alpha'}{K}_{\alpha\beta\gamma:\lambda}\right)\nonumber \\ &\quad +{D_{\alpha\beta\gamma:\alpha'}} + E_{\alpha\beta\gamma:\alpha'}, + F_{\alpha\beta\gamma:\alpha'} + G_{\alpha\beta\gamma:\alpha'}. \end{align}
We approximated the matrix $$\left \langle \Phi (t)\left |a^\dagger _{\alpha '}a^\dagger _{\beta '}a_{\delta }a_{\gamma }a_{\beta }a_{\alpha }\right |\Phi (t)\right \rangle$$ by antisymmetrized product combinations of $$n_{\alpha \alpha '}$$, $$\kappa _{\alpha \beta }$$, $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$, and $$K_{\alpha \beta \gamma :\alpha '}$$. The terms in $$D_{\alpha \beta \gamma :\alpha '}$$ and $$E_{\alpha \beta \gamma :\alpha '}$$ describe the coupling to the pairing tensor and to the product of three pairing tensors, respectively. The terms in $$F_{\alpha \beta \gamma :\alpha '}$$ describe correlations involving $$K_{\alpha \beta \gamma :\alpha '}$$. The coupling to $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ is contained in $$G_{\alpha \beta \gamma :\alpha '}$$. The matrices in Eq. (16) are given in Appendix B. Equations (13) and (15) may be written in matrix form as in TDHFB:
(17)
$$i\dot{\mathcal R}-\left[{\mathcal H},{\mathcal R}\right]=\left[{\mathcal V},{\mathcal K}\right],$$
where, in obvious notation,
(18)
\begin{align} {\mathcal R}&=\left( \begin{matrix} n& \kappa\\ -\kappa^*& 1-n^* \end{matrix} \right), \end{align}

(19)
\begin{align} {\mathcal H}&=\left( \begin{matrix} \epsilon& \Delta\\ -\Delta^*& -\epsilon^* \end{matrix}\right), \end{align}

(20)
\begin{align} {\mathcal K}&=\left(\begin{matrix} {\mathcal C}& K\\ -K^*& -{\mathcal C}^* \end{matrix}\right), \end{align}

(21)
\begin{align} {\mathcal V}&=\left( \begin{matrix} v& 0\\ 0& -v^* \end{matrix}\right). \end{align}
The ETDHFB equation Eq. (13) conserves on average the total number of particles $$N=\sum _\alpha n_{\alpha \alpha }$$, as is easily shown by taking the trace of Eq. (13). Since our formalism only deals with small number fluctuations around the HFB ground state, ETDHFB cannot recover number symmetry broken by HFB and TDHFB though higher-order effects are included, as will be shown below.

The total energy $$E_{\rm tot}$$ is expressed as

(22)
$$E_{\rm tot}=\sum_{\alpha}\epsilon_\alpha n_{\alpha\alpha} +\frac{1}{2}\sum_{\alpha\beta\alpha'\beta'}\left\langle\alpha\beta|v|\alpha'\beta'\right\rangle \rho_{\alpha'\beta'\alpha\beta},$$
which can be divided into the mean-field energy $$E_{\rm MF}$$, the pairing energy $$E_{\rm pair}$$, and the correlation energy $$E_{\rm cor}$$, given by
(23)
\begin{align} E_{\rm MF}&=\sum_{\alpha}\epsilon_\alpha n_{\alpha\alpha} +\frac{1}{2}\sum_{\alpha\beta\alpha'\beta'}\left\langle\alpha\beta|v|\alpha'\beta'\right\rangle_A n_{\alpha'\alpha}n_{\beta'\beta}, \end{align}

(24)
\begin{align} E_{\rm pair}&=\frac{1}{2}\sum_{\alpha\beta}\Delta_{\alpha\beta}\kappa^*_{\alpha\beta}, \end{align}

(25)
\begin{align} E_{\rm cor}&=\frac{1}{2}\sum_{\alpha\beta\alpha'\beta'}\left\langle\alpha\beta|v|\alpha'\beta'\right\rangle {\mathcal C}_{\alpha'\beta'\alpha\beta}. \end{align}
To conserve $$E_{\rm tot}$$, we need all ETDHFB equations: Eqs. (13), (14), (15), and (16).

In principle the ground state in ETDHFB is given as a stationary solution of the ETDHFB equations. We obtain an approximate ground state using an adiabatic method [20]: Starting from the HFB ground state, we solve the ETDHFB equations with gradually increasing time-dependent residual interaction such that $$v({\boldsymbol r},t)=v({\boldsymbol r})\times t/T$$. This method is motivated by the Gell-Mann–Low theorem [21], which has often been used as a practical way to obtain an approximately stationary solution of time-dependent equations [17]. To suppress spurious oscillations, we need to use large $$T$$. We found that $$T=4\times 2\pi /\omega$$ is sufficiently large, and this value of $$T$$ is used in the following numerical calculations. Of course, the adiabatic method cannot be applied to the cases where $$\omega \approx 0$$ or the interaction is strong. Since the chemical potential $$\mu$$ is not adjusted at each time step to be consistent with $$n_{\alpha \alpha '}$$, $$\kappa _{\alpha \beta }$$ and $$K_{\alpha \beta \gamma :\alpha '}$$ gain a phase given by $$\exp [\Delta \mu t]$$, where $$\Delta \mu$$ is the deviation of the chemical potential from the HFB value.

### Perturbative expressions

To better understand the correlations included in ETDHFB, we present the perturbative expressions for $${\mathcal C_{\alpha \beta \alpha '\beta '}}$$. In Eq. (14) the terms in $$P_{\alpha \beta \alpha '\beta '}$$ and $$H_{\alpha \beta \alpha '\beta '}$$ contain $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$, and $$T_{\alpha \beta \alpha '\beta '}$$ includes $$K_{\alpha \beta \gamma :\alpha '}$$. Therefore, the lowest-order corrections are from $$B_{\alpha \beta \alpha '\beta '}$$ and $$S_{\alpha \beta \alpha '\beta '}$$. The perturbative expression for $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ obtained from $$B_{\alpha \beta \alpha '\beta '}$$ is written as

(26)
$${\mathcal C}_{1\alpha\beta\alpha'\beta'}=-\frac{\left\langle\alpha\beta|v|\alpha'\beta'\right\rangle_A} {\epsilon_\alpha+\epsilon_\beta-\epsilon_{\alpha'}-\epsilon_{\beta'}} \left(\bar{n}^0_\alpha\bar{n}^0_\beta n^0_{\alpha'}n^0_{\beta'}-{n}^0_\alpha{n}^0_\beta \bar{n}^0_{\alpha'}\bar{n}^0_{\beta'}\right).$$
Here, $$n^0_{\alpha }$$ is the occupation probability in HFB and $$\bar {n}^0_\alpha =1-n^0_{\alpha }$$. For simplicity we assume that the single-particle energy $$\epsilon _{\alpha \alpha '}$$ and the occupation matrix are diagonal. $${\mathcal C}_{1\alpha \beta \alpha '\beta '}$$ describes the 2p–2h and 2h–2p excitations. The perturbative expression for $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ obtained from $$S_{\alpha \beta \alpha '\beta '}$$ in Eq. (14) is given by
(27)
\begin{align} {\mathcal C}_{2\alpha\beta\alpha'\beta'}&=\frac{1} {\epsilon_\alpha+\epsilon_\beta-\epsilon_{\alpha'}-\epsilon_{\beta'}} \sum_{\lambda\lambda'}\left[ \left\langle\alpha\lambda'|v|\alpha'\lambda\right\rangle_A\kappa^0_{\beta\lambda}\kappa^{0*}_{\beta'\lambda'}\left(n^0_\alpha-n^0_{\alpha'}\right)\right.\nonumber \\ &\quad +\left\langle\beta{\lambda'}|v|\beta'\lambda\right\rangle_A\kappa^0_{\alpha\lambda}\kappa^{0*}_{\alpha'\lambda'}\left(n^0_\beta-n^0_{\beta'}\right)- \left\langle\alpha\lambda'|v|\beta'\lambda\right\rangle_A\kappa^0_{\beta\lambda}\kappa^{0*}_{\alpha'\lambda'}\left(n^0_\alpha-n^0_{\beta'}\right)\nonumber \\ &\quad - \left.\left\langle\beta\lambda'|v|\alpha'\lambda\right\rangle_A\kappa^0_{\alpha\lambda}\kappa^{0*}_{\beta'\lambda'}\left(n^0_\beta-n^0_{\alpha'}\right)\right]. \end{align}
The correlation energy obtained from $${\mathcal C}_{1\alpha \beta \alpha '\beta '}$$ is interpreted as a correction to the mean-field energy $$E_{\rm MF}$$, whereas that from $${\mathcal C}_{2\alpha \beta \alpha '\beta '}$$ is a correction to the pairing energy $$E_{\rm pair}$$.

## Results

Since our purpose is to make a qualitative analysis of HFB and ETDHFB by comparing with EDA, we consider a small system of $$N=6$$ and use the minimum number of single-particle states which are necessary to calculate higher-order effects: We use the $$1s$$, $$1p$$, $$2s$$, and $$1d$$ states for the HFB and ETDHFB calculations and also for the exact diagonalization. The non-interacting configuration consists of the fully occupied $$1s$$ state and partially occupied $$1p$$ states.

### Contact interaction

We first consider a non-dipolar case where atoms interact with an attractive contact interaction $$g\delta ^3\left ({\boldsymbol r}_1-{\boldsymbol r}_2\right )$$. The ground-state energy $$E_{\rm tot}$$ and pairing energy $$E_{\rm pair}$$ calculated in HFB are shown in Figs. 1 and 2 with the dashed line as a function of $$C=|g|/\omega \xi ^3$$, where $$\xi$$ is the oscillator length $$\xi =\sqrt {1/m\omega }$$. In the case of nuclei for which $$\omega \approx 10$$ MeV is applied, $$C=5$$ corresponds to $$|g|\approx 400$$ MeV fm$$^3$$, which is similar to the strength of nuclear pairing interactions commonly used for a small single-particle space [22]. The results in ETDHFB are shown with the squares. The dot-dashed lines in Figs. 1 and 2 depict the results in EDA. Though HFB slightly underestimates the correlation effects, the ground-state energies calculated in HFB and ETDHF agree well with those in EDA, while the pairing energies in HFB and, to lesser extent, those in ETDHFB deviate from the values in EDA. This indicates that there is a cancellation [22] between the errors in the mean-field energy and $$E_{\rm pair}$$ in HFB or $$E_{\rm cor}+E_{\rm pair}$$ in ETDHFB. The accuracy of HFB seems to depend on the systems considered. For example, the authors of Ref. [22] show that in the case of tin isotopes HFB underestimates the correlation energy by a few tens percent. As shown in Fig. 2, $$|E_{\rm pair}|$$ in ETDHFB is always smaller than that in HFB. This may be interpreted as the pairing correlation being screened by higher-order effects. It may also imply that the number symmetry broken by HFB is partially recovered in ETDHFB: for example, at $$C=6$$ the variance $$\sigma ^2=\langle \hat {N}^2\rangle -\langle \hat {N}\rangle ^2$$, which is expressed by $$\rho _{\alpha \beta \alpha '\beta '}$$ [Eq. (11)], is $$3.0$$ in HFB and $$2.4$$ in ETDHFB.

To investigate the intrinsic structure of the ground state, we consider the pair correlation function $$C\left ({\boldsymbol r}s{\boldsymbol r}'s':{\boldsymbol r}s{\boldsymbol r}'s'\right )$$ given by

(28)
$$C\left({\boldsymbol r}s{\boldsymbol r}'s':{\boldsymbol r}s{\boldsymbol r}'s'\right) =\sum_{\alpha\beta\alpha'\beta'}C_{\alpha(s)\beta(s')\alpha'(s)\beta'(s)} \phi_{\alpha'}^*\left({\boldsymbol r}\right)\phi^*_{\beta'}\left({\boldsymbol r}'\right)\phi_{\alpha}\left({\boldsymbol r}\right)\phi_{\beta}\left({\boldsymbol r}'\right).$$
Here, $$s$$ is the $$z$$ component of a spin coordinate and $$\phi _\alpha$$ the harmonic oscillator wave function. In HFB the two-body correlation matrix is given by
(29)
$$C_{\alpha(s)\beta(s')\alpha'(s)\beta'(s)}=\kappa_{\alpha(s)\beta(s')}\kappa^*_{\alpha'(s)\beta'(s')},$$
and $$C\left ({\boldsymbol r}s{\boldsymbol r}'s':{\boldsymbol r}s{\boldsymbol r}'s'\right )$$ is written as
(30)
$$C\left({\boldsymbol r}s{\boldsymbol r}'s':{\boldsymbol r}s{\boldsymbol r}'s'\right)=\left|\kappa\left({\boldsymbol r}s{\boldsymbol r}'s'\right)\right|^2,$$
where $$\kappa \left ({\boldsymbol r}s{\boldsymbol r}'s'\right )$$ is defined as
(31)
$$\kappa({\boldsymbol r}s{\boldsymbol r}'s') =\sum_{\alpha\beta}\kappa_{\alpha(s)\beta(s')}\phi_\alpha({\boldsymbol r})\phi_\beta({\boldsymbol r}').$$
In ETDHFB the correlated part of the two-body density matrix is given as
(32)
$$C_{\alpha(s)\beta(s')\alpha'(s)\beta'(s)}=\kappa_{\alpha(s)\beta(s')}\kappa^*_{\alpha'(s)\beta'(s')} +{\mathcal C}_{\alpha(s)\beta(s')\alpha'(s)\beta'(s)}.$$

Fig. 1.

Ground-state energy $$E_{\rm tot}$$ in HFB (dashed line) for $$N=6$$ as a function of $$C=|g|/\omega \xi ^3$$. The results in ETDHFB are shown with the squares. The dot-dashed line depicts the results in EDA.

Fig. 1.

Ground-state energy $$E_{\rm tot}$$ in HFB (dashed line) for $$N=6$$ as a function of $$C=|g|/\omega \xi ^3$$. The results in ETDHFB are shown with the squares. The dot-dashed line depicts the results in EDA.

Fig. 2.

Pairing energy $$E_{\rm pair}$$ in HFB (dashed line) for $$N=6$$ as a function of $$C=|g|/\omega \xi ^3$$. The ETDHFB results for $$E_{\rm pair}$$ and $$E_{\rm cor}+E_{\rm pair}$$ are shown with the circles and squares, respectively. The dot-dashed line depicts $$E_{\rm tot}-E_{\rm MF}$$ in EDA.

Fig. 2.

Pairing energy $$E_{\rm pair}$$ in HFB (dashed line) for $$N=6$$ as a function of $$C=|g|/\omega \xi ^3$$. The ETDHFB results for $$E_{\rm pair}$$ and $$E_{\rm cor}+E_{\rm pair}$$ are shown with the circles and squares, respectively. The dot-dashed line depicts $$E_{\rm tot}-E_{\rm MF}$$ in EDA.

Fig. 3.

Contour plot of the pair correlation function $$C({\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow :{\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow )$$ (in arbitrary units) in the $$xy$$ plane calculated in HFB for $$N=6$$ and $$C=10$$, where $${\boldsymbol r}'/\xi =(1,0,0)$$.

Fig. 3.

Contour plot of the pair correlation function $$C({\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow :{\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow )$$ (in arbitrary units) in the $$xy$$ plane calculated in HFB for $$N=6$$ and $$C=10$$, where $${\boldsymbol r}'/\xi =(1,0,0)$$.

The pair correlation function gives the distribution of an atom with spin $$s$$ at $${\boldsymbol r}$$ when the other atom with spin $$s'$$ is located at $${\boldsymbol r}'$$. It describes the non-trivial two-body correlations after the trivial exchange correlations are subtracted. The contour plot of $$C\left ({\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow :{\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow \right )$$ in the $$xy$$ plane calculated in HFB for $$N=6$$ and $$C=10$$ is shown in Fig. 3. The position of $${\boldsymbol r}'$$ is chosen at $${\boldsymbol r}'/\xi =(1,0,0)$$. The result in ETDHFB is similar to that in HFB and is not shown. The two atoms with opposite spin directions favor the same position in the case of the contact interaction. This is because only the $$L=0$$ component (the $$s$$ wave) of the relative motion is allowed in the contact interaction. The result in EDA is shown in Fig. 4. The two distributions agree well except that the result in EDA has regions where the correlation function becomes slightly negative. As shown above, in the case of the attractive contact interaction, HFB gives a good description of the results in EDA and ETDHFB slightly improves the HFB results.

Fig. 4.

As Fig. 3 but for the result in EDA.

Fig. 4.

As Fig. 3 but for the result in EDA.

### Dipole interaction

We now consider the case of the dipole–dipole interaction. The ground-state energy $$E_{\rm tot}$$ and pairing energy $$E_{\rm pair}$$ calculated in HFB are shown in Figs. 5 and 6 with the dashed line as a function of $$\chi =d^2/\omega \xi ^3$$. The results in ETDHFB are shown with squares. The dot-dashed lines in Figs. 5 and 6 depict the results in EDA. The range of $$\chi$$ considered in Figs. 5 and 6 may be rather large for the current experimental situations [1]. For example, $$\chi$$ for a gas of $$^{161}$$Dy trapped in a harmonic potential with $$\omega =2\pi \times 500$$ Hz is about 0.2. A large value of $$\chi$$ may be realized for a dipolar gas confined in a lattice [1, 4].

Fig. 5.

Ground-state energy $$E_{\rm tot}$$ in HFB (dashed line) for $$N=6$$ as a function of $$\chi =d^2/\omega \xi ^3$$. The results in ETDHFB are shown with the squares. The dot-dashed line depicts the results in EDA.

Fig. 5.

Ground-state energy $$E_{\rm tot}$$ in HFB (dashed line) for $$N=6$$ as a function of $$\chi =d^2/\omega \xi ^3$$. The results in ETDHFB are shown with the squares. The dot-dashed line depicts the results in EDA.

Fig. 6.

Pairing energy $$E_{\rm pair}$$ in HFB (dashed line) for $$N=6$$ as a function of $$\chi =d^2/\omega \xi ^3$$. The ETDHFB results for $$E_{\rm pair}$$ and $$E_{\rm cor}+E_{\rm pair}$$ are shown with the circles and squares, respectively. The dot-dashed line depicts $$E_{\rm tot}-E_{\rm MF}$$ in EDA.

Fig. 6.

Pairing energy $$E_{\rm pair}$$ in HFB (dashed line) for $$N=6$$ as a function of $$\chi =d^2/\omega \xi ^3$$. The ETDHFB results for $$E_{\rm pair}$$ and $$E_{\rm cor}+E_{\rm pair}$$ are shown with the circles and squares, respectively. The dot-dashed line depicts $$E_{\rm tot}-E_{\rm MF}$$ in EDA.

The ground-state energies in HFB are larger than those in EDA and the deviation becomes larger with increasing $$\chi$$, while the results in ETDHFB agree well with the EDA results. The situation for $$E_{\rm pair}$$ and $$E_{\rm cor}+E_{\rm pair}$$ is similar, as shown in Fig. 6. It is clear from Figs. 5 and 6 that HFB cannot give sufficient correlations in the case of the dipole–dipole interaction. The largest $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ has combinations of indices such as $$\alpha =1s$$, $$\beta =2d$$, $$\alpha '=1p$$, and $$\beta '=1p$$, which cannot be expressed by the pairing tensors because $$\kappa _{\alpha \beta }$$ does not have the combination $$\alpha =1s$$ and $$\beta =2d$$. This is the reason why HFB underestimates atom–atom correlations. Small $$\left |E_{\rm pair}\right |$$ in ETDHFB may again be interpreted as the screening or partial recovery of the number symmetry broken by HFB. For example, at $$\chi =0.6$$ the variance $$\sigma ^2=\left \langle \hat {N}^2\right \rangle -\left \langle \hat {N}\right \rangle ^2$$ is $$1.9$$ in HFB and $$0.6$$ in ETDHFB. The above conclusion that the $$p$$-wave correlations induced by the dipole–dipole interaction cannot be fully expressed by the pairing tensor in HFB would be applicable for larger-$$N$$ systems. To show this, we performed the ETDHFB calculations for the $$N=16$$ system using the $$1p$$, $$2s$$, $$1d$$, and $$1f$$ states to solve the ETDHFB equations. The non-interacting configuration consists of the fully occupied $$1s$$ and $$1p$$ states and partially occupied $$2s$$ and $$1d$$ states. The results obtained for $$E_{\rm tot}$$ and $$E_{\rm cor}+E_{\rm pair}$$ are compared with the results in HFB in Figs. 7 and 8. The results in EDA are not shown because they are not easily available for this system. As is the case for the $$N=6$$ system, the atom–atom correlations in HFB are significantly smaller than those in ETDHFB.

Fig. 7.

Ground-state energy $$E_{\rm tot}$$ in HFB (triangles) for $$N=16$$ as a function of $$\chi =d^2/\omega \xi ^3$$. The results in ETDHFB are shown with the squares.

Fig. 7.

Ground-state energy $$E_{\rm tot}$$ in HFB (triangles) for $$N=16$$ as a function of $$\chi =d^2/\omega \xi ^3$$. The results in ETDHFB are shown with the squares.

Fig. 8.

Pairing energy $$E_{\rm pair}$$ in HFB (triangles) for $$N=16$$ as a function of $$\chi =d^2/\omega \xi ^3$$. The ETDHFB results for $$E_{\rm pair}$$ and $$E_{\rm cor}+E_{\rm pair}$$ are shown with the circles and squares, respectively.

Fig. 8.

Pairing energy $$E_{\rm pair}$$ in HFB (triangles) for $$N=16$$ as a function of $$\chi =d^2/\omega \xi ^3$$. The ETDHFB results for $$E_{\rm pair}$$ and $$E_{\rm cor}+E_{\rm pair}$$ are shown with the circles and squares, respectively.

Fig. 9.

Contour plot of the two-body correlation function $$C\left ({\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow :{\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow \right )$$ (in arbitrary units) in the $$xy$$ plane calculated in HFB for $$N=6$$ and $$\chi =1$$, where $${{\boldsymbol r}'/\xi }=(1,0,0)$$.

Fig. 9.

Contour plot of the two-body correlation function $$C\left ({\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow :{\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow \right )$$ (in arbitrary units) in the $$xy$$ plane calculated in HFB for $$N=6$$ and $$\chi =1$$, where $${{\boldsymbol r}'/\xi }=(1,0,0)$$.

Fig. 10.

As Fig. 9 but for ETDHFB.

Fig. 10.

As Fig. 9 but for ETDHFB.

The contour plot of $$C\left ({\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow :{\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow \right )$$ in the $$xy$$ plane calculated in HFB for $$N=6$$ and $$\chi =1$$ is shown in Fig. 9. The position of $${\boldsymbol r}'$$ is chosen at $${\boldsymbol r}'/\xi =(1,0,0)$$. The results in EDA and ETDHFB are shown in Figs. 10 and 11, respectively. In the case of the dipole–diple interaction there is no $$s$$-wave contribution and the dominant contribution is from the $$L=1$$ component (the $$p$$ wave) of the relative angular momentum. This is the reason that the spin-up atom is separated from the spin-down atom, as shown in Figs. 9–11. In HFB the dominant contribution is from the pairing tensors with the combinations $$(\alpha ,\beta )=\left (1p_{m=1 \uparrow },1p_{m=-1\downarrow }\right )$$ and $$(\alpha ,\beta )=\left (1p_{m=-1 \uparrow },1p_{m=1 \downarrow }\right )$$, where $$m$$ is the azimuthal quantum number and the arrows indicate spin direction. Since these pairing tensors have opposite sign, the angular distribution looks like $$\sin ^2\phi$$, where $$\phi =\arctan (y/x)$$. The correlation function in EDA becomes strongly negative at $${{\boldsymbol r}'/\xi }=(1,0,0)$$ where an atom with spin down is placed. As a consequence, the two-body density matrix $$\rho \left ({\boldsymbol r}'\uparrow {\boldsymbol r}'\downarrow :{\boldsymbol r}'\uparrow {\boldsymbol r}'\downarrow \right )$$, which is obtained from Eq. (11), becomes approximately 0. The distribution in ETDHFB also becomes negative around $${\boldsymbol r}'$$, though it is not so strong as the EDA result, whereas HFB cannot describe such a behavior of the pair-correlation function around $${\boldsymbol r}'$$. This is consistent with the fact that HFB underestimates the correlation energy. The pair correlation function calculated using the perturbative expression for $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ [Eqs. (26) and (27)] is shown in Fig. 12. The diagonal assumption of $$\epsilon _{\alpha \alpha '}$$ and $$n_{\alpha \alpha '}$$ is not well justified in the case of the dipole–dipole interaction. Therefore, the distribution in Fig. 12 has only a qualitative meaning. Both $${\mathcal C}_{1\alpha \beta \alpha '\beta '}$$ and $${\mathcal C}_{2\alpha \beta \alpha '\beta '}$$ give a quadrupole distribution to $$C\left ({\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow :{\boldsymbol r}\uparrow {\boldsymbol r}'\downarrow \right )$$ in the $$xy$$ plane. The quadrupole distribution of the pair correlation function is determined by $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ with $$\alpha =1p_{m=1}\uparrow$$, $$\beta =1p_{m=-1}\downarrow$$, $$\alpha '=1p_{m=-1}\uparrow$$, and $$\beta '=1p_{m=1}\downarrow$$. Since the difference in the azimuthal quantum number between $$\beta$$ and $$\beta '$$ is 2, the pair correlation function shows the quadrupole distribution. Such a correlation between $$\beta$$ and $$\beta '$$ cannot be expressed by Eq. (30).

Fig. 11.

As Fig. 9 but for EDA.

Fig. 11.

As Fig. 9 but for EDA.

Fig. 12.

As Fig. 9 but for the result obtained from the perturbative expression for $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ [Eqs. (26) and (27)].

Fig. 12.

As Fig. 9 but for the result obtained from the perturbative expression for $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ [Eqs. (26) and (27)].

## Summary

The pairing properties of a trapped fermion gas were studied using the Hartree–Fock–Bogoliubov theory (HFB) and the results were compared with those obtained from the exact diagonalization approach (EDA). It was found that in the case of an attractive contact interaction, the results in HFB agree well with those in EDA, though it is usually anticipated that HFB is not so appropriate for small systems. In the case of a magnetic dipole–dipole interaction, HFB was found to underestimate correlation energy and give quite different atom–atom correlations than the results in EDA. This indicates that the $$p$$-wave correlations induced by the dipole–dipole interaction cannot be fully expressed by the pairing tensor in HFB. This conclusion would be applicable for larger-$$N$$ systems. The extended time-dependent HFB (ETDHFB), which is derived from the truncation of equations of motion for reduced density matrices, was also applied and it was found that the difficulties of HFB can be remedied by ETDHFB. Since the pairing correlations mainly affect the occupation of the single-particle states near the Fermi level, ETDHFB may also be applicable for large-$$N$$ systems by using a limited number of single-particle states near the Fermi level to solve the ETDHFB equations.

### Terms in Eq. (14)

We present the terms in the equations of motion for $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$ and $$K_{\alpha \beta \gamma :\alpha '}$$. Since decomposition of higher-order density matrices to lower-order ones involves various combinations due to the fact that the total wavefunction is not an eigenstate of the number operator, these equations contain many terms. We try to explain the meanings of each term as clearly as possible.

The terms in Eq. (14) are given below. $$B_{\alpha \beta \alpha '\beta '}$$ describes the 2p–2h and 2h–2p excitations as in TDDM [16].

\begin{align*} B_{\alpha\beta\alpha'\beta'}&=\sum_{\lambda_1\lambda_2\lambda_3\lambda_4} \left\langle\lambda_1\lambda_2|v|\lambda_3\lambda_4\right\rangle_A \left[\left(\delta_{\alpha\lambda_1}-n_{\alpha\lambda_1}\right)\left(\delta_{\beta\lambda_2}-n_{\beta\lambda_2}\right) n_{\lambda_3\alpha'}n_{\lambda_4\beta'}\right.\\ &\quad{} - \left. n_{\alpha\lambda_1}n_{\beta\lambda_2}\left(\delta_{\lambda_3\alpha'}-n_{\lambda_3\alpha'}\right) \left(\delta_{\lambda_4\beta'}-n_{\lambda_4\beta'}\right)\right]. \end{align*}
Particle–particle and h–h correlations which are not included in the pairing tensor are described by $$P_{\alpha \beta \alpha '\beta '}$$:
\begin{align*} P_{\alpha\beta\alpha'\beta'}&=\sum_{\lambda_1\lambda_2\lambda_3\lambda_4} \left\langle\lambda_1\lambda_2|v|\lambda_3\lambda_4\right\rangle \left[\left(\delta_{\alpha\lambda_1}\delta_{\beta\lambda_2} -\delta_{\alpha\lambda_1}n_{\beta\lambda_2} -n_{\alpha\lambda_1}\delta_{\beta\lambda_2}\right) {\mathcal C}_{\lambda_3\lambda_4\alpha'\beta'}\right. \\ &\quad-\left.\left(\delta_{\lambda_3\alpha'}\delta_{\lambda_4\beta'} -\delta_{\lambda_3\alpha'}n_{\lambda_4\beta'} -n_{\lambda_3\alpha'}\delta_{\lambda_4\beta'}\right) {\mathcal C}_{\alpha\beta\lambda_1\lambda_2}\right]. \end{align*}
$$H_{\alpha \beta \alpha '\beta '}$$ describes p–h correlations:
\begin{align*} H_{\alpha\beta\alpha'\beta'}&=\sum_{\lambda_1\lambda_2\lambda_3\lambda_4} \left\langle\lambda_1\lambda_2|v|\lambda_3\lambda_4\right\rangle_A \left[\delta_{\alpha\lambda_1}\left(n_{\lambda_3\alpha'}{\mathcal C}_{\lambda_4\beta\lambda_2\beta'} -n_{\lambda_3\beta'}{\mathcal C}_{\lambda_4\beta\lambda_2\alpha'}\right)\right. \nonumber \\ &\quad +\delta_{\beta\lambda_2}\left(n_{\lambda_4\beta'}{\mathcal C}_{\lambda_3\alpha\lambda_1\alpha'} -n_{\lambda_4\alpha'}{\mathcal C}_{\lambda_3\alpha\lambda_1\beta'}\right) -\delta_{\alpha'\lambda_3}\left(n_{\alpha\lambda_1}{\mathcal C}_{\lambda_4\beta\lambda_2\beta'} -n_{\beta\lambda_1}{\mathcal C}_{\lambda_4\alpha\lambda_2\beta'}\right) \nonumber \\ &\quad -\left.\delta_{\beta'\lambda_4}\left(n_{\beta\lambda_2}{\mathcal C}_{\lambda_3\alpha\lambda_1\alpha'} -n_{\alpha\lambda_2}{\mathcal C}_{\lambda_3\beta\lambda_1\alpha'}\right)\right]. \end{align*}
The coupling to the pairing tensor is given by $$S_{\alpha \beta \alpha '\beta '}$$:
\begin{align*} S_{\alpha\beta\alpha'\beta'}&=\sum_{\lambda_1\lambda_2\lambda_3\lambda_4} \left\langle\lambda_1\lambda_2|v|\lambda_3\lambda_4\right\rangle_A \left[\delta_{\alpha\lambda_1}\left(n_{\lambda_3\alpha'}\kappa_{\lambda_4\beta}\kappa^*_{\lambda_2\beta'} -n_{\lambda_3\beta'}\kappa_{\lambda_4\beta}\kappa^*_{\lambda_2\alpha'}\right)\right.\nonumber \\ &\quad+\delta_{\beta\lambda_2}\left(n_{\lambda_4\beta'}\kappa_{\lambda_3\alpha}\kappa^*_{\lambda_1\alpha'} -n_{\lambda_4\alpha'}\kappa_{\lambda_3\alpha}\kappa^*_{\lambda_1\beta'}\right) -\delta_{\alpha'\lambda_3}\left(n_{\alpha\lambda_1}\kappa_{\lambda_4\beta}\kappa^*_{\lambda_2\beta'} -n_{\beta\lambda_1}\kappa_{\lambda_4\alpha}\kappa^*_{\lambda_2\beta'}\right)\nonumber \\ &\quad -\left.\delta_{\beta'\lambda_4}\left(n_{\beta\lambda_2}\kappa_{\lambda_3\alpha}\kappa^*_{\lambda_1\alpha'} -n_{\alpha\lambda_2}\kappa_{\lambda_3\beta}\kappa^*_{\lambda_1\alpha'}\right)\right]. \end{align*}
$$T_{\alpha \beta \alpha '\beta '}$$ expresses the coupling to $$K_{\alpha \beta \gamma :\alpha '}$$:
\begin{align*} T_{\alpha\beta\alpha'\beta'}&=\sum_{\lambda} \left(\Delta_{\alpha\lambda}K^*_{\lambda\beta'\alpha':\beta}-\Delta_{\beta\lambda}K^*_{\lambda\beta'\alpha':\alpha} -\Delta^*_{\alpha'\lambda}K_{\alpha\beta\lambda:\beta'}+\Delta^*_{\beta'\lambda}K_{\alpha\beta\lambda:\alpha'}\right)\\ &\quad +\frac{1}{2}\sum_{\lambda_1\lambda_2\lambda_3\lambda_4}\left\langle\lambda_1\lambda_2|v|\lambda_3\lambda_4\right\rangle_A \left[\delta_{\alpha\lambda_1}\left(2\kappa_{\beta\lambda_4}K^*_{\beta'\lambda_2\alpha':\lambda_3} +\kappa^*_{\beta'\lambda_2}K_{\beta\lambda_4\lambda_3:\alpha'} -\kappa^*_{\alpha'\lambda_2}K_{\beta\lambda_4\lambda_3:\beta'}\right)\right. \nonumber \\ &\quad -\delta_{\beta\lambda_1}\left(2\kappa_{\alpha\lambda_4}K^*_{\beta'\lambda_2\alpha':\lambda_3} +\kappa^*_{\beta'\lambda_2}K_{\alpha\lambda_4\lambda_3:\alpha'} -\kappa^*_{\alpha'\lambda_2}K_{\alpha\lambda_4\lambda_3:\beta'}\right)\nonumber \\ &\quad -\delta_{\alpha'\lambda_3}\left(2\kappa^*_{\lambda_2\beta'}K_{\alpha\lambda_4\beta:\lambda_1} +\kappa_{\alpha\lambda_4}K^*_{\lambda_1\lambda_2\beta':\beta} -\kappa_{\beta\lambda_4}K^*_{\lambda_1\lambda_2\beta':\alpha}\right)\nonumber \\ &\quad +\left.\delta_{\beta'\lambda_3}\left(2\kappa^*_{\lambda_2\alpha'}K_{\alpha\lambda_4\beta:\lambda_1} +\kappa_{\alpha\lambda_4}K^*_{\lambda_1\lambda_2\alpha':\beta} -\kappa_{\beta\lambda_4}K^*_{\lambda_1\lambda_2\alpha':\alpha}\right)\right]. \end{align*}
The terms in the first sum describe the coupling to the pairing potential. The terms in the second sum include both p–p (and h–h) and p–h correlations in infinite order. In the derivation of Eq. (14) we neglected the genuine three-body density matrix $$\left \langle \Phi (t)\left |a^\dagger _{\alpha '}a^\dagger _{\beta '}a^\dagger _{\gamma '}a_{\gamma }a_{\beta }a_{\alpha }\right |\Phi (t)\right \rangle$$ as in TDDM.

#### Terms in Eq. (16)

The terms in Eq. (16) are given below. $$D_{\alpha \beta \gamma :\alpha '}$$ describes the coupling to one pairing tensor:

\begin{align*} D_{\alpha\beta\gamma:\alpha'}&=-\sum_{\lambda_1\lambda_2}\left(\left\langle\alpha\beta|v|\lambda_1\lambda_2\rangle_A{\kappa_{\gamma\lambda_2}}+\langle\beta\gamma|v|\lambda_1\lambda_2\right\rangle_A{\kappa_{\alpha\lambda_2}} -\left\langle\alpha\gamma|v|\lambda_1\lambda_2\right\rangle_A{\kappa_{\beta\lambda_2}}\right){n_{\lambda_1\alpha'}}\nonumber \\ &\quad +\sum_{\lambda_1\lambda_2\lambda_3}\left[\left\langle\alpha\lambda_1|v|\lambda_2\lambda_3\right\rangle_A\left(n_{\beta\lambda_1}{\kappa_{\gamma\lambda_3}} -n_{\gamma\lambda_1}\kappa_{\beta\lambda_3}\right) +\left\langle\beta\lambda_1|v|\lambda_2\lambda_3\right\rangle_A\left(n_{\gamma\lambda_1}{\kappa_{\alpha\lambda_3}} -n_{\alpha\lambda_1}\kappa_{\gamma\lambda_3}\right)\right.\nonumber \\ &\quad +\left.\left\langle\gamma\lambda_1|v|\lambda_2\lambda_3\right\rangle_A\left(n_{\alpha\lambda_1}{\kappa_{\beta\lambda_3}} -n_{\beta\lambda_1}\kappa_{\alpha\lambda_3}\right)\right]{n_{\lambda_2\alpha'}}\nonumber \\ &\quad +\sum_{\lambda_1\lambda_2\lambda_3}\left\langle\lambda_1\lambda_2|v|\alpha'\lambda_3\right\rangle_A \left(n_{\alpha\lambda_1}n_{\gamma\lambda_2}\kappa_{\beta\lambda_3} -n_{\beta\lambda_1}n_{\gamma\lambda_2}\kappa_{\alpha\lambda_3} -n_{\alpha\lambda_1}n_{\beta\lambda_2}\kappa_{\gamma\lambda_3}\right). \end{align*}

The coupling to three pairing tensors is given by $$E_{\alpha \beta \gamma :\alpha '}$$:

\begin{align*} E_{\alpha\beta\gamma:\alpha'}&=-\sum_{\lambda_1\lambda_2\lambda_3}\left(\left\langle\alpha\lambda_1|v|\lambda_2\lambda_3\right\rangle_A\kappa_{\beta\lambda_2}\kappa_{\gamma\lambda_3} -\left\langle\beta\lambda_1|v|\lambda_2\lambda_3\right\rangle_A \kappa_{\alpha\lambda_2}\kappa_{\gamma\lambda_3}\right.\nonumber \\ &\quad +\left.\left\langle\gamma\lambda_1|v|\lambda_2\lambda_3\right\rangle_A\kappa_{\alpha\lambda_2}\kappa_{\beta\lambda_3}\right)\kappa^*_{\alpha'\lambda_1}. \end{align*}
These terms express the modification of the two-particle propagator due to the pairing correlations with other particles. The terms in $$F_{\alpha \beta \gamma :\alpha '}$$ describe correlations among $$K_{\alpha \beta \gamma :\alpha '}$$:
The terms in the first sum describe p–p (and h–h) correlations while those in the second sum are p–h correlations. Since these terms contain $$K_{\alpha \beta \gamma :\alpha '}$$, they describe the coupling to collective p–p, h–h, and p–h correlations. The terms in $$G_{\alpha \beta \gamma :\alpha '}$$ are:
\begin{align*} G_{\alpha\beta\gamma:\alpha'}&= \sum_{\lambda}{\left(\Delta_{\alpha\lambda}{\mathcal C}_{\beta\gamma\alpha'\lambda}-\Delta_{\beta\lambda}{\mathcal C}_{\alpha\gamma\alpha'\lambda}+\Delta_{\gamma\lambda}{\mathcal C}_{\alpha\beta\alpha'\lambda}\right)}\nonumber \\ &\quad -\sum_{\lambda_1\lambda_2\lambda_3\lambda_4}\left\langle\lambda_1\lambda_2|v|\lambda_3\lambda_4\right\rangle_A \left[\left(\delta_{\alpha\lambda_1}\kappa_{\beta\lambda_3}-\delta_{\beta\lambda_1}\kappa_{\alpha\lambda_3}\right){\mathcal C}_{\gamma\lambda_4\alpha'\lambda_2}\right) \nonumber \\ &\quad +\left.\left.\left.\left(\delta_{\beta\lambda_1}\kappa_{\gamma\lambda_3}-\delta_{\gamma\lambda_1}\kappa_{\beta\lambda_3}\right){\mathcal C}_{\alpha\lambda_4\alpha'\lambda_2}\right) -\left(\delta_{\alpha\lambda_1}\kappa_{\gamma\lambda_3}-\delta_{\gamma\lambda_1}\kappa_{\alpha\lambda_3}\right){\mathcal C}_{\beta\lambda_4\alpha'\lambda_2}\right)\right]. \end{align*}
These terms describe the coupling to $${\mathcal C}_{\alpha \beta \alpha '\beta '}$$. In the above derivation of Eq. (16), the genuine correlated matrices $$\left \langle \Phi (t)\left |a^\dagger _{\alpha '}a^\dagger _{\beta '}a_{\delta }a_{\gamma }a_{\beta }a_{\alpha }\right |\Phi (t)\right \rangle$$ and $$\left \langle \Phi (t)\left |a_{\delta }a_{\gamma }a_{\beta }a_{\alpha }\right |\Phi (t)\right \rangle$$ are neglected.

## References

1
Lu
M.
,
Burdick
N. Q.
, and
Lev
B. L.
,
Phys. Rev. Lett.

108
,
215301
(
2012
).
()
2
Aikawa
K.
,
Frisch
A.
,
Mark
M.
,
Baier
S.
,
Grimm
R.
, and
Ferlaino
F.
,
Phys. Rev. Lett.

112
,
010404
(
2014
).
()
3
Lu
M.
,
Burdick
N. Q.
,
Youn
S. H.
, and
Lev
B. L.
,
Phys. Rev. Lett.

107
,
190401
(
2012
).
()
4
Li
Y.
and
Wu
C.
,
J. Phys.: Condens. Matter

26
,
493203
(
2014
).
()
5
Ring
P.
and
Schuck
P
,
The Nuclear Many-Body Problem
,
(Springer-Verlag, Berlin, 1980)
.
6
You
L.
and
Marinescu
M.
,
Phys. Rev. A

60
,
2324
(
1999
).
()
7
Li
Y.
and
Wu
C.
,
Sci. Rep.

2
,
392
(
2012
).
()
8
Tamagaki
R.
,
Prog. Theor. Phys.

44
,
905
(
1970
).
()
9
Baldo
M.
,
Ø.
Elgarøy
,
Engvik
L.
,
Hjorth-Jensen
M.
, and
Schulze
H.-J.
,
Phys. Rev. C

58
,
1921
(
1998
).
()
10
Zuo
W.
,
Cui
C. X.
,
Lombardo
U.
, and
Schulze
H.-J.
,
Phys. Rev. C

78
,
015805
(
2008
).
()
11
Osterloh
K.
, N. Barberán
, and
Lewenstein
M.
,
Phys. Rev. Lett.

99
,
160403
(
2007
).
()
12
N.
Barberán
,
Lewenstein
M.
,
Osterloh
K.
, and
Dagnino
D.
,
Phys. Rev. A

73
,
063623
(
2006
).
()
13
Popp
M.
,
Paredes
B.
, and
Cirac
J. I.
,
Phys. Rev. A

70
,
053612
(
2004
).
()
14
Tohyama
M.
and
Schuck
P.
,
Phys. Rev. C

91
,
034316
(
2015
).
()
15
Wang
S. J.
and
Cassing
W.
,
Ann. Phys.

159
,
328
(
1985
).
()
16
Gong
M.
and
Tohyama
M.
,
Z. Phys. A

335
,
153
(
1990
).
()
17
Pfitzner
A.
,
Cassing
W.
, and
Peter
A.
,
Nucl. Phys. A

577
,
753
(
1994
).
()
18
Tohyama
M.
and
Schuck
P.
,
Eur. Phys. J. A

50
,
77
(
2014
).
()
19
Barger
V. D.
and
Olsson
M. G
,
Classical Electricity and Magnetism
,
(Allyn and Bacon, Boston, 1987)
.
20
Tohyama
M.
,
Phys. Rev. A

71
,
043613
(
2005
).
()
21
Gell-Mann
M.
and
Low
F.
,
Phys. Rev.

84
,
350
(
1951
).
()
22
Sandulescu
N.
and
Bertsch
G. F.
,
Phys. Rev. C

78
,
064318
(
2008
).
()
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.