## 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

### 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

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

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

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

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

The equation for the pairing tensor is also extended such that

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

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

## 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

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.

### 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].

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.

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

## 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].

#### Terms in Eq. (16)

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

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