Coupled pair approach for strongly-interacting trapped fermionic atoms
We present a coupled pair approach for studying few-body physics in harmonically trapped ultracold gases. The method is applied to a two-component Fermi system of particles. A stochastically variational gaussian expansion method is applied, focusing on optimization of the two-body correlations present in the strongly interacting, or unitary, limit. The groundstate energy of the four-, six- and eight-body problem with equal spin populations is calculated with high accuracy and minimal computational effort. We also calculate the structural properties of these systems and discuss their implication for the many-body ultracold gas and other few-body calculations.
pacs:03.75.Ss, 05.30.Fk 31.15.ac, 34.50.-s
Ultracold atomic two-component Fermi gases under harmonic confinement have become an important field of study for fundamental quantum phenomena. The tunability of the interspecies -wave scattering length – the dominant interaction channel – makes these systems ideal for exploring the strongly interacting, or unitary, regime at the BEC-BCS crossover where the scattering length diverges and becomes the dominant length-scale in the system Shin et al. (2008); Chin et al. (2010); Braaten and Hammer (2006); Zwierlein et al. (2006); Blume (2012). Many-body calculations based on perturbative methods fail at unitarity due to the divergence of the scattering length Bloch et al. (2008); Giorgini et al. (2008). Alternative techniques involving Monte Carlo integration and effective interactions Gilbreth and Alhassid (2012); Blume et al. (2007); Mukherjee and Alhassid (2013) require accurate knowledge of the nodal surfaces to be used as references for antisymmetric wavefunctions for fermionic systems. Density functional theory requires an accurate energy functional for the study of many-body systems Bulgac and Yu (2002); Xianlong et al. (2006); Papenbrock (2005). Studies of few-body systems provide benchmarks for optimization and refinement of these many-body calculations and future experiments Werner and Castin (2006a, b); Kestner and Duan (2007); Daily and Blume (2010); Liu et al. (2009, 2010); Stöferle et al. (2006).
Few-body calculations can also be directly applied as atom traps become more sophisticated and offer the possibility of trapping only a few atoms in one trap Serwane et al. (2011) or a few atoms on each site of an optical lattice Bakr et al. (2009); Stöferle et al. (2006). Two-body correlations have been observed to play an important role in these systems Zürn et al. (2013). Extensions to more complex systems involving the application of external fields expands the known set of universal relations Mulkerin et al. (2012a, b).
The in-principle exact calculation of harmonically trapped few-body systems with zero-range -wave interactions has been greatly extended in recent years. The exact wavefunction and energy spectrum for two unlike atoms in a trap was found by Busch et al. Busch et al. (1998). Knowledge of the two-body system has spurred calculations of the three-body problem using the adiabatic hyperspherical method Werner and Castin (2006a); Kestner and Duan (2007); Liu et al. (2010). Exact diagonalization using the stochastic variation of a correlated gaussian basis has enabled calculation of energetics and structural properties of up to six trapped fermions for finite range interactions Daily and Blume (2010); Blume and Daily (2011). However, as the number of particles increases the Hilbert space grows exponentially and exact diagonalization becomes intractable. The challenge is to extend these calculations as far as possible beyond the two- and three- body problem towards the many-body regime.
In this work we present an improved methodology based on the stochastic variation of a gaussian basis Mitroy et al. (2013); Suzuki and Varga (1998); Hiyama et al. (2003) that allows calculation of energy levels and structural properties of a two-component Fermi system with up to atoms. The atomic hyperfine states which form the two components are treated as two arbitrary spin- states and the associated statistics only allows interactions between unlike particles. In this work we restrict ourselves to the case of even with equal population in each spin state, i.e. , as proof of concept of the approach. The groundstate of this system has zero total orbital angular momentum and spin, simplifying the calculation. The extension to more general cases is straightforward but requires greater computational effort.
The key idea of the coupled pair approach is to consider only the essential correlations. Interactions only occur between two unlike fermions and all other correlations are captured in the non-interacting correlations between dimers whose behavior is governed by the trap. This problem can be solved using variational methods but the optimization procedure is only applied at the two-body level simplifying the calculation for larger where exact diagonalization is normally intractable. With these considerations we use a gaussian expansion of the relative wavefunction and stochastic optimization to calculate the groundstate energy and structural properties of up to fermionic atoms in a harmonic trap.
Section II outlines the basic formalism for trapped fermions. In Section III we discuss the problem in detail and introduce the coordinate channels and their importance for the few-body problem. Section IV extends the method to and and outlines the details of the coupled pair approach. Then, in Section V we calculate the single-particle density and pair correlation functions and discuss their importance to the many-body system. Finally, in Section VI we summarize and discuss the extension of the approach to more general cases.
Ii -body problem and general approach
The Hamiltonian for harmonically trapped atoms of equal mass is
where is the position of particle , is the trapping frequency and is the interparticle interaction potential. The sum is restricted to interactions between unlike fermions and interactions between more than two particles are omitted. Here we only consider the equal mass case but the following discussion can be generalized to arbitrary masses. Firstly, the center-of-mass motion is decoupled from Eq. (1) and solved separately, whereupon we may assume that it is in the groundstate with energy and wavefunction
where is the center-of-mass coordinate and is the harmonic oscillator length.
All the interparticle interactions are then contained in the relative Hamiltonian and the relative wavefunction is an antisymmetrized function of relative coordinate vectors. The choice of the relative coordinates is not unique and can be chosen to take advantage of the symmetry of the system and reduce the complexity of the problem. Different sets of coordinates can represent different channels in which the particles are correlated, and we can include multiple channels in our ansatz for the relative wavefunction
where is an unsymmetrized wavefunction of the channel and the set of coordinates for each channel is expressed as a supervector . The operator projects these channel wavefunctions onto the correct antisymmetric space.
The problem of two trapped fermions has been solved analytically for the zero-ranged Fermi pseudo-potential and reproduces the cusp in the relative wavefunction Busch et al. (1998), where is the -wave scattering length. Unfortunately, these analytic solutions are unsuitable for use in the -body problem, due to the difficulty of treating the many cusps that occur and because the non-interacting harmonic oscillator states have poor convergence at unitarity. Instead we use a gaussian basis to expand the of each channel
where the are expansion coefficients, and the are the widths of the gaussians, which serve as variational parameters for the model. Each of the relative vectors represents a correlation within the problem and are treated independently. Hence, each basis term is separable in the vectors, as represented by the product in Eq. (4). In general, this product can have more terms corresponding to any number of correlations, but as explained below, we choose the same number of correlations as relative coordinates. The basis size is determined by the number of terms in the sums and can be expanded to improve the solution in accordance with the variational principle Mitroy et al. (2013). The groundstate of the equal spin component problem has zero relative angular momentum so our ansatz needs no explicit angular component. The gaussian basis functions are not orthogonal but they are simple to manipulate and are effective at replicating correlations at any length scale, including the short-ranged interparticle interactions, while also being efficient enough to scale to larger systems.
In the unitary limit the -wave scattering length diverges and becomes the only important length-scale associated with the interparticle interaction.The details of the short-ranged interparticle interaction potential are unimportant, provided it can support a single bound state. We chose a gaussian basis for its flexibility and ability to access all length scales so we also choose a gaussian potential
where, for any width there is a depth such that the potential supports a single resonant bound state and has a divergent scattering length corresponding to the unitary limit. In the limit the gaussian potential tends towards a regularized contact interaction and the bound state is well behaved. Moreover, it simplifies the calculation of matrix elements when using a gaussian basis and the appropriate values of and can be found with elementary scattering theory. Universal properties only emerge in the true zero-range limit but the values of considered here are sufficiently small for the properties of the system to be considered very close to the true groundstate, for . There are possible interacting pairs for the equal spin component system and we include all of them in the Hamiltonian (1).
The product of gaussians in Eq. (4) can be conveniently expressed with an symmetric matrix as
where the superscript ‘T’ denotes matrix transposition. The matrix elements of all terms in the Hamiltonian (1) with the gaussian potential Eq. (5) can be found from these correlation matrices Suzuki and Varga (1998).
By diagonalizing the relative Hamiltonian we not only obtain the energy spectrum but the relative wavefunction. Combining this with the center-of-mass wavefunction we obtain the total wavefunction for the -body problem. From we can calculate a general structural property
where (and ) is a coordinate describing the property of interest and is normalized to unity. Here is a general set of coordinates such as the center-of-mass plus relative coordinates as defined above or the single-particle coordinates. These quantities are related to the density matrices of the system and are calculated in a similar way Suzuki and Varga (1998); Blume and Daily (2011). In particular we calculate the single-particle reduced density , with in Eq. (7) and the (scaled) pair correlation function , with in Eq. (7). is the density of either spin species and is the probability of finding a pair of opposite spin fermions of size .
As seen by the choice of basis, our approach identifies interparticle correlations in the -body problem with the relative coordinates. These ‘channels’ serve the dual purpose of being a separable set of coordinates in which to perform the integration of the Schrödinger equation while also directly corresponding to the most important correlations in the problem. Other correlations can be included as off-diagonal elements in . A common approach is to include all two-particle correlations Daily and Blume (2010); Rittenhouse et al. (2011); Mitroy et al. (2013), but having more variational parameters quickly becomes intractable for . Furthermore, correlations between like fermions that cannot interact are much less significant at unitarity, and all the important correlations are naturally included in the choice of coordinates. By only including parameters from each channel the calculation of individual matrix elements is more efficient, since the are diagonal, and the basis we use focuses more on the most important two-particle correlations. Therefore, the choice of coordinates is the key to not only encapsulating the important properties of the -body problem but to make it tractable as increases. We first illustrate these ideas in the case of .
Iii Four-body problem:
iii.1 Coordinate channels
For the four-body problem the relative coordinates can be constructed in two ways, K-type and H-type. K-type coordinates are constructed by iteratively defining the relative vector between the center-of-mass of a subgroup of particles and one extra particle; they are the canonical Jacobi coordinates. Physically, in the four-body problem the K-type channels represent the correlations between a pair and two free particles. H-type coordinates begin by defining two pairs and then the relative vector between the centers-of-mass of the pairs. These channels represent the correlations within two interacting dimers, and then the dimer-dimer correlation.
There are ways to construct any coordinate channel depending on the order in which the particles are correlated and so the four-body problem has 48 possible coordinate channels Hiyama et al. (2003). However, we do not need to include all channels since the symmetry of the problem will make many of them redundant. The antisymmetrizing operator for two spin-up and two spin-down particles is , where permutes the and particles. Throughout this work we have adopted the notation that odd and even indices label the two different spin components. Under the action of the permutation operators in and considering the interaction terms included in the Hamiltonian Eq. (1) the 48 K-type and H-type channels in Eq. (3) for the four-body problem reduce to only three linearly independent channels, shown in Fig. 1 and given explicitly in Table 1. These channels contain all the correlations required for the four-body problem. The reduced mass is the same for all relative coordinates and preserves the volume element in the transformation from single-particle coordinates. The antisymmetrizer also has the effect of ensuring that all possible interacting pairs of particles produce a cusp in the wavefunction, without having to explicitly include them all in the basis.
iii.2 Groundstate energy for
For excited states or weaker interactions the correlations of free particles is important and the contributions from the K-type channels must be included. However, for the groundstate at unitarity we intuitively expect that the H-type channel is sufficient. In each channel the basis is chosen using the stochastic variational method whereby a set of gaussian widths is chosen semi-stochastically and the relative Hamiltonian is constructed and diagonalized. This process is iterated until the lowest eigenenergy converges. These calculations are repeated for to and show a linear trend that is extrapolated to the limit . Using all three channels from Fig. 1 we achieve a groundstate energy of , very consistent with previous calculations Daily and Blume (2010); Rittenhouse et al. (2011). The uncertainty is in the last digit. With this calculation as a benchmark, we can compare to a calculation using only the H-type channel, for which we obtain , or a difference of from the multiple channel calculation. These calculations are summarized in Table 2.
The H-type channel by itself gives a very good approximation to the true groundstate energy, but as a further test we examine the structural properties of the system at unitarity. In Fig. 2 we plot the scaled pair correlation function of the four-body groundstate at unitarity using different combinations of channels. The calculation is performed using only the H-type channel (blue, solid), only the K-type channels (red, dashed), and all three (black, dotted). As shown by the inset, the calculation using only the H-type channel is very close to the calculation using all three channels. The minimal effect of the adding the K-type channels to the H-type channel is to make the pair sizes slightly smaller, reflecting higher-order pair correlations. The K-type channels by themselves are very different from the full calculation, especially at small . This is due to the fact that they do not allow for as many pairs of unlike fermions with a small separation.
These results confirm that the groundstate at unitarity is well-represented by only the H-type channel. That is, the most important correlations in the system are those between two particles interacting via an -wave contact interaction and then between the dimers whereas correlations involving free particles contribute only a small perturbation. Interactions involving more than two particles are also negligible. Therefore, when we consider scaling to larger in the next section, we seek to discard the insignificant terms and devote computational power to the interacting pairs.
Iv Extension to higher
iv.1 Coordinate channels and optimization
In the four-body problem, the different correlations contained in the K-type and H-type channels leads to a clear distinction between the results. For larger there arise more possible sets or ‘shapes’ of coordinate channels, including hybrids of the generalized K-type and H-type channels introduced earlier and each set allows permutations before considering the antisymmetry of the wavefunction. The different coordinate channels can be characterized by the types of correlations they naturally represent and any channel can be viewed as a set of subsystems. We identify two types of correlations, interacting-pair correlations (IPCs) between two fermions in different spin states and non-interacting correlations (NICs) involving more than two particles. The latter includes correlations between dimers, and correlations between a subcluster (two or more particles) and a single particle. The distinction is important because the IPCs in the relative wavefunction must reproduce the cusp between any two interacting fermions. Although correlations involving more than two particles can be treated via effective interactions Petrov et al. (2004), cusps in the NICs are either suppressed by Fermi statistics or are much weaker than the two-body interaction and the correlations are on a larger length scale. These higher-order effects have not been considered explicitly in this work but their effects can be considered to be incorporated in the NICs which are governed primarily by correlations of order . For higher H-type refers only to those channels which have the maximum IPCs. That is, they are constructed by first pairing up all particles then building the correlations between pairs. All other channels are referred to as (generalized) K-type, even if they contain several IPCs.
The antisymmetrizer for general is , where is the symmetric group. This amounts to permutation terms with associated minus signs for each exchange of two identical particles. By applying this operator we find that there is only one H-type channel for and two H-type channels for ; these are shown in Fig. 3. For , Fig. 3(a) shows the H-type channel representing the correlation between a tetramer and an additional pair. Figures 3(b) and 3(c) show, respectively, the H-type channels with either the correlation between two H-type tetramers, or the correlation between a six-body subcluster and an additional pair. All three include the internal correlations of the smaller subclusters. Below we show that these channels are sufficient for calculating the groundstate of the six- and eight-fermion problem to high accuracy.
Even with only variational parameters in one or two channels it is still not practical with available computational resources to variationally optimize the entire problem when . However, another advantage to identifying coordinates with correlations is that the basis does not require optimization with respect to the -body problem. Instead, each NIC and IPC is optimized as an independent subsystem of reduced mass . The energy spectrum without interactions and at unitarity is known exactly for Busch et al. (1998), and it is simple to optimize these pairs with a gaussian basis and a gaussian potential, Eq. (5). This not only makes the optimization procedure extremely fast, but only a very small number of basis states is needed in each correlation to reproduce the first few energy levels to high accuracy. This allows us to use a very small basis focused on the most important length scales, and , and extend the approach to larger . Specifically, the NICs are only of order the trap size so require a smaller basis than the IPCs which also need to access the interparticle potential, i.e. contain terms with . One drawback to using such a small basis in the subsystems is that it limits how small we can choose and thus how accurately we can extrapolate . In practice however, the formal requirement that is satisfied by with only a small error compared to the true zero-range limit.
The method is variational even though it does not seek to find convergence of the solution at the -body level. The flexibility and accuracy of this method comes from including all the important correlations in the choice of coordinates. We still use as large a basis size as possible but the terms are distributed more efficiently. The gaussian basis functions are not orthogonal so care must be taken to avoid linear dependence between the differnet subsystems. However, because the optimization is performed first at the level, we only need to construct and diagonalize the -body Hamiltonian once, rather than at each step of a full optimization procedure. This leads to a large reduction in computational time especially as the number of permutations grows with .
iv.2 Groundstate energy for and
Similar to the calculation presented in Section III for , for we first perform a large calculation to compare against results involving a restricted basis size and choice of channels. Initially, we consider five channels, including several K-type channels, to better incorporate single-particle excitations and the number of basis functions in each correlation is maximized up to the limit of computational resources while maintaining sufficient linear independence. Although still large, the basis size for the ICPs and NICs for is smaller than for . Hence, we consider larger values of : . However, we find that the linear extrapolation as is still valid. In the zero-range limit we obtain a groundstate energy of , which is lower than the calculation using all two-particle correlations 111Blume and Daily Blume and Daily (2011) report a lower figure of by also extrapolating the basis size to infinity, however this is not strictly in accordance with the variational principle. and so is a lower upper bound on the true groundstate energy. This is principally because although the basis sizes are comparable we have distributed the basis among only the most important correlations so it accesses a large part of the relevant Hilbert space. Secondly, we repeat the calculation using only the H-type channel and obtain in the zero-range limit. This differs from the larger calculation by . This demonstrates that the correlations included in the H-type channels are dominant in the groundstate of the system. Details of these results are also summarized in Table 2.
These results demonstrate that the restriction to only H-type channels using a basis of optimized pairs is valid, so we now extend this method to . There are two linearly independent H-type channels [see Fig. 3(b,c)], which must both be included to access all the two-particle IPCs, and there is a much larger number of permutations. The calculations are performed for larger : , and the extrapolation to gives . This result is a lower upper bound to the groundstate energy from Monte Carlo calculations, Blume et al. (2007), or an effective interaction method, Mukherjee and Alhassid (2013).
Although the total basis size for is maximized, the large number of correlations in the problem means that, in comparison to and , we can only use a small number of gaussians in each subsystem due to limited computational resources. This also means we cannot compare to a calculation that includes K-type channels. The basis is ‘reduced’ at the level in that it is restricted to the bare minimum needed to reproduce the cusp behavior of the IPCs and the NICs are taken to be of similar size to the trap length, i.e. . The ‘reduced’ basis also requires larger values of in order to maintain similar accuracy in the optimization of the IPCs and NICs. Previous calculations for and used a ‘full’ basis size, which essentially allowed arbitrary accuracy in the two-body problems used for optimization. To justify using a ‘reduced’ basis size we repeated the calculations for and using the same number of gaussians in each IPC and NIC of the H-type channel as was used in the calculation. The results and details of all groundstate calculations are summarized in Table 2. Using the reduced basis size for and gives groundstate energies, in the zero-range limit, of and , respectively. The results are within and , respectively, of the ‘full’ basis results. We expect that the error in in the zero-range limit is of similar order.
V Structural properties
The coupled pair approach enables the calculation of the structural properties for up to . In Fig. 4(a) we plot the normalized single-spin species reduced density at unitarity, for using only the H-type channels with the largest possible basis and smallest possible for each case. We see the flattening of the small- density and emergence of a small peak at non-zero for the six and eight-body, which may indicate formation of shell structure in the system. The more prominent peak in Ref. Blume and Daily (2011) may be due to the inclusion of single-particle correlations, but we note that it was performed at a larger value of .
In Fig. 4(b) we plot the (scaled) pair correlation function at unitarity for the same parameters and basis sizes as Fig. 4(a). As increases the main peak at is enhanced and becomes the main feature. This is expected as the formation of large weakly-bound dimers due to unitary -wave scattering becomes more likely. The value of in the limit is related to the number of pairs of fermions of opposite spin separated by a small distance , thereby allowing them to interact Tan (2008); Werner and Castin (2012); Daily and Blume (2009). In the thermodynamic limit this provides a strong link between the physics of few- and many-body systems but for the values of considered here the shell structure plays a larger role.
The curves given in Fig. 4 are obtained using the best calculation possible for each in terms of basis size and the value of . In particular, for the value of is significantly larger than for and although it is still considered small. The near vertical lines in Fig. 4(b) are due to the small finite range effect of non-zero . Decreasing produces a small change in the peak values of and , but this increase is bounded in the limit. Increasing the basis size has the same effect. This effect can be seen by calculating the structural properties of the system for and for a range of and for the ‘reduced’ basis as determined by the system.
We have accurately calculated the groundstate energy and structural properties of a two-component few-fermion system at unitarity with an even number of particles up to . The coupled pair approach minimizes computational effort by retaining only the most important two-body correlations in the unitary limit, allowing the extension to a higher number of particles. For we present a lower upper bound on the groundstate of . The two-particle correlation function has two peaks at and , representing the number of particles coming within the range of the interparticle interaction and the size of the weakly bound dimers, respectively. As expected for a Fermi system at unitarity, the second peak increases with due to the dominance of the strong -wave interactions.
All calculations were performed on desktop computers where the principal limiting factor on the computational time for increasing is the number of permutations required to antisymmetrize the wavefunction. The size of the basis used for each two-body correlation is limited by available memory but even for , in which at most five gaussians were used for each correlation, the zero-range ground state energy obtained is lower than other techniques. Given these computational limitations, we found that the extrapolation to the zero-range limit was more significant in calculating a lower groundstate energy. Performing the calculation for as small a value of as possible, while still maintaining high accuracy in the optimization of the IPCs, led to a better linear fit and lower zero-range groundstate energy than increasing the basis size at larger .
The ground state of the -fermion problem with equal spin components has zero total angular momentum allowing us to ignore the angular component of the wavefunction in the ansatz of Eq. (6). The coupled pair approach can be extended to excited states with different symmetry and even other systems with unequal spin populations, odd or unequal masses. These problems would require the addition of K-type channels and explicit angular functions since single-particle excitations become more important but otherwise the same principles described in this work apply. Advances in the trapping of atomic gases in confined dimensions or with different interparticle interactions Bloch et al. (2008) open up many other avenues in which gaussian expansion and coupled pair methods may apply.
In the unitary limit perturbative many-body techniques fail and new techniques are required. The coupled pair approach is a more efficient method of calculating energetic and structural properties of ultracold Fermi systems with a few atoms. By shifting the computational problem to the most important subsystems of the problem we have highlighted the significance of two-body-correlations as well as pushed the calculation to higher . These results can provide more accurate benchmarks for experiments and calculations in ultracold few-body physics to bridge the gap to the many-body gas.
Acknowledgements.H.M.Q. acknowledges the support of the ARC Centre of Excellence for Coherent X-ray Science and the ARC Centre of Excellence for Advanced Molecular Imaging.
- Shin et al. (2008) Y.-i. Shin, C. H. Schunck, A. Schirotzek, and W. Ketterle, Nature 451, 689 (2008).
- Chin et al. (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
- Braaten and Hammer (2006) E. Braaten and H.-W. Hammer, Physics Reports 428, 259 (2006).
- Zwierlein et al. (2006) M. W. Zwierlein, C. H. Schunck, A. Schirotzek, and W. Ketterle, Nature 442, 54 (2006).
- Blume (2012) D. Blume, Rep. Prog. Phys. 75, 046401 (2012).
- Bloch et al. (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
- Giorgini et al. (2008) S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
- Gilbreth and Alhassid (2012) C. N. Gilbreth and Y. Alhassid, Phys. Rev. A 85, 033621 (2012).
- Blume et al. (2007) D. Blume, J. von Stecher, and C. H. Greene, Phys. Rev. Lett. 99, 233201 (2007).
- Mukherjee and Alhassid (2013) A. Mukherjee and Y. Alhassid, Phys. Rev. A 88, 053622 (2013).
- Bulgac and Yu (2002) A. Bulgac and Y. Yu, Phys. Rev. Lett. 88, 042504 (2002).
- Xianlong et al. (2006) G. Xianlong, M. Polini, R. Asgari, and M. P. Tosi, Phys. Rev. A 73, 033609 (2006).
- Papenbrock (2005) T. Papenbrock, Phys. Rev. A 72, 041603 (2005).
- Werner and Castin (2006a) F. Werner and Y. Castin, Phys. Rev. Lett. 97, 150401 (2006a).
- Werner and Castin (2006b) F. Werner and Y. Castin, Phys. Rev. A 74, 053604 (2006b).
- Kestner and Duan (2007) J. P. Kestner and L.-M. Duan, Phys. Rev. A 76, 033611 (2007).
- Daily and Blume (2010) K. M. Daily and D. Blume, Phys. Rev. A 81, 053615 (2010).
- Liu et al. (2009) X.-J. Liu, H. Hu, and P. D. Drummond, Phys. Rev. Lett. 102, 160401 (2009).
- Liu et al. (2010) X.-J. Liu, H. Hu, and P. D. Drummond, Phys. Rev. A 82, 023619 (2010).
- Stöferle et al. (2006) T. Stöferle, H. Moritz, K. Günter, M. Köhl, and T. Esslinger, Phys. Rev. Lett. 96, 030401 (2006).
- Serwane et al. (2011) F. Serwane, G. Zürn, T. Lompe, T. B. Ottenstein, A. N. Wenz, and S. Jochim, Science 332, 336 (2011).
- Bakr et al. (2009) W. S. Bakr, J. I. Gillen, A. Peng, S. Folling, and M. Greiner, Nature 462, 74 (2009).
- Zürn et al. (2013) G. Zürn, A. N. Wenz, S. Murmann, A. Bergschneider, T. Lompe, and S. Jochim, Phys. Rev. Lett. 111, 175302 (2013).
- Mulkerin et al. (2012a) B. C. Mulkerin, C. J. Bradly, H. M. Quiney, and A. M. Martin, Phys. Rev. A 85, 053636 (2012a).
- Mulkerin et al. (2012b) B. C. Mulkerin, C. J. Bradly, H. M. Quiney, and A. M. Martin, Phys. Rev. A 86, 053631 (2012b).
- Busch et al. (1998) T. Busch, B.-G. Englert, K. Rzazewski, and M. Wilkens, Found. Phys. 28, 549 (1998).
- Blume and Daily (2011) D. Blume and K. M. Daily, C. R. Phys. 12, 86 (2011).
- Mitroy et al. (2013) J. Mitroy, S. Bubin, W. Horiuchi, Y. Suzuki, L. Adamowicz, W. Cencek, K. Szalewicz, J. Komasa, D. Blume, and K. Varga, Rev. Mod. Phys. 85, 693 (2013).
- Suzuki and Varga (1998) Y. Suzuki and K. Varga, Stochastic Variational Approach to Quantum-Mechanical Few-Body Problems (Springer-Verlag, 1998).
- Hiyama et al. (2003) E. Hiyama, Y. Kino, and M. Kamimura, Progress in Particle and Nuclear Physics 51, 223 (2003).
- Rittenhouse et al. (2011) S. T. Rittenhouse, J. von Stecher, J. P. D’Incao, N. P. Mehta, and C. H. Greene, J. Phys. B: At., Mol. Opt. Phys. 44, 172001 (2011).
- Petrov et al. (2004) D. S. Petrov, C. Salomon, and G. V. Shlyapnikov, Phys. Rev. Lett. 93, 090404 (2004).
- (33) Blume and Daily Blume and Daily (2011) report a lower figure of by also extrapolating the basis size to infinity, however this is not strictly in accordance with the variational principle.
- Tan (2008) S. Tan, Ann. Phys. 323, 29712986 (2008).
- Werner and Castin (2012) F. Werner and Y. Castin, Phys. Rev. A 86, 013626 (2012).
- Daily and Blume (2009) K. M. Daily and D. Blume, Phys. Rev. A 80, 053626 (2009).