{'fraction_non_alphanumeric': 0.046806685255151706, 'fraction_numerical': 0.015211760241540315, 'mean_word_length': 4.150577838205303, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum theory of bright matter wave solitons in harmonic confinement

This paper investigates bright quantum-matter-wave solitons beyond the Gross-Pitaevskii equation (GPE). As proposals for interferometry and creating nonlocal quantum superpositions have been formed, it has become necessary to investigate effects not present in mean-field models. We investigate the effect of harmonic confinement on the internal degrees of freedom, as the ratio of zero-point harmonic oscillator length to classical soliton length, for different numbers of atoms. We derive a first-order energy correction for the addition of a harmonic potential to the many-body wave function and use this to create a variational technique based on energy minimization of this wave function for an arbitrary number of atoms, and include numerics based on diagonalization of the Hamiltonian in a basis of harmonic oscillator Fock states. Finally we compare agreement between a Hartree product ground state and the Bethe ansatz solution with a Gaussian envelope localizing the center of mass and show a region of good agreement.


I. INTRODUCTION
Since the experimental realization of Bose-Einstein condensation (BEC) with dilute atomic gases [1,2], much progress has been made in the degree of control possible in terms of external potentials and control over interactions via magnetic and optical Feshbach resonances. Recently, it has become possible to create condensates of atomic species with scattering lengths that can be tuned to be negative [3][4][5], for example hyperfine levels of 85 Rb, 7 Li and 133 Cs [6,7]. Chromium is also shown to have negatively tunable scattering lengths, but also has significant long-range dipole forces [8].
Such attractive condensates have the remarkable property of 'self trapping', i.e. being localized (at least in terms of pair correlations) on a length scale shorter than would be expected for a non-interacting condensate. In the limit in which there is one dimension where the non-interacting ground state would be infinitely wide, the Gross Pitaevskii equation (GPE) predicts the ground state would still be localized to a finite size. Therefore with the GPE, these condensates behave as solitary matter waves, or in the quasi 1D case, as classical solitons. The transitional region between 1D and 3D has been investigated using variational methods [9]. Experiments measured systems with a lifetime of several seconds, both for the case of a single wavepacket [3] (a ground state) and multiple smaller wavepackets [4,10], referred to as soliton trains. In addition to experimental results, theoretical results relating to the interaction of such systems with potential barriers predicted effects such as enhanced reflection and transmission [11,12].
In this situation, the GPE predicts an infinite number of conserved quantities within the system and thus complete integrability [13]. As a result the inverse scatting transform can be used to obtain solutions [14] that are a combination of bright solitons, which do not change shape as the quantum pressure is exactly balanced by the nonlinear interaction, and radiation, which does. In fact any initial condition for the GPE * d.i.h.holdaway@dur.ac.uk equation can be broken up into these components [14]. In the case of systems of a multiple soliton system, individual solitons can collide with other solitons without a transfer of energy between them, resulting in only an asymptotic position and phase shift. Such localized matter waves (which are typically of the order of a few micrometers in width) could also theoretically be split coherently into multiple parts [15] and prove useful for interferometry [16,17], studies of quantum reflection [12,18] or probing surface potentials [19].
Under certain circumstances such solitons behave like classical particles with finite range interactions [20] even in the presence of harmonic confinement. However, in reality, such objects should behave as quantum particles (i.e. with no substructure, but with the location determined by a wave-function obeying quantum mechanical laws rather than a specific position). Despite the GPEs success in describing many phenomena in BEC, even for very small numbers of atoms [21], this quantum mechanical center-of-mass behavior is totally lost under the approximation of a product state wave-function. Therefore we consider a full many-body quantum description, making use of the usual pseudo-potential approximation.
The dynamics of the center of mass of an interacting gas in a harmonic potential are independent of the interactions, giving rise to the so called "Kohn mode" [22]. More generally, any potential which looks locally harmonic on the length scales dictated by the internal degrees of freedom (in this case, the classical soliton length) can be considered to only weakly couple the center-of-mass to other degrees of freedom. As a result, this behavior needs to be considered separately. In this weak coupling approximation the center-of-mass behaves like a non-interacting particle of mass N M, which can be localized on scales far wider than a classical soliton length or even not at all. This delocalisation is not present in the mean field model, however the Kohn mode is still present for translational motion of the center of mass.
Exact results exist for Bose gases in free space with periodic boundary conditions for both repulsive [23] and attractive [24,25] interactions; these have the advantage of explicitly separating the center-of-mass component of the wave-function and being accurate even for very small numbers of atoms and fragmented states, for which the GPE is not. It has however been shown that in certain situations, in the limit N → ∞, g → 0 with Ng = constant, the GPE functional is an exact description of the system [26]. Modern numerical approaches are available, such as the multi-configurationaltime-dependent Hartree method, which have been used to study bright matter wave solitons [27]. The separation properties of the center-of-mass in free space and harmonic traps make available the possibilities of creating non local superpositions [27,28] (cf. [29]) with sufficient numbers of atoms to make them detectable, in ways not predicted in the classical field description of the GPE.
In the following we study the Lieb-Liniger(-McGuire) gas [23,24], a 1D system of identical Bosons with attractive contact interactions, with the addition of a harmonic trapping potential. We have present a series of analytic and numerical techniques to study many-body effects, making use of the separability of the many-body wave-function into a center-ofmass component and a relative component. This paper focuses on the ground state of the system and energy corrections to the relative components of the wave-function as trapping is increased, along with estimates of the overlap with the freespace relative ground state.
We consider a unit system with the Hartree soliton length set to unity and √ γ a dimensionless ratio between this length and the harmonic oscillator length in the axial direction, such that γ → 0 recovers the free case solved by the Bethe Ansatz. We use a numerical method based on exact diagonalization of the Hamiltonian over a basis set of Hermite functions, truncated up to a maximum energy, in order to determine the many body ground state energies, combined with a variational method for low γ.
This paper is organized as follows: Section II introduces the exact results in one dimensional infinite systems (using the Lieb-Liniger model [23]) and the unit rescaling used to keep the mean field soliton length constant throughout the paper. Also included is the separability of the many-body Hamiltonian and the existence of the Kohn mode as well as the exact eigenstates for two interacting bosons in a harmonic potential. Section III derives a perturbative energy correction to the relative ground state energy from the introduction of a harmonic trapping potential, along with a variational procedure to estimate the ground state in the limit of weak trapping. Section IV introduces the numerical method used to perform calculations in the many-body system for varying 1D harmonic trapping potential, using a basis set of harmonic oscillator eigenstates, which are projected to a center-of-mass excitation basis. Section V numerically investigates changes to relative component (i.e. having excluded the center of mass) of the ground state as the trapping potential is increased. These calculations are performed for different numbers of atoms and compared with predictions based on the GPE. Section V C examines quantitatively the overlap between the mean field approximation of the ground state and the many-body solutions, along with finding a regime of agreement where the difference between the two models is small in every respect. Section VI summarizes and comments on the results.

II. PRELIMINARIES
A. System overview and rescaling

The many-body Hamiltonian in first and second quantization
We consider a system of identical bosonic atoms within a cylindrically symmetric prolate (the radial frequency ω r is greater than the axial frequency ω x ) harmonic trapping potential V(x, y, z) = M[ω 2 x x 2 +ω 2 r (y 2 +z 2 )]/2, where M is the atomic mass. We further consider the system to be sufficiently lowtemperature for the atomic interactions to be pure s-wave and contact-like, and assume that we are in an appropriate parameter regime so that the radial modes can be considered "frozen out" for low-energy states, taking the Gaussian form of radial harmonic oscillator ground states [30][31][32]. Finally, we assume the interactions to be attractive.
Integrating out the radial degrees of freedom, we obtain the 1D Hamiltonian, in second-quantized form, and hence for an arbitrary number of atoms, as follows: where g 1D = 2 ω r |a s |, and a s is the (assumed negative) s-wave scattering length [33]. The coordinate-space representation for the corresponding first-quantized form of this Hamiltonian for N atoms is then (2) where we use x as a shorthand for {x 1 , x 2 , x 3 , . . . , x N }, the N individual particle coordinates. In the absence of any trapping potential (i.e., ω x = 0), the system is formally integrable [23], and the first exact results for the case of attractive interactions were obtained in 1964 by McGuire [24].

Hartree factorization: The Gross-Pitaevskii equation
In this approximation one assumes the many-body wavefunction ψ( x) to be factorizable into product form, such that each individual atom is described by the same singleparticle wavefunction φ(x k ). Hence, the Hartree wavefunc- of such a stationary state with respect to variations in φ(x k ) leads to [34] with φ(x) normalized to unity, and µ a Lagrange multiplier, given by Equation (3) is the one-dimensional time-independent Gross-Pitaevskii equation (GPE) [33], which formally tends to an exact description as N → ∞ while g 1D N is held constant [35,36]. In this limit g 1D (N − 1) ≈ g 1D N, and it is more typical for the coefficient of the nonlinearity in Eq. (3) to be set proportional to N. As we will also consider small particle numbers, in what follows we choose to retain the proportionality to (N − 1).

Rescaling to dimensionless form
It is convenient to rescale our description of the system in terms of an effective = M = g 1D (N − 1) = 1 unit system, referred to as "soliton units" [20,30]. Space, time and energy scales are then given in units of 2 /Mg 1D (N − 1) (the classical soliton length [30]), 3 We work within this system of units from this point onwards. Equation (1) then simplifies tô  (6) and the GPE [Eq. (3)] becomes We have introduced the dimensionless parameter γ, which is the square of the ratio of the classical soliton length to the harmonic length √ /Mω x [30], i.e., Within our chosen system of units γ appears in the rescaled Hamiltonian and GPE as a dimensionless effective trap frequency. This also reveals γ to be the only free parameter in the GPE, which as a description of the system is effectively a classical field limit, and the particle number N appears as an additional free parameter in the fully quantal Hamiltonian.
B. Known exact results

Exact results in free space
In the case where there is no axial trapping potential, i.e., γ = 0, the many-body eigenstates [24,37] and GPE stationary states [14,30,34] are known. The stationary solutions to the GPE [Eq. (7)] minimizing the energy are classical bright solitons [20,30] where the value of x 0 is arbitrary (as is that of an irrelevant global phase). The Hartree approximation to the many-body wavefunction is thus and is localized around x 0 , i.e., ψ H ( x) → 0 as |x k − x 0 | → ∞. Exact solutions also exist in box and periodic boundary conditions, given by Jacobi elliptic functions [38]. The exact ground state wavefunction for Eq. (6) with γ = 0 (see Appendices A and B 2) is an N-particle bound state, proportional to [24,37] As the Hamiltonian [Eq. (6)] has no external potential, all its eigenfunctions are independent of the center-of-mass coordinate (see section II B 2) and there exists a continuum of moving N-particle bound state eigenfunctions ψ(P, x) = e iPx C ψ G ( x)/ √ 2π. The normalization convention is then such that d xψ * (P, x)ψ(P , x) = δ(P − P ), and the ground state is written as 2π. The exact ground state is thus completely delocalized in the center-of-mass, and is only localized in the sense that ψ G ( x) → 0 as |x k − x j | → ∞, for any k, j.
The localized Hartree solution ψ H ( x) violates the translational symmetry requirement imposed by the absence of external potentials in the Hamiltonian due to the tacit assumption that the minimizing wavefunction should vanish as x k → ±∞ [34]. Note, however, that the particle densities about a specified value R of the center-of-mass coordinate [given by the expectation value of δ(R− x C ) N k=1 δ(x− x k )] corresponding to ψ H ( x) [Eq. (10)] and ψ G ( x) [Eq. (11)] agree to order 1/N [39], and hence are identical in the limit N → ∞. The energies E H and E G corresponding to the wavefunctions ψ H ( x) and ψ G ( x) are given by .
As one would expect, the exact eigenenergy E G is less than E H , and the difference in energies per particle (E H − E G )/N = 1/12(N − 1) vanishes as N → ∞.

Separation of the center-of-mass coordinate x C
In the case of any external potential being either harmonic or nonexistent, the center-of-mass dynamics separate and are independent of any two-body interactions. Consequently the center-of-mass eigenstates are simple harmonic oscillator eigenstates or plane waves, respectively, in the former case this is referred to as the Kohn mode. This may be readily seen by expressing the first-quantized form of the Hamiltonian Eq. (6) in terms of the Jacobi coordinates, i.e., x C [Eq. (12)] together with for k ∈ {2, 3, 4, . . . , N}. The Hamiltonian can then be phrased ξ is a shorthand for {ξ 2 , ξ 3 , ξ 4 , . . . , ξ N }, and we have used the identity x k − x j = ξ k + k−1 = j+1 ξ / − [( j − 1)/ j]ξ j (with b > a and ξ 1 ≡ x C ). In cases where the upper limit of a sum is less than its lower limit, the sum is taken = 0.
Hence, the normalized ground state of H C is exactly with eigenenergy = γ/2.

Two interacting bosons in a harmonic potential
The case of two identical bosons in a harmonic potential with contact (δ-function) interactions is also exactly solvable [40,41]. In this case the eigenfunctions of H R (ξ 2 ), defined through H R (ξ 2 )φ n (ξ 2 ) = E R,n φ n (ξ 2 ), are given by φ n (ξ 2 ) = N n U(−ν n , 1/2, γξ 2 2 /2)e −γξ 2 2 /4 , where U(a, b, z) is the Tricomi confluent hypergeometric function [42], and N n is a normalization constant. The ν n are implicit solutions of and set the eigenvalues of H R (ξ 2 ) through Attractive interactions must reduce E R,0 from the noninteracting case, so that E R,0 < γ/2 ⇒ ν 0 < 0. As outlined in appendix C, it then follows that in the limit γ → 0 (interaction dominated regime) E R,0 → −1/4 + O(γ 2 ). This is in agreement with the total ground state energy E G for the case of two attractively interacting bosons in free space [Eq. (14)], as one would expect due to the center-of-mass energy of the free space ground state being = 0. In the opposite limit of γ −1 → 0 (trap dominated regime) harmonic oscillator eigenvalues and eigenfunctions must result, i.e., E n → (2n + 1/2)γ and U(−ν n , 1/2, γξ 2 2 /2) → H 2n ( where the H 2n are even Hermite polynomials. 1

A. Interaction dominated limit in a harmonic potential
In the case where γ 1, we may consider the effect of the trap to be dominated by the effect of the interactions, and therefore negligible in H R . As there are no interactions present in H C , the effect of the trap is in this case always significant, even in the interaction dominated regime.
We may therefore consider a limiting case Hamiltonian H 0 , composed of H R [Eq. (17)] with γ = 0, plus H C [Eq. (16)]. Written in terms of conventional single particle coordinates, and the correctly normalized ground state ψ 0 can be put together from Eq. (11) multiplied by Eq. (18), i.e., ψ 0 ≡ ψ C ψ G , with the sum of the corresponding eigenvalues determining the overall energy E 0 . Hence, in terms of single particle coordinates and, from Eq. (14) plus γ/2 (the harmonic oscillator zeropoint energy), This interaction dominated limit does not correspond to any physical system but is a useful starting point for perturbation theory.

B. Perturbation results
To proceed from the approximated Hamiltonian (23), we include the effect of the harmonic trap on the relative degrees of freedom via Rayleigh-Schrödinger perturbation theory. The full Hamiltonian (2) can be written as As ∆H(x) ∝ γ 2 , we expect perturbation theory to yield particularly good results in the limit of small γ. For the first-order energy correction to the ground state, which also serves as a definition of the bra-ket notation, we find (Appendix B 4): The sum in Eq. (27) is simply the second Harmonic number, for which the asymptotic behavior in the N 1 limit is given by [43] Thus, asymptotically the energy correction goes as For large N, this coincides with the result obtained using the free space Hartree solution, given in Eq. (10), as an approximation for the ground state (Appendix D) These results are displayed in Fig. 1; for small N there is a large difference between the result predicted by the Hartree product state [Eq. (27)] and the result predicted by the exact many-body ground state [Eq. (30)] of the approximate Hamiltonian (22). There is also a weak number dependence from the Harmonic series in Eq. (27). However as N 1 both methods give the same energy correction per atom, π 2 γ 2 /6. For N = 1000 the relative difference (E (1) H − E (1) )/E (1) ≈ 0.0016 is already small. First-order energy correction per atom for a many-body soliton with the center-of-mass in the lowest eigenstate of a harmonic oscillator given in Eq. (23). The exact solution given in Eq. (27) and its expansion up to next-to-leading order, given in Eq. (29), begin to agree for N 10 with the latter always underestimating the true value. Both curves approach the approximate result predicted by the Hartree approximation [Eq. (30)]. The relative difference between different predictions lies below 1% for N 165.

C. Variational minimization
In order to improve the value for the ground state energy beyond the first-order-perturbation-theory result (27), we use the (normalized) variational ansatz with γ > 0 and the constant N the same as Eq. (11) 2 , which is calculated in appendix B 2. Since the center-of-mass wavefunction is unchanged, we will only have a correction to the relative energies, these are calculated in Appendices B 3 and B 4. The total energy for this wavefunction is for the expectation values of each section of the relative Hamiltonian, with E G being the (negative) ground state energy of the free soliton in soliton units given in Eq. (14) and E (1) the first-order correction given by Eq. (27). In order to calculate the energy minimum, the derivative of Eq. (33) with respect to λ has to be zero: which (for λ 0) is equivalent to a 4th order polynomial in λ where the constant κ is defined as the ratio of the ground state energy and first-order correction For fixed N, κ ∝ γ 2 , the value of this prefactor is an increasing function of N with a minimum of κ = 2γ 2 at N = 2 with an asymptotic limit of [cf. Eq. (29)]: Thus, κ is small for γ 1. Equation (35) has four roots, only one of which is real and positive, which is the root of interest. The exact analytic solution is given in Appendix E.
If κ 1, this solution is approximately which leads to the minimum in the energy of As the variational Ansatz (31) does not affect the centerof-mass part of the wavefunction, calculating the overlap between this variational Ansatz and the state Eq. (23) (i.e. the λ = 1 state) is an interesting physical quantity: its modulus squared is the fraction of the relative wavefunction which is projected to the relative ground state if the trapping potential was turned off quasi-instantaneously (cf. [44]). The overlap is given by (see Appendix B 3): Using the approximation (38), the overlap (40) approximately is [1+κ 2 /8+O(κ 3 )] (1−N) and we thus expect the overlap to vanish in the limit N → ∞ for κ > 0. Rather than investigating the total wavefunction overlap (40), the Nth root of Eq. (40), an effective single particle overlap, is a more suitable value in the limit N 1 as it tends to a constant as N → ∞ and is related to comparing two GPE orbitals. Note that for two Hartree-product wavefunctions, the effective single particle overlap would be independent of N, but the Nth root of Eq. (40) still is N-dependent due to the N-dependence of λ 0 [Eq. (38), cf. Appendix E]. Figure 2 (a) shows the overlap (40) as a function of γ for various particle numbers. For λ 0 , the exact value given in Appendix E was used. As expected, the N-dependence is quite strong. Figure 2 (b) shows the Nth root of the overlap (40), i.e. the effective single-particle overlap. The effective singleparticle overlap is larger than 0.99 for γ 0.15 for all N, indicating that ψ 0 ( x) from (23) is still a good description in this parameter regime and the trap has had little effect on the internal degrees of freedom. The limit N → ∞ is nearly reached for particle numbers as low as N = 100 [note that in panel (a), the limit N → ∞ would lie on the coordinate axes].  40), of the variationally obtained solutions for different rescaled trap frequencies γ with the free space ground state solution (γ = 0) with a Gaussian envelope for the center of mass. Effective singleparticle overlap is treated as the Nth root of the total overlap as for two different product states this is independent of number and equal to the overlap between the single particle wavefunctions. Bottom to top the solid lines on both graphs correspond to N = 100, 10, 6, 3, 2, the dashed line corresponds to the N → ∞ limit of the variational many-body solution [using κ from Eq. (37)] and is very close to the N = 100 line.

A. Overview
While the focus of the previous section lies on the case of small γ, the numerical methods introduced in this section work well for γ 0.16.

B. Computation procedure
We expand the field operator over the set of Hermite functions where H k are the Hermite polynomials, givingΨ(x) = k ϕ k (x)â k . The Hamiltonian (5) in this basis can be split into three separate parts: the kinetic Hamiltonian, the potential Hamiltonian, and the interaction Hamiltonian. The factor of f k mn is the integral of four Hermite functions (with W set to unity) over all space, i.e.
Here, the functions are also real, so there is no need to take the complex conjugates. This can be calculated exactly in terms of gamma functions Γ(x) and a standard hypergeometric function 3 F 2 evaluated at unity [45] Without interactions, the ideal gas Hamiltonian is given bŷ For W = √ γ, the basis states (23) are eigenstates of the noninteracting Hamiltonian. The total Hamiltonian can therefore be expressed aŝ and we refer the ground state of this as |ψ g (γ) and the ground state energy as

C. Truncation and projection to center-of-mass excitation basis
In order to do computations we must only use a finite basis set, which will introduce the inaccuracy. This is discussed in Appendix F. Essentially all possible states for which the eigenenergy related to the Hamiltonian (47) lies below an energy cut-off E cut are included, we refer to this set as the 'truncated basis'.
The brute force procedure would now be to calculate all the matrix elements of the interaction Hamiltonian (44) using this truncated basis and add the matrix of energies of the kinetic and potential Hamiltonians and diagonalize this to get the eigenstates and energies. However we can reduce the size of the truncated basis set used in this computation by recalling from Eq. (16) that the center-of-mass Hamiltonian commutes with the relative Hamiltonian and thus they have separate eigenstates.
Inspired by the ladder operator treatment of a single particle in a harmonic oscillator, we define Noting we can express the center of mass Hamiltonian as , as is the case of the single particle ladder functions. Moving to second quantization, we can equivalent operators in terms of creation and annihilation operators in our basis (for W = √ γ) via [46] whereÂ + = (Â − ) † . These satisfy [Ĥ 0 ,Â ± ] = ±γÂ ± and thus energy levels spaced in units of γ; they also commute with the interaction Hamiltonian [Ĥ I ,Â ± ] = 0 as required. Note that for N = 1,Â + acting on the ground state can be used to construct all the eigenstates of the system. We then can express the center-of-mass Hamiltonian aŝ In order to use this property to reduce the basis, we must first transform our truncated basis (E ≤ E cut ) to basis states which are eigenstates of the center-of-mass Hamiltonian, this procedure is detailed in Appendix F 3 along with the reduction in basis size it achieves. We keep only eigenstates ofĤ C with eigenvalue γ/2 (center-of-mass ground state) and project the Hamiltonian (48) from the truncated basis to this new and reduced basis. For high N ≥ E cut /γ − N/2, the reduction asymptotes to π/ 6(E cut /γ − N/2) for E cut /γ − N/2 1. For details see Appendix F 4.

D. Using different-width Hermite functions
Using functions with a W = √ γ, such that they are eigenstates ofĤ ideal , is not desirable in the γ → 0 limit because the basis will consist of states much wider than the wavefunction we are using them to construct. For an infinite basis, the ground state should be independent of the basis used to describe the system (in our case, it should be independent of the value of W). For numerical calculations, the basis will be finite and thus some choices of W are better than others. In Sec. V, we will calculate the ground state for γ = 0 in order to determine the optimal value for W to be used in the calculations.
For arbitrary W, the Hamiltonian now reads: which includes extra mixing terms in the ideal gas Hamiltonian (47). This causes a fairly significant issue in that it is no longer possible to exactly separate center-of-mass eigenstates in this basis, meaning the full basis would need to be used in order to achieve the center-of-mass ground state. Using this method with just a truncated Hilbert space and no projection to the sub space with zero center-of-mass excitation would make achieving convergence painfully slow. The solution to this is therefore to reduce the basis in the same way as before, but accept that the center-of-mass wavefunction we end up with is given by which is not an eigenstate and has energy E C = (W 2 + γ 2 W −2 )/4 rather than the true γ/2, thus we know the true ground state is the wavefunction we obtained, . This approach has the huge advantage that, if W is kept constant, the occupation of the basis states for the ground state should change very little as γ → 0 where they will tend to the solutions on the infinite line.

E. Numerical ground states within the GPE approximation
Within the GPE approximation, we can obtain the ground state by solving Eq. (7) as the ground state is the only stationary state of the system. The method used here is to again expand over a finite basis set of Hermite functions of arbitrary width scaling W then to produce a set of η + 1 nonlinear equations in the coefficient set c (which will be real), by integrating Eq. (7) multiplied by ϕ k (x) over all space, for k = {0, . . . , η}, giving There is also an (η + 2)th equation, relating to the normalization k |c k | 2 = 1. Denoting the vector with an equation at each position as F(c), we wish to solve F = 0. We use Newton's method (as in [46]) to iteratively solve for c, via where J is the η + 1 by η + 2 Jacobian matrix associated with F. η is increased until convergence is achieved.

V. EFFECTS OF HARMONIC CONFINEMENT
A. Ground state energy Using the methods from the previous two sections, we investigate the effect an external potential has on the relative component of the ground state |ψ g (γ) [cf. Eq. (49)]. This is important to quantify how soliton-like the state is, along with what excitations can be expected if the state is released quasi-instantaneously from the potential.Such dynamics have already been considered using the GPE in [44]. Figure 3 shows ∆E/N, the energy difference per atom between the numerically calculated ground state energy E g (γ) and the ground state energy of the artificial Hamiltonian (23) given by Eq. (24), for a range of γ and N values. It is produced by calculating the ground state energy via the three numerical methods, namely exact diagonalisation in a basis of hermite functions with either optimized widths for weak trapping (shown in table I), or widths which are eigenstates of the non interacting problem, and variational minimization, for a range of γ and taking the smallest value. This is because, due to the variational principle, all of these techniques produce only values greater than or equal to the ground state energy, hence the lowest is the best estimate.

B. Universal behaviour
If we consider instead a rescalingγ = γ(N − 1) 2 /N 2 and ∆Ẽ = ∆E(N − 1)/N (i.e. converting to a unit system in which g 1D N = 1 as opposed to g 1D (N − 1) = 1), a more universal behaviour is present in ∆Ẽ, with little number dependence as shown in fig. 3 (b). To see this analytically, we note that for γ 1, our variational result of Eq. (39) for the energy is basis size   3  84  2  631  6  38  1  3009  10  28  0.5  2534  100  24  0.5  1575   TABLE I. This table shows the parameters used in calculating the graph in Fig. 3, where N is atom number, η is the cut-off and W is the width taken for the fixed width calculations (coarsely chosen to minimize the ground-state energy at γ = 0). The reduced basis size is the number of states (with zero center-of-mass excitation) used in the exact diagonalisation of the Hamiltonian. The basis are chosen to be a reasonable computational size, however the numerics are less reliable for small γ applicable. ∆E is obtained by subtracting the factor of E G + γ 2 /2, then converting to our rescaled units we have the N dependent factor of orderγ 2 (which is the rescaled first order energy correction) is 2 for N = 2 and decreases monotonically to π 2 /6 ≈ 1.6 as N → ∞, hence for very smallγ, the N = 2 line is largest, but the difference is very small. The orderγ 4 term has negligible number dependence and so is unlikely effect the ordering of these lines within for the range of variational models validity. On the other end of the scale, as γ → ∞ (the trap dominated system) we can neglect interactions in Eq. (5), giving a ground state of a product of Gaussians of width 1/γ, subtracting the center-of-mass energy gives ∆E → γ(N − 1)/2N + O( √ γ) or, in our rescaled units, ∆Ẽ →γ/2+O( √γ ) and so to leading order, the N dependence vanishes.

C. The classical soliton limit
As shown in Fig. 3, for low γ the variational ansatz (31) gives the best estimate for the ground state energy of all the methods used in this paper. For low enough γ, this variational ansatz in turn is very close to the product of the free manyparticle solution with a Gaussian center-of-mass wavefunction (23) (see Fig. 2). In the limit of small γ, the integral can thus be used to investigate deviations of the Hartreeproduct wavefunction (10) from the true many-particle ground state (which is well approximated by ψ 0 ( x) for γ 1). As was the case for Fig. 2, the effective single particle overlap B 1/N will also be considered as we are interested to see how well the wavefunction is described by a product state, and  when comparing two product states with different single particle wavefunctions, this quantity is constant with changes to number.
In order to make an educated guess about what range of γ will give a large overlap, we look at the expectation value of the square of the center-of-mass location over the Hartreeproduct wavefunction (Appendix D): This value is identical to the variance of the center of mass, as both the many body and Hartree states are centered about x = 0. A variance calculation can also be performed for Eq. (23), this is particularly simple as the center-of-mass is explicitly separate and is given by As the Hartree product state is uncorrelated, for a large enough N the distribution associated with the center-of-mass location will therefore tend to a Gaussian (with variance of π 2 /3N) via the central limit theorem. We therefore consider an effective center-of-mass wavefunction that is the square root of this dis-tribution, yielding the overlap integral which reaches its maximum I (γ max ) = 1 for As the above analysis focuses on the center-of-mass part of the wavefunction, and thus Eq. (62) is likely to overestimate the overlap B as defined in Eq. (59), Eq. (62) also predicts a γ → 0 behaviour of the form I(γ) ∼ cγ 1/4 with c a constant. Figure 4 shows a numerical calculation of B for a range of γ and N, the integration is performed via Monte Carlo methods, i.e. weighted sampling using random variables with a sech(x/2) 2 /4 distribution (obtained via the ziggurat algorithm [47]) until a standard error of < 10 −4 was obtained. It can be seen from Fig. 4 (a) and (b), that maximum overlap occurs just slightly above γ = 0.16 [close to the analytic estimate (63)] and improves as N increases. Based on our previous discussion of effective center-of-mass width, the top value should relate to the overlap of the relative degrees of freedom, although this is not well defined. Graphs (c) and (d) suggest that most of the γ dependence in B is due to the effective "center-of-mass width" of the Hartree product solution, since the shape of each overlap curve is similar, besides a small offset, to the dash-dotted line. This indicates that a Hartree soliton is a very good approximation to the many body solution if the center-of-mass wavefunction is also localized.
show the Nth root of (a) and (b), effectively overlap at the level of single particles, which tends extremely rapidly to unity as N increases for any γ over the range shown. A useful point that this high overlap implies is that the many-body state Eq. (23), is extremely well approximated by the Hartree product state if the center-of-mass envelope squared is approximately the statistical distribution that would arise from taking the mean of the N independent probability distributions |φ| 2 [with φ given in Eq. (9)] associated to single atom positions in the product state. For this reason the Hartree product state would be expected to well approximate the ground state of the system, even at a many-body level, if the center-of-mass envelope is localized to the size of this distribution (i.e. γ ≈ 3/2π 2 ) by the potential. The converse to this is also true, an initial condition that is given by Eq. (10) is well approximated by Eq. (23) with γ ≈ 3/2π 2 . This could be used to estimate center-of-mass position uncertainty of a state, initially given by a Hartree product wavefunction, as it evolves in time, using known results for the spreading of Gaussian wavepackets.
We also consider how these many-body effects would affect experimental observations. From a measurement the atomic density, one could use the mean of this signal to determine the center-of-mass of the system. If the state of the system is well approximated by Eq. (23), the observed location would vary shot to shot with a probability distribution given by |ψ cm (x C )| 2 ∝ exp(−Nγx 2 C ) (combined with any experimental uncertainties associated with density measurement). For γ < γ max , this distribution would be wider than one would expect using the product approximation, most notably if the center-of-mass wavefunction is wider than a classical soliton width 1/ √ γN 1, this jumping effect would be most clearly visible. Non zero temperature would further increase this effect, by introducing a statistical mixture of excited states of the center of mass. As a purely mechanical analogy, one could think of taking a photo of a swinging pendulum at a random time, the shape always looks the same but its position appears random.

VI. CONCLUSIONS
We study a 1D system of identical Bosons with attractive contact interactions, a Lieb-Liniger(-McGuire) gas, in the presence of a harmonic trapping potential. We present variational and numerical many-body calculations, in both cases making use of the separability of the center-of-mass Hamiltonian to split the problem into relative and center-of-mass degrees of freedom. We use a unit system such that the Hartree soliton length is set to unity ( = m = g(N − 1) = 1), leaving two parameters, the number of atoms, N, and √ γ, the dimensionless ratio between the Hartree soliton and harmonic oscillator lengths.
Our key results are firstly that we have derived a first order energy correction to the ground state of the relative degrees of freedom from the introduction of a harmonic oscillator potential [given in Eq. (27)], which is used in a variational minimization technique. This is proportional to γ 2 and the correction per atom tends to the mean field prediction from below, the relative difference is less than 1% for N > 165.
Secondly we have determined the validity range of γ of our many body ansatz, consisting of the free many-body ground state with a Gaussian envelope as given in Eq. (23). Essentially as the trapped ground state deviates from this it becomes less "soliton like", we quantify this with the "effective single particle overlap", given by the Nth root of the overlap between a variationally obtained ground state and our ansatz. For N large, this overlap is greater than 0.99 for γ 0.16. Numerical calculations of energy in the strongly trapped region, γ > 1, indicate energies are still considerably lower than the non interacting case.
Thirdly we show, via a numerical investigation of overlap between the free Hartree product solution and the free many-body ground state with a Gaussian envelope [given in Eq. (23)] describing the center-of-mass wavefunction, that the two wavefunctions can have high agreement, even at a manybody level. This high overlap occurs when the modulus square of the envelope function matches the probability distribution, associated with the Hartree product, for the center-of-mass position, which occurs when γ ≈ 0.16. However, current experiments with bright matter-wave solitons are such that the center of mass is localized to much less than a soliton width, indicating this is unlikely to be an observable effect.
In addition to these physical results, we outlined a numerical method for computing many body eigenstates, this uses a basis set of harmonic oscillator eigenstates, truncated at a particular energy. This is then project a into a subspace of states with the center-of-mass wavefunction in a specific state, using the ladder operator for center-of-mass excitation. This makes use of the separability and achieves a reduction in the size of the basis set required by a factor of π/ √ 6E cut (where E cut is a cut-off energy) or better, greatly improving speed of the diagonalization and allowing us to investigate the internal degrees of freedom separately.

ACKNOWLEDGMENTS
We would like to thank T.P. Billam and J. Brand for discussions, and the UK EPSRC for funding (Grant No. EP/G056781/1) and the Jack Dodd Centre (S.A.G.) for support.
Appendix A: The 1D free system and the full eigenspectrum via the Bethe Ansatz We briefly recapitulate aspects of the treatment of the attractively interacting Lieb-Liniger gas [24,37] in order to set notation within our chosen system of units.
In order to find the solution to the ground state of Eq. (6) with γ = 0, we note that the wavefunction in the region x 1 < x 2 < . . . < x N is solved by where a sum over all permutations of the set P = {1, . . . , N} is performed to make it symmetric, the energy eigenvalue is thus simply equal to E = k p 2 k /2. Each permutation has a coefficient associated with it that is linked to the boundary conditions when x k = x k+1 , for an interacting system they can be determined by the equation [37] where P is the permutation swapping the kth and (k + 1)th indices, the coefficient of the identity permutation is determined by the normalization condition. The center-of-mass motion is independent of the interactions, and so will have eigenstates of plane waves. In the case of attractive interactions. These momenta can also have very specific imaginary components corresponding to bound clusters of atoms. The ground state of relative motion occurs for in which all the permutation coefficients apart from one (the identity permutation) are equal to zero. Higher eigenstates can have multiple bound state clusters or strings, each with an associated real momentum P m and imaginary components (that must sum to zero) which are spaced in units of 1/(N − 1). If we have η clusters, each of size n m , we have momenta associated with the mth cluster given by This state would normally be denoted as |n 1 , p 1 , n 2 , p 2 , . . . , n η , p η , the total energy of the state, E = k p 2 k /2, scales as though these are η isolated single soliton states and thus the energy eigenvalue is The total number of different combinations of clusters scales p(N), the number of ways to partition N with integers (cf.

Appendix F 4 and [48])
Appendix B: Normalization, energy and overlap using the variational state

Preamble
In order to make use of our variational state given in Eq. (31), we must calculate the normalization constant and expectation value of energy. Calculations for the energy and normalization constants for all the eigenstates in free space (γ = 0) can readily be found in literature [25], it is also the case that the center-of-mass component of the Hamiltonian can be considered separately and so taking a finite center-ofmass component does not significantly alter the calculations. However, the choice of normalization condition for a non local system is somewhat arbitrary and conventions vary between papers, we choose a normalization that means both the relative and center-of-mass parts are normalized to unity with respect to Jacobi coordinates. Also most derivations of the energy rely on the fact that the gradient discontinuity at the points x k = x j in the wave-function, exactly cancel the interaction terms (essentially from the condition of being an eigenstate), and thus these terms can simply be ignored. Because of our variation of λ, this will no longer be the case and thus we are forced to make a more explicit calculation of the kinetic energy. In addition to this we derive a first order energy correction to the relative degrees of freedom, which is a new result.

Normalization
In order to normalize Eq. (11), we could insist that two states with different center-of-mass momenta are orthonormal, i.e. p , N|p, N = δ(p − p) such as was calculated in [37] or consider wave-function to be trapped in a box which we allow to grow infinitely large [44]. However we are interested in the normalization of the free space solution with a Gaussian center-of-mass envelope and freedom to tune a variational parameter, denoted ψ (λ) var in Eq. (31). This result and technique will also be used in appendix B 4 and follows the method of [37]. We consider a Fourier decomposition of the wave-function with σ = 1/(N − 1) corresponding to Eq. (11), however for greater generality we allow this parameter to be free in order to use these results for variational calculations where σ → λ/(N − 1), which would correspond to Eq. (31). Also one may wish to consider instead units in which the harmonic oscillator frequency and length are set to unity and the interaction constant rescaled tog, in which case the replacement σ → λ|g| would be used instead, or indeed in S.I. units σ → λM|g 1d |/ 2 . Essentially this term serves to allow easy conversion between unit systems and making variational manipulation easier.
In the form of Eq. (B1), it is far simpler to perform the integrals of the coordinate variables. Calculating ψ (λ) var |ψ (λ) var in coordinate space will require integration over N spatial integrals and two momentum integrals, however we only need to integrate over the simplex region x 1 ≤ x 2 . . . ≤ x N as by Bose symmetry any integration over any such region will be identical, hence we multiply by factor of N! to include all possibilities for such a regions construction. Within this simplex region, all arguments in the absolute value signs are positive and the wave-function is given by with β(k) = (N + 1 − 2k)σ, and using the notation we can now express the inner product as transformation of variables p = (p 1 + p 2 )/2 and p = p 1 − p 2 with Jacobian unity then allows us to perform the integral over p leaving To perform the remaining integrals we note that y −∞ dx exp(ax + by) = exp[(a + b)y]/a, denoting and noting a(N) = 0, we can recursively use the previous result to perform all but one of the spatial integrals and obtain Integrating the final term gives 2πδ(p √ N) and the momentum integral is then trivial, noting that A(N, 0) = 1/(N − 1)! 2 σ (N−1) gives us the final result for the normalization factor As was mentioned before, the normalization factor for Eq. (31) in our units is obtained by letting σ = λ/(N − 1), in the case of λ = 1 where this relates the ground state in infinitesimal trapping Eq. (11), we refer to this constant simply as N. It is also worth noting that both the com wave-function and the relative are both chosen to be normalized to unity with respect to Jacobi coordinates, hence the (γ/Nπ) 1/4 relates to the center-of-mass part and the rest to the relative component.

Kinetic and interaction energy
We wish calculate the kinetic energy and potential energy of the variational state, i.e. the expectation of Eq. (6) with λ set to zero on Eq. (31), which we will denoteĤ free . Due to the separability of the wave-function and Hamiltonian, it is sufficient to consider only the relative part of the wave-function and note that the center-of-mass kinetic energy is given by γ/4. We first denote ϕ(x 1 , .., x N ) = exp −σ N k=2 k−1 j=1 |x k − x j |/2 , being the relative part of the variational wave-function (up to a normalization factor) and calculate the second derivative with respect to some coordinate x − 1 2 The first term in Eq. (B9) can be split up into terms of the form sgn 2 (x − x b ) = 1, of which there are (N − 1) and terms of the form sgn(x − x a )sgn(x − x b ) with a b, of which there are (N − 1)(N − 2). The former will evaluate to unity by normalization of the wave-function, however the latter terms will equal +1 when x < x a < x b or x < x b < x a and when ; the wave-function must be identical in all these 6 simplicies due to Bose symmetry and so the expected value of these terms will equal 1/3. When the sum over all is performed, these terms will total to N(N − 1)σ 2 (1 + (N − 2)/3)/4. The latter terms of the form δ(x − x k ) can then be combined with those from the interaction part of the Hamiltonian, noting that there are twice as many terms but δ(a − b) = δ(b − a).
all that remains now to calculate the value of the expectation value of the delta function terms. Following the method in Appendix B 2 we integrate over a simplex region −∞ < x 1 < x 2 ... < x N < ∞, as a result of this we need only consider the N − 1 terms of the form δ(x k − x k+1 ) as the rest will be zero. Each integral will be the same as in Appendix B 2 except missing a factor of 2/a(k) for each term δ(x k − x k+1 ), hence the result will equal N−1 k=1 a(k)/2 (using the result y −∞ dx f (x, y)δ(x − y) = f (y, y)/2). Hence finally, substituting in this result into Eq. (B10) we have where E 0 = −N(N + 1)/24(N − 1). This discussion has not mentioned the harmonic envelope of the center-of-mass function, however due to the separability of the Hamiltonian this will only add factor of the center-of-mass kinetic energy, which will be independent of λ regardless of what the centerof-mass wave-function is. This energy term is combined with the potential energy calculation derived in Appendix B 4 to form the basis for a variational principle.

Derivation of the first-order perturbation energy
This section is related to the calculation of ψ (λ) var |V(x)|ψ (λ) var , whereV(x) = γ 2 N k=1 x 2 k /2, note that Eq. (27) is given by this quantity minus the center-of-mass energy. It is again easier not to perform this integral in Jacobi coordinates but to Fourier transform out the center-of-mass, we will also again replace the factor 1/(N − 1) with σ to generalize the results for our variational principle. This calculation is similar to, although more complicated than, the calculation of the normalization factor; to that end we can start the calculation from Eq. (B4) (as no spatial integrals are yet performed) adding in the potential factor V(x) giving The same recursive integral procedure can be applied here except we need two additional results, true for real(k) > 0 Let us consider only the latter part of Eq. (B13) omitting the constant √ Nγ 2 /4πA(N, 0), taking the integrals in order from x 1 to x N , the integral over x will be over a function of the form with k 0 (1) = k 1 (1) = 0 and k 2 (1) = 1 and A( , p ) defined in Eq. (B7). The common prefactor of A( , p ) is the equivalent of k from Eq. (B14). Besides this, each integral will increase the factor in front of the x 2 term by one each time and hence k 2 (l) = l. Contributions to k 1 ( + 1) come from k 1 ( ) and k 2 ( ) and as such Eq. (B14) implies k 1 ( + 1) = k 1 ( ) − 2k 2 ( )(a( ) + ip / √ N) −1 , given that k 1 (1) = 0 this implies Applying this same induction logic down to k 0 gives simply evaluating these at N and performing the final integration over x N then yields we then insert this expression back into Eq. (B13) giving The integral over the δ(p ) term can be performed immediately and gives γ 2 k 0 (N)/2. Considering next the integral over δ (p ); since exp(−p 2 /4γ) has zero gradient at the origin it will not contribute, however the terms will contribute to Eq. (B19), giving Finally for the δ (p ) term we must include A (N, p )| p =0 = −A(N, 0) N−1 k, =1 (1 + δ kl )k /Na(k)a( ) and the differential of a Gaussian, hence we have Summing these three terms together, and substituting k 2 (N) = N, we are left with The first term in this expression is equal to γ/4, which is simply the potential energy of the center-of-mass component.It can be proved via induction [49] that the double sum is equal to thus reinstating σ = λ/(N − 1), the remaining terms simplified down to which is used in Sec. III. Equation Eq. (27) is the first order energy correction to the free soliton with Gaussian center-ofmass envelope and is obtained by subtracting the center-ofmass energy and setting λ = 1 a. Energy correction from potentials of higher powers of x An energy correction for general power law potentials can be derived in the mean field case. For Re(m) with ζ(m) the Riemann zeta function. A similar result would be desirable to calculate energy correction from an anharmonic potential for a quantum soliton, although potentials with m 2, 0 will couple the center-of-mass and relative degrees of freedom together (possibly only very weakly) and so [Ĥ cm ,Ĥ rel ] 0 and this is only of limited use.

Overlap of the relative components of the variational wavefunctions
Finally we consider the overlap between the relative parts of the variational wavefunction with λ > 1 and the ground state in infinitesimal trapping γ = 0, λ = 1, given by This calculation can be achieved by performing the calculations in Appendix B 2 with σ → (1 + λ)/2(N − 1), the resulting factor will not equal unity and instead will be equal to N 1 N λ /N 2 (1+λ)/2 . Therefore that the overlap is given by which is used in Sec III C.
All the cross terms of the form x k x j will evaluate to zero as sech(x) is an even function, thus leaving only the power terms.
By Bose symmetry f (x k ) = f (x j ) and thus the value of all the terms in the first sum will be identical to the single particle correction and we have The relation between energy level degeneracy in systems of identical particles and number partitioning has been investigated in [51][52][53] and references therein. Within a one dimensional Harmonic oscillator, the degeneracy for N distinguishable particles scales the same as the degeneracy for one particle in an N dimensional spherically symmetric potential. This is not the case for indistinguishable particles, to calculate these we must use introduce integer partition functions. We introduce the notation p([a, b], m) being the number of ways to partition an integer m using only integers a ≤ z ≤ b, in order to compute these for a given b, we use the recursion relation This works by noting that we can divide a partition into two distinct sets, partitions which uses only numbers larger than a, being p([a + 1, b], m), and partitions which uses a at least once in the partitions, p([a, b], m − a). Also p([a, b], 0) = 0 by convention.
Using the usual Fock space representation of these harmonic oscillator states |N 0 , N 1 , . . . with k N k = N and defin-ingẼ as the energy of the state (with no interactions) minus the ground state energy divided by γ Given that each occupancy of the kth mode raises the energy by k it can be seen that the degeneracy of the energy levelẼ is given by the number of ways to partitionẼ using N nonnegative integers. Denoting Φ(Ẽ, ) as the ways to partitionẼ in numbers we have It is also known that this sum is equal to the number of ways to partition an integer 'Ẽ' using only numbers less than or equal to N i.e. g(Ẽ, N) = p([1, N],Ẽ).

Truncation with an energy cut-off
In order to make a basis computationally manageable, it must truncated to be made to be finite. This is achieved by only taking states with energy less than an arbitrary cut-off η, note that this also implies that N k = 0 if k > η. The size of this truncated Fock state basis is given by The reason an energy cut-off is chosen rather than a mode cutoff at η (although as mentioned before this is implicit in an energy cut-off method) is two fold. Firstly in order to project into the center-of-mass and relative excitation basis we require all the states with a given energyẼ [the Hamiltonian (52) is block diagonal], if we do not have all those states the projection is not possible. Secondly having just a mode cut-off would include the state |0, 0, . . . , N withẼ = Nη, but not the state |N − 1, 0, . . . 0, 1 (one occupancy in the η + 1th mode) withẼ = η + 1, as long as harmonic oscillator energy remains a non negligible quantity, the former state will have almost no mixing to the ground state, making it a very inefficient truncation.

Deriving the projector to the center-of-mass basis
As we have expressed the Hamiltonian (48) in terms ofâ † k andâ k , the creation and annihilation operator for bosons in mode k, it is far simpler to compute the matrix elements in terms of basis states in the |N 0 , N 1 , . . . , occupation notation. Therefore we wish to calculate the elements and then project into eigenstates of the center-of-mass HamiltonianĤ cm given in Eq. (52). It is therefore sufficient to diagonalizeÂ +Â− , [given by Eq. (51)] as this is the only operator dependence in H cm , using basis states of the form |N 0 , N 1 , . . . . This gives a square matrixP of eigenvectors of center-of-mass, which can project the truncated Fock state basis into this new basis, and a vector of eigenvalues. This is computationally simple asÂ +Â− cannot mix states of different energies and therefore is block diagonal when states are ordered by energy and each block can be diagonalized separately. By removing all the columns ofP with associated eigenvalues not equal to zero (meaning they have excitations in the center-of-mass mode) we are left with a rectangular matrixP which projects into this ground state of center-of-mass excitation subspace that we call the 'reduced basis'. UsingP results in a far smaller basis set (discussed in the next subsection) without changing any of the relative dynamics, however it is not immediately clear what states in this new basis relate to. Given that each partition ofẼ into N positive integers has the interpretation that each integer k represents a single occupancy in the kth mode, one may ask what the relation to quantum numbers is of partitions in terms of integers less than or equal to N, for instanceẼ = 2 can be partitioned by 1 + 1 and 2. Given that we know the ladder operator associated with the center-of-mass modeÂ ± , satisfies [Ĥ 0 ,Â ± ] = ±γÂ ± and is thus spaced in steps of unity times γ, we can associate all the 1's in a given partition with a quanta in this mode. Assuming we have quanta in the center-ofmass mode, this leaves all the numbers 2 ≤ z ≤ N as ways to partitionẼ − , which must then relate to some relative excitation modes. Going back toẼ = 2 the partition, 2 = 1 + 1 is two quanta in the center-of-mass mode i.e.Â +Â+ |N, 0, . . . and the partition 2 = 2 is one quanta in the first relative mode.
In order to help understand this we examine the N = 2 case in first quantization, using Jacobi coordinates [Eq. (15)]. The Hamiltonian can be expressed in two commuting parts For distinguishable atoms, these would each have normal harmonic oscillator eigenstates (up to a scaling factor), which can be multiplied together to create a many-body eigenstate. However, we require Bose symmetry of the many-body wavefunction: ψ(x 1 , x 2 ) = ψ(x 2 , x 1 ); in terms of Jacobi coordinates this implies no conditions on x C but that ψ(x C , ξ 2 ) = ψ(x C , −ξ 2 ) and hence odd eigenstates for H rel are disallowed and relative energy levels are spaced in units of 2.

Basis size reduction
As mentioned in Appendix F 3, the center-of-mass mode ladder operatorÂ ± of Eq. (51) has an energy spacing of unity, implying relative excitation modes are spaced in units of 2, 3, .., N.
In our reduced basis, the subset with the center-of-mass in the ground state, we can no longer partitionẼ using the number 1. Therefore the energy degeneracyg(Ẽ) of levelẼ in the reduced basis, is the number of ways to partitionẼ using integers z satisfying 2 ≤ z ≤ N, i.e.g(Ẽ, N) = p([2, N],Ẽ). Therefore the number of basis states in the reduced basis relative to the occupation number basis with cut ofẼ is given by .

(F8)
Essentially this property can be seen from projecting the set of kets with energy η into the center-of-mass excitation basis, this set will contain all the relative excited states with energy less than or equal to η, but with additional center-of-mass excitation. The basis reduction for N = 2 can be calculated by noting there are k/2 + 1 ways to partition k using 1 and 2 (the notation k means round k down to an integer), thus the reduced if η even 1 + η + (η 2 − 1)/4 if η odd. (F9) To leading order the reduction ∆(η, N) goes as 2/η. Such simple analytic expressions are not known for general N, however we have the following expression by Ramanujan [48] p([1, N ≥ η], η) ∼ 1 4η √ 3 exp this can be used to get an asymptotic estimate of the basis reduction by replacing the sum in Eq. (F8) with an integral, giving which will be our best estimate for the reduction achieved for large N, note that this improves slower than the ∝ 1/η reduction for the N = 2 case. This asymptotic estimate is included in Fig. 5, along with the reduction for intermediate values of N.