Are smooth pseudopotentials a good choice for representing short-range interactions?
AAre smooth pseudopotentials a good choice for representing short-range interactions?
Péter Jeszenszki,
1, 2, 3, ∗ Ali Alavi,
3, 4, † and Joachim Brand
1, 2, 3, ‡ Dodd-Walls Centre for Photonics and Quantum Technology, PO Box 56, Dunedin 9056, New Zealand New Zealand Institute for Advanced Study, and Centre for Theoretical Chemistry and Physics,Massey University, Private Bag 102904 North Shore, Auckland 0745, New Zealand Max Planck Institute for Solid State Research, Heisenbergstraße 1, 70569 Stuttgart, Germany Department of Chemistry, University of Cambridge,Lensfield Road, Cambridge, CB2 1EW, United Kingdom (Dated: March 7, 2019)When seeking a numerical representation of a quantum-mechanical multiparticle problem it istempting to replace a singular short-range interaction by a smooth finite-range pseudopotential.Finite basis set expansions, e.g., in Fock space, are then guaranteed to converge exponentially. Theneed to faithfully represent the artificial length scale of the pseudopotential, however, places a costlyburden on the basis set. Here we discuss scaling relations for the required size of the basis set anddemonstrate the basis set convergence on the example of a two-dimensional system of few fermionswith short-range s -wave interactions in a harmonic trapping potential. In particular we show that thenumber of harmonic-oscillator basis functions needed to reach a regime of exponential convergencefor a Gaussian pseudopotential scales with the fourth power of the pseudopotential length scale,which can be improved to quadratic scaling when the basis functions are rescaled appropriately.Numerical examples for three fermions with up to a few hundred single-particle basis functionsare presented and implications for the feasibility of accurate numerical multiparticle simulations ofinteracting ultracold-atom systems are discussed. Keywords: Gauss potential, harmonic trap, exact diagonalization, exponential convergence
I. INTRODUCTION
There is increasing interest in the theoretical descrip-tion of multiparticle systems of interacting ultracoldatoms thanks to the recent progress in experimental re-alizations [1–5]. In particular we may expect excitingdevelopments in microtraps [6–8] with tens of particleswhere accessing strongly correlated regimes of quantum-Hall-like physics seems feasible [9–12].The theoretical description of atom-atom interactionsis significantly simplified at ultracold temperatures wheredetails of the interaction potentials can be neglected infavor of a single constant, the s -wave scattering length a s ,to define a physical model with contact interactions [13].Despite these simplifications, the complexity of many-particle quantum mechanics still makes it a very difficultproblem to solve, where exact solutions are only availablein special cases in one spatial dimension [14–18] or for upto three particles in a harmonic trap [19–21].A straightforward and generalist approach to repre-senting the many-body problem for computational treat-ment is to introduce a discrete and necessarily finite ba-sis of smooth single-particle wave functions from which afinite but still potentially very large Fock-space is con-structed to represent the many-body Hamiltonian asa matrix. Finding eigenstates and eigenvalues of thefull matrix is known as exact diagonalization or full ∗ [email protected] † [email protected] ‡ [email protected] configuration-interaction [22–26], but many different ap-proximation schemes have also been followed [27]. In par-ticular, standard approaches of ab initio quantum chem-istry or nuclear physics like the coupled-cluster [28] ormulti configurational self-consistent field theory [29] allcan be formulated in this language as they rely on anunderlying single-particle basis. Also Monte Carlo (orother) approaches that rely on a lattice discretization ofcontinuous space fall into the same category [30], as theunderlying single-particle space can be represented as adiscrete set of plane waves.One of the complications in the numerical treatmentof contact interactions with basis set expansions stemsfrom the nonanalytic behavior of the wave function atthe point of particle coalescence. At this point, the ap-propriate Bethe-Peierls boundary conditions demand acusp in one spatial dimension, i.e. a point of nondiffer-entiability [14], in two dimensions a logarithmic diver-gence, and in three dimension a /r divergence of thewave function [31–33]. While in one dimension a Dirac δ function pseudopotential provides a well-defined modelfor contact interactions, the convergence of basis set ex-pansions is algebraic and painfully slow [28, 34]. Basisset expansions in two and three dimensions diverge forbare contact interaction [35–37] and basis-set-dependentrenormalization procedures have to be used in order toobtain convergent and correct results [37–40]. In the bestcase, renormalized contact interactions will lead to alge-braic convergence in the size of the finite single-particlebasis set [30, 33, 34]. Some of us have recently described atranscorrelated method where the singular nature of thecontact interaction is reduced by means of a similaritytransformation of the Hamiltonian [34] (see also Ref. [41] a r X i v : . [ c ond - m a t . qu a n t - g a s ] M a r for a related idea). This promising approach results inan improved power-law scaling but the convergence stillremains algebraic.Here, we consider a different approach where the con-tact interaction is replaced by a smooth finite-range pseu-dopotential. This has the advantage that basis set expan-sions will converge exponentially for appropriately chosensingle-particle basis functions [42]. Specifically we con-sider what the requirements are for the basis set to reachthe regime of exponential convergence and whether thisapproach is feasible for multiparticle simulations. Ex-amples of finite-range pseudopotentials used in the lit-erature are the Troullier-Martins [43, 44], Pöschl-Teller[45, 46], and Gaussian potentials [26, 36, 47–57], in addi-tion to the square well popular in diffusion Monte Carlosimulations[48].When using finite range pseudopotentials to representshort-range interactions, an interpolation in the widthof the pseudopotential should be done to the zero-rangelimit [58]. In order to approach this limit, the length scaleof the pseudopotential should be significantly smallerthan other physically relevant length scales of the prob-lem, in particular the mean particle separation and lengthscales imposed by external potentials. In order to reach aregime where the basis set expansion converges exponen-tially, however, the basis set needs to resolve the smallestlength scale of the pseudopotential. At the same time,the large length scales of the problem, i.e., the (Thomas-Fermi) size of the cold atom cloud, or the size of thecontainer, also have to be represented by the basis set.This hierarchy of length scales, typically spanning at leastone but possibly several orders of magnitude, presents achallenge for accurate numerical simulations. While thesize of the single-particle basis (quantified by the num-ber of single-particle functions M ) is determined by thishierarchy of length scales, the size of the full many-bodyproblem also depends strongly on the number of particles N . Specifically for spinless bosons, the size of the rele-vant part of Fock space is (cid:0) N + M − M (cid:1) , whereas for spin- fermions the total dimension is (cid:0) MN ↑ (cid:1)(cid:0) MN ↓ (cid:1) , where N ↑ and N ↓ are the numbers of up- and down-spin particles, re-spectively.In this work, we specifically consider ultracoldfermionic atoms in a harmonic oscillator trapping po-tential where the potential in one of the three trappingdirections is so tight that the problem can be consid-ered two-dimensional. We furthermore choose a Gaussianpseudopotential to model attractive s -wave interactionsbetween spin-up and spin-down particles [59]. For theunderlying single-particle basis we consider two cases:(1) a basis that is defined by the single-particle eigen-states of the isotropic two-dimensional harmonic trap-ping potential, and (2) the same set of basis functionswith scaled spatial coordinates by a scaling factor γ . Us-ing the known properties of the harmonic oscillator eigen-functions we show that the basis set size M required toresolve the chosen length scale of the pseudopotential l res scales as ( l/l res ) where l = (cid:112) (cid:126) /mω is the harmonic os- cillator length scale of the trapping potential for case (1).Allowing the basis functions to be scaled by γ as incase (2) leads to an improved scaling of ( l/l res ) whilestill faithfully resolving the small length scale l res and afixed large length scale that is determined by the particle-number and interaction strength. We provide estimatesfor these length scales and the required scaling param-eters. Numerical examples show the convergence of theground-state energy for three fermions obtained by exactdiagonalization with single-particle basis sets of up to 231Fock-Darwin orbitals. In order to compute the matrixelements of the Gaussian interaction potential with theFock-Darwin basis of this size, a careful algorithm basedon recursion formulas had to be developed in order toavoid an excessive accumulation of round-off errors. Thisalgorithm is described in Appendix B.This paper is organized as follows: After defining theHamiltonian in Sec. II and introducing the methodologyin Sec. III, we discuss the main results of the paper inSec. IV. Examples for the numerical convergence witha harmonic oscillator basis are presented in Sec. IV Abefore deriving analytical formulas for the required min-imum basis set size in Sec. IV B for the unscaled andin Sec. IV C for a scaled harmonic oscillator basis. Therequired scaling factor is estimated in Sec. IV D wherealso numerical results for the scaled basis are presented.Implications of our findings for the feasibility of accuratecomputations of larger multiparticle problems are dis-cussed in Sec. V. Two appendixes define the Fock-Darwinorbital basis used (Appendix A) and detail the explicitformulas and the algorithm used to compute the matrixelements (Appendix B). II. HAMILTONIAN
We consider ultracold fermions in a two-dimensionalharmonic trap, H = H osc + N ↑ (cid:88) i N ↓ (cid:88) j V ( r i ↑ , r j ↓ ) , (1) H osc = (cid:88) σ N σ (cid:88) i (cid:18) − (cid:126) m ∇ iσ + mω r iσ (cid:19) , (2)where m is the mass of the fermions, ω is the harmonicoscillator strength, r iσ is the position of the i th particlewith the spin σ , and V ( r i ↑ , r j ↓ ) is the interaction poten-tial between the fermions. Dividing the operator H osc with (cid:126) ω , Eq. (2) takes the form H osc (cid:126) ω = (cid:88) σ N σ (cid:88) i (cid:18) − l ∇ iσ + 12 l r iσ (cid:19) , where l = (cid:113) (cid:126) mω is the harmonic oscillator length scaleof the trapping potential. The interaction between theparticles is described with a Gaussian pseudopotential V ( r i ↑ , r j ↓ ) = − V R e − ( r i ↑− r j ↓ ) R . (3)The parameters V and R control the strength and widthof the interaction potential, respectively. These param-eters can be converted to the s -wave scattering lengthusing a simple approximate formula [59], a s R = √ (cid:32) − γ E (cid:126) V m + n (cid:88) i α i V V − W i (cid:33) , (4)where a s is the s -wave scattering length in two dimen-sions, γ E is the Euler-Mascheroni constant with the ap-proximate value of γ E ≈ . , and α i and W i areparameters fitted to direct numerical calculations [60].Accurate numerical values of the parameters are given inTable III. in Ref. [59] for i ≤ . Alternatively, more com-plicated numerical approaches can be applied [59, 61]. III. BASIS-SET EXPANSION
For our numerical approach, we compute the ground-state energy of a multiple-fermion system following theexact diagonalization approach. Starting from a finitesingle-particle basis of size M (i.e., with M spin orbitals),the multiparticle wave function is expanded as a linearcombination | Ψ (cid:105) = (cid:88) n C n | Φ n (cid:105) , (5)of states in the Fock basis | Φ n (cid:105) = M (cid:89) i =1 (ˆ c † i ) n i | vac (cid:105) , (6)where n = ( n , . . . , n i , . . . , n M ) and n i is the occupationnumber of the n th single-particle basis function (spinrorbital). The fermionic Fock states | Φ n (cid:105) (often referredto as Slater determinants) are constructed from the com-plete set of index vectors n for fixed particle number N ,with N = M (cid:88) i =1 n i . (7)The exact diagonalization approach (also referred to asfull configuration interaction) refers to considering themultiparticle Hamiltonian (1) with the chosen particle-number content projected onto the basis of states (6) asa matrix and finding its eigenvalues and eigenvectors.In this paper we use a single-particle basis constructedfrom the spinful eigenstates of the two-dimensional har-monic oscillator (2). The explicit form of the basis func-tions used in the numerical procedure is presented in Ap-pendix A. Even though the one-body and two-body in-tegrals needed for the relevant matrix elements can be expressed analytically (see, e.g., Ref. [26]), obtaining ac-curate numerical values is challenging due to the prolifer-ation of rounding errors during floating-point arithmetic.We have therefore developed an iterative algorithm forthe evaluation of the two-body integrals, which alleviatesthis problem. The details are presented in Appendix B.For the numerical procedure we determine the ground-state energy E = min C n (cid:104) Ψ | H | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) ≡ (cid:104) H (cid:105) , (8)with a matrix-free approach: Using a variant of the powermethod [62], we iteratively rotate an initial state ontothe ground state vector without having to construct thematrix explicitly. Numerical computations are done us-ing the NECI [63] software in deterministic mode. Evenlarger Hilbert spaces could be explored using stochasticalgorithms for exact diagonalization such as Full Con-figuration Interaction Quantum Monte Carlo (FCIQMC)[64]. IV. RESOLVING THE GAUSSIANPSEUDOPOTENTIALA. Unscaled harmonic oscillator basis
Figure 1 shows the ground state energy of three in-teracting fermions (two spin-up and one spin-down) ac-cording to the Hamiltonian (1) after full diagonalizationwith a finite Fock basis. We are using the Fock-Darwinform (A1) of the (unscaled) harmonic oscillator eigen-functions of the single-particle Hamiltonian (2) up toshell n = 20 , which yields up to M = 231 single-particlebasis functions. The maximal dimension of the computa-tional Hilbert space for the three fermions with the zerototal angular momentum slot is ≈ . × , which isalready a significant size for the deterministic diagonal-ization with available computational resources.It is clearly seen in Fig. 1 that these basis set sizesare not sufficient to enter a regime of exponential conver-gence except for Figs. 1(a) and 1(b), where this regimeis reached for the last few data points as seen from theinsets. In these cases the width of the pseudopotential R ≈ l is close to the length scale of the trapping potential.Since the mean particle separation in the harmonic trapwill also be of the same order l , or even smaller, this pseu-dopotential does not provide a useful approximation forthe zero-range contact interactions that are relevant formodeling experiments with ultracold neutral atoms. Asthe lowest energy values reached for each pseudopotentialwidth R change significantly between the different panelsof Fig. 1, it is also apparent that R (cid:28) l is a necessarycondition for a useful, convergent approximation of thezero-range limit. Even without more sophisticated anal-ysis, it is apparent from the results of Fig. 1 that the nec-essary extrapolations the infinite basis set ( M → ∞ ) andzero-range ( R → ) limit will be challenging to achieve. -4.2-4.1-4-3.9-3.8-3.7-3.6-3.5 0 50 100 150 200 250 (a) E / − h ω M -6-5-4-3-2-1 0 0 50 100 150 200 250 l n [( E - E c ) / − h ω ] M -8-7.5-7-6.5-6-5.5 0 50 100 150 200 250 (b) E / − h ω M -3-2.5-2-1.5-1-0.5 0 0.5 1 0 50 100 150 200 250 l n [( E - E c ) / − h ω ] M -45-40-35-30-25-20-15-10 0 50 100 150 200 250 (c) E / − h ω M -40-35-30-25-20-15-10 0 50 100 150 200 250 (d) E / − h ω M Figure 1. Convergence of the ground-state energy for three fermions with attractive Gaussian pseudopotential interactions in theunscaled harmonic oscillator basis. The energy from exact diagonalization in the finite multiparticle basis of Eq. (6) is plotted vsthe size M of the single-particle basis for Gaussian pseudopotentials of different widths R : (a) R = l , (b) R = 0 . l , (c) R = 0 . l ,and (d) R = 0 . l , where l = (cid:112) (cid:126) /mω is the length scale of the harmonic trapping potential. The insets in Figs. 1(a) and 1(b)show logarithmic plots of the same data as the main graph and demonstrate that a regime of exponential convergence is reached.The extrapolated limiting values of the energy E c are obtained by nonlinear fitting of the exponential function Ae − BM + E c to the last three data points: (a) E c = − . (cid:126) ω and (b) E c = − . (cid:126) ω . The interaction strength ln( l/a s ) = 3 . is keptconstant for all panels and the corresponding amplitude parameters V are determined numerically following the proceduredescribed in Ref. [59]: (a) V = 19 . (cid:126) ω , (b) V = 19 . (cid:126) ω , (c) V = 18 . (cid:126) ω , and (d) V = 13 . (cid:126) ω . B. Length scale resolution
Since we are using a smooth Gaussian pseudopotential,we should expect that, for a sufficiently large basis set,the energy will converge exponentially to a limiting valuewith the size of the single-particle basis set. Indeed, itis well known that the Fourier transform of a Gaussianfunction yields again a Gaussian, which decays, in fact,faster than exponential in the tails. Sampling a Gaus-sian potential function in momentum space, should thuslead to at least exponential convergence, once the basisset is large enough to sample the tails of the Fourier-transformed Gaussian in momentum space. The neces-sary condition to reach this regime is that the basis set can resolve length scales that are smaller than the lengthscale R of the Gaussian.We now use this argument as a motivation to considerthe length scale resolution of a two-dimensional harmonicoscillator basis. In order to keep the basis set indepen-dent from the Hamiltonian of Eq. (2), we consider ba-sis functions that are eigenfunctions of a harmonic os-cillator with frequency ˜ ω and corresponding length scale ˜ l = (cid:112) (cid:126) /m ˜ ω . For a one-dimensional harmonic oscillator,the p th excited state has a spatial extent that can beestimated by the classical turning point x t : ( p + 12 ) (cid:126) ˜ ω = 12 m ˜ ω x t , (9)or x t = √ p + 1 ˜ l . The set of M = p + 1 harmonicoscillator functions up to the p th excited states providesan approximately homogeneous sampling of the interval [ − x t , x t ] at a length scale l res = 2 x t M = 2 √ p + 1 p + 1 ˜ l. (10)In order to connect this result to the number M of two-dimensional harmonic oscillator basis functions, we con-struct the latter as a product basis (of Hermite functions)with an energy cutoff. This yields M = i + j ≤ p (cid:88) i,j =0 p + 1)( p + 2)2 . (11)We can thus relate the resolution length scale l res to thesize of the basis and obtain l res = 4 (cid:112) √ M + 1 − √ M + 1 − l, (12)which can be solved for M to yield M ≈ (cid:32) ˜ ll res (cid:33) , (13)where lower order terms were neglected assuming ˜ l (cid:29) l res .Equation (13) provides an estimate for the size of thesingle-particle basis needed to resolve a length scale l res .For the situation of Sec. IV A where ˜ l = l we can estimatethe minimum size of the basis set to be able to resolvethe pseudopotential length scale R as M min ≈ (cid:18) lR (cid:19) , (14)i.e. the required basis set size increases rapidly when thelength scale R (and thus the range) of the pseudopoten-tial is decreased.Specifically, it means that reaching an exponentiallyconvergent regime should be quite achievable when thepseudopotential width is of the same order of magnitudeas the oscillator length l . For R = l and R = 0 . l wewould require a minimum of 32 and 78 single-particlebasis functions, respectively, which means about 16,000and 200,000 Fock states for three fermions. This is con-sistent with the numerical results of Fig. 1 obtained withup to M = 230 single-particle basis functions.In order to explore the physics of short-range interac-tions, however, we may need to use narrower pseudopo-tentials. With modest choices of R = 0 . l and R = 0 . l the number of required single-particle basis functions al-ready increases to about 4,000 and 300,000, respectively,corresponding to about × and × multiparti-cle basis functions. In this case, even the storage of thewave function is very expensive as it would require 0.2and 100,000 terabytes of memory, respectively. Althoughexploiting symmetries, using approximation schemes, orstochastic sampling techniques can reduce these require-ments [47, 64], the quartic scaling in Eq. (13) does notseem pleasant. C. Scaled harmonic oscillator basis
The required size of the single-particle basis can bedecreased by introducing an appropriate scaling of thebasis functions. The main idea is the following: Increas-ing the number of basis functions not only improves theresolution length scale l res (by making it smaller) butalso increases the largest length scale that can be de-scribed by the basis, which is given by x t = 2 √ p + 1 ˜ l ,i.e., through the classical turning point. This means thatthe basis functions can be scaled according to the basisset size while still resolving the largest length scale ofthe problem (e.g., twice the Thomas-Fermi radius of acold atomic cloud). Let us denote this largest relevantlength scale as γl , where γ is the dimensionless form ofthis length scale in units of the length scale l of the trap-ping potential. Note that γ is determined by the physicalproperties of the system (i.e., particle number, interac-tion strength, etc.) and is thus independent of the basisset size. For the required basis set length scale we thenobtain ˜ l = γ √ p + 1 l . (15)With Eq. (10) this yields l res = γl/ ( p + 1) and usingEq. (11) to solve for M we obtain M = γ (cid:18) ll res (cid:19) + γ ll res . (16)Applying this result to the problem of resolving aGaussian pseudopotential with length scale l res = R andtaking the leading term for l (cid:29) R , we obtain the revisedrelation for the minimum size of the single-particle basis, M min ≈ γ (cid:18) lR (cid:19) , (17)when the length scale of the single-particle basis is opti-mally scaled with the size of the basis set. Compared tothe result (14) for the unscaled basis, the power-law scal-ing of the required basis set size with the pseudopotentiallength scale R in Eq. (17) is improved by two orders. D. Estimating the interaction-dependent prefactor As γl represents the largest length scale that has tobe resolved by the single-particle basis, i.e., the maximalspatial extent of the system, the factor γ depends on thedetails of the Hamiltonian, which, for our example, arethe particle number content and the strength of the con-tact interaction. Although the exact value of γ is difficultto calculate we can estimate its value and present upperand lower bounds from limiting cases that are simple toanalyze.Specifically for fermions with attractive short-range in-teractions as per Eq. (1), the size of the trapped nonin-teracting gas cloud will provide an upper bound, while alower bound can be obtained from the strongly interact-ing limit where fermion pairs form deeply bound com-posite bosons while excess fermions remain with littleresidual interactions.Starting with the upper bound, we consider a nonin-teracting Fermi gas with a possibly unequal populationof spin-up and spin-down fermions. The largest lengthscale is determined by the Fermi pressure of the majoritycomponent with the quantum number p majority along asingle spatial dimension γ upper = 2 (cid:112) p majority + 1 . (18)The parameter p majority can be expressed in terms of par-ticle number N majority by considering Eq. (11) and usingthe information that only one particle can occupy onespatial orbital from the same spin component to yield γ upper = 2 (cid:118)(cid:117)(cid:117)(cid:116) (cid:38) (cid:112) N majority + 1 − (cid:39) + 1 , (19)where (cid:100) x (cid:101) is the ceiling function [65].For strong short-range attractive interactions, pairs ofspin-up and spin-down particles form point like effectivebosons [66, 67]. The interactions between the effectivebosons are repulsive but vanish in the limit of strong at-traction between fermions. In this limit we thus obtaina noninteracting system where all of the bosons occupythe lowest single-particle orbital. In the spin-balancedsystem, the largest length scale will thus be given by thelength scale of the harmonic oscillator trapping potential.This length scale provides a lower bound of the largestlength scale of the system as the effective repulsive in-teractions can only increase the average particle-particledistance. We thus obtain a lower bound of γ sblower = 2 , (20)where the index ‘ sb ’ stands for the spin balanced case.In the spin-imbalanced case, excess fermions from themajority component that are not locked-up in effectivebosons who will maintain Fermi pressure. Indeed, the ex-cess fermions have weak repulsive interactions with theeffective bosons that also vanish in the strongly inter-acting limit (regarding the original interactions betweenfermions) [66, 67]. Hence, the lower bound for the largestlength scale can be tightened by considering a noninter-acting Fermi gas of the excess fermions following the sameargument as above. We obtain γ silower = 2 (cid:118)(cid:117)(cid:117)(cid:116) (cid:38) (cid:112) | N ↑ − N ↓ | + 1 − (cid:39) + 1 , (21)where the index " si " stands for the spin-imbalanced case.We can see from Eqs. (19) and (21) that the largestrelevant length scale γl increases with particle number.To leading order the bounds become γ upper ≈ (cid:112) N majority , (22) γ silower ≈ (cid:113) | N ↑ − N ↓ | . (23) Comparing with expression (17) for the minimum sizeof the single-particle basis, we see that M min increaseswith the square root of the particle number (that is N majority and | N ↑ − N ↓ | , respectively), i.e., requiring alarger single-particle basis for a larger particle number.We also see that the largest interaction dependence of M min can be expected in the spin-balanced case, whereas M min will approximately remain independent of interac-tions for large spin polarization (e.g., polaron physics),where N majority ≈ | N ↑ − N ↓ | .For large particle numbers the large length scale γl is well approximated by Thomas-Fermi theory as γl =2 R TF , where R TF is the Thomas-Fermi radius. InThomas-Fermi theory, the single-particle density n ( r ) isfound from solving µ = V ext ( r ) + µ [ n ( r )] , (24)where V ext ( r ) = mω r is the external potential, µ [ n ( r )] represents the chemical potential at local equilibrium[from the equation of state µ ( n ) of the homogeneous gas],and the constant µ is the chemical potential of the fi-nite system [13, 68, 69]. The Thomas-Fermi radius R TF is then the value of r where n ( r ) drops to zero.For our case of a two-dimensional BCS mean-field the-ory yields a Thomas-Fermi radius R TF that is indepen-dent of particle interactions [70], i.e., the result (22). Inthe strongly interacting limit of the balanced Fermi gas,we may use the known asymptotic form of the equation ofstate for a two-dimensional gas of bosonic dimers [71, 72] µ ( n ) = − π (cid:126) m n ln ( a dd n ) , (25)where the scattering length of bosonic dimers a dd is ap-proximately a dd ≈ . a s [67, 73]. The Thomas-Fermiradius can be determined from Eq. (24). After some al-gebraic manipulations we obtain R TF ≈ (cid:115) a dd N d + 4 πl a dd ln (cid:18) πl e γ E +1 N d a dd (cid:19) , where N d is the number of the bosonic dimers and γ E isthe Euler-Mascheroni constant. The resulting approxi-mation for the scaling factor γ is γ = 2 R TF /l .In order to test the prediction of Eq. (17) for reducedrequirements for the size of a scaled single-particle basiswe consider again the example of two spin-up particlesand one spin-down particle with attractive interactions.The upper and lower bounds for γ can be easily cal-culated from Eqs. (19) and (21) to give γ upper = 2 √ and γ lower = 2 . The required minimum number of basisfunctions M min from Eq. (17) then yields 20 and 60, re-spectively, for R = 0 . l [reduced from M min ≈ R = 0 . l the number decreases from 300,000 to 200–600basis functions. In Fig. 2 we show the ground-state en-ergy R = 0 . l with different scaling factors γ . The resultsdemonstrate that the regime of exponential convergence -65-60-55-50-45-40-35-30-25-20-15-10 0 50 100 150 200 250 (a) E / − h ω M l~= l γ =4 γ =2 √ γ =2 √ γ =2-62.4-62.2-62-61.8 120 160 200 240 -4-3-2-1 0 1 2 3 4 0 40 80 120 160 200 240 (b) l n [( E - E c ) / − h ω ] M l ~= l γ =4 γ =2 √ γ =2 √ γ =2 Figure 2. Energy convergence in a scaled harmonic oscillatorbasis at R = 0 . l with different values of the factor γ . Forreference, we also show results for the unscaled basis (greencrosses) identical to Fig. 1(c). (b) A logarithmic represen-tation of the same data as in (a). The extrapolated limit-ing value of the energy E c is obtained by nonlinear fittingof the exponential function Ae − BM + E c to the last threedata points of the γ = 2 √ data with E c = − . (cid:126) ω . V = 18 . (cid:126) ωl corresponding to ln( l/a s ) = 3 . . can be reached with the scaled basis, even though it wasunattainable with the unscaled basis with the availablecomputational resources. The smallest scaling factors of γ = 2 and γ = 2 √ are seen to result in an (unphysical)energy bias. This can be understood by the fact that thelower bound γ = 2 and γ = 2 √ underestimate the sys-tem size and thus the scaled basis set does not cover thewhole area occupied by finite particle density. Increasingthe value of γ from the lower bound eliminates the biasbut also eventually leads to a reduced convergence rate.At the upper bound ( γ = 2 √ ) the computation is seem-ingly free of bias but convergence of the energy is greatlyimproved compared to the unscaled basis set. Hence,we find that the upper bound can safely be used for theaccurate determination of the ground-state energy. Extrapolating our results to larger particle numbers,we need to consider the following issues: First, addingmore particles at constant M increases the size of Hilbertspace (and thus computational complexity) due to the bi-nomial scaling with N and M . Whether the size M of thesingle-particle basis has to be changed depends on howthe relevant length scales change. The smallest lengthscale only depends on the choice of the pseudopotentialand is invariant with particle number. Thus the required M for the unscaled harmonic oscillator basis is indepen-dent of particle number (as long as M is larger than themajority-component particle number). For the scaled-basis approach we have to consider that the relevantlargest length scale may change. For example, addinganother particle in the minority component (say goingfrom three to four fermions where two are spin-up andtwo are spin-down) will compress the wave function to-wards the center of the trap for attractive interactionsand thus reduce the size of the relative wave function.Also the center-of-mass wave function will shrink due tothe larger total mass. Therefore, an even smaller valueof γ can be applied leading to smaller M min . However,if the number of the particles in the majority componentexceeds three particles, one of the particles in the nonin-teracting case occupies the next shell, which increases thesize of the system and the upper bound for γ increases to γ upper = 2 √ . In this case, the number of single-particlebasis functions increases to M min ≈ , which is stillmuch less than the required number of the unscaled ba-sis functions, M min ≈ . The increment of M min isless severe for larger particle numbers, again due to theeffect of attractive interactions. Using M min ≈ scaledsingle-particle basis functions, we can describe about 15particles in the majority component, which can be 30particles in the spin-balanced case. For repulsive inter-action, the relevant largest length scale of the system, theThomas-Fermi radius, will generally increase with inter-action strength leading to a larger scaling factor γ .Finally, we would like to illustrate that accessing theregime of exponential convergence is significantly easierfor the total energy than for other physical quantitiesthat more sensitively probe the short-range correlationsin the wave function. Let us examine the energy partition E = E osc + E int , where E osc = (cid:104) H osc (cid:105) collects the single-particle contributions of kinetic and potential energy [seeEq. (1)] , and the interaction energy E int probes the two-particle density matrix. While both quantities reach ex-ponential convergence for R = l and R = 0 . l in a similarway as the total energy (data not shown), they do notconverge for the narrower Gaussians within the availablerange of M . Figure 3 demonstrates this behavior for R = 0 . l where the total energy seemingly converges anddoes not significantly change for M > on the scaleof the plot, but the single-particle and interaction energyparts still show a strong M -dependence within the acces-sible range of basis set sizes. While E osc and E int serveas useful and sensitive measures of low-order correlationsfor finite-size Gaussian pseudopotentials, please note that -200-150-100-50 0 50 100 150 0 50 100 150 200 250 E n e r g y / − h ω M EE osc E int Figure 3. Convergence of energy components with respect toincreasing size of the single-particle basis set M . While thetotal energy E = E osc + E int converges nicely with increasing M at R = 0 . l and γ = 2 √ (data from Fig. 2), the single-particle part of the energy E osc and the interaction energy E int , which more sensitively probe short-range correlations inthe wave function, have not yet separately converged withinthe accessible range of M values. they lose meaning in the zero range limit R → wherethey separately diverge. Only the sum (i.e., the total en-ergy) remains meaningful, which is a known feature ofzero-range pseudopotentials in two and three dimensions[19]. V. CONCLUSION AND OUTLOOK
In this work we considered the question whethersmooth pseudopotentials are a good choice for represent-ing short-range interactions in numerical approaches thatrely on a Fock-state basis constructed from a finite setof single-particle functions. The combination of smoothpseudopotentials and an asymptotically complete set of(smooth) single-particle basis functions promises expo-nential convergence. This regime of exponential conver-gence can only be reached, however, if the basis set islarge enough to resolve the relevant the length scales ofthe problem.In order to isolate the effects of the single-particle ba-sis we have used an exact diagonalization procedure tocapture the many-particle quantum physics. For an ex-ample system of experimental interest (a two-dimensionalharmonic trapping potential with attractively-interactingfermions) we have derived simple expressions for the re-quired minimum size of the single-particle basis in orderto resolve a given pseudopotential length scale. The al-gebraic scaling of l / can be improved to l / by scalingthe basis set but remains algebraic with the required res-olution length scale. An additional algebraic scaling ofthe required basis set size is found to apply to the par- ticle number. With numerical example calculations forthree particles we could demonstrate that the exponen-tial convergence regime could indeed be reached, albeitin a regime where the pseudopotential length scale is notmuch smaller than the particle separation of length scaleof the trapping potential.In order to faithfully represent the physics of short-range interacting particles, as is relevant for ultracoldquantum gases of neutral atoms, it would be necessary toreduce the length scale of the pseudopotential much fur-ther and extrapolate to the zero range limit. Althoughsuch extrapolation has been demonstrated in few-bodysystems using different computational approaches [58], itdoes not seem feasible for the current approach. Look-ing towards the treatment of larger particle numbers andsimultaneous extrapolation to the zero-range limit, onehas to consider that the benefits of asymptotically expo-nential convergence of the single-particle basis are coun-teracted by the algebraically scaling requirements for theminimum size of the single-particle basis, both with par-ticle number and with the ratio of extremal length scalesthat have to be resolved.It would be interesting to examine the basis set conver-gence in Jacobi coordinates, where one, typically large,length scale can be removed by separating out the center-of-mass motion. The center-of-mass length scale can be-come the largest length scale in the system at strongattractive interaction and weak trapping (e.g., for low-mass bright solitons or droplets), and hence its elimina-tion from the numerical calculations should accelerate theconvergence properties according to the arguments dis-cussed in Sec. IV. In addition, the overall computationalcomplexity is reduced by eliminating the center-of-massdegrees of freedom. Calculations in Jacobi coordinateshave been performed, e.g., for three particles in one di-mension [74, 75]. Extensions to higher dimensions andlarger particle numbers are complicated by the fact thatfermionic or bosonic permutation symmetries manifestthemselves by more complex symmetries in Jacobi coor-dinates that are more difficult to treat [74].A remaining alternative approach is to replace thesmooth pseudopotential by a renormalized contact in-teraction. This has the advantage of removing the ar-tificial length scale of the pseudopotential, while at thesame time the property of exponential convergence is lostand replaced by algebraic convergence [33, 34]. Extrap-olation to the limit of zero range interactions has beensuccessfully demonstrated for up to 66 fermions in threedimensions with an auxiliary-field quantum Monte Carloapproach [30]. Recently developed approaches like thetranscorrelated method for short-range interactions canfurther improve the power-law scaling [34]. Extendingthe applicability of this method to trapped systems is apromising avenue for future work. VI. ACKNOWLEDGEMENT
We thank Tal Levy for initial work on this project,which indicated the inadequacy of conventional basisset expansions and motivated the analysis presented inthis paper. This work was supported by the MarsdenFund of New Zealand (Contract No. MAU1604), fromgovernment funding managed by the Royal Society TeAp¯arangi.
Appendix A: Fock-Darwin orbitals
In the main text of the paper, we discussed the con-vergence properties from an analytic point of view, wherea product basis of one-dimensional basis functions pro-vides an intuitive picture for the analysis. For numericalcalculations it is more advantageous to apply a set of or-bitals that satisfy the symmetries of the system. Thishelps to restrict the problem to a single irreducible rep-resentation of the symmetry operator, which reduces therequired number of the many-body basis functions andthus the requirements for computer memory and CPUtime.We here use simultaneous eigenfunctions of the har-monic oscillator and the angular momentum operatorknown as Fock-Darwin orbitals
L ϕ n(cid:96) = (cid:126) (cid:96) ϕ n(cid:96) , where L is the angular momentum operator, ϕ n(cid:96) is thesingle-particle eigenfunction function, and n and (cid:96) arequantum numbers with non-negative integer and integervalues. The eigenfunction ϕ n(cid:96) can be easily given in polarcoordinates, ϕ n,(cid:96) ( r,ϑ ) = (cid:115) n !˜ l π ( n + | (cid:96) | )! (cid:18) r ˜ l (cid:19) | (cid:96) | × e − ( r/ ˜ l ) e i(cid:96)ϑ L | (cid:96) | n (cid:18) r ˜ l (cid:19) , (A1)where L | (cid:96) | n ( x ) is the associated Laguerre polynomial.In the numerical calculation the finite single-particlebasis set is chosen according to the total quantum number ¯ n = 2 n + (cid:96) , representing a “shell.” All single-particleorbitals are selected where ¯ n is smaller than or equal toa maximal value ¯ n max .The number M of spatial single-particle orbitals canbe expressed with ¯ n max , M = (¯ n max + 2) (¯ n max + 1)2 . (A2)In this paper the largest ¯ n max is 20, which corresponds to231 spatial orbitals. The total number of the many-bodybasis functions for the three particles is equal to around × . Considering only those many-body states withprojected angular momentum of (cid:126) , the computationalspace can be reduced by about an order to . × basis functions. Appendix B: Evaluation of the matrix elements
The matrix elements of the Hamiltonian can be evalu-ated according to the Slater-Condon rules [76–78], whichexpress them as a linear combination of one-particle in-tegrals (cid:104) ϕ n (cid:96) | H osc | ϕ n (cid:96) (cid:105) = (cid:90) d r ϕ ∗ n (cid:96) ( r ) H osc ϕ n (cid:96) ( r ) , and two-particle integrals (cid:104) ϕ n (cid:96) ϕ n (cid:96) | V | ϕ n (cid:96) ϕ n (cid:96) (cid:105) = (cid:90) d r d r ϕ ∗ n (cid:96) ( r ) ϕ ∗ n (cid:96) ( r ) V ( r − r ) ϕ n (cid:96) ( r ) ϕ n (cid:96) ( r ) . The integrals are calculated with the single-particle ba-sis described in Appendix A. The explicit expressions aregiven in the following sections.
1. Evaluation of one-particle integrals
The one-particle integral can be evaluated analyticallyproviding an easily implementable formula (cid:68) ϕ n (cid:96) (cid:12)(cid:12)(cid:12) ˆ H osc (cid:12)(cid:12)(cid:12) ϕ n (cid:96) (cid:69) = δ (cid:96) (cid:96) (cid:126) ω (cid:32) (cid:16) ˜ l/l (cid:17) (cid:16) ˜ l/l (cid:17) (2 n + | (cid:96) | + 1) δ n n + 1 − (˜ l/l ) (˜ l/l ) (cid:112) n ( n + | m | ) δ n ,n +1 + 1 − (˜ l/l ) (˜ l/l ) (cid:112) ( n + 1)( n + | m | + 1) δ n ,n − (cid:33) .
2. Evaluation of two-particle integrals