Effective renormalized multi-body interactions of harmonically confined ultracold neutral bosons

We calculate the renormalized effective 2-, 3-, and 4-body interactions for N neutral ultracold bosons in the ground state of an isotropic harmonic trap, assuming 2-body interactions modeled with the combination of a zero-range and energy-dependent pseudopotential. We work to third-order in the scattering length a defined at zero collision energy, which is necessary to obtain both the leading-order effective 4-body interaction and consistently include finite-range corrections for realistic 2-body interactions. The leading-order, effective 3- and 4-body interaction energies are U3 = -(0.85576...)(a/l)^2 + 2.7921(1)(a/l)^3 + O[(a/l)^4] and U4 = +(2.43317...)(a/l)^3 + O[(a\l)^4], where w and l are the harmonic oscillator frequency and length, respectively, and energies are in units of hbar*w. The one-standard deviation error 0.0001 for the third-order coefficient in U3 is due to numerical uncertainty in estimating a slowly converging sum; the other two coefficients are either analytically or numerically exact. The effective 3- and 4-body interactions can play an important role in the dynamics of tightly confined and strongly correlated systems. We also performed numerical simulations for a finite-range boson-boson potential, and it was comparison to the zero-range predictions which revealed that finite-range effects must be taken into account for a realistic third-order treatment. In particular, we show that the energy-dependent pseudopotential accurately captures, through third order, the finite-range physics, and in combination with the multi-body effective interactions gives excellent agreement with the numerical simulations, validating our theoretical analysis and predictions.


I. INTRODUCTION
Effective multi-body interactions arise when quantum fluctuations dress the intrinsic interactions between particles. They play a central role in quantum field theories and exemplify the significant difference between interactions in classical and quantum theories. For example, even for a quantum field that has only intrinsic two-body interactions at high energies, at low-energy scales, after the high-energy degrees of freedom are coarse-grained away, the field will manifest at some level effective nbody interactions. The ability to trap and control systems of ultracold neutral atoms [1,2] has created new opportunities to study this physics in the laboratory. Effective three-body interactions in the limit of large twobody scattering length have in particular received a great deal of attention, motivated both by the predictions of universal behaviors [3][4][5][6][7][8][9] and the ability to use ultracold atoms to study physics ranging from molecular [10] to nuclear scales [11,12]. Recently, attention has focused on Efimov-like states and universal behaviors for four-body systems, again in the limit of large scattering lengths [13][14][15][16].
Here, we focus on the opposite regime of weakly in- * Electronic address: pjohnson@american.edu teracting neutral bosons with small scattering lengths. Even in this limit, effective higher-body interactions can be important, particularly for tightly confined or strongly correlated particles. This is seen dramatically in [17], where a superfluid of bosonic atoms is quenched by suddenly increasing the depth of an optical lattice. After the quench, which creates a non-equilibrium state of strongly correlated bosons, beating effects due to multiple distinct interaction energies, as expected from effective three-and higher-body interactions [18,19], are seen in the collapse and revival oscillations of the first-order coherence. Effective multi-body interactions should also have played a role in previous collapse and revival experiments [20][21][22], although in those cases inhomogeneities may have masked their signature. More recently, effective threeand four-body interactions have been used to demonstrate atom-number sensitive photon-assisted tunneling in optical lattices [23], and their influence has been seen in precision measurements on Mott-insulator states of ultracold atoms [24]. A number of studies also suggest that elastic multi-body interactions can play an interesting role in generating exotic quantum phases in optical lattices or modifying the superfluid to Mott-insulator phase transition [25][26][27][28][29][30][31][32].
In this paper, we use renormalized quantum field theory [33] to calculate the perturbative ground-state energy for N ultracold neutral bosons in a three dimensional isotropic harmonic potential with angular frequency ω, and extract from it the effective m-body interaction en-ergies U 2 (ω), U 3 (ω), and U 4 (ω) as a function of ω. The key purpose of the present paper is to (i) systematically develop a renormalized quantum field theory approach for ultracold trapped bosons including finite-range effects, (ii) determine the leading-order four-body interaction, and (iii) validate the formalism through comparison with numerical results. To obtain effective four-body interaction energies it is necessary to work through third order in the two-body scattering length. We use renormalized perturbation theory (see [33]), which develops an expansion around physical as opposed to bare coupling parameters, to systematically cancel the multiple divergences that arise at higher-orders in quantum field perturbation theory. (In this paper, the physical coupling parameter is defined in terms of the measured scattering length, or alternatively the measured energy shift, for two interacting ultracold boson in a harmonic trap at a specified trap frequency.) Renormalized perturbation theory, which is more commonly used in high-energy physics, in this context naturally describes how the effective interactions depend on trap frequency. An example of the power of renormalized perturbation theory to capture low-energy physics is that we independently reproduce, through third order, the two-body ground-state energies calculated in [34]. More fundamentally, the analysis in this paper provides an explicit example of renormalization physics and running coupling constants that can be directly probed using trapped ultracold bosonic atoms, and used to test central concepts in effective field theory.
To calculate effective interactions for confined bosons, we first assumed that the two-body interactions could be described in the low-energy, s-wave limit by an energyindependent zero-range δ-function pseudopotential. To test our perturbative predictions, we then numerically calculated N -boson ground-state energies using a finiterange two-body Gaussian model potential. Comparison with the numerical results revealed that finite-range effects must also be taken into account for an accurate description of realistically interacting bosons. In this paper, we show that both the finite range effects and effective interactions are accurately captured by the combination of zero-range and energy-dependent δ-function pseudopotentials. Including the finite-range corrections, we are able to validate our analytic and numerical calculations of all perturbation theory coefficients through third order.
The basic idea in our approach is the following: we "integrate out" excited vibrational states thereby trading a multi-orbital theory with intrinsic two-body interactions for a single-orbital theory with effective multi-body interactions. The latter can provide a simple but powerful alternative description of the low-energy few-body physics. The quantum fluctuations to excited states both dress the two-body interactions and generate effective higher-body interactions. The idea is illustrated in Fig. 1. We showed in [18,19] how this approach can be used to approximately incorporate the influence of higher bands via the simple modification of adding higher-body interactions i = 0 (nlm = 000) i > 0 Multi-orbital model Single-orbital model to the single-band Bose-Hubbard model [35,36].
Beyond applying directly to ultracold neutral bosons in an isotropic harmonic potential, our results can give qualitative insight into the effective interactions for other trapping potentials. They can also be used for rough approximations to the effective two-, three-, and fourbody interactions in anisotropic potentials, and for neutral bosons in optical lattices. In the latter case, however, anharmonicities are important. For example, we estimate an approximately 30% anharmonic correction to the three-body interactions for 87 Rb in typical lattices. The role of anharmonicities for collapse-and-revival dynamics in optical lattice systems has been analyzed further in [37]. Inhomogeneities and the effect of a background harmonic potential on lattice collapse-and-revival dynamics has been studied in [38,39].
Tunneling also has an influence on collapse and revival in optical lattices [40][41][42]. In deep (post-quench) lattices the typical tunneling energy is nearly an order of magnitude smaller than the effective three-body interaction energy, making the latter effect dominant. Tunneling should, however, be of comparable importance to the effective four-body interactions. Approaches applying effective interaction methods to tunneling in lattice or multi-well systems include [43][44][45][46], and related methods for analyzing physics involving interactions, correlations, higher bands, and quantum tunneling in lattice systems include [47][48][49][50][51][52]. Fermionic systems and fermionboson mixtures also yield interesting types of effective interactions that have received increasing attention (e.g., [31,37,[53][54][55][56]), as well as three-body interactions of fermions and polar molecules in lattices [25].
For experiments with 87 Rb at typical lattice densities the recombination rate [5,57] is one or more or-ders of magnitude smaller than the frequencies associated with both the effective three-and four-body energies, and therefore the elastic effective interactions described in the present paper are more important than inelastic multi-body interactions driving loss. Roughly, we expect three-body recombination to scale at fourth order in the scattering length [58], and in the future we would like to understand both elastic and inelastic interactions in a unified framework. The role of effective three-body interactions in thermalizing a homogenous 1D Bose gas has also been studied [59], and it would be interesting to investigate this physics in the context of a 3D optical lattice system.
The remainder of this paper is organized as follows. In Sec. II, we provide an overview of our results. Section III compares the perturbation theory predictions to numerical estimates for finite-range interactions. Sections IV and V describe the details of the renormalized perturbation theory used to obtain the effective multi-body interactions. Section IV defines the renormalized Hamiltonian and derives the first-and second-order corrections, while Sec. V derives the two-, three-, and four-body interaction energies through third order. Section VI summarizes our results and conclusions. Finally, the appendices give derivations of a number of technical results used in the paper.

II. OVERVIEW
We find the effective interactions of N ultracold bosons in the ground state of an isotropic harmonic oscillator with pairwise interactions modeled by a zero-range δfunction pseudopotential where r i is the position vector of the i th boson. We assume there are no intrinsic three-or higher-body interactions. The two-body coupling constant g 2 is related to a t (0), at first order in perturbation theory, by , where m A is the boson mass, a t (0) is the physical s-wave scattering length measured in the limit that the trap frequency and collision energy go to zero, and O([a t (0)] 2 ) are terms of order [a t (0)] 2 and higher. At higher orders, the relationship between g 2 and a t (0) is modified, and in Secs. IV and V we generalize the perturbation theory as an expansion around the physical trap scattering length a t (ω 0 ) defined for a harmonic potential with frequency ω 0 . In this overview, we summarize our results to third order in a t (0), i.e., the special case ω 0 = 0. We obtain the ground-state energy of N bosons as an expansion E = n=0 E (n) , where E (n) is proportional to [a t (0)] n . Throughout this paper energies are expressed in units of the harmonic oscillator energy ω. The zeroth-order (one-body) energy is E (0) (ω) = ε 0 N, where ε 0 = 3/2 is the dimensionless single-particle ground-state energy. The n th -order energies for n > 0 can be expanded as where N m is the binomial coefficient. The sum goes up to the minimum of N and n + 1, and the n th -order contributions to the m-body interaction energies (in units of ω) are where the harmonic oscillator length for an isotropic potential with frequency ω is independently reproduce the results in [34], if the exact solution found there is expanded through third order. The coefficient c 2 , in particular, is nontrivial and provides a strong consistency check that the renormalized perturbation theory captures the two-body low-energy interactions correctly.
The analytic value of the three-body coefficient c 3 was previously found in [18]. The coefficient c  Table I combines both analytic and approximate numerical results, and the uncertainty is due to the slow convergence of one of the numerically determined sums (see App. B 2).
We also obtain the coefficient c 4 , which gives the leading-order contribution to the effective four-body energy. The coefficient c have similar magnitudes, and consequently we need to include the effective three-body corrections when effective four-body effects are important or of interest. At the end of Sec. II, we show that the correction from the third-order terms becomes significant for ultracold atoms in trap potentials with relatively tight confinement. The coefficients c In Sec. III, we compare the predictions for zero-range interactions to numerical calculations for a Gaussian boson-boson interaction potential and find significant effects from its finite-range nature. We show that these are accurately modeled by adding to the zero-range pseudopotential V 2 an energy-dependent (higher-derivative) pseudopotential [11] 4,2 are obtained to very high precision, but slow convergence of the expression giving α  which has been symmetrized to make it Hermitian. The operators ← − ∇ ij and − → ∇ ij are gradients with respect to the relative separation r i − r j , acting to the left and right, respectively. The coupling constant is where r eff is the effective range [60]. To first-order in g 2 , the shift to the N -body ground-state energy is with The superscript (1,2) indicates that the term is first order in r eff and second order in a t (0), and d (1,2) 2 is given in Table I.
The potential V 2 is proportional to r eff [a t (0)] 2 /σ(ω) 3 and we consider in this paper a regime where a t (0) ≈ r eff σ(ω), such that V 2 and therefore U (ω) can be treated as if the contribution is third order in a t (0). This approach is supported by the comparison between the perturbative energies and the energies for the Gaussian potential with spatial widths r 0 0.01σ(ω) in Sec. III. Adding the contribution U (ω) to the two-body interaction energy extends our results to more realistic systems, like ultracold atoms that interact through finiterange van der Waals potentials. Equation (2) organizes the N -body energy in powers of the free-space s-wave scattering length a t (0). Alternatively, combining our results, we can reorganize the energy in terms of m-body contributions as where through third order the two-body interaction energy is the three-body interaction energy is and the four-body interaction energy is The four-body interaction energy U 4 (ω), although comparatively small, can lead to qualitatively important effects, particularly for traps with stronger confinement. For example, for N = 4 87 Rb atoms and a t (0)/σ(ω) = 0.05, corresponding to a 10 4 Hz trap frequency, the four-body energy should generate a distinct approximately 60 Hz beating frequency in collapse-andrevival oscillations, using our harmonic trap results to estimate the energy in an optical lattice potential. These effects should be measurable as long as tunneling and trap inhomogeneities are sufficiently reduced [17].
Using the effective interaction energies in Eqs. (10), (11), and (12), we can construct a single-orbital effective Hamiltonian whereâ (â † ) annihilates (creates) a boson in a renormalized single-particle ground state. The effective Hamiltonian can be used to incorporate some higher-band physics, via effective multi-body interactions, into a single-band Bose-Hubbard model [18]. The effective interaction energies can be tuned by changing either the scattering length a t (0), for example with a Feshbach resonance [2], or the trap frequency ω of the confinement [61,62]. For example, for a fixed a t (0), this tuning follows from rewriting the U m (ω) in terms of the characteristic scattering energy ω s = 2 /m A [a t (0)] 2 . That is, we writeŨ m (ω) = U m (ω)(ω/ω s ) such that Figure 2 shows, for the case of a zero-range potential (i.e. r eff = 0), the two-body energiesŨ         3 (ω) shows the second-order three-body result found previously in [18], due to the c    2 (ω) shows the good agreement with the exact two-body results from [34] for the regularized zero-range potential.
It is interesting to directly compare the relative sizes of the second-and third-order corrections for 87 Rb in a trap. For small magnetic field strengths, the 87 Rb scattering length and effective range are approximately 5.3 nm and 7.9 nm, respectively [2]. For a trap frequency of 10 2 Hz, and thus a t (0)/σ(ω) = 0.005 ("weak" confinement), the third-order two-body terms c [a t (0)/σ(ω)] 2 [r eff (0)/σ(ω)] are 1% and 2% of the second-order two-body contribution c Similarly, the third-order three-and four-body terms c 3 are each about 1.5% of the second-order three-body contribution c For a trap frequency of 10 4 Hz, and thus a t (0)/σ(ω) = 0.05 ("strong" confinement), the third-order two-body terms increase giving approximately 10% and 20% corrections compared to the second-order two-body contribution. Similarly, the third-order three-and fourbody terms increase giving approximately 15% corrections compared to the second-order three-body contribution. (Notice, however, that the third-order effective two-body coefficient and the finite-range coefficient have opposite signs, and hence their contributions partially cancel.) In typical optical lattice collapse-and-revival experiments with 87 Rb the confinement is even stronger  2 (ω) gives values using the exact two-body results forŨ2(ω) from [34]. and the ratio a t (0)/σ(ω) is on the order of 0.05 − 0.10 [17,[20][21][22]. In this regime we expect non-perturbative effects to also become increasingly important.

III. COMPARISON OF PERTURBATIVE ENERGIES WITH ENERGIES FOR FINITE-RANGE INTERACTIONS
This section compares the predictions of the perturbative ground-state energies for a zero-range δ-function interaction potential, summarized in Sec. II and derived in Secs. IV and V, and numerically obtained energies for N -boson systems with finite-range interactions. We show that the leading-order contribution of an energy-dependent pseudopotential accurately captures the finite-range effects, and allows us to also validate the analytic and numerical coefficients found from the zero-range perturbation theory.
We use a finite-range interaction model based on a Gaussian two-body potential V g (r) = V 0 exp[−(r/r 0 ) 2 /2] with depth (or height) V 0 and width r 0 [63,64]. For a given width r 0 , we adjust the depth V 0 such that V g (r) produces the physical free-space s-wave scattering length a t (0) at zero collision energy. We restrict ourselves to  figure) and the effective range r eff /σ(ω) (in the inset) as a function of at(0)/σ(ω) for the Gaussian potential with r0 = 0.005σ(ω) and r0 = 0.01σ(ω), respectively.
depths V 0 for which V g supports no two-body s-wave bound state in free-space. This implies that V 0 is positive for a t (0) > 0 and negative for a t (0) < 0.
An energy-dependent free-space scattering length for two particles with relative energy E rel and relative wave number where δ f (k rel ) is the free-space s-wave phase shift. The effect of a finite-range potential on the free-space scattering of two ultracold bosons can be captured by Taylorexpanding δ f (k rel ) [60,65], giving where r eff is the effective range parameter which describes the lowest-order energy dependence of the phase shift [61,62]. Figure III shows the effective range r eff and the "volume" r eff [a t (0)] 2 for two bosons interacting with the Gaussian potential with two different choices of r 0 /σ(ω) 1. (The volume factor here characterizes the leading-order effective-range correction to s-wave scattering.) We extract r eff by fitting the numerically evaluated − tan(δ f (k rel ))/k rel to the right-hand-side of Eq. (18) for small scattering energies. The effective range is positive for negative a t (0), negative for small positive a t (0), and diverges as a t (0) → 0. Importantly, since a t (0) = 0 implies V 0 = 0 (no scattering potential), the volume r eff [a t (0)] 2 also vanishes when a t (0) = 0, as seen in the main part of Fig. III. The divergent behavior of the effective range is also observed for realistic van der Waals potentials [66] and indeed for any potential that falls off faster than 1/r 5 [65], although for these potentials (unlike the Gaussian) r eff [a t (0)] 2 is finite but non-zero in the limit a t (0) → 0.
We determine the ground-state energy of N = 3 and N = 4 bosons interacting through the Gaussian model potential under external spherically symmetric harmonic confinement using a basis set expansion that expresses the relative N -body wave function in terms of explicitly correlated Gaussians [64] The u k denote expansion coefficients, N b is the number of basis functions, and S symmetrizes the wave function under the exchange of any pair of bosons. The ij , chosen stochastically from the interval [r 0 /5, 4σ(ω)], are optimized semistochastically following the scheme outlined in Ref. [64]. In brief, the variational method works as follows. Assume we have a basis set consisting of j −1 basis functions that yields a ground-state energy estimate E j−1 . To add the j th basis function (j ≤ N b ), we generate a few thousand trial functions. For each trial function, we solve for a trial ground-state energy by diagonalizing a j × j dimensional generalized eigenvalue problem. (It is a generalized eigenvalue problem because the basis functions are nonorthogonal.) We choose as the j th basis function the one which makes E j smallest, and repeat this process for the (j + 1) th basis function until j = N b . A key benefit of the explicitly correlated basis functions is that the Hamiltonian and overlap matrix elements have compact analytical expressions [64].
Convergence is analyzed by investigating the dependence of the energies on N b and by performing calculations for different sets of widths v (k) ij . To meaningfully compare numerical three-and four-body energies E FR for the finite-range (FR) interaction potential with perturbative results up to order [a t (0)] 3 , the numerical accuracy of the finite-range energies should be notably better than |a t (0)/σ(ω)| 3 . For example, for |a t (0)| = 0.001σ(ω) and |a t (0)| = 0.01σ(ω), this implies numerical accuracy better than 10 −9 and 10 −6 , respectively. An analysis of the basis set error shows that our N -body energies are sufficiently accurate to test the perturbative predictions up to order [a t (0)] 3 for |a t (0)| 0.1r 0 , using about 100 and 500 basis functions for N = 3 and N = 4, respectively. Our numerical accuracy is insufficient to test the perturbative predictions for smaller |a t (0)|. Figure 4 shows the quantity [E FR − E (0) − E (1) ] × 10 4 versus a t (0)/σ(ω), with the finite-range energies E FR numerically computed using r 0 = 0.01σ(ω) (the blue squares) and r 0 = 0.005σ(ω) (the red circles). We have subtracted the energies E (0) and E (1) obtained from the perturbative theory to better examine the physics beyond first order in a t (0). The solid line is [E (2) + E (3) ] × 10 4 from the perturbative theory with r eff = 0. Panels (a) and (b) give the energies for N = 3 and 4 bosons, respectively. For N = 3, we see that finite-range corrections to the zero-range theory become more significant for increasing r 0 .
In Figs. 5(a) and (b), we multiply the N = 3 and 4 en- The nonperturbative numerical results are for potentials with r 0 = 0.005σ(ω) and r 0 = 0.01σ(ω). The figures show that the scaled numerical results are singular near zero scattering length, and only approach the zerorange perturbative results with increasing |a t (0)|. Moreover, by decreasing r 0 the difference between the perturbative results and the scaled finite-range energies is reduced, and we conclude that the divergences at a t (0) = 0 are due to the finite range of the Gaussian potential. Multiplying the energies by [σ(ω)/a t (0)] 2 has magnified the finite-range corrections, showing that an effective field theory description for finite-range potentials requires corrections to the zero-range δ-function potential.
We can calculate the leading-order influence of a finiterange potential by including the energy-dependent zerorange pseudopotential of Eq. (5). For the N -boson ground state, the pseudopotential gives to first order in g 2 an energy shift E (1,2) [see Eq. (7)]. At this order, the addition of V 2 is equivalent to replacing a t (0) by a f (E rel ), with Eq. (18) evaluated at the relative zero-point energy E rel = 3/2 (in units of ω) of two non-interacting bosons in the trap. The solid lines in Fig. 5 show 2 as a function of a t (0)/σ(ω) for N = 3 and N = 4 trapped bosons, respectively. Combining the perturbative predictions for zero-range contributions E (2) + E (3) and the effective-range correction E (1,2) gives excellent agreement with the nonperturbative finite-range energies. The comparison validates the perturbation theory and predictions derived in this paper for effective interactions including finite-range correc-tions, through third order in a t (0). It also shows that the divergences in Fig. 5 at a t (0) = 0 are due to the divergence of the effective range shown in the inset of Fig. III. Finally, we note that the energy shift is proportional to the volume r eff [a t (0)] 2 and goes to zero at a t (0) = 0, as expected.

A. Hamiltonian and renormalization condition
The numerical results in Sec. III show that finite-range effects are important at third order in perturbation theory for realistic bosons. We incorporate these corrections by modeling the pairwise collisions of ultracold bosons by combining the zero-range pseudopotential where a bare is now identified as the bare scattering length, and the effective-range potential which has the bare coupling constant The interactions of N ultracold neutral bosons can be described in quantum field theory with the Hamiltonian H = H 0 +H I , where H 0 is the single-particle Hamiltonian and The field operatorsψ(r) andψ † (r) respectively annihilate and create a boson at position r. We assume the absence of intrinsic three-or higher-body interactions. The bosonic field is expanded over isotropic harmonic oscillator states with frequency ω aŝ withâ i annihilating a boson in orbital φ i (r). In the following we use the shorthand notation i = {nlm}, denoting the (dimensionless) single-particle energies as ε i = ε nlm = (2n + l + 3/2), where n, l = 0, 1, 2, ..., and i = 0 is the {nlm} = {000} single-particle vibrational ground state. Substituting Eq. (24) into H and dividing by ω, we define the dimensionless Hamiltonian and The matrix elements and are normalized such that K 00;00 = 2/π and K 00;00 = (3/4) 2/π, with the semi-colon separating initial and final states and The factors of [σ(ω)] 3 and [σ(ω)] 5 make the matrix elements dimensionless and ω-independent. As explained in Sec. II, we assume a regime where H I can be treated as third order in perturbation theory.
The noninteracting ground state containing N bosons in the i = 0 (i.e., nlm = 000) vibrational ground state is using N | a † 0 a † 0 a 0 a 0 |N = N (N − 1) and recalling that |N denotes N bosons in the non-interacting vibrational state nlm = 000. The two-body, first-order coefficient is α (1) 2 = 2/π. At higher orders in H I , there are divergences due to the δ-function potential (see e.g. [67,68]). We regulate these by either truncating sums over intermediate states at a high-energy cutoff ω c , or by using an exponential regulator function. The former is more convenient for numerical approximations, while the latter is more convenient for analytic results. In either case, we find at second order that U

diverges as
√ ω c and renormalization is required. Although this can be done using bare perturbation theory, in which infinities are absorbed by appropriately redefining bare parameters, we use the method of renormalized perturbation theory which provides a systematic and self-consistent approach for calculations beyond second order involving multiple divergent terms. Renormalized perturbation theory (e.g., see [33]) reexpresses the bare scattering length as a bare = a t (ω 0 ) + a ct (ω 0 ). We assume a separation of length scales lc at(ω0) σ(ω0), or equivalently a separation of energy scales ωc ωs ω0, where ωs = /mAat(ω0) 2 . The harmonic oscillator lengths σ(ω) and σ(ω0) are of the same order, although σ(ω0) is not necessarily larger than σ(ω). Similarly, we assume that at(ω0), at(0), r eff , and the Gaussian width r0 are of the same order. The order of length scales within a group is arbitrary.
A renormalization condition defines a t (ω 0 ) as the physical scattering length for two bosons in a trap at frequency ω 0 . The cutoff dependent remainder a ct (ω 0 ) is called a counterterm. For brevity, this notation suppresses the dependence of a ct (ω 0 ) on ω c . In the following, we call a t (ω 0 ) the "trap scattering length" at frequency ω 0 , to distinguish it from the energy-dependent free-space scattering length a f (E rel ) defined in Eq. (17). In the limits of zero relative collision energy and ω 0 = 0, the trap and free-space scattering lengths are equal, i.e., a t (0) = a f (0). With the combination of V 2 and V 2 , the trap scattering length a t (ω 0 ) includes both the effects of the ω 0dependent dressing by quantum fluctuations to higher orbitals and finite-range effects. Note that the trap scattering length a t (ω 0 ) does not, in general, equal the freespace scattering length a f (E rel ) defined in Eq. (18) because the latter does not correctly capture the influence of the harmonic confinement on the quantum fluctuations to higher orbitals.
Together with the renormalization condition, the other key ingredient in renormalized perturbation theory is that the leading-order scattering length counterterm a ct (ω 0 ) is proportional to [a t (ω 0 )] 2 ; in other words, it is a second-and higher-order contribution. This, plus the renormalization condition, systematically reorganizes the perturbation theory, order-by-order, so that it is an expansion in the physical value a t (ω 0 ) instead of a bare . Figure 6 summarizes the relationship between the characteristic length and energy scales for our model system of trapped ultracold bosons.
Substituting Eq. (30) into Eqs. (25) and (26) gives (31) where the zero-range and counterterm operators are and the effective-range operator is The renormalized perturbation theory is then organized based on the observation that The single counterterm operator V ct (ω; ω 0 ) cancels all divergences from the operator V (ω; ω 0 ), at all orders in perturbation theory. In contrast, the effectiverange operator V (ω; ω 0 ) leads to a nonrenormalizable field theory with the consequence that new counterterm operators are required at every order in perturbation theory beyond first order in g 2 ; because we are only working to first order in g 2 in this paper, no additional counterterms are needed. Note that the frequency ω 0 at which a t (ω 0 ) is defined and the trap frequency ω for which we want to compute energies are independent. In the overview, we summarized our results for the special case where ω 0 = 0. The general case of arbitrary ω 0 facilitates renormalization of the perturbation theory. More importantly, the renormalized perturbation theory is "calibrated" to a measured value of a t (ω 0 ) at a desired trap frequency ω 0 , and is then used to predict energies for trap frequencies ω not generally equal to ω 0 .
We can now compute the ground-state energy We have used the semi-colon notation in Eqs. (31), (32), (33), (34), and (35) to distinguish between the roles of the frequencies ω and ω 0 . Before renormalization, the interaction energies U m (ω; ω 0 ), found from perturbation theory in H I (ω; ω 0 ), are functions of a t (ω 0 ) and a ct (ω 0 ). The renormalization condition can be expressed as which, in practice, is solved for a ct (ω 0 ) to the desired order in perturbation theory. Another way of describing the renormalization condition is that a ct (ω 0 ) is tuned such that the first-order result is exact and the secondand higher-order corrections to the two-body energy vanish when evaluated for two bosons in a trap with ω = ω 0 . After renormalization, the interaction energies U m (ω; ω 0 ) depend only on a t (ω 0 ) and, moreover, the ω-dependence of the ground-state energy satisfies E(ω; ω 0 ) = E(ω; ω 0 ), for any pair of frequencies ω 0 and ω 0 .

B. Energy at first-order in scattering length
We use renormalized Rayleigh-Schrödinger (RS) perturbation theory to compute the N -boson ground-state We separate the contributions at each order into m-body energies, such that U m (ω; ω 0 ) = n U (n) m (ω; ω 0 ). The zeroth-order term is E (0) (ω) = ε 0 N . The first-order energy shift is . Comparing to Eq. (35), we see that the two-body energy to first-order for any ω and ω 0 is with c (1) For a trap with ω = ω 0 , the renormalization condition says that U 2 (ω 0 ; ω 0 ) = 2/π[a t (ω 0 )/σ(ω 0 )] is the exact two-body energy. For ω = ω 0 , U 2 (ω = ω 0 ; ω 0 ) is the leading order contribution to the full two-body energy U 2 (ω; ω 0 ), but, as shown in the following sections, there are higher-order corrections that become increasingly important the more ω differs from ω 0 .

C. Energy at second-order in scattering length
The second-order energy shift is given by where V ij,kl = ij|V |kl and V ct;ij,kl = ij|V ct |kl . The notation |ij = Z ijâ † iâ † jâ 0â0 |N denotes the state with either one or two particles excited from the non-interacting ground state, ∆ε ij = ε i + ε j − 2ε 0 , Z ij is a normalization factor, and ij = 00 denotes summing over all i, j except i = j = 0. Equation (40) is modified from the usual RS perturbation theory because of the presence of the O([a t (ω 0 )] 2 ) interaction term V ct , which generates the counterterm contribution.
Using Eqs. (32) and (33), we have and the expectation value is with respect to the noninteracting ground state |N − 2 ∝â 0â0 |N . The notation ij = 00, kl indicates that the sum is over all i, j, k, l except i = j = 0. Wick's theorem gives  42) can be understood as leading to an effective three-body interaction, whereas the second term is a correction to the two-body interaction.
The second-order interaction energies U 2 (ω; ω 0 ) and U and The expressions for β 2 (ω) and α 3 are defined in Table II, which also shows the explicit values calculated in Appendices A and B for an isotropic harmonic trap. We use the notation that α 2 (ω) diverges with cutoff as ω c /ω, where ω c /ω is approximately the number of harmonic oscillator levels included in the sum as a function of ω, for a fixed cutoff ω c . The coefficient α (2) 3 is convergent and in the limit ω c /ω → ∞ is independent of ω. In the following, we only indicate the explicit ω dependence for coefficients that remain sensitive to ω in the limit ω c /ω → ∞, e.g., we write β (2) 2 (ω) but α (2) 3 . We use a hard cutoff to numerically evaluate the coefficients α 5 , we find analytic results in the limit ω c /ω → ∞. Finally, using the exponential regulator, we obtain analytic results for the coefficients β   . Intermediate states have one or more excited particles and contribute an energy denominator 1/∆ε ij . For example, the diagram is a graphical representation for the term α . We obtain combinatorial prefactors [e.g. −6 for U 3 (ω, ω 0 )] from Wick's theorem by counting the number of equivalent contractions, dividing by 2 for every factor of a t (ω 0 ) or a ct (ω 0 ), and multiplying by m! for an m-body term.
The renormalization condition through second order is U 2 (ω; ω 0 ) = U Solving for the counterterm gives Substituting into Eq. (43) gives where the function can be used for any ω. (We have used σ(ω 0 )/σ(ω) = ω/ω 0 above to simplify the expressions.) The form of the expression for the coefficient c The renormalization condition is automatically satisfied since c Combining the first-and second-order contributions for the two-body interaction energy gives The coefficient α 3 in Eq. (44) is finite and does not require a regulator. For the three-body interaction energy we obtain where c 3 = −0.85576... This value was previously obtained in [18], and is also calculated in App. B 1.

V. EFFECTIVE INTERACTIONS THROUGH THIRD ORDER
We now extend our analysis to third order in the scattering length a t (ω 0 ). This is necessary to obtain the leading-order effective four-body interaction. Including the counterterm and effective-range interaction, the formula for the third-order energy shift is − V 00,00 V ct;00,ij V ij,00 ∆ε ij + V 00,00 .
The first term on the right-hand-side of Eq. (54) gives Noninteracting ground state (3 bosons) Perturbation theory diagram for this process.
FIG. 7: Sequence of boson-boson interaction induced transitions to higher orbitals. This example generates corrections to the ground state energy that can be viewed as an effective three-body interaction. The process, which involves three interaction vertices, arises at third order in perturbation theory and gives the energy shift α where the expectation value is with respect to the noninteracting ground state with (N − 2) bosons. Applying Wick's theorem, Eq. (55) expands as We find effective two-, three-, and four-body interactions from the terms with four, three, and two contractions, respectively. The zero-contraction term vanishes since :â iâjâ † kâ † lâ i â j â † k â † l : = 0. Comparing to Eq. (35), we see that there is a two-body contribution β 3 , β 3 , etc., are given in Table II, along with the associated diagrams, asymptotic behavior, and explicit forms for an isotropic harmonic oscillator potential (calculated in Appendices A and B). Figure 7 illustrates one of the sequences of virtual transitions giving rise to α

.
We next use Wick's theorem to evaluate the second term on the right-hand-side of Eq. (54), finding where α contribute to several effective multi-body energies. From Table II 5 [a t (ω 0 )/σ(ω)] 3 = look like processes requiring four and five distinct particles, respectively. They also appear to be composed of "disconnected" sub-diagrams. In RS perturbation theory, however, the second term in Eq. (54) can be reinterpreted in terms of particles going "backward" in time (right to left), or alternatively an interpretation can be given in terms of holes. For example, the term −α . Similarly, if only one particle goes back in time, it can collide with a third particle, giving the threebody contribution −6α , which is also connected. In this paper, these related two-, three-, and four-body diagrams have the same numerical value: = = . The third term on the right-hand-side of Eq. (54) gives the two-and three-body counterterm contributions and These counterterm contributions, shown in Table II, cancel the divergences from , , and . The disconnected 5-body contribution from Eq. (57), generated by the term with a single contraction, cancels with a + e −ωc/ω  The coefficients for all interaction processes contributing to the two-, three-, and four-body interaction energies through third-order in perturbation theory in ξt = at(ω0)/σ(ω). The first column shows the diagrams from which the m-body, n th -order coefficients α (n) m and β (n) m (ω) can be reconstructed. The coefficients as multidimensional sums are given in the second column. Sums are over all indices i, j, k, ... except combinations that give a zero energy term in the denominator. The third column gives the asymptotic behavior of the coefficients in terms of the cutoff ωc and a constant a. The last column gives the explicit values for the coefficients for an isotropic harmonic oscillator potential. These values are obtained in the Appendices. The table also shows the counterterm processes, the leading-order effective-range contribution, and other relations needed for the renormalized perturbation theory.
the disconnected five-body term in Eq. (56), and there is no effective five-body interaction at third order. Finally, the last term on the right-hand-side of Eq. (54) gives the effective-range contribution from which we extract the effective-range two-body interaction energy The coefficient α is given in Table II. The special case ω 0 = 0 gives Eq. (7) and Eq. (8).

A. Two-body interaction energy
Adding all two-body contributions through third order, we obtain Note that all diagrams in Eq. (64) are connected, when interpreted in terms of both forward-and backward propagating particles. All coefficients are given in Table II. For brevity, in Eq. (64) and the following we adopt The counterterm, found in the previous section to second order, must now be recalculated using the renormalization condition through third order. This adds a third-order term which cancels the divergence from β (ω; ω 0 ). Solving the renormalization condition U 2 (ω; ω 0 ) = U Diagrammatically, this can be expressed as with all diagrams evaluated at ω = ω 0 . By including U (1,2) 2 (ω; ω 0 ) in the counterterm equation, the renormalization condition for a t (ω 0 ) includes both zero-range and effective-range contributions. If we do not include U (1,2) 2 (ω; ω 0 ) in the renormalization condition, then a t (ω 0 ) is the trap scattering length for zero-range potentials. Substituting the counterterm from Eq. (65) into Eq. (64) and using α 2 (ω)] 2 , which is proven in Appendix B 7, we find after some algebra that where c (ω, ω 0 ) are given in Table III. Recall that in the formula U 2 (ω; ω 0 ) the first argument ω is the trap frequency for which we are interested in predicting the two-body energy, and the second argument ω 0 is the trap frequency at which the two-body trap scattering length a t (ω 0 ) is defined or measured. The coefficients for ω 0 = 0 are given in Table I. If r eff = 0, these values reproduce through third order the exact solution for the ground state of two harmonically trapped bosons with zero-range interactions found in [34]. This agreement between the quantum mechanical and quantum field theory solutions is a nice illustration of how the renormalized effective field theory captures the correct low-energy physics. Interestingly, if r eff = 0, our result for U 2 (ω, 0) still agrees with the solution in [34], if that solution is Taylor expanded in a t (0) and r eff [a t (0)] 2 after making the substitution a f (0) → a f (0) + (1/2)r eff [a f (0)] 2 k 2 rel , providing further evidence of the universality of the higherorder perturbative results derived here.
Another important special case is ω = ω 0 . Since c (ω 0 , ω 0 ) = 0, the predicted two-body energy is reproducing the renormalization condition that a t (ω 0 ) is the physical trap scattering length for two bosons at frequency ω 0 .

Frequency-Dependent Effective Interaction Coefficients
Two-body c

]
Four-body c   (ω, ω0), which determines the leading-order effective-range correction, for neutral bosons in a harmonic potential of frequency ω, in terms of the scattering length at(ω0) defined at trap frequency ω0. The special case ω0 = 0 reduces to the results given in Table 1.

B. Three-body interaction energy
In [18], we obtained the effective three-body interaction energy to second order. We now determine the next-order correction by combining all three-body con-tributions through third order, giving Representing the three-body contributions from α 4,3 and α (3) 5 using reversed (left to right) particle lines, as previously described, we again find that only connected diagrams contribute.
For the three-body energy, it is sufficient to use the second-order counterterm in Eq. (46). In Appendix B 8, we show that β , is finite. After some algebra, we find that where c 3 and c 3 (ω, ω 0 ) are given in Table III. If ω 0 equals zero, we find The error reflects a one standard deviation uncertainty due to the extrapolation of the numerical estimate for α 3 to the limit ω c → ∞ (see App. B 2). Another special case is ω = ω 0 , giving (72)

C. Four-body interaction energy
Finally, we calculate the leading order contribution to the effective four-body interaction energy. We find with coefficient 4,1 + 48α is independent of ω and ω 0 . The leading-order contribution to the four-body energy does not require renormalization, as is true for all leading-order m-body terms. Comparison of c (3) 4 and c (3) 3 reveals, however, that they are of similar magnitude and therefore for a consistent and accurate treatment both corrections need to be included. Because c 3 requires renormalization, we see why the systematic renormalization of divergences is needed even though the leading-order contribution to the four-body interaction energy could be obtained without these considerations.

VI. SUMMARY
We have derived effective two-, three-, and four-body interaction energies for N bosons in an isotropic harmonic trap of frequency ω. These energies are functions of the trap scattering length a t (ω 0 ) and harmonic oscillator length σ(ω), and include both renormalization effects due to quantum fluctuations to higher-orbitals and leading-order finite-range corrections. The frequency ω 0 at which the scattering length is defined plays a role closely analogous to the low-energy scale at which coupling constants are defined in high-energy effective field theories (e.g., see [33]). The formulas for the interaction energies are given in Eqs. (67), (69), and(73), and are expressed in terms of the functions c (n) m (ω; ω 0 ) given in Table III. In turn, these functions require the coefficients α  Table II. The special case when ω 0 = 0 is summarized in Table I and Eqs. (10), (11), and (12). In Sec. III, we showed that these results give excellent agreement to numerical simulations for ultracold bosons interacting through a Gaussian model potential.
We find at third-order in a t (ω 0 ) that the shifts to the effective three-and four-body interaction energies are comparable, showing that the renormalized threebody interaction needs to be taken into account when the leading-order four-body interactions are considered. In the future, we plan to use this formalism to determine the effective multi-body interactions for other potentials, such as anisotropic traps or the anharmonic sites of an optical lattice. The cross-over from the perturbative small scattering length regime to the universal regime of Efimov physics is also very interesting and diagrammatic resummation techniques can be used to study the onset of nonperturbative behaviors. For example, collapse and revival experiments suggest that four-and higherbody interactions may be present in the data [17], but our results also show that in these systems a t (0)/σ(ω 0 ) is large enough for significant nonperturbative effects to be important, and we would like to better understand this physics within our framework. A unified description of elastic and inelastic interactions (e.g. three-body recombination physics [5,57]) would also be useful.
More immediately, the results in this paper can be ap-plied to investigations of finite-range interactions, can be used for precision experiments probing for the possible existence of intrinsic three-and higher-body interactions, and can enable explorations of fundamental concepts in effective field theory including renormalization and energy-dependent (running) coupling constants. For example, the influence of intrinsic higher-body interactions would cause deviations from our predictions, which are based on only intrinsic two-body interactions. Moreover, we can engineer and exploit useful effective interactions using a combination of magnetic Feshbach resonances [2] and the ability to tune the few-body interactions by controlling the trap parameters and shape [53,61]. One of our longer-term goals is to use this physics to develop nonlinear measurement techniques. For example, the nonlinear dynamics seen in collapse-and-revival experiments can lead to better than shot-noise measurements of the m-body interaction energies, or it may be possible to exploit strongly correlated non-equilibrium states in lattices for new types of sensing. In this way, the rich physics of renormalization and nonlinear quantum dynamics could be used to create new types ultra-cold atom simulators, quantum information processors, or quantum sensors.

VII. ACKNOWLEDGEMENTS
PRJ and ET acknowledge support from the U.S. Army Research Office under contract/grant 60661PH. PRJ acknowledges additional support from the Research Corporation for Science Advancement and computing resources provided by the American University High Performance Computing System. ET acknowledges support from a NSF Physical Frontier Center. DB and XYY gratefully acknowledge fruitful discussions with Kevin Daily and support by the NSF through grant PHY-0855332. PRJ thanks Nathan Harshman for helpful discussions.
Appendix A: δ-function boson-boson interaction matrix elements for an isotropic harmonic trap This appendix derives interaction matrix elements for bosons in an isotropic harmonic oscillator trap with frequency ω and zero-range δ-function interactions. Alternative methods for obtaining these matrix elements are given in [69,70].
This expression also gives the matrix element for the transition |0i → |0k .

Matrix elements in relative and center-of-mass particle basis
It is simpler to compute some matrix elements by switching to a basis of states |ĩj = |nlm, N LM , with normalized relative and center-of-mass wavefunctionsφ nlm (r)Φ N LM (R) defined in terms of coordinates r = (r 1 − r 2 )/ √ 2 and R = (r 1 + r 2 )/ √ 2, and (dimensionless) two-particle energy differences ∆ε nlm,N LM ≡ 2n + l + 2N + L. (A11) Working in the |nlm, N LM basis and using the fact that the interactions conserve the center-ofmass motion, the matrix elements for the transitions |kl → |ĩj are Kĩj ;kl = K nlm,n l m ;N LM,N L M = K rel (n, n )δ l,0 δ m,0 δ l ,0 δ m ,0 δ N,N δ L,L δ M,M , where only depends on the principle quantum numbers for the relative motion. Below we use the fact that K rel (n, n ) factors as and K rel (n, 0) = 2 π 3/4 Γ(n + 3/2) Γ(n + 1) .
1. Three-body, second-order coefficient α In the single-particle basis |n 1 l 1 m 1 , n 2 l 2 m 2 , the contribution has the coefficient where the sum i =0 is over all allowed single-particle states excluding the ground state. Due to angular momentum conservation only i = {nlm} with l = m = 0 contribute, and ∆ε i0 = 2n. Evaluating the sum gives the analytic result The terms in the summand become smaller exponentially with n, and α 3 converges to the asymptotic form a + O(e −ωc/ω ). This behavior will be true for all "tree diagrams" which, like , have no closed loops.
Inserting exponential regulators for each energy denominator in Eq. (B18), we obtain The factorization of β 3 in the finite part α 2 (ω) is important for the renormalization of the three-body interaction at third order.