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.


Introduction
Few-fermion systems are the building blocks of matter. Atoms and nuclei are wellknown 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 [3], and the atomic interaction strength may be tuned via the use of Feshbach resonances [4]. The spin-1/2 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 [9]. 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 [8]. 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 [29], 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 N ↑ spin up and N ↓ spin down fermions for the imbalanced case where N ↑ > N ↓ = 1. Zero-range interactions of strength g 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 N ↑ = 2 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 N ↑ and g, but the trapped case has only been solved exactly for N ↑ = 1 [29] and for resonant interaction in three dimensions [37].

Results
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 b = /mω, energies in units of the trap oscillator energy ω, σ = ±, andσ = −σ. The interaction strength, g, becomes dimensionless in units of ωb. 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 x → −x. We concentrate mostly on the first non-trivial case beyond the two-fermion system, which is N ↑ = 2, N ↓ = 1 (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 Appendix 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 Appendix A). The numerical method used in this work therefore represents a significant advance in the description of strongly interacting, finite-size quantum systems. Energy spectrum for the 2+1 system (excluding the center-of-mass contribution). Diverging states on the attractive side are dimmed. The Busch-model for a 1+1 system plus the energy of a spectator particle has been added for comparison.
In Fig. 1 we show the numerically obtained energy spectrum for the 2+1 system; plotted on a cylinder where the fermionization limit |g| → ∞ is on the front while weak-coupling |g| → 0 is on the back. We plot the total intrinsic energy, i.e. the total energy minus ω/2 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 [38]. The first interesting feature is the ground state behavior; it starts as a strongly bound dimer plus a particle for g → −∞, wraps around the cylinder to the non-interacting limit, |g| → 0, and then becomes energetically degenerate with two other states at g → ∞. 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 E = 4 ω, corresponds to having one particle in each of the three lowest HO states. At g = ∞ it becomes degenerate with two interacting states. Note also that we have many molecular branches close to |g| = ∞ for g < 0 (dimmed curves in Fig. 2). Starting from g → 0 + (far left in the figure) and following the odd ground state we see that it makes a jump around g = ∞ before becoming non-interacting at g → 0 − (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 |g| → ∞ 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 g = ∞ 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 [39]. 4 2 0 2 4 Non-interacting state Analytical model A particularly interesting feature of the interacting states, as they cross over to the attractive side of |g| → ∞, is their density distributions that we show in Fig. 3 for the odd ground state. While they are approaching the fermionization limit |g| → ∞ 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 |g| = ∞, 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 Appendix B). Since the states are degenerate at |g| = ∞, 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 twocomponent Fermi systems. These general structures can be experimentally investigated by tunneling [2,40].

Discussion
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 g one can explore the structure on both sides of the resonance [2]. It is possible to go diabatically from the repulsive ground state and onto the repulsive branch on the g < 0 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 g > 0. 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 [41].  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 [42] 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 [28]. 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 g > 0 that becomes exact at g → ∞ and reproduces both degeneracies, densities, and occupation numbers (see Appendix B). Lastly, we note that a recent paper by Gharashi and Blume [55] 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 [55], and both are in agreement with the exact solution for large repulsive interactions presented in reference [56]. However, they do not agree with the results based on symmetry arguments presented in reference [24]. 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 [56] and some explicit examples are given in the supplementary material of reference [55].
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 |n . Each many-body basis state is written as |n 1 . . . n N ↑ ⊗ |n ↓ i.e. a product of a HO antisymmetrized state of the N ↑ spin up particles and a HO single particle state for the spin down particle. The corresponding HO energy is (n 1 + . . . + n N ↑ + n ↓ + N 2 ) where N = N ↑ + N ↓ is the total number of particles. The model space truncation is defined by a total upper limit, i.e. n 1 +. . .+n N ↑ +n ↓ ≤ n tot . Since we are only interested in the intrinsic dynamics of states, a Lawson projection term [57] 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 V eff P is constructed in a truncated two-body space P , defined as the set of two-body relative HO states whose radial quantum number are smaller or equal to a cutoff n P . The effective force V eff P is designed such that its solutions correspond to exact solutions given by the Busch formula [29]. Using a unitary transformation, we construct a two-body effective Hamiltonian H eff P as [33] where E (2) P P is the diagonal matrix formed by the n P + 1 lowest exact energies given by the Busch formula, and U P P is the matrix whose rows are formed by the corresponding eigenvectors projected on P . The effective interaction V eff P is obtained from H eff P by subtracting the HO potential. For each cutoff n P , we diagonalize the many-body Schrödinger equation with V eff P and increase n tot until convergence of the many-body energies is reached [33]. We find that n tot = n P + 2 is sufficient to capture the properties of the effective interaction. With this choice, we can then study the energy convergence as a function of n tot . By construction, this unitary transformation approach will reproduce exact bare Hamiltonian results (both energy spectrum and wave functions) in the n P → ∞ limit.  Figure A1. Relative error (E ntot − E ∞ )/E ∞ in the ground state energy of the 2 + 1 system as a function of the many-body model-space truncation n tot . We show results for several bare interaction strengths (represented by circles, diamonds, and squares), and the effective-interaction results for the strongest case, 1/g = 0.01, (crosses). In this latter case, n tot = n P + 2, where n P is the truncation of the two-body space.
The excellent convergence property of our method is demonstrated in Fig. A1 that shows the relative error (E ntot −E ∞ )/E ∞ in the ground state energy of the 2+1 system as a function of the size of the model space n tot . We show results using the bare interaction for several different strengths, and one sequence of results obtained with the effective interaction V eff P for the strongest case. For each interaction strength, we define E ∞ as the converged result obtained with V eff P . Fast convergence will then be characterized by the relative difference being close to zero already for small n tot . As expected, with the bare interaction, the convergence with n tot is much slower for the strong interactions: for 1/g = 1 the relative error is ∼ 1% for n tot = 64, whereas for 1/g = 0.01, the relative error is still ∼ 4% for the same model space. On the other hand, for 1/g = 0.01, the result obtained with the effective interaction is within 0.01% of the fully converged value already at n tot = 14. Dashed lines correspond to a fit to the bare interaction results using the functional form: E ntot = E ∞ + cn −λ tot , with E ∞ , c, λ as free fit parameters.