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

We employ \textit{ab initio} methods of quantum chemistry to investigate spin-1/2 fermions interacting via a two-body contact potential in a one-dimensional harmonic trap. The convergence of the total energy with the size of the one-particle basis set is analytically investigated for the two-body problem and the same form of the convergence formula is numerically confirmed to be valid for the many-body case. Benchmark calculations for two to six fermions with the full configuration interaction method equivalent to the exact diagonalization approach, and the coupled cluster method including single, double, triple, and quadruple excitations are presented. The convergence of the correlation energy with the level of excitations included in the coupled cluster model is analyzed. The range of the interaction strength for which single-reference coupled cluster methods work is examined. Next, the coupled cluster method restricted to single, double, and noniterative triple excitations, CCSD(T), is employed to study a two-component Fermi gas composed of 6 to 80 atoms in a one-dimensional harmonic trap. The density profiles of trapped atomic clouds are also reported. Finally, a comparison with experimental results for few-fermion systems is presented. Upcoming possible applications and extensions of the presented approach are discussed.


I. INTRODUCTION
Ultracold gases are highly controllable systems ideal for investigating different phenomena of quantum many-body physics [1][2][3][4][5][6]. They can be prepared in a well-defined quantum state, carefully manipulated, and accurately measured, and thus can serve as a perfect tool for quantum simulations [7]. On one hand, they can be used to realize highly nontrivial states of matter. On the other hand, the problems of condensed-matter or other areas of physics can be mapped on and solved by such quantum simulators. Reducing the dimensionality of the trapped gases to one dimension brings a plethora of new interesting possibilities [8]. Experimental studies of the Tonks-Girardeau gas [9][10][11] and of the super Tonks-Girardeau gas [12] are the first and most eminent examples.
Levi Tonks (1897Tonks ( -1971 was the first to consider in 1936 [13] the problem of equation of state simultaneously for one, two, and three-dimensional gases of hard elastic spheres -this has led him to the concept of a (classical) gas of impenetrable particles in 1D. Marvin Girardeau (1930Girardeau ( -2015 was a real pioneer of the studies of the quantum impenetrable gas of bosons, known since then as Tonks-Girardeau gas. In particular, Girardeau investigated in 1960 the relation between systems of impenetrable bosons and fermions in one dimension [14]. This seminal work attained an extreme importance in the cold atom era, and has led to numerous nearly equally seminal generalizations and developments: studies of 1D Coulomb gas [15], studies of trapped bosonic gases in 1D [16], studies of quantum dynamics and quantum solitons [17], general theory of Fermi-Bose [18,19] and anyon-fermion mapping [20], investigations of super-Tonks-Girardeau gas [21,22], of 1D dipolar gases [23], and much more. From the point of view of the present work the most important were more recent generalizations and applications of Girardeau's approach to soluble models of strongly interacting 1D ultracold gas mixtures [24,25], and especially to spinor Fermi gases [26][27][28][29]. Here we investigate precisely the problem of 1D fermionic gas of atoms of spin 1/2, trapped in a harmonic potential and in the presence of strong interactions.
Recently, experiments achieving a deterministic preparation of tunable several-fermion systems in a one-dimensional trap became possible [30]. This opened the new area of ultracold research on systems with a complete control over all degrees of freedom: the particle number, the internal and motional states of the particles, and the strength of the interparticle interactions. The fermionization of two distinguishable fermions [31], the formation of a Fermi sea [32], pairing in few-fermion systems [33], two fermions in a double well potential [34], and antiferromagnetic Heisenberg spin chain of few cold atoms [35] have been investigated experimentally and are a promise of upcoming new and equally fascinating research.
The aforementioned experiments allowed to observe for the first time the transition between the few-fermion limit and the many-fermion limit of trapped atoms at ultralow temperatures. The emergence of the many-body properties of the physical systems is crucial across all areas of research. The Fermi gases in one dimension in the many-body regime have been studied intensively over the years [36][37][38][39][40][41][42][43][44][45]. Recently, several papers used analytic methods and exact diagonalization to investigate in detail the few-body regime [46][47][48][49][50][51][52][53][54][55][56][57]. Nevertheless, new numerical approaches providing additional insight on the experimental findings and having the predictive power of proposing new experiments are always welcome. Especially methods that can cover both regimes of few and many fermions are of great interest.
The goal of the present work is to provide new numerical tools to investigate the Fermi systems of few to many cold atoms. More specifically, we consider spin-1/2 fermions confined in a one-dimensional harmonic trap and interacting via a two-body contact potential. For this aim, we employ a quantum chemistry approach -the coupled cluster (CC) method [58][59][60][61][62][63][64][65][66][67][68]. This method has successfully been applied to study various properties of atoms, molecules, and condensed phases -see, for instance, Refs. [69,70] for applications to high-precision spectroscopy of ultracold molecules, Ref. [71] for simulations of liquid water properties, or Ref. [72] to determine the structure and characteristics of molecular crystals.
In the condensed-matter physics, the coupled cluster method has up to now been successfully applied to ultracold gases of bosonic atoms in traps [73,74] and allowed to describe correlations beyond the mean-field regime in a Bose-Einstein condensate of thousands of atoms. However, the ground state of bosonic systems and bosonic many-body wave functions have much simpler form than the fermionic counterparts, and therefore the advantage of the CC method was not fully pronounced in that case. The CC approach has also found numerous applications in studies of spin-1/2 lattice models, both in one-dimensional chains and in two-dimensional square lattices (see, e.g., Refs. [75,76]). In the following, we will show that the CC method proves to be ideally suited to study the problem of many fermionic atoms in one-dimensional traps.
The plan of our paper is as follows. Section II describes the theoretical framework, including the definition of the many-body Hamiltonian in subsection II A, and the discussion of the exact solution for the two-body case, and convergence of the energy with the size of basis set in subsection II B. This is followed by the summary of employed many-body approaches in subsection II C. Section III presents the results on the convergence with the size of the one-particle basis set in subsection III A, the convergence with the excitation level included in the coupled cluster Ansatz for the wave function in subsection III B, the density profiles in subsection III C, and the comparison with experiments in subsection III D. Section IV summarizes our paper and discusses possible future applications of the developed approach and further developments. Details of the analytic solutions of the two-body case of subsection II B, and the derivations of the corresponding convergence laws with size of the basis set are presented in appendices A-C.

A. The many-body Hamiltonian
The Hamiltonian describing a system of N spin-1/2 fermions (atoms) in a one-dimensional harmonic trap readŝ where x i represents the coordinate of the i-th 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 ω = m = = 1. This amounts to measuring energies E in units of ω, lengths in units of the harmonic oscillator characteristic length a ho = /(mω), and the interaction strength g in units of ωa 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 ↑ +N ↓ describes the total number of atoms (N ↑ atoms with the spin projection 1/2, and N ↓ atom with the spin projection −1/2). We will consider either systems with N ↑ = N ↓ , that is, in the simplest spin singlet S = 0 state, or ones with N ↑ = N ↓ such that the total spin S = |N ↑ − N ↓ |, 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: 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.
B. The two-body case and the convergence with the size of the one-particle basis set 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 where 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 = ǫ + n + 1, while U (a, b, x) stands for the Tricomi function [80]. The energy of the relative motion ǫ is determined from where Γ(.) 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 ǫ n b converging to the exact value ǫ as (see appendix B for details) 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. Ref. [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 > 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) where x = g −1 (ǫ), as given by Eq. (4), x ′ ≡ ∂ ǫ [g −1 (ǫ)], and A = −1/( √ 2πx ′ ). 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 Eq. (5) and Eq. (6) are numerically confirmed to be the same.
C. Summary of the many-body approaches 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 ↑ − N ↓ |, 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 (cf. [84]). We will denote the N functions occupied in the reference Slater determinant by {φ α } N α=1 , while the remaining one-particle functions that are not occupied in the reference determinant will be denoted by {φ ρ } 2n b ρ=N +1 , where n b is the number of the harmonic oscillator functions used to expand the oneparticle solutions of the Hartree-Fock equations. Note that the total number of one-particle functions is 2n b instead of n b since φ α and φ ρ include the dependence on the spin coordinate through the spin functions. Obviously, the relation n b ≫ N must hold. There is a theorem stating that if a set of one-particle functions spans the one-particle Hilbert space L 2 (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 full configuration interaction method the wave function of the many-fermion system is represented as where the CI operatorĈ may be written as a linear combination of l-particle excitation operators, with coefficients C α1...α l ρ1...ρ 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 Eq. (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 {φ α } N α=1 ⊕ {φ ρ } 2n b ρ=N +1 set of the one-particle functions. The l-particle excitation operator e ρ1...ρ l α1...α l is given by the product of one-particle excitation operators e ρ1 α1 · · · e ρ l α l . The one-particle excitation operators are in turn defined in terms of the conventional fermionic creation and annihilation operators e ρ α = a † ρ a α , where the creation and annihilation operators are defined with the respect to the physical vacuum |0 , but applied to the Fermi vacuum Φ [86]. Note that the one-particle functions φ α and φ ρ 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 φ α 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 α1...α l ρ1...ρ l and the energy are obtained by diagonalizing the Hamiltonian matrix constructed from the matrix elements H lm = e ρ1...ρ l α1...α l Φ|Ĥ|e ρ1...ρm α1...αm Φ . 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 (2n b ) N +2 for N ≪ 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 Eq. (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 coupled cluster (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 where the cluster operatorT is given byT Note that by comparison with Eq. (7) the following relation between the FCI and cluster operator must hold,T = ln(1 +Ĉ), but since the algebra of the excitation operators is nilpotent, the logarithm, as well as the exponential function in Eq. (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 while the cluster amplitudes t α1...α l ρ1...ρ l are obtained by solving the following system of nonlinear algebraic equations 0 = e ρ1...ρ l α1...α l Φ|e −TĤ eT Φ , l = 1, ..., N.
Note that by virtue of the Baker-Campbell-Hausdorff formula the exponential factor e −TĤ eT reduces to multiple commutators ofĤ andT , 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 coupled cluster 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, Eq. (9), the CC method truncated to single and double excitations effectively includes triply, quadruply, and higher excited determinants through the products, e.g.T 1T2 orT 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][92][93]. The computational cost of these methods scales as n 6 b , n 8 b , and n 10 b , 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 7 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. 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 many-body perturbation theory (MBPT) [63], or Møller-Plesset perturbation theory (MPPT) [95][96][97], it is useful to rewrite the equations for the cluster operatorT in the following operator form [98,99] T whereŴ is the correlation (fluctuation) potential in the Møller-Plesset partitioning of the Hamiltonian F is the Fock operator, and [Ŵ ,T ] n denotes the n-fold nested commutator: The nested commutator expansion in Eq. (13) terminates after the quadruple commutators since in our case the operator W contains only two-particle interactions. The resolvent superoperatorR is defined for arbitrary operator X as [98,99] where and ǫ κ 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 superoperatorsR n are always well defined. Note that for a givenX the operator Y =R n (X) can be viewed as a formal solution of the equation Obviously, this solution is unique if we assume thatŶ belongs to the linear span of the excitation operators e α1...αn ρ1...ρn . The many-body perturbation expansion of the energy is obtained by introducing the following parametrization of the correlation operatorŴ ,Ŵ where the complex parameter λ is introduced to derive the perturbation expansions of Eqs. (11) and (13). Obviously, the physical value of λ is equal to one. By substituting the parametrized Eq. (18) into Eq. (14) the cluster operator became dependent on λ, and can be expanded as a power series in λ Substituting the above expansion in the expression for the energy (11) and collecting all terms at λ n gives the expression for the nth-order correction to the energy in the many-body perturbation theory. The cluster operators necessary to evaluate this expression are obtained by substituting Eq. (19) and again collecting all terms at λ n . It follows immediately from this short sketch of the derivation that the MBPT and CC theories are strongly connected. The analysis reported in Ref. [99] shows that the CCSD, CCSDT, and CCSDTQ are valid through the third, fourth, and fifth order of the many-body perturbation theory. 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 0 it may be obtained as the expectation value of the operator N i=1 δ(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 φ i (x 0 )φ j (x 0 ) added to the ij-th element of the one-particle Hamiltonian matrix.
The full configuration interaction and coupled cluster 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 ↑ = N ↓ case, whereas orbitals obtained with the unrestricted Hartree-Fock method are employed for the N ↑ = N ↓ case.

III. NUMERICAL RESULTS AND DISCUSSION
A. Convergence with the size of the one-particle basis set 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 Eq. (5) and Eq. (6) for the two-body problem hold, we compare their predictions with the exact results in Fig. 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 [Eq. (B5) in appendix B] and variational energies equivalent to exact diagonalization results. The former one should be described by the asymptotic expansion of Eq. (5), whereas the latter ones should coincide with the asymptotic expansion of Eq. (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.
To evaluate the adequateness of the functional formula given by Eqs. (5) and (6) for extrapolation of the variational energies to the complete basis set limit, Fig. 1 presents also the fit of the variational results to the formula E ∞ + An . This extrapolation formula behaves extremely well [see also the inset of Fig. 1]. For instance, the extrapolated energies in the complete basis set limit E ∞ for g = 1 and g = 5 are 1.306744 and 1.726780, whilst the exact values equal to 1.306745 and 1.726771, respectively. Note that in the biggest basis set available, n b = 200, the calculated energies are 1.310545 and 1.745325, respectively, so that the accuracy gain of three order of magnitude is impressive and very important.  5) and (6). Results with inclusion of the leading-order term only are given as the short-dashed lines, results with inclusion of the first two terms are given as the dashed lines, and with inclusion of all three terms are given as the dot-dashed lines.
The solid black line is the fit of the variational results to the formula, E∞ + bn and the inset shows the performance of this fit over full set of data.
Absolute error with respect to the complete basis set limit of the ground state energy of the 2+2 (a) and 3+3 (b) systems as a function of the one-particle basis set size for several interaction strengths obtained with the FCI method. Lines correspond to fits to the extrapolation function given by Eq. (20).
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 (c.f. 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 where E ∞ is the extrapolated energy to the infinite number of the basis functions, E n b 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 Fig. 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 (N ↑ + N ↓ means N ↑ atoms with the spin projection 1/2, and N ↓ 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 Eq. (20). The agreement between the computed points and the analytical fits to the extrapolation function is excellent, supporting the correctness of our extrapolation scheme. 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 Fig. 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 Ref. [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 Eq. (20), especially in the limit of intermediate and large interaction strength. The apparent convergence of the results observed by the Authors of Ref. [48] is solely due to the pathological convergence pattern as n −1/2 b . 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.
As seen in Figs. 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 Tonks-Girardeau (TG) gas (when the interaction strength is infinite, g = ∞), 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 Ref. [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 = ∞. 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 ω 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 Tonks-Girardeau 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.
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 Fig. 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 Fig. 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 n b 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 Refs. [107,108].
B. Convergence with the excitation level included in the coupled cluster Ansatz for the wave function 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, Eq. (10), on the correlation and total energies. We start the discussion of our results with Fig. 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 Fig. 6 shows that the basis set dependence of the energy is relatively weak in the weakly correlated repulsive regime (g > 0), and very pronounced in the strongly correlated attractive case (g < 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 Tonks-Girardeau 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.
The correlation energy presented in Fig. 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 coupled cluster method), which start from the mean-field solution and tend to recover correlation energy.
We now turn to the applicability of the many-body perturbation theory and truncated configuration interaction expansions for the energy calculations. In Fig. 7 we report the percentage of the ground state correlation energy reproduced with the second-order and fourth-order many-body perturbation theory, 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 TABLE I. The percentage of the ground state correlation energy reproduced at the CCSD, CCSD(T), CCSDT, and CCSDTQ levels of the coupled cluster 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. 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 many-body perturbation theory is divergent [109][110][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 4 b for MBPT2 vs. n 7 b for MBPT4. Finally, we note that the variational CISD results do not offer a good advantage over the MBPT results.
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 coupled cluster theory will improve the situation. This is indeed the case. Fig. 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 coupled cluster 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 Fig. 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 Refs. [103,114]. Comparison of Figs. 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.
The discussion of the last two paragraphs is well summarized in Table I, 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 coupled cluster theory, as opposed to the erratic behavior of the many-body perturbation theory, and the poor performance of the configuration interaction method limited to single and double excitations.
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 coupled cluster 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 many-body perturbation theory. 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 coupled cluster wave function does not introduce any significant errors. The energy spectrum of the investigated systems for strong repulsive interactions becomes quasi-degenerate (cf. Figs. 3 and 4) and the correlation energy diverges for both strongly attractive and repulsive interactions (cf. Fig. 6). These two effects lead to the convergence problem of the standard single-reference coupled cluster 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 coupled cluster theory or the ferromagnetic reference state for calculations in the limit of the strong repulsive interaction.

C. Density profiles
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 Ref. [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 Tonks-Girardeau 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 LL (x) = 2 N/2−1 i=0 |φ i (x)| 2 for the Lieb-Liniger gas of hard-core bosonic dimers at g = −∞,  |ϕ 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. 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 Ref. [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 = −∞) and strong repulsion (g = ∞) is clearly visible in Fig. 11(b). Finally, the CC method allows one to address much larger systems than the FCI approach.

D. Comparison with the available experimental data
All the results reported thus far strongly suggest that the approximate variants of the coupled cluster 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 shortrange 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 Fig. 12. We decided to compare the best of our results, i.e. the FCI results, but any variant of the coupled cluster 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 ↑ + 1 atoms consisting of a single impurity interacting with a Fermi sea of N ↑ 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 ∆ E conservatively estimated from the extrapolated energies and the energies computed with the largest n b = 200 basis functions, ∆ E = E n b =200 − E ∞ . 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.
Less satisfactory is the agreement between theory and experiment for the separation energies ∆E(N ) = µ(N ) − µ * (N ), where µ(N ) = E(N ) − E(N − 1) is the chemical potential, E(N ) is the extrapolated energy for the Nbody system in the weakly attractive regime, and µ * (N ) is the chemical potential of the noniteracting system. This is illustrated in the panel (c) of Fig. 12, where we compare our theoretical results with the experimental data of Ref. [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 Ref. [33] is different than in Ref. [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 Ref. [54].

IV. SUMMARY AND CONCLUSIONS
In this paper we have reported the first application of the ab initio methods of quantum chemistry to systems of many interacting fermions in a one-dimensional harmonic trap. Our results can be summarized as follows: • The behavior of the two-body energy as a function of the number of single-particle functions n b included in the calculations for a fixed interaction strength has thoroughly been analyzed and an extrapolation formula has been derived. The convergence with n b is pathologically slow and is reported in the leading order as n • The convergence of the interaction strength as a function of the number of single-particle functions n b for calculations at fixed energy in the two-body case have also been analyzed and an extrapolation formula which does not depend on the energy and on the number of particles in the lowest order has been derived.
• The convergence with the number of basis set functions n b in the many-body case has been analyzed and a suitable three-point extrapolation formula has been found.
• The importance of the use of converged results to correctly describe the physics of the trapped Fermi systems has been pointed out. Shortcomings of the approaches lacking a proper analysis of the convergence issues have been shown.
• Comparison of several quantum many-body methods has been reported. In particular, a careful consideration has been given to the level of the excitation present in the wavefunction in the coupled cluster method. An optimal method has been chosen and used throughout the rest of the study.
• The limitations of the adopted coupled cluster implementation with the standard single-reference state have been discussed together with the analysis of the behavior of the decomposition of the total energy into the mean-field and correlation parts in the limits of strong repulsive and attractive interactions.
• The coupled cluster method restricted to single, double, and noniterative triple excitations CCSD(T) has been applied to describe systems of up to 80 fermionic atoms.
• Density profiles have been computed and analyzed in the full spectrum of the interaction parameter and the number of atoms. The transition between the Lieb-Liniger gas of hard-core bosonic dimers and the Tonks-Girardeau gas has been observed with the FCI method.
• Comparison with the available experimental data including estimates of the computational uncertainties has been provided. Very good agreement between the theory and experiment has been pointed out for the precision measurements based on the RF spectroscopy whereas the perturbed character of the 1D harmonic shape of the trap during measurements based on the observation of the tunneling of particles out of the trap has been confirmed by confrontation with our accurate ab initio calculations with estimated error bars.
The presented numerical approach allows us to get an important insight into the ongoing experiments on the trapped fermions. By extending the number of atoms from less than ten within the exact diagonalization approach (see, e.g., Ref. [56]) to many tens within the coupled cluster method, it is becoming possible to investigate the transition between few-and many-body systems. This is particularly important for a good understanding of the emergence of bulk matter properties, crucial across all areas of physics.
Recently, we have applied the numerical approach developed here to investigate the properties of a two-component Fermi gas trapped in one dimension [115]. We have addressed the question of how the observables, such as the energy, the chemical potential, the pairing gap, and the density profile, evolve as the number of particles increases from very few to many tens. We have found that the energy converges surprisingly rapidly to the many-body result for all interaction strengths between minus and plus infinity, covering the whole range from the molecular bosonic Tonks gas to the atomic (fermionic) one. On the other hand, we observed the emergence of a non-analytic behavior of the pairing gap only when a substantially larger number of particles is present in the trap.
We believe that the results presented here and in Ref. [115] on several fermions in a harmonic trap will be followed by many applications of the proposed numerical approach to investigate interesting physics in different systems, geometries, and dimensions. We foresee that the coupled cluster method will describe atoms trapped in other potentials, such as double wells, arrays of microtraps, or even optical lattices equally well. The method should work also for more complex two-body interaction potentials, e.g., including van der Waals or long range 6dipole-dipole interactions. What is more, the method can be generalized to handle atoms trapped in two or three dimensions. The pathologically slow convergence with the number of single-particle basis functions included in the calculations can be overcome by using explicitly correlated basis functions in analogy with methods developed in quantum chemistry [116]. The use of the explicitly correlated methods can be crucial for an accurate description of systems in a higher number of dimensions. The range of interaction strengths which can be used in calculations with the single-reference coupled cluster method is restricted to intermediate values. This can be overcome by using a ferromagnetic (high-spin) reference state instead of the antiferromagnetic (low-spin) one which is typically used in the standard coupled cluster method. This would allow to describe many tens of interacting atoms in the strongly repulsive regime, and therefore would be ideal for describing atomic clouds in the Tonks-Girardeau limit. Finally, one can employ the coupled cluster method based on the high-spin reference state to describe many strongly-interacting bosonic atoms and investigate in detail the fermionization process in such a system.
In these variables the exact wavefunction, Ψ(x, X), separates into a product of functions φ and ξ dependent only on x and X, respectively. These functions obey the following differential equations where primes denote the usual differentiation over the corresponding variables. The equation for the centre-of-mass motion can immediately be solved, as it coincides with the Schrödinger equation for the quantum harmonic oscillator. Further in the text we shall use the shorthand notation,Ĥ x = − 1 2 ∂ 2 x + 1 2 x 2 +ḡδ(x) − 1/2, so thatĤ x φ(x) = ǫφ(x). The Hamiltonian,Ĥ x , is invariant with respect to the spatial inversion, i.e., x → −x. Therefore, its eigenfunctions possess a definite parity (even or odd). Eigenfunctions that are of odd parity must vanish at x = 0.
A substitution φ(x) = e −x 2 /2 f (x) in the first equation of (A2) gives the corresponding differential equation for f (x) First, we solve the homogeneous differential equation, i.e., without the Dirac delta term on the right-hand-side, which is well-known and the general solution can be written as a linear combination of two functions where C U and C M are some constants necessary to satisfy the initial conditions, and M and U are the Kummer and Tricomi hypergeometric functions [80], respectively. Returning to the original (inhomogeneous) equation, we consider first the odd eigenfunctions. As noted beforehand, they vanish at the origin (x = 0) and are not affected by the presence of the Dirac delta source. Therefore, for the odd states the solution of the quantum harmonic oscillator is simply obtained. The problem of even states is more pronounced as the ground state is nodeless and even. A detailed investigation of the differential equation (A3) reveals that the even solutions are also written in terms of M and U , but the initial conditions need to be chosen properly to account for the presence of the Dirac delta distribution. To find the required initial condition one integrates both sides of Eq. (A3) on a small interval around the origin, i.e., (−ε, +ε). Subsequently, all terms that cannot contribute in the small ε limit are dropped. One can easily estimate which terms are significant by noting that both solutions of the homogeneous equation are regular around x = 0. Finally, the following result is obtained Let us recall that around x = 0 the solutions of the homogeneous equation behave for x > 0 like: The functions M are unable to satisfy the above initial condition. Therefore, we must pick up the solution expressed through the Tricomi functions U and from Eq. (A5) we obtain The above expression cannot be solved with elementary methods. Nonetheless, for a givenḡ the solution is straightforward to obtain numerically. Finally, the exact wavefunction of the Hamiltonian (A1) is given by where C n and C ǫ are the normalization constants for the x-and X-dependent portions, and the corresponding total energy of the system is simply E = ǫ + n + 1.
Having the exact wavefunction of the system, we can analyze the possible cusp-like conditions at the particles coalescence points (x = 0). The wavefunction around x = 0 behaves as which is a direct consequence of the properties of U . The above expression can considerably be simplified with the help of Eq. (A10) which leads to Strikingly, the behavior of the exact wavefunction around x = 0 depends only on the value ofḡ. This is also the counterpart of the Kato's electron-electron cusp condition. In analogy, Eq. (A13) can be written as where the equality sign strictly holds for any X. Due to the requirement of even parity, the exact wavefunction possesses a derivative discontinuity at x = 0 and a "cusp" at the origin. This discontinuity can be a considerable difficulty from the practical point of view. When the wavefunction is modeled with smooth basis set functions, the cusp may be very difficult to reproduce.
Appendix B: Energy convergence with the basis set size -truncated exact wave function In the two-body case the expansion of the exact wavefunction in one-dimensional harmonic oscillator eigenfunctions can be written as By recalling the form of the exact wave function, Eq. (3), one finds that the X-dependent part of the wavefunction is automatically described exactly. Therefore, the actual task in the two-body calculations is to approximate the x-dependent part of the wavefunction, φ(x), by the solutions of the quantum harmonic oscillator eigenproblem, i.e., Due to the orthonormality of ϕ n (x) the coefficients c n obey the formal relationship Note that since the exact wavefunction is even, the coefficients vanish for odd n by symmetry.
As the basis set, Eq. (2), is complete in the second Sobolev space, we can construct systematic approximations to the exact wavefunction by terminating the expansion (B2) at some n b . This gives a family of approximants which are not normalized to the unity. The energy connected with φ (n b ) (x) for each n b is given by Since the exact wavefunction, φ(x), is normalized, the exact energy of the relative motion is simply ǫ = φ|Ĥ x |φ . Let us also define the complementary function, φ (n b ) c (x), given for each n b by the expression We can now precisely state that we are interested in finding the asymptotic expansion of ǫ n b − ǫ as n b → ∞. At this point we must stress that the energy ǫ n b is not strictly equal to the energy calculated in the same basis set. In fact, in the actual finite basis set calculations, the coefficients c n are determined variationally rather than by projection onto the exact wavefunction. Therefore, the variationally determined coefficients have an opportunity to relax and accommodate to the incompleteness of the basis set. However, for large n b this relaxation effect is small and virtually limited to the last coefficient, c n b −1 , as shown by Carroll for the helium atom [117]. Therefore, the coefficients c n and the energies ǫ (n b ) are expected to be very close to the corresponding variational values. Let us recall Eq. (B3) and insert the explicit form of ϕ n (x) To evaluate the integrals v n we need to recall a particular integral representation of U (a, b, z) which is valid for a > 0. Therefore, we confine ourselves here to the regime ǫ < 0 which roughly corresponds toḡ < 0 (attractive potential). However, one can show that the final result derived here remains valid also for ǫ > 0. By inserting the above expression into the definition of v n and exchanging the order of integrations, one finds Let us consider only the inner integral for a moment, denoted by I. The change of variables u = x √ 1 + t gives To simplify the Hermite function under the integral sign one can make use of the following relation where γ = (1+t) − 1 2 . Upon inserting into the integral I and making use of the orthogonality of the Hermite polynomials one arrives at By returning to the definition of v n and slightly rearranging the following formula is obtained and the remaining integral is elementary, so that Therefore, the expression for c n reads c n = C ǫ π 1/4 (−1) n/2 Γ(− 1 2 ǫ) and the coefficients vanish for odd n, i.e., c 2n+1 = 0. The above expression is fairly difficult to handle because of the presence of the factorials. Fortunately, we require only their asymptotic forms, i.e. approximately valid for large n.
With the help of the Stirling formula, n! → √ 2πn n n e −n , one arrives at The first ingredient required for the presented derivation is the asymptotic formula for the square of the norm, φ (n b ) |φ (n b ) , for large value n b . Taking advantage of the orthonormality of ϕ n and normalization of the exact wavefunction one finds Since the second term of the above expression vanishes for large n b , one can write that φ (n b ) |φ (n b ) → 1 as n b → ∞ which is entirely sufficient for the present purposes. One can also verify that a somewhat more accurate formula, accounting for the next term in the asymptotic expansion, would include a term proportional to 1/n 3/2 b . However, this term does not contribute to the final results and is omitted for brevity.
Because of the approximation derived for the denominator in Eq. (B5) one can rewrite it for large n b as By recalling the definition (B6) and rearranging one gets The first of the above matrix elements can be expanded to give where we have additionally used the relationship ϕ 2n (0) → (−1) n √ π n −1/4 as n → ∞, easily obtained from Eq. (2) by using the Stirling approximation. Finally, the above sums (taking into consideration all relevant terms in the Stirling series) can be evaluated with the help of the Euler-MacLaurin resummation formula which leads to Let us also mention that with the same methodology as presented above one can derive further terms in the asymptotic expansion of ǫ n b − ǫ. For example, the expression including additionally terms proportional to n b −3/2 reads Note that the convergence proportional to 1/ √ n b in the leading order is extremely slow. The latter convergence formula may be employed to describe numerical results in the weak and intermediate coupling regime. In the next appendix we will instead outline a different approach, which yields a convergence formula valid for all coupling strengths.
Appendix C: Energy convergence with the basis set size -variational wave function In Appendix B we have considered the problem of the convergence of the calculations with the increasing basis set size for the two-particle case by calculating the energy given by expectation value (B5) with the truncated exact wave function (B4) for a finite interaction strength g. Now, we will present an alternative derivation, allowing for a variational relaxation of the coefficients c n in the truncated wave function (B4). Obviously, since the functional that we minimize is quadratic, the resulting wave function corresponds to the exact ground state wave function of the truncated (projected) Hamiltonian. In the following, we shall first consider the numerical convergence of the inverse coupling constant versus n b , at fixed energy ǫ, and finally derive the convergence of the energy at fixed coupling.
By minimizing the expectation value of the HamiltonianĤ x on the wavefunction (B4) with respect to its coefficients c n , one finds where 2 1/4 f n = k ϕ n | k = ϕ n (0) is the n-th harmonic oscillator wavefunction for the relative coordinate, evaluated at x = 0. We have f 2n+1 = 0, and (f 2n ) 2 = 1 √ 2π (2n − 1)!! 2n!! = 1 √ 2π (2n)! 2 2n (n!) 2 . (C2) By solving Eq. (C1) for g in a procedure similar to that of Busch et al. [78], one finds the result presented in Eq. (4) of the main text, Let us now define x(ǫ) ≡ g −1 (ǫ). In numerical calculations, one has to truncate the basis to a certain maximum number of states n b , and therefore obtains only an approximate value, The convergence rate of this formula at a fixed energy ǫ, as a function of the number of states included in the summation, n b , is given by as may be found using Stirling's formula, and then expanding the fraction inside Eq. (C3) for 2n ≫ |ǫ|.
Let us now consider the numerical error ∆ǫ n b =ǫ n b − ǫ obtained by computing the energy at a fixed g with a basis set containing up to and including n b functions. The error is given by the solution of the implicit equation which can be rewritten as On the left-hand-side of the above equation we can now use Eq. (C5). Assuming that ∆ǫ n b can be expanded in powers of 1/ √ n b and equating the coefficients on both sides, we obtaiñ with A = −1/( √ 2πx ′ ), and x ′ ≡ ∂ ǫ [g −1 (ǫ)]. In the vicinity of the Tonks-Girardeau point where ǫ ≈ 1 we have ǫ = 1 − g −1 8/π + . . ., so that A = 2/π 3/2 . In the weak coupling limit ǫ ≈ 0 instead, we have ǫ = g/ √ 2π + . . . so that A = g 2 /(2π 3/2 ). Both these limits coincide with the analytical results of Ref. [53].