Fermionization of two-component few-fermion systems in a one-dimensional harmonic trap
The nature of strongly interacting Fermi gases and magnetism is one of the most important and studied topics in condensed-matter physics. Still, there are many open questions. A central issue is under what circumstances strong short-range repulsive interactions are enough to drive magnetic correlations. Recent progress in the field of cold atomic gases allows to address this question in very clean systems where both particle numbers, interactions and dimensionality can be tuned. Here we study fermionic few-body systems in a one dimensional harmonic trap using a new rapidly converging effective-interaction technique, plus a novel analytical approach. This allows us to calculate the properties of a single spin-down atom interacting with a number of spin-up particles, a case of much recent experimental interest. Our findings indicate that, in the strongly interacting limit, spin-up and spin-down particles want to separate in the trap, which we interpret as a microscopic precursor of one-dimensional ferromagnetism in imbalanced systems. Our predictions are directly addressable in current experiments on ultracold atomic few-body systems.
Few-fermion systems are the building blocks of matter. Atoms and nuclei are well-known examples, but also systems such as quantum dots, superconducting grains, and other nanoscale structures are of great interest. The key to understanding such structures is first and foremost the relation between the discrete level structure, due to the finite size, and the strength and nature of interparticle interactions. An exciting recent development in atomic physics is the experimental realization of few-body Fermi systems with ultracold atoms [1, 2]. These setups are extremely versatile as the potential that traps the atoms can produce lattices and/or low-dimensional geometries , and the atomic interaction strength may be tuned via the use of Feshbach resonances . The spin- nature of electrons or nucleons is addressable by populating two hyperfine states in the atoms and we thus have a direct mapping from the atomic setup to ordinary matter. We will refer to these two components as spin up and spin down.
A seminal contribution of ultracold atomic gas research is the realization of strongly interacting quantum gases [5, 6, 7, 8] using confinement-induced resonances . The Tonks-Girardeau (TG) gas [10, 11, 12] of strongly repulsive bosons that displays fermionic behavior [13, 14] is one such example [5, 6, 7]. The so-called super-TG (sTG) limit of very strong attractive interactions has also been addressed both theoretically [15, 16, 17, 18, 19, 20, 21, 22, 23] and experimentally . Most recently, the TG and sTG states have been explored in fermionic systems [24, 25, 26, 27, 28]. While the two-body system in a harmonic trap has a well-known exact solution for any interaction strength, known as the Busch model , two-component fermionic few-body systems with more than two particles have not been solved exactly. Although a number of numerical studies have been performed (see discussion below), many questions still remain related to the main difficulty in the handling of very strong interactions in the vicinity of the fermionization limit.
In this article we face this challenge and consider the experimentally accessible situation of a harmonically trapped few-body system in one dimension with spin up and spin down fermions for the imbalanced case where . Zero-range interactions of strength are employed between different spin components, while the identical spin particles remain non-interacting by the Pauli exclusion principle. This is a few-body analog of the fermionic polaron, which is currently under intense study [30, 31, 32]. We solve the few-body problem for various interaction strengths using a numerical technique inspired by developments in nuclear physics [33, 34]. This method is exact for any interaction strength in the limit of infinite model space; it yields both the energy spectrum and energy eigenstate wave functions and indeed shows excellent convergence properties. In addition, we present an exact solution for the case in the fermionization limit of infinite interaction strength, using an analytical model. The classic work of McGuire [35, 36] solved the untrapped case with periodic boundary conditions for arbitrary and , but the trapped case has only been solved exactly for  and for resonant interaction in three dimensions .
We are interested in obtaining numerically exact eigensolutions for the system of trapped atoms with arbitrary strong, zero-range interaction between different spin components. We use harmonic-oscillator (HO) units, in which the Hamiltonian is
where lengths are in units of the oscillator length , energies in units of the trap oscillator energy , , and . The interaction strength, , becomes dimensionless in units of . These units are used for all quantities, unless we explicitly state otherwise. The Hamiltonian is parity invariant and one can classify states as either even or odd under . We concentrate mostly on the first non-trivial case beyond the two-fermion system, which is , (denoted 2+1). The many-body problem is solved using an effective-interaction approach that uses the exact analytical solution of the two-body problem as input. As we discuss below in A, this speeds up the convergence tremendously and allows us to obtain very accurate results for mesoscopic samples with particle numbers of order ten. We stress that our approach is far superior to exact diagonalization with the bare zero-range interaction, which has a very slow convergence (see A). The numerical method used in this work therefore represents a significant advance in the description of strongly interacting, finite-size quantum systems.
In Fig. 1 we show the numerically obtained energy spectrum for the 2+1 system; plotted on a cylinder where the fermionization limit is on the front while weak-coupling is on the back. We plot the total intrinsic energy, i.e. the total energy minus from the center of mass motion in the trap. This way of plotting the spectrum emphasizes the spectral flow and the connection to the Zeldovich rearrangement effect . The first interesting feature is the ground state behavior; it starts as a strongly bound dimer plus a particle for , wraps around the cylinder to the non-interacting limit, , and then becomes energetically degenerate with two other states at . Another representation of the 2+1 spectrum for odd and even parity states is shown in Fig. 2. The horizontal lines correspond to totally antisymmetric states, which are non-interacting in the case of zero-range interactions. For example, the lowest such state, at , corresponds to having one particle in each of the three lowest HO states. At it becomes degenerate with two interacting states. Note also that we have many molecular branches close to for (dimmed curves in Fig. 2). Starting from (far left in the figure) and following the odd ground state we see that it makes a jump around before becoming non-interacting at (far right). This is an analog of the so-called repulsive branch for untrapped polarons [31, 32]. Repulsive branch means that excited states are pushed up on the attractive side of the resonance in constrast to the lower-lying molecular branches that become strongly bound, as shown in Figs. 1 and 2. However, the jump endured by the odd and even states that become degenerate at is quite different. For comparison, we plot the two-body Busch results shifted by the energy of a free spectator particle (dashed blue line in Fig. 2), which turns out to be almost identical to the even parity state at low energy. This even parity state therefore has an atom-dimer structure, with almost no interaction between atom and dimer. This has also been observed in three-dimensional traps .
A particularly interesting feature of the interacting states, as they cross over to the attractive side of , is their density distributions that we show in Fig. 3 for the odd ground state. While they are approaching the fermionization limit for the total density, the spin-resolved densities demonstrate a distinct separation in the trap. This we interpret as a precursor of ferromagnetic behavior in a one-dimensional few-body context for imbalanced systems. In the vicinity of , the ground, first excited, and non-interacting states all have completely different spin-resolved densities; the ground state has the impurity at the center and the first excited state has the impurity at the edge while the non-interacting state yields a three-hump profile independent of spin. This has all been verified using our analytical model (see B). Since the states are degenerate at , this clearly demonstrates that the behavior is not due to energetics but to different correlations in the wave functions. It also shows that the fermionization limit is very different for two-component fermions as compared to fermionization of bosons.
The approach that we have presented can be applied to larger systems. In Fig. 4 we show the spin-resolved densities for the ground states of the 3+1, 6+1, and 9+1 systems as a function of the (repulsive) interaction strength. These spin densities show the same spin separation behavior in the limit of fermionization as the 2+1 case in Fig. 3. Our results imply that this is a more general feature of one-dimensional two-component Fermi systems. These general structures can be experimentally investigated by tunneling [2, 40].
Current experiments on few-fermion systems [1, 2] can study the structures that we predict by performing tunneling measurements that map out the occupancies of the few-body wave function. By varying one can explore the structure on both sides of the resonance . It is possible to go diabatically from the repulsive ground state and onto the repulsive branch on the side since the overlap with the atom+dimer molecular branch is small. It is thus possible to investigate a large part of the parameter space. In Fig. 5 we show the occupancies of different single-particle levels in the trap. Note how well our analytical model reproduces the numerical results for . By selective ejection of the minority particle it is possible to measure the majority occupation number. A preliminary comparison to experimental data shows agreement with our predictions .
Our numerical and analytical findings show that around fermionization the two spins tend to separate in the trap. We interpret this as a few-body precursor of Stoner ferromagnetism  in one dimensional imbalanced systems. Ferromagnetism is hotly pursued topic at the moment both in balanced and imbalanced Fermi systems [43, 44, 45, 46, 47, 48, 49, 50]. As discussed above, this separation of species should be directly observable in current experiments [2, 51, 52]. Other methods have been employed recently to similar systems to study energy spectra [53, 54] and also the density profiles . Here we have presented a complementary method that converges extremely fast for multi-particle systems. We have also provided new analytical insights into the problem by obtaining a wave function for that becomes exact at and reproduces both degeneracies, densities, and occupation numbers (see B). Lastly, we note that a recent paper by Gharashi and Blume  has studied some of the same systems using a different method. Our results for similar systems (2+1 and 3+1) are in agreement with reference , and both are in agreement with the exact solution for large repulsive interactions presented in reference . However, they do not agree with the results based on symmetry arguments presented in reference . The reasons that symmetry and group theoretical arguments and spin algebra cannot be used to determine the eigenstates for large but finite interaction strengths are discussed in the appendix of reference  and some explicit examples are given in the supplementary material of reference .
Magnetism is often considered a bulk property of a system while magnetic correlations such as super-exchange, etc., are typically discussed in the context of just a few particles. Many studies of magnetism in one dimension are conducted in a flat potential and employing periodic boundary conditions starting from the work of McGuire [35, 36]. Dispensing of the periodic boundaries (which are not suitable for the few-body systems studied here), we may consider how our results would change if we had replaced the harmonic trap with a hard-wall box potential. Since the degeneracy in the strongly interacting limit is a result of short-range correlations in the wave function (nodal structures), we do not expect anything to change there. However, since the approach to the strongly interacting regime certainly depends on the single-particle wave functions (as clearly seen in the analytical model presented here) the spectrum will change quantitatively under the constraint that the degeneracies at infinite coupling strength are preserved.
The effective-interaction approach used in this work is key to the quality of our numerical results and to our conclusions. In the construction of these effective interactions we benefit from having access to the exact two-body solutions. However, we stress that, using numerical two-body solutions, this approach can be generalized to study many-body systems in higher dimensions, with finite-range interactions, and in any trapping potential.
Appendix A Effective interaction approach
We solve numerically the many-body Schrödinger equation with the Hamiltonian (1) in a finite basis constructed from the HO single-particle states . Each many-body basis state is written as i.e. a product of a HO antisymmetrized state of the spin up particles and a HO single particle state for the spin down particle. The corresponding HO energy is where is the total number of particles. The model space truncation is defined by a total upper limit, i.e. . Since we are only interested in the intrinsic dynamics of states, a Lawson projection term  is used to push away many-body solutions corresponding to excitations of the center of mass motion.
Instead of the bare zero-range interaction in (1), we consider an effective two-body interaction in order to speed up the convergence of the eigenenergies with respect to the size of the many-body basis. The effective potential is constructed in a truncated two-body space , defined as the set of two-body relative HO states whose radial quantum number are smaller or equal to a cutoff . The effective force is designed such that its solutions correspond to exact solutions given by the Busch formula . Using a unitary transformation, we construct a two-body effective Hamiltonian as 
where is the diagonal matrix formed by the lowest exact energies given by the Busch formula, and is the matrix whose rows are formed by the corresponding eigenvectors projected on . The effective interaction is obtained from by subtracting the HO potential. For each cutoff , we diagonalize the many-body Schrödinger equation with and increase until convergence of the many-body energies is reached . We find that is sufficient to capture the properties of the effective interaction. With this choice, we can then study the energy convergence as a function of . By construction, this unitary transformation approach will reproduce exact bare Hamiltonian results (both energy spectrum and wave functions) in the limit.
The excellent convergence property of our method is demonstrated in Fig. 6 that shows the relative error in the ground state energy of the system as a function of the size of the model space . We show results using the bare interaction for several different strengths, and one sequence of results obtained with the effective interaction for the strongest case. For each interaction strength, we define as the converged result obtained with . Fast convergence will then be characterized by the relative difference being close to zero already for small . As expected, with the bare interaction, the convergence with is much slower for the strong interactions: for the relative error is 1% for , whereas for , the relative error is still 4% for the same model space. On the other hand, for , the result obtained with the effective interaction is within 0.01% of the fully converged value already at . Dashed lines correspond to a fit to the bare interaction results using the functional form: , with as free fit parameters.
Appendix B Analytical model
We now present an analytical approach that captures the behavior of the wave function exactly at . We first define Jacobi coordinates, and , and use these to obtain the spherical variables, and . For , the (unnormalized) eigenstates are and , where and are non-negative integers and is the associated Laguerre polynomial. The corresponding energies are . We need only consider the wave function in the interval since one can use the Pauli principle and parity invariance to extend to . At opposite spins overlap. The full solution can thus be obtained by matching the wave function and its derivative on the line . We have
where and are solutions for and respectively. For these equations are complicated to solve, but by introducing an ad hoc rescaled strength parameter we decouple the equations and can write for , where . This rescaled model becomes exact when . The eigenfunctions and eigenvalue equations can now be obtained by using the free angular solutions and the conditions in Eqs. (3) and (4). The nature of the three-fold degeneracy at fermionization seen in Fig. 1 and Fig. 2 comes from the odd and even solutions of these equations, while the non-interacting has the structure . In the case of , the angular wavefunction for odd parity become and , and for even parity and , while the energies can be obtained from the algebraic equations for odd, , and even solutions. and , are normalization factors. The important point is that and are non-zero in the limit and thus the wave functions are non-zero.
- Serwane F 2011 et al. Science 332 336
- Zürn G et al. 2012 Phys. Rev. Lett. 108 075303
- Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys. 80 885
- Chin C, Grimm R, Julienne P and Tiesinga E 2010 Rev. Mod. Phys. 82 1225
- Paredes B et al. 2004 Nature 429 277
- Kinoshita T, Wenger T and Weiss D S 2004 Science 305 1125
- Kinoshita T, Wenger T and Weiss D S 2005 Phys. Rev. Lett. 95 190406
- Haller E et al. 2009 Science 325 1224
- Olshanii M 1998 Phys. Rev. Lett. 81 938
- Tonks L W 1936 Phys. Rev. 50 955
- Girardeau M D 1960 J. Math. Phys. 1 516
- Lieb E H and Liniger W 1963 Phys. Rev. 130 1605
- Deuretzbacher F et al. 2008 Phys. Rev. Lett. 100 160405
- Christensson J, Forssén C, Åberg S and Reimann S M 2009 Phys. Rev. A 79 012707
- Astrakharchik G E, Blume D, Giorgini S and Granger B E 2004 Phys. Rev. Lett. 92 030402
- Astrakharchik G E, Boronat J, Casulleras J and Giorgini S 2005 Phys. Rev. Lett. 95 190407
- Batchelor M T, Bortz M, Guan X-W and Oelkers N 2005 J. Stat. Mech. 2005 L10001
- Tempfli E, Zöllner S and Schmelcher P 2009 New J. Phys. 11 073015
- Girardeau M D and Astrakharchik G E 2010 Phys. Rev. A 81 061601(R)
- Girardeau M D 2011 Phys. Rev. A 83 011601(R)
- Valiente M 2012 Europhys. Lett. 98 10010
- Brouzos I and Schmelcher P 2012 Phys. Rev. Lett. 108 045301
- Girardeau M D and Astrakharchik G E 2012 Phys. Rev. Lett. 109 235305
- Guan L, Chen S, Wang Y and Ma Z-Q 2009 Phys. Rev. Lett. 102 160402
- Girardeau M D 2010 Phys. Rev. A 82 011607(R)
- Chen S, Guan X W, Yin X, Guan L and Batchelor M T 2010 Phys. Rev. A 81 031608(R)
- Guan L and Chen S 2010 Phys. Rev. Lett. 105 175301
- Brouzos I and Schmelcher P 2013 Phys. Rev. A 87 023605 (2013). Note that the densities provided in this reference are different from those found in the current study.
- Busch T, Englert B-G, Rza̧żewski K and Wilkens M 1998 Found. Phys. 28 549
- Schirotzek A, Wu C-H, Sommer A and Zwierlein M W 2009 Phys. Rev. Lett. 102 230402
- Kohstall C et al. 2012 Nature 485 615
- Koschorreck M et al. 2012 Nature 485 619
- Rotureau J 2013 Eur. Phys. J D 67 153
- Lisetskiy A F, Barrett B R, Kruse M K G, Navrátil P, Stetcu I and Vary J P 2008 Phys.Rev. C 78 044302
- McGuire J B 1965 J. Math. Phys. 6 432
- McGuire J B 1966 J. Math. Phys. 7 123
- Werner F and Castin Y 2006 Phys. Rev. Lett. 97 150401
- Farrell A, MacDonald Z and van Zyl B P 2012 J. Phys. A: Math. Theor. 45 045303
- Daily K M and Blume D 2010 Phys. Rev. A 81 053615
- Rontani M 2012 Phys. Rev. Lett. 108 115302
- We thank the group of Selim Jochim for a comparison of our results to their unpublished data for .
- Stoner E 1933 Philos. Mag. 15 1018
- Jo G-B et al. 2009 Science 325 1521
- Cui X and Zhai H 2010 Phys. Rev. A 81 041602(R)
- Pekker D et al. 2011 Phys. Rev. Lett. 106 050402
- Zhang S and Ho T-L 2011 New J. Phys. 13 055003
- Massignan P and Bruun G M 2011 Eur. Phys. J D 65 83
- Sanner C et al. 2012 Phys. Rev. Lett. 108 240404
- Cui X and Ho T-L 2013 Phys. Rev. Lett. 110 165302
- Massignan P, Yu Z and Bruun G M 2013 Phys. Rev. Lett. 110 230401
- Zürn G PhD thesis, Ruprecht-Karls Universität Heidelberg (2012)
- Sala S et al. 2013 Phys. Rev. Lett. 110 203202
- Harshman N L 2012 Phys. Rev. A 86 052122
- Gharashi S E, Daily K M and Blume D 2012 Phys. Rev. A 86 042702
- Gharashi S E and Blume D 2013 Phys. Rev. Lett. 111 045302
- Volosniev A G et al. 2013 Preprint arXiv:1306.4610v2
- Gloeckner D H and Lawson R D 1974 Phys. Lett. B 53 313
- Zinner N T et al. 2013 Preprint arXiv:1309.7219