PHYSICAL REVIEW A 82, 053635 (2010) Integro-differential equation for Bose-Einstein condensates R. M. Adam 1 and S. A. Sofianos 2 1 South African Nuclear Energy Corporation, P.O. Box 582, Pretoria 0001, South Africa 2 Physics Department, University of South Africa, P.O. Box 392, Pretoria 0001, South Africa (Received 3 May 2010; published 30 November 2010) We use the assumption that the potential for the A-boson system can be written as a sum of pairwise acting forces to decompose the wave function into Faddeev components that fulfill a Faddeev type equation. Expanding these components in terms of potential harmonic (PH) polynomials and projecting on the potential basis for a specific pair of particles results in a two-variable integro-differential equations suitable for A-boson bound-state studies. The solution of the equation requires the evaluation of Jacobi polynomials P α, β K (x ) and of the weight function W (z) which give severe numerical problems for very large A. However, using appropriate limits for A →∞ we obtain a variant equation which depends only on the input two-body interaction, and the kernel in the integral part has a simple analytic form. This equation can be readily applied to a variety of bosonic systems such as microclusters of noble gasses. We employ it to obtain results for A ∈ (10–100) 87 Rb atoms interacting via interatomic interactions and confined by an externally applied trapping potential V trap (r ). Our results are in excellent agreement with those previously obtained using the potential harmonic expansion method (PHEM) and the diffusion Monte Carlo (DMC) method. DOI: 10.1103/PhysRevA.82.053635 PACS number(s): 03.75.Hh, 24.10.Cn I. INTRODUCTION Broadly speaking, the A-particle bound-state problem can be solved using two families of approaches. The first is based on the assumption that the potential can be written as a sum of pairwise acting forces resulting in the wave function of the system being written as a sum of amplitudes for the pairs and fulfilling a Faddeev type equation; in the second family one may use correlation functions and employ one of the numerous variational approaches to solve the problem. In both types of approaches three-body correlations can also be included. The Faddeev approach, introduced in the early 1960s [1], can be applied to systems up to A = 4 and it has been extensively used during the last few decades to study bosonic as well as fermionic systems in a rigorous way. Going beyond the A = 4 system, however, is not at present practical within the Faddeev scheme, as the resulting equations are too complicated to solve and therefore one has no option but to consider instead several assumptions and simplification such as clustering and effective interactions. An alternative to the Faddeev scheme is the use of the hyperspherical harmonics expansion method (HHEM) [2,3], which converts the Schr¨ odinger equation into an infinite set of coupled second-order differential equations. This method is, up to a certain extent, variational since one has to truncate the expansion for numerical purposes and thus one limits the Hilbert space of the wave function, leading to an underestimate of the binding energy. It is therefore clear that any alternative method that includes the two-body (and three-body, if needed) correlations into account is welcome. One such method in which the two-body correlations are taken into account exactly is the integro-differential equation approach (IDEA), valid for A-body systems suggested by Fabre de la Ripelle and collaborators [4–6]. It is based on the expansion of the Faddeev amplitudes into potential harmonics (PH) [7,8] and it can be used in a straightforward manner to calculate bound states of bosonic systems. The same procedure to take correlations into account can be used in fermionic systems, in which case spin (and isospin) is taken into account. This, however, results in some modifications stemming from spin-isospin projections (see, for example, Refs. [9,10]). The IDEA method has been successfully applied in few-body calculations [6,9], in realistic fermionic systems [10], in unequal mass particle systems [11–14], as well as in model calculations for the A = 16 system [15]. In all applications, the binding energies obtained are in good agreement with other results in the literature obtained by other methods. When, however, the number of particles increases, the number of degrees of freedom also increases and the numerical complexity becomes intractable and one has no alternative but to seek methods or simplifications of existing ones suitable for handling many-body systems. The typical number of atoms involved in Bose-Einstein condensation (BEC), for example, is 10 3 –10 6 [16], and consequently studies of the BEC phenomenon are naturally based on quantum Monte Carlo type methods, such as the diffusion Monte Carlo (DMC) [17,18], the variational Monte Carlo (VMC) [19], and the practically exact Green’s-function Monte Carlo (GFMC) [20] methods. A different approach to Monte Carlo methods in studies of BEC is the one employed by Das and collaborators [21–23], and it is based directly on the PH expansion of the Faddeev components resulting in a large system of differential equa- tions. The scheme has been used to study the BEC phenomenon for 87 Rb atoms using repulsive interboson interactions. The method requires the evaluation of Jacobi polynomials P α, β K (z), where α = (D − 5)/2, β = 1/2 + ℓ. D is the dimensionality of the A-boson system, D = 3(A − 1), and z is an angular variable. Furthermore, it requires the use of the so-called weight function W (z) ≡ (1 − z) α (1 + z) β . It is clear that the accuracy in calculating the relevant quantities suffers with increasing A as the P α, β K (z) becomes highly oscillatory while the W (z) has a spike similar to a δ function for z ∼−1, which is difficult to control numerically. As a result, the calculations become cumbersome and practically uncontrollable beyond certain A. A relevant discussion on this problem and how this can, up to a certain extent, be addressed is given in Ref. [23]. 1050-2947/2010/82(5)/053635(7) 053635-1 ©2010 The American Physical Society