Many interacting fermions in a one-dimensional harmonic trap: a quantum-chemical treatment

The Hamiltonian describing a system of N spin-1/2 fermions (atoms) in a one-dimensional harmonic trap reads

$$\begin{eqnarray}&&\hat{H}=-\displaystyle \frac{{{\hslash }}^{2}}{2m}\displaystyle \sum _{i=1}^{N}{\partial }_{{x}_{i}}^{2}+\displaystyle \frac{1}{2}\;m{\omega }^{2}\displaystyle \sum _{i=1}^{N}{x}_{i}^{2}+g\displaystyle \sum _{i\lt j}\delta ({x}_{i}-{x}_{j}),\end{eqnarray} \tag{ 1 }$$

where x_i represents the coordinate of the ith atom, m is the atom mass, ω is the frequency of the trap, and g is the strength of the two-body contact interaction. Throughout the paper we use units of energy and interaction strength that correspond to $\omega =m={\hslash }=1.$ This amounts to measuring energies E in units of ${\hslash }\omega ,$ lengths in units of the harmonic oscillator characteristic length ${a}_{\mathrm{ho}}=\sqrt{{\hslash }/(m\omega )},$ and the interaction strength g in units of ${\hslash }\omega {a}_{{\rm{ho}}}.$ Obviously, the Hamiltonian is symmetric with respect to the transposition of two arbitrary space coordinates, implying that the eigenfunctions must transform according to an irreducible representation of the permutation group S_N. In this paper we limit ourselves to the spin-1/2 atoms, so that the wave function will be fully antisymmetric with respect to the simultaneous transposition of any space and spin coordinates of two atoms. From now on we assume that $N={N}_{\uparrow }+{N}_{\downarrow }$ describes the total number of atoms (${N}_{\uparrow }$ atoms with the spin projection 1/2, and ${N}_{\downarrow }$ atom with the spin projection $-1/2$). We will consider either systems with ${N}_{\uparrow }={N}_{\downarrow },$ that is, in the simplest spin singlet S = 0 state, or ones with ${N}_{\uparrow }\ne {N}_{\downarrow }$ such that the total spin $S=| {N}_{\uparrow }-{N}_{\downarrow }| ,$ that is, in a spin-stretched (high-spin) configuration.

In almost all many-body theories it is customary to use the so-called algebraic approximation [77], i.e. to use a finite set of one-particle functions. These one-particle functions are usually expanded in terms of some known basis functions. In our case the most natural choice of the basis functions is obviously given by the eigenfunctions of the one-dimensional harmonic oscillator:

$$\begin{eqnarray}&&{\varphi }_{n}(x)=\displaystyle \frac{1}{\sqrt{{2}^{n}n!}}\displaystyle \frac{1}{{\pi }^{1/4}}{{\rm{e}}}^{-\frac{{x}^{2}}{2}}{H}_{n}\left(x\right),\end{eqnarray} \tag{ 2 }$$

where H_n(x) are the Hermite polynomials. The integrals of the one-particle part of the Hamiltonian are diagonal, and given by the eigenenergies of the harmonic oscillator, while the two-particle integrals may be calculated numerically using the exact Gauss–Hermite quadratures.

Obviously, eigenfunctions of the harmonic oscillator provide the basis set that is the simplest for numerical applications. There are many other choices, that require more numerical effort, but should assure better numerical precision and convergence: an obvious example in the case of non-harmonic potential is to use the corresponding orbitals that are exact eigenfunctions of the one-particle Hamiltonian. Explorations of these choices go, however, beyond the scope of this paper.

To investigate in detail the properties of the considered system it is useful to solve analytically the two-body case first. The exact solution of the two-body problem in a one-dimensional harmonic trap was found by Busch et al[78] and by Franke-Arnold et al [79]. For convenience and to introduce the notation we decided to include a detailed description of the solution of the aforementioned Hamiltonian in appendix A. The ground state wave function is

$$\begin{eqnarray}&&{\rm{\Psi }}(X,x)={C}_{n}{C}_{\epsilon }\;{{\rm{e}}}^{-{X}^{2}/2-{x}^{2}/2}{H}_{n}(X)U(-\epsilon /2,1/2,{x}^{2}),\end{eqnarray} \tag{ 3 }$$

where $X=({x}_{1}+{x}_{2})/\sqrt{2}$ and $x=({x}_{1}-{x}_{2})/\sqrt{2}$ are the center of mass and relative coordinates, and C_n and C are the normalization constants for the X- and x-dependent portions. The corresponding total energy of the system (including explicitly the zero-point energies of the two particles) is then $E=\epsilon +n+1,$ while $U(a,b,x)$ stands for the Tricomi function [80]. The energy of the relative motion is determined from

$$\begin{eqnarray}&&\displaystyle \frac{1}{g}=\displaystyle \frac{{\rm{\Gamma }}\left(\frac{2-\epsilon }{2}\right)}{\sqrt{2}\epsilon {\rm{\Gamma }}\left(\frac{1-\epsilon }{2}\right)},\end{eqnarray} \tag{ 4 }$$

where ${\rm{\Gamma }}(.)$ is the Gamma function. See appendix A for more details.

To address the problem of the energy convergence with increasing number n_b of one-particle functions in the basis set, we have pursued two complementary strategies. The first follows the one developed by Hill [81] for the case of the helium atom, where the exact wavefunction, expanded in terms of an infinite sum over basis functions, is truncated at a fixed number n_b of basis set functions. This leads to an approximated energy ${\epsilon }_{{n}_{b}}$ converging to the exact value as (see appendix B for details)

$$\begin{eqnarray}&&{\epsilon }_{{n}_{b}}-\epsilon =\displaystyle \frac{{2}^{5/2}{C}_{\epsilon }^{2}}{{\rm{\Gamma }}{(-\frac{1}{2}\epsilon )}^{2}}\left(\displaystyle \frac{1}{\sqrt{{n}_{b}}}+\displaystyle \frac{g}{\pi {n}_{b}}\right.\left.+\displaystyle \frac{15+8\epsilon }{48\;{n}_{b}^{3/2}}\right)+{\mathcal{O}}\left({n}_{b}^{-2}\right).\end{eqnarray} \tag{ 5 }$$

Note that the formula found by Hill, although derived explicitly for the two-electron case, has been successfully applied across a wide range of atomic and molecular calculations (see e.g. [82, 83]). Our aim here is similar, in that the above functional form of the two-body extrapolation formula will prove crucial for obtaining accurate results for systems with $N\gt 2,$ as discussed below.

Since the approach outlined above does not allow for the optimization of the coefficients in a smaller Hilbert space, we also followed a second strategy, based on an actual variational minimization of the energy in the space spanned by the states contained in the basis set. The result reads (see appendix C for details)

$$\begin{eqnarray}&&{\tilde{\epsilon }}_{{n}_{b}}-\epsilon ={\mathcal{A}}\left\{\displaystyle \frac{1}{\sqrt{{n}_{b}}}-\displaystyle \frac{{\mathcal{A}}x^{\prime\prime} }{2x^{\prime} }\displaystyle \frac{1}{{n}_{b}}+\right.\left.\left[\displaystyle \frac{5+4\epsilon }{24}+\displaystyle \frac{{{\mathcal{A}}}^{2}}{2}{\left(\displaystyle \frac{x^{\prime\prime} }{x^{\prime} }\right)}^{2}-\displaystyle \frac{x\prime \prime \prime }{3x^{\prime} }\right]\displaystyle \frac{1}{{n}_{b}^{3/2}}\right\}+{\mathcal{O}}({n}_{b}^{-2}),\end{eqnarray} \tag{ 6 }$$

where $x={g}^{-1}(\epsilon ),$ as given by equation (4), $x^{\prime} \equiv {\partial }_{\epsilon }[{g}^{-1}(\epsilon )],$ and ${\mathcal{A}}=-1/(\sqrt{2}\pi x^{\prime} ).$

Note that both approaches yield the same functional form of the convergence formula, differing in coefficients that are specific to the two-body case, and leading terms in equations (5) and (6) are numerically confirmed to be the same.

As already stated above, almost all many-body theories are based on the algebraic approximation [77], i.e., on the parametrization of the wave function by expansion in a finite set of basis functions. Since we deal with states with total spin S = 0 or high-spin (spin stretched) states with the total spin $S=| {N}_{\uparrow }-{N}_{\downarrow }| ,$ the first approximation to the exact wave function will be the Slater determinant Φ built of one-particle functions obtained by solving the Hartree–Fock equations of the mean-field theory (see [84]). We will denote the N functions occupied in the reference Slater determinant by ${\{{\phi }_{\alpha }\}}_{\alpha =1}^{N},$ while the remaining one-particle functions that are not occupied in the reference determinant will be denoted by ${\{{\phi }_{\rho }\}}_{\rho =N+1}^{2{n}_{b}},$ where n_b is the number of the harmonic oscillator functions used to expand the one-particle solutions of the Hartree–Fock equations. Note that the total number of one-particle functions is 2n_b instead of n_b since ${\phi }_{\alpha }$ and ${\phi }_{\rho }$ include the dependence on the spin coordinate through the spin functions. Obviously, the relation ${n}_{b}\gg N$ must hold. There is a theorem stating that if a set of one-particle functions spans the one-particle Hilbert space ${L}^{2}({{\mathbb{R}}}^{3}),$ then the set of all N-particle determinantal wave functions constructed from this one-particle set will span the antisymmetric part of the N particle Hilbert space [85]. In other words, the N particle wave function can rigorously be expanded in terms of the Slater determinants built from the complete set of one-particle functions. The same theorem holds in the algebraic approximation, and it is the basis of exact diagonalization, otherwise referred to as the full configuration interaction (FCI) method. The difference between the mean-field and FCI energies is called the correlation energy.

In the FCI method the wave function of the many-fermion system is represented as

$$\begin{eqnarray}&&{\rm{\Psi }}=(1+\hat{C}){\rm{\Phi }},\end{eqnarray} \tag{ 7 }$$

where the CI operator $\hat{C}$ may be written as a linear combination of l-particle excitation operators,

$$\begin{eqnarray}&&\hat{C}=\displaystyle \sum _{l=1}^{N}\hat{{C}_{l}},\ \ \ \ \ \hat{{C}_{l}}=\displaystyle \frac{1}{{(l!)}^{2}}{C}_{{\rho }_{1}...{\rho }_{l}}^{{\alpha }_{1}...{\alpha }_{l}}{e}_{{\alpha }_{1}...{\alpha }_{l}}^{{\rho }_{1}...{\rho }_{l}},\end{eqnarray} \tag{ 8 }$$

with coefficients ${C}_{{\rho }_{1}...{\rho }_{l}}^{{\alpha }_{1}...{\alpha }_{l}}.$ Summation over the repeated lower and upper indices is assumed. Note that the indices α and ρ refer to the one-particle functions that are occupied or empty in the reference Slater determinant Φ, so equation (7) can be viewed as an operator representation of the expansion of the exact wave function in terms of the all possible N-particle determinantal wave functions constructed from the ${\{{\phi }_{\alpha }\}}_{\alpha =1}^{N}\oplus {\{{\phi }_{\rho }\}}_{\rho =N+1}^{2{n}_{b}}$ set of the one-particle functions. The l-particle excitation operator ${e}_{{\alpha }_{1}...{\alpha }_{l}}^{{\rho }_{1}...{\rho }_{l}}$ is given by the product of one-particle excitation operators ${e}_{{\alpha }_{1}}^{{\rho }_{1}}\cdots {e}_{{\alpha }_{l}}^{{\rho }_{l}}.$ The one-particle excitation operators are in turn defined in terms of the conventional fermionic creation and annihilation operators ${e}_{\alpha }^{\rho }={a}_{\rho }^{\dagger }{a}_{\alpha },$ where the creation and annihilation operators are defined with the respect to the physical vacuum $| 0\rangle ,$ but applied to the Fermi vacuum Φ [86]. Note that the one-particle functions ${\phi }_{\alpha }$ and ${\phi }_{\rho }$ are always orthogonal, so that the excitation operators commute. Note also that our approach is limited to the N particle Hilbert space (layer of the Fock space with the fixed number of spin-1/2 fermions). This means that the algebra of the excitation operators is not only commutative, but also nilpotent. The latter property easily follows from the fact that we can replace at most N one-particle functions ${\phi }_{\alpha }$ in the reference Slater determinant, and any further action will give zero as the result. An important consequence of this property is that any transcendental function of the excitation operators that has a well defined Taylor expansion will reduce to a finite polynomial. The CI coefficients ${C}_{{\rho }_{1}...{\rho }_{l}}^{{\alpha }_{1}...{\alpha }_{l}}$ and the energy are obtained by diagonalizing the Hamiltonian matrix constructed from the matrix elements ${{\mathbb{H}}}_{{lm}}=\langle {e}_{{\alpha }_{1}...{\alpha }_{l}}^{{\rho }_{1}...{\rho }_{l}}{\rm{\Phi }}| \hat{H}| {e}_{{\alpha }_{1}...{\alpha }_{m}}^{{\rho }_{1}...{\rho }_{m}}{\rm{\Phi }}\rangle .$ This explains why the FCI method is also referred to as the exact diagonalization method.

The computational cost of the FCI method is prohibitive, as it scales as ${N}^{2}{(2{n}_{b})}^{N+2}$ for $N\ll {n}_{b}$ [87], which limits this method to small systems (small N) or small number of basis functions n_b. One possible way to cure this problem is to limit the number of excitations included in equation (8), e.g., to single and double excitations. Such a truncation may lead to considerable savings of computer time, but unfortunately has a serious drawback. Any truncated CI method is no longer size-extensive, which means the energy of two non-interacting systems is not the sum of the energies of these two systems.

To overcome the size-extensivity problem of the limited CI expansions, the CC method was introduced, first in the nuclear physics [58] and slightly later in quantum chemistry [59, 60] and in the electron gas theory [88]. In this method the wave function is given by the following exponential ansatz

$$\begin{eqnarray}&&{\rm{\Psi }}={{\rm{e}}}^{\hat{T}}{\rm{\Phi }},\end{eqnarray} \tag{ 9 }$$

where the cluster operator $\hat{T}$ is given by

$$\begin{eqnarray}&&\hat{T}=\displaystyle \sum _{l=1}^{N}\hat{{T}_{l}},\ \ \ \ \ \hat{{T}_{l}}=\displaystyle \frac{1}{{(l!)}^{2}}{t}_{{\rho }_{1}...{\rho }_{l}}^{{\alpha }_{1}...{\alpha }_{l}}{e}_{{\alpha }_{1}...{\alpha }_{l}}^{{\rho }_{1}...{\rho }_{l}}.\end{eqnarray} \tag{ 10 }$$

Note that by comparison with equation (7) the following relation between the FCI and cluster operator must hold, $\hat{T}=\mathrm{ln}(1+\hat{C}),$ but since the algebra of the excitation operators is nilpotent, the logarithm, as well as the exponential function in equation (9), reduce to finite polynomials. The CC method is non-variational, in the sense that it does not involve any optimization over a set of variational parameters. The energy is given by the expression

$$\begin{eqnarray}&&E=\langle {\rm{\Phi }}| {{\rm{e}}}^{-\hat{T}}\hat{H}{{\rm{e}}}^{\hat{T}}{\rm{\Phi }}\rangle ,\end{eqnarray} \tag{ 11 }$$

while the cluster amplitudes ${t}_{{\rho }_{1}...{\rho }_{l}}^{{\alpha }_{1}...{\alpha }_{l}}$ are obtained by solving the following system of nonlinear algebraic equations

$$\begin{eqnarray}&&0=\langle {e}_{{\alpha }_{1}...{\alpha }_{l}}^{{\rho }_{1}...{\rho }_{l}}{\rm{\Phi }}| {{\rm{e}}}^{-\hat{T}}\hat{H}{{\rm{e}}}^{\hat{T}}{\rm{\Phi }}\rangle ,\ \ \ \ \ l=1,\ldots ,N.\end{eqnarray} \tag{ 12 }$$

Note that by virtue of the Baker–Campbell–Hausdorff formula the exponential factor ${{\rm{e}}}^{-\hat{T}}\hat{H}{{\rm{e}}}^{\hat{T}}$ reduces to multiple commutators of $\hat{H}$ and $\hat{T},$ and one can prove that this commutator expansion is finite and contains at most four-fold commutators if only two-particle interactions are present in the Hamiltonian [58].

Obviously, the CC method including all excitations is fully equivalent to the FCI method, and its computational cost is as high. However, the truncated CC methods are much less time consuming, and unlike the truncated CI methods, remain size-extensive. Due to the exponential form of the ansatz, equation (9), the CC method truncated to single and double excitations effectively includes triply, quadruply, and higher excited determinants through the products, e.g. ${\hat{T}}_{1}{\hat{T}}_{2}$ or ${\hat{T}}_{2}^{2}.$ In the present work we will use the series of approximations with the cluster operator limited to single and double (CCSD) [89], single, double, and triple (CCSDT) [90], and single, double, triple, and quadruple excitations (CCSDTQ) [91–93]. The computational cost of these methods scales as ${n}_{b}^{6},$ ${n}_{b}^{8},$ and ${n}_{b}^{10},$ respectively. For a system consisting of a dozen of atoms methods including triple excitations cannot reasonably be used, while the experience gained for many-electron systems suggests that the triple excitations have an important contribution to the correlation energy. To reduce the computational cost of the CC method including triple excitations, the CCSD(T) method with the n_b⁷ scaling was introduced [94]. In this method the CCSD equations are solved iteratively, and the contribution of the triple excitations to the correlation energy is evaluated from the expression based on the many-body perturbation theory (MBPT). It turned out that the CCSD(T) method is very accurate for many properties of atoms and molecules and, as for now, it is considered as the golden standard of quantum chemistry. As we show in the next section, this method can also be applied with success to the system of ultracold atoms in a one-dimensional harmonic trap interacting via a short-range contact type potential.

In order to introduce the perturbation expansion of the correlation energy, often referred to as MBPT [63], or Møller–Plesset perturbation theory [95–97], it is useful to rewrite the equations for the cluster operator $\hat{T}$ in the following operator form [98, 99]

$$\begin{eqnarray}&&\hat{T}=\hat{{\mathcal{R}}}\left({{\rm{e}}}^{-\hat{T}}\hat{W}{{\rm{e}}}^{\hat{T}}\right)=\hat{{\mathcal{R}}}\left(\hat{W}+[\hat{W},\hat{T}]+\displaystyle \frac{1}{2}{[\hat{W},\hat{T}]}_{2}+\displaystyle \frac{1}{6}{[\hat{W},\hat{T}]}_{3}+\displaystyle \frac{1}{24}{[\hat{W},\hat{T}]}_{4}\right),\end{eqnarray} \tag{ 13 }$$

where $\hat{W}$ is the correlation (fluctuation) potential in the Møller–Plesset partitioning of the Hamiltonian

$$\begin{eqnarray}&&\hat{H}=\hat{F}+\hat{W},\end{eqnarray} \tag{ 14 }$$

$\hat{F}$ is the Fock operator, and ${[\hat{W},\hat{T}]}_{n}$ denotes the n-fold nested commutator: ${[\hat{W},\hat{T}]}_{0}=\hat{W},$ ${[\hat{W},\hat{\hat{T}}]}_{n+1}=[{[\hat{W},\hat{T}]}_{n},\hat{T}].$ The nested commutator expansion in equation (13) terminates after the quadruple commutators since in our case the operator W contains only two-particle interactions. The resolvent superoperator $\hat{{\mathcal{R}}}$ is defined for arbitrary operator $\hat{X}$ as [98, 99]

$$\begin{eqnarray}&&\hat{{\mathcal{R}}}=\displaystyle \sum _{n=1}^{N}{\hat{{\mathcal{R}}}}_{n},\ \ \ \ \ {\hat{{\mathcal{R}}}}_{n}(\hat{X})={\left(\displaystyle \frac{1}{n!}\right)}^{2}\langle {e}_{{\rho }_{1}...{\rho }_{n}}^{{\alpha }_{1}...{\alpha }_{n}}{\rm{\Phi }}| \hat{X}{\rm{\Phi }}\rangle \displaystyle \frac{{e}_{{\alpha }_{1}...{\alpha }_{n}}^{{\rho }_{1}...{\rho }_{n}}}{{\epsilon }_{{\rho }_{1}...{\rho }_{n}}^{{\alpha }_{1}...{\alpha }_{n}}},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{\epsilon }_{{\rho }_{1}...{\rho }_{n}}^{{\alpha }_{1}...{\alpha }_{n}}={\epsilon }_{{\alpha }_{1}}+\cdots +{\epsilon }_{{\alpha }_{n}}-{\epsilon }_{{\rho }_{1}}-\cdots -{\epsilon }_{{\rho }_{n}},\end{eqnarray} \tag{ 16 }$$

and ${\epsilon }_{\kappa }$ denotes the one-particle energy associated with the one-particle function labeled by the index κ. We assume that the energy of the highest occupied one-particle function is smaller than the energy of the lowest unoccupied in the reference determinant Φ, so the superoperators ${\hat{{\mathcal{R}}}}_{n}$ are always well defined. Note that for a given $\hat{X}$ the operator $\hat{Y}={\hat{{\mathcal{R}}}}_{n}(\hat{X})$ can be viewed as a formal solution of the equation

$$\begin{eqnarray}&&\langle {e}_{{\alpha }_{1}...{\alpha }_{n}}^{{\rho }_{1}...{\rho }_{n}}{\rm{\Phi }}| ([\hat{F},\hat{Y}]+\hat{X}){\rm{\Phi }}\rangle =0.\end{eqnarray} \tag{ 17 }$$

Obviously, this solution is unique if we assume that $\hat{Y}$ belongs to the linear span of the excitation operators ${e}_{{\rho }_{1}...{\rho }_{n}}^{{\alpha }_{1}...{\alpha }_{n}}.$

The many-body perturbation expansion of the energy is obtained by introducing the following parametrization of the correlation operator $\hat{W},$

$$\begin{eqnarray}&&\hat{W}\to \lambda \hat{W},\end{eqnarray} \tag{ 18 }$$

where the complex parameter λ is introduced to derive the perturbation expansions of equations (11) and (13). Obviously, the physical value of λ is equal to one. By substituting the parametrized equation (18) into (14) the cluster operator became dependent on λ, and can be expanded as a power series in λ

$$\begin{eqnarray}&&\hat{T}=\displaystyle \sum _{k=1}^{\infty }{\lambda }^{k}{\hat{T}}^{(k)}.\end{eqnarray} \tag{ 19 }$$

Substituting the above expansion in the expression for the energy (11) and collecting all terms at ${\lambda }^{n}$ gives the expression for the nth-order correction to the energy in the MBPT. The cluster operators necessary to evaluate this expression are obtained by substituting equation (19) and again collecting all terms at ${\lambda }^{n}.$

It follows immediately from this short sketch of the derivation that the MBPT and CC theories are strongly connected. The analysis reported in [99] shows that the CCSD, CCSDT, and CCSDTQ are valid through the third, fourth, and fifth order of the MBPT. The golden standard of quantum chemistry, the CCSD(T) method, is also valid through the fourth order of MBPT, although it is computationally less demanding than the full CCSDT theory.

One of the observables that can be obtained with the presented method is the density profile of the trapped atoms. At a point x₀ it may be obtained as the expectation value of the operator ${\displaystyle \sum }_{i=1}^{N}\delta ({x}_{i}-{x}_{0})$ within the FCI calculations, and from the Hellmann–Feynman theorem as the first derivative of the energy in the presence of the perturbation given by the same operator in the CC calculations. In practice, we compute this derivative by using the finite-field approach with a perturbation of the form ${\phi }_{i}({x}_{0}){\phi }_{j}({x}_{0})$ added to the ijth element of the one-particle Hamiltonian matrix.

The FCI and CC calculations were performed with the customized versions of the HECTOR [100] and ACESII codes [101], respectively. The CC and truncated CI calculations use the reference state Φ built of orbitals obtained with the restricted Hartree–Fock method for the ${N}_{\uparrow }={N}_{\downarrow }$ case, whereas orbitals obtained with the unrestricted Hartree–Fock method are employed for the ${N}_{\uparrow }\ne {N}_{\downarrow }$ case.

Examples of a slow convergence of the results with respect to the basis set size are well-known in the literature. The conventional exact diagonalization calculations for solving the electronic Schrödinger equation converge as ${L}^{-3},$ where L is the highest angular momentum present in the one-electron basis set [81, 102]. Even slower convergence (${L}^{-1}$) was observed for some relativistic corrections arising from the perturbative approach based on the Breit–Pauli Hamiltonian [103].

In order to check that the convergence formulas of equations (5) and (6) for the two-body problem hold, we compare their predictions with the exact results in figure 1. Two sets of calculated energies relative to the exact value are presented as a function of the one-particle basis set size: non-variational energies obtained as the expectation value with the truncated exact wavefunction (equation (B.5) in appendix B) and variational energies equivalent to exact diagonalization results. The former one should be described by the asymptotic expansion of equation (5), whereas the latter ones should coincide with the asymptotic expansion of equation (6). The asymptotic formulas with increasing number of terms are presented. The first terms of the asymptotic expansions give a reasonable estimate of the convergence rate and upon inclusion of the second and third terms the differences between the analytical formulas and calculated values are greatly reduced. Therefore, the formulas (5) and (6) are valid and can serve as a guide for further investigations.

Zoom In Zoom Out Reset image size

To evaluate the adequateness of the functional formula given by equations (5) and (6) for extrapolation of the variational energies to the complete basis set limit, figure 1 presents also the fit of the variational results to the formula ${E}_{\infty }+{{An}}_{b}^{-1/2}+{{Bn}}_{b}^{-1}+{{Cn}}_{b}^{-3/2}.$ This extrapolation formula behaves extremely well (see also the inset of figure 1). For instance, the extrapolated energies in the complete basis set limit ${E}_{\infty }$ for g = 1 and g = 5 are 1.306 744 and 1.726 780, while the exact values equal to 1.306 745 and 1.726 771, respectively. Note that in the biggest basis set available, ${n}_{b}=200,$ the calculated energies are 1.310 545 and 1.745 325, respectively, so that the accuracy gain of three order of magnitude is impressive and very important.

Given the very slow convergence rate of the two-body energy with the number of the basis functions n_b, it is very important to consider this convergence in the many-body case. One can expect that the convergence pattern obtained for the two-body case is almost certainly valid for the many-body case. This fact was observed in the electronic structure calculations on atoms and molecules, and is related to the fact that the correlated pair functions of two electrons have the same behavior around the electrons coalescence points as in the discussed two-body example. Clearly, this analytic behavior determines the convergence rate of the results towards the exact value. A very similar situation is found for the partial wave expansion of the exact solution of the Schrödinger equation. The well-known ${L}^{-3}$ convergence pattern has rigorously been derived only for the helium atom [81, 102], but was successfully applied also for many-electron atoms/molecules [104, 105] and is widely accepted to be universal (see the discussion given by King [106]).

Following the discussion above and guided by the expressions (5) and (6), we decided to adopt the following three-term extrapolation formula

$$\begin{eqnarray}&&{E}_{{n}_{b}}-{E}_{\infty }=\displaystyle \frac{A}{\sqrt{{n}_{b}}}+\displaystyle \frac{B}{{n}_{b}},\end{eqnarray} \tag{ 20 }$$

where ${E}_{\infty }$ is the extrapolated energy to the infinite number of the basis functions, E_nb is the energy computed with n_b basis functions, and A and B are the fit parameters. To check the correctness of the above expression we performed FCI calculations for the ground state of a balanced two-component Fermi gas. We have found that the convergence pattern with the number of the basis functions exactly follows our extrapolation formula. The accuracy of our formula is illustrated in figure 2. There we show the absolute error with respect to the complete basis set limit in the ground state energy for the 2 + 2 and 3 + 3 systems (${\text{}}{N}_{\uparrow }\;+\;{\text{}}{N}_{\downarrow }$ means ${N}_{\uparrow }$ atoms with the spin projection 1/2, and ${N}_{\downarrow }$ atom with the spin projection $-1/2$) as a function of the number of one-particle basis functions for several interaction strengths obtained with the FCI method, as well as the error plotted as a function of n_b according to equation (20). The agreement between the computed points and the analytical fits to the extrapolation function is excellent, supporting the correctness of our extrapolation scheme.

Zoom In Zoom Out Reset image size

In order to get saturated results, careful studies of the convergence are necessary. As stated above, the system of many identical fermions with spin-1/2 in a 1D trap interacting via the contact potential has an even worse convergence than the one observed for the electronic Schrödinger equation [81, 102, 103], and apparent convergence may be observed. This is illustrated in figure 3 where we show the energy spectra for a few atoms in a trap obtained with a small number of the basis functions (n_b = 20 as in [48]) and in the complete basis set limit. An inspection of this figure shows that the results obtained with 20 basis functions are far from the complete basis set limit obtained from equation (20), especially in the limit of intermediate and large interaction strength. The apparent convergence of the results observed by the authors of [48] is solely due to the pathological convergence pattern as ${n}_{b}^{-1/2}.$ As will be shown in the next section, a quantitative picture of the physics, and a quantitative agreement with various experimental data, can only be obtained if the extrapolation to the complete basis set limit is properly done.

Zoom In Zoom Out Reset image size

As seen in figures 2 and 3, both the absolute and relative errors due to the finite basis set size are increasing with the increasing interaction strength g. In the limit of very strong interactions, particularly approaching the unitary limit of TG gas (when the interaction strength is infinite, $g=\infty$), this error can become larger than the separation between the ground and excited states and will lead to artefacts in the description of physical phenomena both at zero and finite temperatures. The authors of [48] investigated with the exact diagonalization method the few-fermion physics in the TG limit when the spectrum is quasi-degenerate, but the careful convergence analysis reveals that the results obtained with the basis set as small as n_b = 20 give only a qualitative picture in this regime of interaction strength.

Figure 4 shows the energy spectrum for the example of the 3 + 1 system as a function of $1/g$ obtained with the FCI method in the limit of the strong interaction in the basis of n_b = 20 one-particle functions and in the complete basis set limit ${n}_{b}=\infty .$ The complete basis set limit results were obtained by extrapolation from the numerical results in basis of 30, 40, 50, and 60 functions. In the TG limit the four states become degenerate. Unfortunately, any calculation with a finite number of one-particle basis set functions will always locate a crossing of these states at a finite g, potentially leading to incorrect conclusions. For n_b = 20 the crossing appears at g = 17.1 and corresponds to the large error of $0.4\;\omega$ in the ground-state energy. Therefore, the extrapolation to the complete basis set limit is not only needed to improve accuracy, but is necessary to have a quantitatively correct picture of the physics. The numerical results in the complete basis set limit agree perfectly well with the analytical results obtained in the vicinity of the TG limit with the perturbation theory [50, 53]. The slopes of the curves fitted to the numerical data are $\{-7.14,-3.58,-1.18,0\},$ in a very good agreement with the exact contact coefficients $\{-7.08,-3.58,-1.18,0\}.$ This confirms that the presented extrapolation scheme works also very well in the regime of strong interactions allowing to approach the TG limit with the accurate finite basis set calculations. Interestingly, the contact coefficients corresponding to non-converged results, even in basis set as small as ${n}_{b}=20,$ are surprisingly close to the correct values for the TG limit.

Zoom In Zoom Out Reset image size

As we have shown above, the energy spectrum in the TG limit obtained in the calculations with a single one-particle basis set has always an artificial crossing of the states at a finite value of the interaction strength g. In figure 5 we plot the values of $1/g$ for which this crossing occurs as a function of the one-particle basis set size for a few systems. Interestingly, the interaction strength at which the states in the ground manifold become degenerate depends solely on the number of the basis functions and not on the number of atoms, neither on their state. This observation agrees with our prediction on the convergence of $1/g$ for a given energy (in this case the energy of the ferromagnetic state) with the size of the one-particle basis set to be independent in the leading order on both interaction strength and energy. The leading term of the convergence formula for the two-body problem, derived analytically in appendix C, is shown as a solid black line in figure 5. This observation suggests a second approach to get complete basis set limit results and to reach the TG limit with accurate finite basis set calculations, that is, instead of extrapolating the energy, one can extrapolate $1/g$ for fixed energies by using the universal convergence formula. Exactly the same convergence coefficient valid for all presented few-body cases and observed universality suggest even a third approach to get accurately converged results with finite basis set calculations at a given interaction strength g. The idea is to use in the numerical calculations a renormalized coupling constant g_nb explicitly dependent on the high-energy cut-off set by n_b, which reproduces exact two-body result in a smaller Hilbert space spanned by n_b basis functions. In this way, the inaccurate description of the short-distance two-body physics introduced by the high-energy cutoff may be cured order by order. A similar idea was discussed, e.g., in [107, 108].

Zoom In Zoom Out Reset image size

Having solved the problem with the convergence of the results with the number of basis functions, we now turn to the effect of the truncation of the excitation in the cluster operator, equation (10), on the correlation and total energies. We start the discussion of our results with figure 6 where we report the correlation energy for the 2 + 2 and 3 + 3 systems as a function of the interaction strength for two one-particle basis set sizes and the complete basis set limit. An inspection of figure 6 shows that the basis set dependence of the energy is relatively weak in the weakly correlated repulsive regime ($g\gt 0$), and very pronounced in the strongly correlated attractive case ($g\lt 0$). These significantly different behaviors result from the fact that in the limit of the strong repulsion even distinguishable fermions tend to occupy different one-particle states and the total wave function approaches the structure of the ferromagnetic state (the TG gas limit). On the other hand, the fermions with the opposite spin projection tend to pair in the case of the attractive interaction and become tightly bound (hard-core) bosonic dimers in the limit of the strong attraction (the Lieb–Liniger gas limit). The description of these tightly bound pairs is challenging for calculations in the one-particle functions basis sets and explicitly correlated methods could potentially overcome the problem and significantly accelerate the convergence. It is important to stress at this point that the mean-field energy converges very fast with the number of one-particle functions, so the convergence problems are related to the basis saturation of the correlation energy.

Zoom In Zoom Out Reset image size

The correlation energy presented in figure 6 diverges linearly with the interaction strength g for both negative and positive values. One should note that in the weak interaction regime the distinction between the mean-field and correlation energies is reasonable. For intermediate and especially strong interaction regimes the mean-field description fails completely and diverging correlation energy compensates unphysically diverging mean-field energy, so no physical meaning should be attributed to this divergence. However, this behavior of the mean-field and correlation energies affects the performance of post-mean-field methods in the strong interaction regime (including the standard CC method), which start from the mean-field solution and tend to recover correlation energy.

We now turn to the applicability of the MBPT and truncated configuration interaction expansions for the energy calculations. In figure 7 we report the percentage of the ground state correlation energy reproduced with the second-order and fourth-order MBPT, MBPT2 and MBPT4, and the configuration interaction method limited to single and double excitations, CISD, for the 2 + 1, 2 + 2, 3 + 2, and 3 + 3 systems as functions of the interaction strength. As expected from the electronic structure calculations on atoms and molecules the performance of the MBPT2 and MBPT4 methods is not very good, and quite erratic as a function of the interaction strength. Given the fact that for atoms and molecules, the MBPT is divergent [109–111] as perturbation theories applied in a different physical context [112, 113], it is not surprising that the MBPT4 theory is not a big improvement over the MBPT2 approach, despite a much higher theoretical complexity and computational time scaling with the size of the basis, n_b⁴ for MBPT2 versus n_b⁷ for MBPT4. Finally, we note that the variational CISD results do not offer a good advantage over the MBPT results.

Zoom In Zoom Out Reset image size

Since the MBPT and limited CI theory perform badly, one may ask if a selective infinite-order summation of some MBPT diagrams with different variants of the CC theory will improve the situation. This is indeed the case. Figure 8 shows the percentage of the ground state correlation energy reproduced with the CCSD, CCSD(T), CCSDT, and CCSDTQ methods for the 2 + 1, 2 + 2, 3 + 2, and 3 + 3 systems as functions of the interaction strength. First of all note that for the 2 + 1 system CCSDT is equivalent to the FCI theory, and for the 2 + 2 system the same statement is valid for CCSDTQ. This means that for these systems these particular methods will reproduce 100% of the energy. It may be surprising at first glance that some of the approximate variants of the CC theory overestimate the total energy of the system. This is due to the fact that the truncated CC methods are not variational, and one can obtain e.g. 101% of the energy. An inspection of figure 8 shows that all CC methods perform very well, although there is no clear convergence pattern with the excitation operators included in the CC wave function. Indeed, the CCSD method tends to slightly underestimate the energy, while methods including triple (and possibly quadrupole) excitations tend to slightly overestimate. In general, the percentage errors are of the order of a few percent. Note also that the CCSDT method does not offer big advantage over the CCSD(T) approach, despite of much higher computational requirements. This is in line with the results obtained for atoms and molecules, see for instance [103, 114]. Comparison of figures 7 and 8 shows that the cluster expansion of the wave function and the summation of the diagrams involving the single and double excitation in the MBPT theory to infinite order is crucial for an accurate description of many-atom systems in a one-dimensional harmonic trap.

Zoom In Zoom Out Reset image size

The discussion of the last two paragraphs is well summarized in table 1, where we show the CC, MBPT, and CISD results for selected values of the interaction strength. An analysis of the numerical results reported in this table confirms a very good performance of various variants of the CC theory, as opposed to the erratic behavior of the MBPT, and the poor performance of the configuration interaction method limited to single and double excitations.

Table 1.  The percentage of the ground state correlation energy reproduced at the CCSD, CCSD(T), CCSDT, and CCSDTQ levels of the CC theory and the many-body perturbation theory MBPT2, MBPT4, and the truncated configuration interaction methods CISD in the 2 + 1, 2 + 2, 3 + 2, and 3 + 3 systems for the different interaction strengths. N/D (no data) indicates that calculation was not feasible with used software.

System	CCSD	CCSD(T)	CCSDT	CCSDTQ	MBPT2	MBPT4	CISD
g = 2
$N=2+1$	99.60	100.07	100	—	117.33	100.31	99.42
$N=2+2$	99.30	100.10	99.95	100	115.26	101.49	98.70
$N=3+2$	99.35	100.07	99.97	N/D	113.48	100.18	97.40
$N=3+3$	99.23	100.06	99.94	99.999	120.07	101.65	96.17
g = 4
$N=2+1$	97.82	99.69	100	—	92.42	98.28	95.31
$N=2+2$	98.53	100.72	99.90	100	153.09	127.21	95.39
$N=3+2$	96.89	100.18	99.86	N/D	97.54	95.05	87.18
$N=3+3$	98.14	100.56	99.84	99.997	144.33	116.83	90.64
$g=-4$
$N=2+1$	98.33	99.46	100	—	76.51	100.08	98.24
$N=2+2$	99.53	101.97	102.64	100	68.49	96.08	87.60
$N=3+2$	97.96	100.10	100.35	N/D	97.54	95.05	87.18
$N=3+3$	99.79	102.35	102.99	100.003	144.33	116.83	90.64

Finally, we note that the performance discussed above of the truncated CC methods shows a weak dependence on the number of one-particle basis functions. Similarly, we do not observe any increase of the relative error with the number of atoms, although our comparison with the FCI results is restricted to the maximum of six atoms. For a larger number of atoms we can evaluate the contributions of the noniterative triple excitations in the CCSD(T) method. Figure 9(a) shows the interaction energy in the ground state of the 10 + 10, 20 + 20, 30 + 30, and 40 + 40 systems obtained at the CCSD and CCSD(T) levels of the CC theory. Figure 9(b) presents the percentage of the ground state interaction energy accounted by the difference between energy obtained with the CCSD(T) and CCSD methods, i.e. by the connected triples contribution calculated noniteratively within the MBPT. Interestingly, for a given interaction strength the contribution of the excitations higher than double to the ground-state interaction energy is becoming less important with increasing number of atoms in the system. Based on this fact we predict that the performance of the CCSD(T) method for larger systems is as good as for the investigated small systems and the lack of higher excitations in the CC wave function does not introduce any significant errors.

Zoom In Zoom Out Reset image size

The energy spectrum of the investigated systems for strong repulsive interactions becomes quasi-degenerate (see figures 3 and 4) and the correlation energy diverges for both strongly attractive and repulsive interactions (see figure 6). These two effects lead to the convergence problem of the standard single-reference CC method starting from the antiferromagnetic reference state as used in the present study, restricting the interaction strength g that can effectively be used in actual calculations to intermediate values between −6 and 6. The possible solution to overcome this problem is the use of the multireference version of the CC theory or the ferromagnetic reference state for calculations in the limit of the strong repulsive interaction.

The density profiles of the trapped clouds are important observables that can be measured experimentally and used to monitor the evolution of the interatomic interactions and resulting states of a many-body system.

In [115] we have analyzed the density profiles for the 3 + 3 system with the FCI method and for the 15 + 15 system with the CC approach. We have shown that both methods perfectly reproduce the density profiles of the non-interacting gas, and that the FCI calculations allow to approach both the limit of strong attraction (the Lieb–Liniger gas) when the fermionic atoms of opposite spin projection pair into hard-core bosonic dimers, and the limit of strong repulsion (the TG gas) when even distinguishable fermions must occupy different one-particle levels. The CC method allows to describe the evolution of the density profiles in the range of intermediate values of the interaction strength. The densities obtained with two methods agree with predictions of the local density approximation applied to the solution of the Gaudin–Yang integral equations describing a homogeneous gas [38] providing a further confirmation of the validity of the local density approximation for investigating a trapped 1D gas. Here, we present other examples for a few fermion systems obtained with the FCI method and for the many fermion systems calculated with the CC approach.

Figure 10 shows the density profiles for the 2 + 2 and 3 + 1 systems obtained with the FCI method for three interactions strengths: strongly repulsive (g = 50), moderately repulsive (g = 5) and strongly attractive ($g=-10$), together with the analytic results. In the balanced case, it is straightforward to find analytical expressions for the limiting cases of strong attraction (${n}_{{\rm{LL}}}(x)=2{\displaystyle \sum }_{i=0}^{N/2-1}| {\tilde{\varphi }}_{i}(x){| }^{2}$ for the Lieb–Liniger gas of hard-core bosonic dimers at $g=-\infty ,$ with ${\tilde{\varphi }}_{i}(x)$ the ith eigenfunction of the harmonic oscillator for a dimer of mass $2m$), strong repulsion (${n}_{{\rm{TG}}}(x)={\displaystyle \sum }_{i=0}^{N-1}| {\varphi }_{i}(x){| }^{2}$ for the TG gas of 'fermionized' fermions at $g=\infty$), and no interaction (${n}_{0}(x)=2{\displaystyle \sum }_{i=0}^{N/2-1}| {\varphi }_{i}(x){| }^{2}$ for g = 0). In the case of strongly repulsive interaction, the FCI results are indistinguishable from the analytic result confirming perfect performance of the FCI method in this limit of the interaction strength. The regime of the strong attraction is much harder to describe due to the presence of strong two-body correlations when the hard-core bosonic dimers are formed. The presented results for $g=-10$ approach the analytic results but obtaining fully converged results for much more attractive interaction strengths, even with the extrapolation used within the paper, is very challenging. A possible solution to overcome the very slow convergence in the limit of strong attractive interaction can be to use the explicitly correlated method or to include the formation of the hard-core bosonic dimers directly in the structure of the wave function.

Zoom In Zoom Out Reset image size

Figure 11(a) shows the density profiles for the 5 + 5, 10 + 10, and 20 + 20 systems obtained with the CCSD(T) method for the interaction strength g = 1. With an increasing number of atoms the evolution of the size of the cloud can be observed. As we have shown in [115], the overall shape of the profiles can be described by the typical Thomas–Fermi profile of an inverted parabola. The relative size of the density oscillations are smaller and smaller with the increasing number of atoms, and the density profile approaches the exact shape given by the local density approximation in the thermodynamic limit. Figure 11(b) shows the density profiles for the 5 + 5 system obtained with the CCSD(T) method for the three interactions strengths $g=-4,$ g = 2, and g = 6. The present CC calculations are restricted to intermediate values of the interactions strength, which however extend far beyond the mean field regime [115]. Moreover, the accessible range is sufficient to observe significant modifications of the shape and the size of the cloud. The evolution of the cloud's shape towards the analytic results of the limiting cases of strong attraction ($g=-\infty$) and strong repulsion ($g=\infty$) is clearly visible in figure 11(b). Finally, the CC method allows one to address much larger systems than the FCI approach.

Zoom In Zoom Out Reset image size

All the results reported thus far strongly suggest that the approximate variants of the CC method combined with the efficient extrapolation schemes based on the rigorous analytical formulas derived for the two-body case perform very well for a few-body and many-body systems of identical spin-1/2 atoms interacting via the short-range contact potential in a 1D harmonic trap. However, the most stringent test of the accuracy of any theoretical model is the comparison with precision experiments.

We report such a comparison in figure 12. We decided to compare the best of our results, i.e. the FCI results, but any variant of the CC method that would include triple excitation would lead to the same results within $1\%,$ indistinguishable from the FCI values on the scale of the plot. Panel (a) of this figure shows the comparison between the measured [32] and computed energies (with respect to the ground-state energy of the noninteracting system E_F) of a system of $N={N}_{\uparrow }+1$ atoms consisting of a single impurity interacting with a Fermi sea of ${N}_{\uparrow }$ identical fermions for various repulsive interactions. An inspection of this figure shows that the agreement between theory and experiment is indeed excellent. Note parenthetically that we include on this figure both the experimental error bars and the computational uncertainties ${{\rm{\Delta }}}_{E}$ conservatively estimated from the extrapolated energies and the energies computed with the largest n_b = 200 basis functions, ${{\rm{\Delta }}}_{E}={E}_{{n}_{b}=200}-{E}_{\infty }.$ On the scale of this graph the experimental and theoretical error bars are indistinguishable, so in the panel (b) of this figure we report the differences between the theoretical and experimental energies with the respective error bars. The agreement viewed in this way is also very good, and the energy difference always lies within the combined theoretical and experimental uncertainties.

Zoom In Zoom Out Reset image size

Less satisfactory is the agreement between theory and experiment for the separation energies ${\rm{\Delta }}E(N)=\mu (N)-{\mu }^{*}(N),$ where $\mu (N)=E(N)-E(N-1)$ is the chemical potential, E(N) is the extrapolated energy for the N-body system in the weakly attractive regime, and ${\mu }^{*}(N)$ is the chemical potential of the noniteracting system. This is illustrated in the panel (c) of figure 12, where we compare our theoretical results with the experimental data of [33]. An inspection of this figure shows that a relatively good agreement (but not within the experimental error bars) is observed for odd values of N, and important disagreement by a factor of roughly 3/2 for even N. Note that the experimental technique used in [33] is different than in [32]. The latter one is based on an accurate non-destructive measurement of the RF spectrum of the systems, whereas in the former one the system is probed by deforming the trapping potential and by observing the tunneling of particles out of the trap. We believe that the main source of disagreement in the second case comes from the perturbed character of the 1D harmonic shape of the trap during measurements and an approximate determination of the interaction strength in the experiment [33]. Indeed, the results of the calculations with a smaller (in the absolute value) interaction strength g are in a much better agreement for both even and odd values of N, as it has also recently been shown by the authors of [54].