An interpolatory ansatz captures the physics of one-dimensional confined Fermi systems

Interacting one-dimensional quantum systems play a pivotal role in physics. Exact solutions can be obtained for the homogeneous case using the Bethe ansatz and bosonisation techniques. However, these approaches are not applicable when external confinement is present. Recent theoretical advances beyond the Bethe ansatz and bosonisation allow us to predict the behaviour of one-dimensional confined systems with strong short-range interactions, and new experiments with cold atomic Fermi gases have already confirmed these theories. Here we demonstrate that a simple linear combination of the strongly interacting solution with the well-known solution in the limit of vanishing interactions provides a simple and accurate description of the system for all values of the interaction strength. This indicates that one can indeed capture the physics of confined one-dimensional systems by knowledge of the limits using wave functions that are much easier to handle than the output of typical numerical approaches. We demonstrate our scheme for experimentally relevant systems with up to six particles. Moreover, we show that our method works also in the case of mixed systems of particles with different masses. This is an important feature because these systems are known to be non-integrable and thus not solvable by the Bethe ansatz technique.

compared to entirely numerical approaches, which typically have the solution represented on a large basis set with many non-zero contributions.
To describe the main idea we will focus on two-component Fermi systems of N ↑ particles with spin projection up and N ↓ particles with spin projection down. The Hamiltonian for the N = N ↑ + N ↓ system is where m is the mass and the sums run over the number of particles, i = 1, … , N, x i is the coordinate, and p i the momentum of the ith particle. The external confinement, V ext (x), is assumed to be the same for all particles. Furthermore, we will assume that V ext (x) is parity invariant and has at least N bound states. The interaction strength g is positive for repulsive interactions and negative for attractive interactions. We note that since the parity operator commutes with the Hamiltonian, parity is conserved when g changes. While we only discuss the two-component Fermi system below, the case of bosons or Bose-Fermi mixtures is similar in spirit and in formalism and we return briefly to this extension in the outlook. For the majority of the discussion we assume all particles have equal mass, m, but we also discuss the important extension to mass-imbalanced systems. Notice that in the case of two-component Fermi systems, particles of the same spin projection will not interact due to the antisymmetry required by the Pauli principle which implies that the zero-range interaction will vanish for similar spin projections. For this reason we may use the general form of the Hamiltonian in Eq. (1) with the sum over i < j also for fermionic systems.
Consider now an N-body system and assume that we have solved the problem for g = 0 with energy eigenstate |γ 0 〉 and for 1/g = 0 with eigenstate |γ ∞ 〉 . Now form the linear combination γ α γ α γ = + .
∞ ∞ (2) 0 0 This interpolatory ansatz is motivated by the intuitive idea that the wave function of a system with intermediate-strength interactions contains a mixture of qualities from the wave functions with weak and strong interactions. We may compute 〈 γ|H|γ〉 as a function of α 0 and α ∞ and look for optimum values. As we shall demonstrate in this paper, |γ〉 provides a simple yet accurate description of the system for any value of g. The only exception is deeply bound states to which we return below.
We will demonstrate that the interpolatory ansatz in Eq. (2) can capture the qualitative features of the eigenstates of the Hamiltonian in Eq. (1) and is quantitatively accurate at the level of a few percent. Furthermore, we show that the expression for the optimum energy of |γ〉 may be modified slightly to make it perturbatively correct in both limits of the interaction strength. With this modification, the interpolatory ansatz provides an approximation to the eigenenergy with an accuracy that is comparable to state-of-the-art numerical methods, though the ansatz is far simpler.
We provide a proof of principle by considering some important examples from the few-body limit that are experimentally relevant at the moment. For this purpose we restrict to a harmonic confinement, i.e. ω = V x m x ( ) ext 1 2 2 2 , throughout (here ω is the angular frequency of the oscillator).
The results we present here are: • Analytical expressions for the ansatz parameters α 0 and α ∞ that depend exclusively on two matrix elements in |γ 0 〉 and |γ ∞ 〉 (as well as the eigenenergies of these limiting states). • The case N = 2 where one has the exact solution available 31 . Here we show that our method is accurate for both repulsive and attractive interactions to less than a few percent. • The N = 3 case where no analytical solutions are known for general g. Here our ansatz provides very accurate results for all g. The results are compared to exact numerical diagonalisation utilising a unitary transformation of the interaction Hamiltonian. • Energies for the impurity limit with N ↓ = 1 and N ↑ = 1-5 which we compare to experiments and find excellent agreement. We also discuss the Anderson orthogonality catastrophe for this system, which is related to the coefficient α 0 . • An extension of the method to systems with particles of different mass. The examples we discuss are threebody systems and we compare to exact numerical results based on the correlated Gaussian method. This is the first application of the correlated Gaussian method to mass-imbalanced systems in 1D that we are aware of.
The ansatz we propose can be used to get very simple expressions for different observables as one needs to compute only a few matrix elements between the g = 0 and 1/g = 0 states of interest. Our method is directly extendable to bosonic systems or mixed systems, as long as one has access to the two limiting wave functions, and while we have focused on harmonically confined systems it can be straightforwardly extended to any other form of confinement. Furthermore, one may systematically improve the ansatz by adding more states.

Results and Discussion
Let |γ 0 〉 be an energy eigenstate in the non-interaction limit (that is, g = 0) with eigenenergy E 0 . We now adiabatically change the interaction strength from g = 0 to |g| → ∞ and in turn we adiabatically change |γ 0 〉 into a new state denoted |γ ∞ 〉 with energy E ∞ . The wave function of the state |γ ∞ 〉 vanishes whenever the position coordinates of any two particles coincide, and |γ ∞ 〉 is thus unaffected by the interaction potential Scientific RepoRts | 6:28362 | DOI: 10.1038/srep28362 As our ansatz we now construct the trial state where α 0 and α ∞ are real parameters. Assuming |γ 0 〉 and |γ ∞ 〉 are normalised, the energy of the trial state is We use a variational approach to solving the Schrödinger equation by identifying stationary points of the trial state energy functional Eq. (5). We thus select the values of α 0 and α ∞ that optimise the energy of the trial state for a given value of g. This will yield energies and eigenstates that, although approximate, turn out to be extremely accurate as discussed below.
Before presenting the results of the variational calculation, we shall briefly examine Eq. (5) in the limiting cases of the interaction strength: If we require that α ∞ /α 0 → 0 for g → 0 such that the trial state approaches |γ 0 〉 in the non-interacting limit, Eq. (5) gives a first-order term in g of the form 〈 γ 0 |V|γ 0 〉 , in agreement with first order perturbation theory. We caution, however, that requiring α 0 /α ∞ → 0 (i.e. |γ〉 → |γ ∞ 〉 ) for 1/g → 0, does not automatically ensure that the first order expansion of Eq. (5) in 1/g is equal to that of the exact eigenstate. We will come back to this point later on.
The coefficients that yield stationary points of Eq. (5) are given by (see the Supplementary Materials for details) This gives the energy The denominator in Eq. (7) is positive because |〈 γ 0 |γ ∞ 〉 | ≤ 1 by the Cauchy-Schwarz inequality. Hence, + E opt ( ) is the energy maximum while − E opt ( ) is the energy minimum. To ensure the correct energies in the limits g = 0 and 1/g = 0, we use the energy minimum, − E opt ( ) , to approximate the eigenenergy whenever g > 0, while for g < 0 we use the energy maximum, + E opt ( ) . While the interpolatory ansatz is extremely simple, it has a shortcoming in the 1/g → 0 limit where it does not reproduce the slope of the energy. As is shown in the following, we may, however, modify the ansatz slightly to correct this behaviour. Letting q ≡ − 1/g, the first-order expansion where K 0 = 〈 γ 0 |V|γ 0 〉 /g is the corresponding slope of the energy curve in the limit of vanishing interactions. This demonstrates that the important quantity in the slope is the overlap 〈 γ 0 |γ ∞ 〉 . In the original philosophy of the ansatz, we exploit that we know both wave functions in this overlap exactly and thus also the overlap itself; this leaves no unfixed parameters. Realising that this yields a discrepancy we have explored how to modify this assumption in order to improve the approximation.
To this end, we note that the derivation of Eq. (7) actually does not depend on |γ ∞ 〉 being an energy eigenstate. It must be a state with energy E ∞ (with respect to the non-interacting Hamiltonian), but beyond that the only requirement of |γ ∞ 〉 is that V|γ ∞ 〉 = 0. Furthermore, E opt only depends on |γ ∞ 〉 through the squared wave-function overlap 〈 γ 0 |γ ∞ 〉 2 . Thus, if we substitute 〈 γ 0 |γ ∞ 〉 2 in Eqs (7 and 8) with some parameter λ, we may regard E opt and ∞ K opt as functions of λ. We can then select λ such that E opt (λ) becomes perturbatively correct, that is λ = where ∞ K exact is the slope of the true eigenenergy curve at q = 0 (which is known exactly using the formalism of A.G.V. et al. 22 ). We shall refer to this perturbatively correct modification, E opt (λ), as the modified ansatz. Note that this modification breaks variational bounds since we have no a priori knowledge of any trial state whose energy is E opt (λ). By going from the original interpolatory ansatz to the modified ansatz, we lose information about the wave function, but gain the correct slope of the energy at 1/g → 0. Below we will see that the modified ansatz increases the accuracy significantly as compared to the original ansatz.
We now proceed to discuss examples of the ansatz. Throughout our discussion, the external potential is taken to be harmonic ω = V x mx ( ) ext 1 2 2 2 and the same for all particles. This is the most widely studied and experimen-tally relevant case at the moment so it will be our focus here. Henceforth, we use natural units  = m = ω = 1 such that energies are given in units of ω, lengths in units of , and interaction strengths g in units of σ ω. Two particles. The N = 2 case is special as analytical results for any g are available due to the seminal work of Busch et al. 31 . It is therefore an important benchmark case for our approach. First we note that we will only be interested in the relative energy as the center of mass decouples in the harmonic trap. Note that this decoupling is not an essential assumption of our method and is merely a convenience. The details on how to construct the ansatz states for the two-particle case are given in the Methods section below. The energy spectrum using the interpolatory ansatz of Eq. (7) is shown in Fig. 1. Only the even-parity solutions are shown on the plot as odd-parity states are unchanged by the zero-range interaction. The figure also includes experimental measurements of the ground-state energy. The experiment has been conducted by S. Jochim's ultra-cold atoms group in Heidelberg, and the data originates from Wenz et al. 19 . This data has been corrected for the imperfections of the trap as described in Wenz et al. 19 . As seen, the experimental data agrees with the interpolatory ansatz well within the experimental uncertainties.
If we expand the energy in terms of g, the first-order term agrees with the result of ordinary non-degenerate perturbation theory. Similarly, in the limit 1/g = 0, the energy curve (as a function of q) of the true eigenstate has the same slope as that of our interpolatory ansatz.
The inset in Fig. 1 shows a zoom of the energy spectrum on the ground state and compares the energy predicted by Eq. (7) with the exact solution. For q > 0 the energy of the interpolatory ansatz is within 0.05 of the exact energy in the vicinity of q ~ 0.4 and even less elsewhere. The error decreases as we move up in the spectrum. For q < 0 the deviation is less than 0.006 for the ground state; again greatest around q ~ − 0.4. On this side, the error increases for the excited states, but we find that it is bounded by about 0.03.
We conclude that the ansatz is extremely accurate for the N = 2 case where we can compare to analytical results 31 . One may also compare the wave functions and again find extremely good agreement (see Supplementary  Fig. S1).
For attractive interactions, a deeply bound molecular state exists that we have so far ignored. However, it turns out that the ansatz of Eq. (7) can be extended to also give extremely accurate results for the deeply bound state. As is shown in the Supplementary Materials, this can be done by including an additional state in the ansatz that has the correct asymptotic behaviour as g → − ∞ . We stress that this is in fact in complete agreement with the universal philosophy of the ansatz method, i.e. interpolation between (known) extremes. Thus to address deeply bound states one needs a state in the extreme limit of large negative energy. This yields a very precise approximation also for the deeply bound state. This highlights the universal nature of our approach. We will not pursue the deeply bound states any further in this paper.
Three particles. The simplest non-trivial example of a three-body two-component Fermi system has N ↑ = 2 and N ↓ = 1. The interaction potential is Energy spectrum for the two-particle system. Relative energy of two distinguishable fermions according to the interpolatory ansatz, showing the low-energy part of the spectrum. Only states with even parity are plotted as odd-parity states are not influenced by the zero-range interaction. The system also has a deeply bound state, but this is not plotted. Squares are experimental data points (note that the error bars are smaller than the data points) 19 . The inset is a zoom of the ground-state energy comparing the interpolatory ansatz (solid) with the exact eigenenergy (dashed) for the ground state on the repulsive side as it crosses to the attractive side.
Scientific RepoRts | 6:28362 | DOI: 10.1038/srep28362 Again the third interaction term will vanish for identical fermions but we keep it for generality. Eigenstates of the harmonic Hamiltonian are described by two quantum numbers, ν ≥ 0 and μ ≥ 1, which we shall call the radial quantum number and the angular quantum number, respectively. The eigenenergies are The quantum numbers ν and μ can be used to describe the energy eigenstates in both the non-interacting limit and the infinite-interaction limit.
Constructing the ansatz. Recall that |γ 0 〉 and |γ ∞ 〉 denote eigenstates in the non-interacting limit and the infinite-interaction limit, respectively, and that we are looking for states that are adiabatically connected. As states with different radial quantum number ν are orthogonal, we assume that the adiabatically connected states have the same ν values. Parity p is exactly conserved and thus also the same for two states in the ansatz. The quantity that changes with g is therefore the μ quantum number, and we call the limiting values μ 0 and μ ∞ , respectively. The angular quantum numbers are related by For repulsive interactions, μ 0 < μ ∞ , and for attractive interactions, μ 0 > μ ∞ , cf. Eq. (11). The optimum energy of the trial state can now be determined using Eq. (7) with E 0 and E ∞ given by Eq. (11).
Energy spectrum. Figure 2 shows the energy spectrum for E < 8 (deeply bound states are not considered) both as calculated using the interpolatory ansatz and by exact numerical diagonalisation. The interpolatory ansatz clearly describes the qualitative features of the spectrum well; both for repulsive and attractive interactions. In this energy-range, the ν = 0 states also give a good quantitative match to the numerical results. However, the error seems to grow with ν.
We see in Fig. 2 that trial states with ν = 0 and Δ μ = 1 offer a particularly good approximation to the corresponding eigenstates in the repulsive region (see e.g., the first excited state). This is because the slope of the energy curve, that is = ∂ ∂ , at q = 0 is the same as that of the exact eigenenergy. For the trial states with Δ μ = 2, this is generally not the case. This result suggests that ∞ K opt is an important quantity to reproduce correctly in an attempt such as the present to describe the physics of our problem through a simple ansatz. Note also, that the slope of the energy curve has a discontinuity at 1/g = 0, which contradicts the expectation that the states go smoothly through this region.
If we now enforce the correct slope by a modification of the interpolatory ansatz as proposed in the discussion following Eq. (8), we arrive at the spectrum shown in Fig. 3. We see that the modified ansatz agrees better with the Figure 2. Energy spectrum of the interpolatory ansatz for the three-particle system. Red curves are energies of trial states with odd parity, and blue curves are those with even parity. Solid curves represent states with ν = 0, dashed curves ν = 1, and dotted curves ν = 2. Circles are numerical calculations. Deeply bound states are excluded from the plot. Squares are experimental data-points 19 . The error bars on the experimental data points are smaller than the squares and are therefore not shown.
numerical results than the original ansatz; especially for states with ν > 0. There are, however, still deviations on the attractive side of the spectrum, and for high energies, also on the repulsive side. We shall give a quantitative discussion of the quality of the approximation when discussing the impurity system below. We have included the experimental measurements of the ground-state energy 19 in the figure and see that both ansatz and modified ansatz agrees very nicely with experiment, although for large g the modified ansatz naturally does better.
Mass-imbalanced systems. In the case of different masses, one typically uses another length scale given by In Fig. 4 we show the energy spectrum obtained by the interpolatory ansatz for M/m = 1/2 and M/m = 2 and compare this with numerically calculated results using the correlated Gaussian approach (see the Supplementary Materials for details on the numerical methods). The agreement between the ansatz and the numerical results is striking for the low-energy part of the spectrum considered here, and we see that the ansatz can be extended also to mass-imbalanced systems. . Energy spectrum for the three-particle system-modified ansatz. Energy spectrum for the threeparticle system as in Fig. 2, but with the interpolatory ansatz modified to be perturbatively correct. Impurity systems. We now consider a system of N fermions among which one particle (x 1 ) is spin-down and the N ↑ = N − 1 remaining particles (x 2 , … , x N ) are spin-up. Taking into account that interactions between identical fermions vanish, we can write the interaction term as We restrict the discussion to the ground state with repulsive interactions as this has been a focus of recent experimental attention 19 . The considerations can be extended to obtain more states in the spectrum, to the attractive side, and to deeply bound states in the same manner as in the previous examples.
We denote by |γ A 〉 the Slater determinant of the single-particle harmonic eigenstates 〈 x k |n〉 = ψ n (x k ) for k = 1, … , N and n = 0, … , N − 1. Here ψ n (x) is the single-particle eigenstate of the harmonic oscillator Hamiltonian in coordinate space with quantum number n and argument x. The state |γ A 〉 is antisymmetric with respect to interchange of any two coordinates, so 〈 γ A |V|γ A 〉 = 0. We can write the antisymmetric state as , S N is the symmetric group of order N, and Π (σ ) indicates the integration region x σ(1) < … < x σ(N) . In each region, Π (σ), the wave function of the ground state in the infinite-interaction limit is proportional to 〈 x|γ A 〉 . The (normalised) ground state in the infinite-interaction limit is then A for some coefficients a = (a 1 , … , a N ) to be determined by the method of A.G.V. et al. 22 (see the Supplementary Materials). The ground state in the non-interaction limit is the single-particle harmonic eigenstate ψ 0 (x 1 ) multiplied by the Slater determinant for the remaining particles in the states ψ n (x k ) for k = 2, … , N and n = 0, … , N − 2.
For systems of four and five particles, the interpolatory ansatz with these |γ 0 〉 and |γ ∞ 〉 yields integrals that can be evaluated analytically. For larger systems, the integrals are readily evaluated numerically. The resulting energies for N = 4-6 are plotted as dashed lines in Fig. 5. Here we see a very good agreement between the ansatz and the numerically exact results, and in turn excellent agreement with experimental measurements 19 . This indicates that the energetics of the system is captured by the interpolatory ansatz with high accuracy.
However, as we discussed briefly above, the ansatz does not generally reproduce the correct first-order energy term at 1/g = 0. Assuming that the non-interacting state should remain untouched, this prompts us to investigate whether another state can be found in the 1/g = 0 limit that can replace |γ ∞ 〉 in Eq. (4) and in turn give the exact result for the slope of the energy at 1/g = 0. As noted above, Eq. (7) remains valid if we substitute the eigenstate |γ ∞ 〉 with any other state with energy E ∞ that obeys V|γ ∞ 〉 = 0. In particular, any choice of a in Eq. (16) would work, and E opt only depends on a through the wave-function overlap 〈 γ 0 |γ ∞ 〉 . Hence, because E opt is monotonic in Figure 5. Energy spectrum of the impurity system. Energy of the ground state as predicted by the interpolatory ansatz with |γ ∞ 〉 being the energy eigenstate in the q = 0 limit (dotted curve) and the modified ansatz (solid curve) for N = 4, 5, 6. This is compared with exact numerical calculations (circles) and experimental data (squares) 19 . Note that the error bars on the experimental data points are smaller than the squares and are only discernible for the points at small g.
Scientific RepoRts | 6:28362 | DOI: 10.1038/srep28362 〈 γ 0 |γ ∞ 〉 2 , we may optimise the energy with respect to a for all values of g by maximising 〈 γ 0 |γ ∞ (a)〉 2 with respect to a. Leaving the details to the Supplementary Materials, the optimum is For the ground state of the N = 3-6 systems, however, 〈 γ 0 |γ ∞ (a max )〉 2 is very close to the known exact value of 〈 γ 0 |γ ∞ 〉 2 , and is not large enough to make the slope of the energy correct in the strongly-interacting limit, that is, 〈 γ 0 |γ ∞ (a max )〉 2 < λ with λ given by Eq. (9). This indicates that we cannot find a state in the infinite-interaction limit that both has the correct zeroth-order energy, E ∞ , and satisfies the delta-boundary conditions, V|γ ∞ 〉 = 0. However, we caution that this is under the assumptions that the wave function of |γ ∞ 〉 is continuous and has discontinuous derivatives that satisfy the boundary conditions imposed by the zero-range interaction.
We see that the interpolatory ansatz-with whatever choice of |γ ∞ 〉 -does not reproduce the correct energy slope. The modified ansatz, however, does. It is worthwhile to note that the inability to find a state, |γ〉 = α 0 |γ 0 〉 + α ∞ |γ ∞ 〉 , whose expectation value is that predicted by the modified ansatz, does not imply that the modified ansatz cannot be used. The two methods rely on different information about the system: The interpolatory ansatz requires the knowledge of the wave function at g = 0 and 1/g = 0, whereas the modified ansatz uses the knowledge of the energy behaviour. Therefore, one can use the modified anzatz to estimate the energies, and the interpolatory ansatz to approximate the wave functions. Figure 5 compares the ground-state energy of the modified ansatz with results from an exact numerical diagonalisation as well as experimental data 19 . Here we see an improved agreement for large g. The deviations of the modified ansatz from exact numerical diagonalisation are shown in greater detail in Fig. 6. The errors are more than an order of magnitude smaller than those of the unmodified ansatz (shown in Supplementary Fig. S4).
We plot the error scaled against N ↑ instead of E 0 or E, because E 0 scales as ↑ N 2 while Δ E = N ↑ . As seen in Fig. 5, the maximum in the error moves to larger interaction strength for higher N ↑ , i.e. with system size. There seems to be only very little increase in the magnitude of the error for larger N ↑ for the system sizes we have studied. Notice that the error is slightly negative for 1/g ≃ 0 in the case N = 3-6. Since we know the exact energies and slopes around 1/g = 0, the likely cause of this is that we are pushing the accuracy of the exact numerical diagonalisation method here.
Our first observation is that the approximation of the ground-state energy offered by the modified ansatz is so good that it can compete and even in some cases beat state-of-the-art numerical exact methods. The drawback seems quite evident, i.e. we cannot be sure that any state exists in the infinite-interaction limit that when used in Eq. (5) would reproduce the energy of the modified ansatz. Thus, the modified ansatz does not immediately give us any information about the wave function in spite of its near perfect approximation of the energy. We also note that the modified ansatz breaks the variational bound on the ground state and can in principle have lower energy as compared to the exact result. However, we have clearly demonstrated that the slope of the energy at 1/g = 0 is an extremely important quantity for these systems as it appears crucial to reproduce in order to capture the energetics. This highlights the important role played by the recently developed theory of strongly interacting confined systems 22 .
Anderson overlap. Finally, we discuss the so-called Anderson overlap, which is the wave-function overlap between the non-interacting eigenstate, |γ 0 〉 , and the interacting state, |γ〉 , for some value of the interaction strength, g. This quantity is related to the Anderson orthogonality catastrophe 32 , which states that the Anderson overlap is zero in the thermodynamic limit; that is 〈 γ 0 |γ〉 → 0 for N → ∞ , and in particular 〈 γ 0 |γ ∞ 〉 → 0 for N → ∞ .
In Fig. 7 we illustrate how the overlap 〈 γ 0 |γ〉 2 -with |γ〉 given by the interpolatory ansatz-decreases as a function of g. Due to the fact that we only consider a finite-sized system, the overlap does not approach zero as g → ∞ , but the plot clearly shows that the overlap from our interpolatory ansatz tends to decrease as expected. The fact that the ansatz gives a very accurate approximation for the energy of the system does not immediately imply that this is also the case for the wave function. We leave this question for future studies, and the overlaps presented here are thus predictions based on the ansatz. It should be compared either to elaborate exact numerical calculations or to experimental measurements.

Outlook
We have proposed a simple interpolatory ansatz for approximating the energy eigenstate of a confined, one-dimensional system of interacting particles. The ansatz is a linear combination of known eigenstates in the extreme limits of the interaction strength, g → 0 and 1/g → 0, respectively. Thanks to recent advances in the description of the eigenstates in the 1/g = 0 limit, both these wave functions are now available. An analytical expression for the optimum energy of this ansatz is presented which is an elementary function of only two matrix elements; the interaction energy of the eigenstate at g = 0, and the wave-function overlap between the eigenstates in the two limits of g → 0 and 1/g → 0. By focusing on harmonically trapped impurity systems of fermions, we have demonstrated that the ansatz is able to capture the physics of such a system. It gives us a highly accurate approximation for the energy and it also gives us a very simple expression for the wave function.
For both the two-and three-particle systems, we have been able to reproduce the entire energy spectrum with the interpolatory ansatz, save for deeply bound molecular states. The ansatz can be extended to describe deeply bound states as well; this has been shown specifically for the two-particle system. Taking the three-particle system as an example, we have also demonstrated that the ansatz works equally well for mass-imbalanced systems. A future extension of this study might investigate mass-imbalanced systems of four or more particles. It should be noted that the bottleneck here is that there are generally quite few known results about mass-imbalanced systems in the 1/g → 0 limit 33 . It is an open problem to find a general method that yields exact eigenstates in the strongly interacting regime for mass-imbalanced systems.
A drawback of our ansatz is that it is in general not perturbatively correct in the strongly-interacting regime. More precisely, if we take the first order derivative of the energy with respect to 1/g it deviates from the known exact result. We may, however, modify the expression for the energy of the interpolatory ansatz slightly such that it is perturbatively correct to linear order in g for g → 0 and to linear order in 1/g for 1/g → 0. The modified ansatz has great simplicity and accuracy at a level that is competitive with state-of-the-art numerical methods for obtaining the energy of the ground state for arbitrary g. Due to its simplicity it should provide a very useful tool.
We note that although the results presented here assume that the particles are trapped in an external harmonic trap, the formalism is completely general and can be applied for arbitrary external potentials with at least N bound single-particle states for N-body systems. As we have shown, the relative deviation of the energy obtained from the ansatz grows only very slowly with N, and there should be no problem in extending the technique to even larger systems than considered here. The decisive quantity is the overlap of the g = 0 and 1/g = 0 wave functions. This overlap may be computed using the same methods that have recently been used to compute spin chain models for strongly interacting fermions 34 and thus scaling to larger N of order 30 or 40 is certainly within reach.
A future direction would be to consider a generalised version of the interpolatory ansatz where one systematically adds more states at g = 0 and 1/g = 0 in order to gradually improve the comparison. In addition, it is relatively straightforward to apply the interpolatory ansatz in the case of strongly interacting bosons 35,36 , or mixed systems 22,28 . The requirements are knowledge of states in the two limits and their overlaps so that the interpolation can be performed. The formulas for the interpolated energy given here still apply. An example could be an impurity interacting strongly with a Tonks-Girardeau gas of hard-core bosons, which is a topic of great recent interest [37][38][39] .

Methods
In the following, we provide the details of the methods used in applying the interpolatory ansatz to two-, threeand many-particle systems.
Details of the two-particle system. We consider a system of two distinguishable fermions and we define In the absence of interactions, the energy eigenstates are the harmonic eigenstates denoted |n〉 with integer n ≥ 0 (we are only concerned with the motion relative to the center of mass).
Exact solutions of the two-particle problem are available for arbitrary values of g 31,40 . The energy of an exact energy eigenstate is given indirectly by 40 and the wave function of the state is Constructing the ansatz. Let |γ 0 〉 = |n 1 〉 be an energy eigenstate in the non-interaction limit. Here |n 1 〉 is a single-particle eigenstate of a harmonic oscillator Hamiltonian (in the relative coordinate x) with quantum number n 1 . Furthermore, let |γ ∞ 〉 be the corresponding eigenstate in the infinite-interaction limit-that is, |γ 0 〉 → |γ ∞ 〉 as the interaction strength is changed adiabatically from g = 0 to 1/g = 0. By conservation of parity, the states |γ 0 〉 and |γ ∞ 〉 have the same parity. If |γ 0 〉 is odd, 〈 x = 0|γ 0 〉 = 0 and thus 〈 γ 0 |V|γ 0 〉 = 0. Hence, odd harmonic eigenstates do not change as we introduce a non-zero interaction. The even harmonic eigenstates do, however, change; the correct even eigenstate in the infinite-interaction limit is 2 with n 2 = n 1 + 1 for repulsive interactions (g > 0) and n 2 = n 1 − 1 for attractive interactions (g < 0).
Details of the three-particle system. Before we employ the interpolatory ansatz, we first separate out the center-of-mass motion using hyperspherical coordinates 35,41 . This is done merely for convenience and is not in any way essential for the approach. Defining , the hyperradius is given by ρ = + x y 2 2 and the hyperangle is defined by tanφ = y/x. The Hamiltonian of the relative motion becomes where ∇ 2 is the Laplacian in polar coordinates (ρ, φ) and is the interaction.
Limiting cases. The quantum numbers ν and μ can be used to describe the energy eigenstates in both the non-interacting limit and the infinite-interaction limit. The eigenstate wave function in both limits has the general form 41 Scientific RepoRts | 6:28362 | DOI: 10.1038/srep28362 , 2 / 2 2 where ν µ L is a generalised Laguerre polynomial. In the non-interacting limit, the (normalised) angular part of the wave function is where p is the parity of the wave function. In this limit p = (− 1) μ . In the infinite-interaction limit, the wave function vanishes at the lines φ = − 5π/6, − π/2, − π/6, π/6, π/2, 5π/6. In the regions between these lines, the wave function solves the Schrödinger wave equation for the harmonic oscillator without interactions.
The eigenenergies in the limit 1/g = 0 is given by Eq. (11) with μ = 3, 6, 9, 12, … . Each allowed eigenenergy is three-fold degenerate (not counting states with a different radial quantum number ν): One of the eigenstates in each energy triplet is a harmonic eigenstate with μ = 3, 6, 9, 12, … which is unaffected by the interactions. Between the two remaining eigenstates in the triplet, one is odd and the other is even.
The μ eigenvalue can be found by using parity symmetry, the Pauli principle and the delta-boundary conditions of the interaction potential. Using these conditions, one can setup an equation from which μ can be obtained in the limits g = 0 and 1/g = 0: For solutions with odd parity, μ solves the equation