Ultracold few fermionic atoms in needle-shaped double wells: spin chains and resonating spin clusters from microscopic Hamiltonians emulated via antiferromagnetic Heisenberg and t-J models

Advances with trapped ultracold atoms intensified interest in simulating complex physical phenomena, including quantum magnetism and transitions from itinerant to non-itinerant behavior. Here we show formation of antiferromagnetic ground states of few ultracold fermionic atoms in single and double well (DW) traps, through microscopic Hamiltonian exact diagonalization for two DW arrangements: (i) two linearly oriented one-dimensional, 1D, wells, and (ii) two coupled parallel wells, forming a trap of two-dimensional, 2D, nature. The spectra and spin-resolved conditional probabilities reveal for both cases, under strong repulsion, atomic spatial localization at extemporaneously created sites, forming quantum molecular magnetic structures with non-itinerant character. These findings usher future theoretical and experimental explorations into the highly-correlated behavior of ultracold strongly-repelling fermionic atoms in higher dimensions, beyond the fermionization physics that is strictly applicable only in the 1D case. The results for four atoms are well described with finite Heisenberg spin-chain and cluster models. The numerical simulations of three fermionic atoms in symmetric double wells reveal the emergent appearance of coupled resonating 2D Heisenberg clusters, whose emulation requires the use of a t-J-like model, akin to that used in investigations of high T$_c$ superconductivity. The highly entangled states discovered in the microscopic and model calculations of controllably detuned, asymmetric, double wells suggest three-cold-atom DW quantum computing qubits.


Introduction
The unparalleled experimental advances and control achieved in the field of ultracold atoms have rekindled an intense interest in emulating magnetic behavior using ultracold atoms in optical traps [1,2]. Quantum magnetism and spintronics in both extended [3,4,5,6,7] and finite-size [8,9,10,11,12] systems have a long history. Recently antiferromagnetism (AFM) without the assistance of an external periodic ordering potential has been demonstrated experimentally for N = 3 and N = 4 ultracold 6 Li atoms confined in a single-well one-dimensional optical trap [13]. Moreover, progress aiming at bottom-up approaches to fermionic many-body systems, addressing entanglement, quantum information, and quantum magnetism in particular, is predicated on experimental developments of which the recently created double-well (DW) ultracold atom traps [14,15,16] are the first steps.
To date the theoretical studies of magnetism of a few ultracold atoms have mainly addressed [17,18,19,20,21,22] strictly one-dimensional systems trapped within a single well (see, however, Refs. [19,20] for double-well configurations), where the fermionization theory [23,24,25] (applicable to 1D systems in the limit of infinite strength of the contact interaction) can assist in inventing analytic forms for the correlated many-body wave functions. However, the recently demonstrated ability to create needle-shaped double well traps [15], and the anticipated near-future further development of small arrays of such needle-like quasi-1D traps in a parallel arrangement (PA, see schematics in Fig. 1), whose corresponding physics incorporates certain twodimensional (2D) aspects [26], enjoin the development of additional conceptual and computational theoretical methodologies.
We remark that the physics of ultracold atoms in 1D and 3D single traps has been investigated also away from the fermionization point using a Lippman-Schwinger equation approach, see Refs. [27,28] and [29], respectively. Similarly, states of ultracold fermions in a single strictly-1D trap, away from the fermionization limit, using an exact diagonalization of the many-body Hamiltonian have been reported [30].
In this paper, using large-scale configuration-interaction (CI) calculations as means for exact diagonalization of the microscopic Hamiltonian, we report that for N = 4 (even number) strongly interacting ultracold fermions in a double-well trap with parallel arrangement [DWPA [26]; see Figs. 1(II,III)] the many-body problem can be reduced to that of a 2D rectangular AFM Heisenberg ring. The associated mapping between the many-body wave function and the spin eigenfunctions [31] for N = 3 and N = 4 electrons confined in single and double semiconductor quantum dots has been predicted in previous studies [11,12] to occur through the formation of quantum molecular structures in the regime of strong long-range Coulombic repulsion. Such molecular structures are usually referred to as Wigner molecules (WMs) [32]. For N = 3 (odd number) ultracold fermions, few-body quantum magnetism requires introduction of a more complex t-J-type [5,6] model; here the t-J model consists of two coupled and resonating triangular 2D Heisenberg clusters. In all cases, we find AFM ordered ground states.
We remark that the emergence of a resonance associated with the symmetrization of the many-body wave function in two-center/three-electron bonded systems is well known [33,34,35] in theoretical chemistry and in particular in the valence-bond treatment of the three-electron bond which controls the formation of molecules like He + 2 , NO, and F − 2 . The concept of the three-electon resonant bond and its significance were introduced in 1931 in a seminal paper by Linus Pauling [36,37].
The emergence of the simple, as well as the resonating, Heisenberg clusters is a consequence of the spatial localization of the strongly-interacting, highly correlated fermionic atoms and the formation [26] of quantum 1D and 2D molecule-like structures, referred to as ultracold Wigner molecules (UCWMs). The name of Wigner is used here in the context of ultracold atoms in order to emphasize the universal aspects that are present in the few-fermion molecular structures irrespective of the nature of the repelling two-body interaction, i.e., contact versus long-range Coulomb. In this way the concept of UCWM extends and incorporates [38] the fermionization physics [23,24,25] beyond the restricted 1D case.
We note that due to the quantum character [32] of the WMs, the spatial localization of fermions is not necessarily pointlike as in the classical electronic Wigner crystal [39,40,41]. However, depending on the strength of the Coulomb repulsion compared to the quantal kinetic energy, the regime of high-degree localization can be reached also in Wigner molecules formed by electrons confined in quantum dots [32,42,43]. In contrast, fermions interacting via a a repulsive contact interaction cannot attain a similar degree of strong spatial localization; as a result, the UCWMs retain their full quantal character even in the limit of infinite repulsion.
The resonant coupling of magnetic configurations through the tunneling of electrons between occupied and empty sites has long been studied. Two well-known relevant fields are: (i) the socalled direct exchange mechanism [3,4] (related to ferromagnetism in the mixed-valency manganites of perovskite structure), and (ii) the t-J model [5,44] which modifies (away from the half filling) the antiferromagnetic Heisenberg Hamiltonian associated with the Mott insulator at half-filling. The t-J model has attracted much attention, because it has been proposed for explaining the high-T c superconductivity arising in the case of underdoped insulators [44]. Due to the antiferromagnetic aspect, our resonating model Hamiltonian for N = 3 fermions (see Section 4.2 below) represents a finite variant of the t-J model. The emergence of the t-J model in this work suggests future investigations of fundamental aspects associated with the physics of high-T c Figure 1. Energy versus −1/g spectra, SPDs (green surfaces), and spin-resolved CPDs (red surfaces) of N = 4 strongly-repelling 6 Li atoms in a double-well confinement with a parallel arrangement of the two 1D traps, as a function of the interwell separation d and/or interwell barrier V b . Schematic (I) shows a SW in the y-direction, and (II) displays a symmetric DWPA, where interwell tunneling in the x-direction occurs along the entire y-range of the trap, conferring 2D aspects to the trap. Insert (III) shows the sites in the 2D Heisenberg-ring spin model. (a,b,c) d = 0 (single well). (d) d = 2.5 µm and V b = 6.08 kHz (lower barrier). (e,f,g,h,i) d = 2.5 µm and V b = 11.14 kHz (high barrier). In all cases, the confinement frequencies of the 1D traps are ω x = 6.6 kHz and ω y = 1 kHz. The SPDs and CPDs in (b,c) and (f,g) correspond to the S = 0, S z = 0 CI ground states [gs, brown curves in the associated spectra (a) and (e)] at the point (marked by a star) −1/g = −0.1( √ 2l 0y ω y ) −1 . The SPD and CPD in (h,i) correspond to the S = 0, S z = 0 CI excited state [orange curve in the associated spectrum (e)]. g here is the 1D contact-interaction strength along the y direction [26]. All three CPDs display the distributions of the two down spins when the fixed spin-up fermion (see the black arrow) is placed at r 0 = (0, −0.8 µm) in (c), r 0 = (−1.3 µm, 1 µm) in (g), and r 0 = (−1.28 µm, 0.94 µm) in (i). l 0y(x) = [ /(M ω y(x) )] 1/2 is the harmonic-oscillator length; M = 9.99 × 10 −27 kg is the mass of 6 Li. The zero of energy in the spectra corresponds to the ground-state total energy of the corresponding noninteracting system, that is, to 4 ω y + 2 ω x = 17.20 ω y in (a), 13.45 ω y in (d), and 14.91 ω y in (e).
superconductivity via studies utilizing the ability to prepare and measure trapped ultracold fermionic atom systems with precise control over the number of atoms and the strength of interatomic interactions. We complement our investigations by further highlighting the differences arising from the different geometries of the traps in both the parallel and linear (LI) arrangement; in analogy to the DWPA designation defined above, a double-well trap with linear arrangement of the two needle-like wells will be denoted as DWLI; see schematics in Figs. 2(I,II). Our theoretical predictions can be directly confirmed using the recently developed experimental techniques. We stress again that the regime of ultracold Wigner-molecule (UCWM) formation and of the associated simple-Heisenbergchain and t-J resonating-spin-chains magnetism appears for strong interparticle interactions and contrasts sharply with the regime of itinerant magnetism [45,46], which appears for weaker interactions. The microscopic treatment of itinerant magnetism (weaker interactions) can be handled within mean-field approaches (e.g. Hartree-Fock), whereas the regime of spin chains (strong repulsion) considered here entails conservation of the total spin and requires more sophisticated approaches like the full CI.
Finally, we note that the related three-electron system in semiconductor double quantum dots has recently attracted a major attention in conjunction with the fabrication and implementation of pulse-gated fast hybrid qubits for solid-state-based quantum computing [47,48]. These advances and the fascinating physics of double-welltrapped three ultracold fermionic atoms that we uncover, and in particular the high degree of entanglement predicted by us for strong interatomic repulsion (see Sections 3, 4, and 5) and the very slow decoherence in such traps, suggest future exploration of this system as a robust ulracold 3-atom DW qubit.
Before leaving the Introduction, we wish to clarify that the term antiferromagnetic is used by us to characterize finite systems having a ground state with the minimum possible value of the total spin, i.e., S = 1/2 for N = 3 fermions and S = 0 for N = 4 fermions.
The plan of the paper is as follows: A statement of the many-body Hamiltonian, including a description of the doublewell employed by us, is given in Section 2.
In Section 3, we describe investigations concerning four ultracold 6 Li atoms in double-well traps with both parallel and linear arrangement of the individual needlelike wells. A comparison with the case of four fermions in a quasi-1D single well is also included in order to appreciate the rich additional magnetic behaviors associated with a double well. Section 3 is divided into two subsections, with Sec. 3.1 describing results of purely microscopic CI calculations, and Sec. 3.2 establishing the mapping onto the Heisenberg 4-fermion phenomenological Hamiltonian.
Section 4 presents our studies concerning the case of three ultracold 6 Li atoms in DWPA and DWLI traps, as well as the comparison with the corresponding case of a single well. Sec. 4.1 describes CI results for both the symmetric and tilted cases; for the tilted case, this section establishes the mapping onto a 3-fermion Heisenberg model. Going beyond the Heisenberg model, Sec. 4.2 introduces the t−J model and establishes, in analogy with the CI results, its validity for describing the case of symmetric double wells.
Section 5 describes the entanglement properties of the CI many-body wave functions.
The appendices provide detailed information concerning the mathematical formalism associated with the spin eigenfunctions and the finite Heisenberg and t − J models. In particular, Appendix A provides a brief description of the branching diagram that describes the multiplicities (number of degenerate spin states) of a given total spin S. This Appendix also presents in the Ising basis the general formulas that describe a spin eigenfunction (i) with S = 0 and S z = 0 for four fermions and (ii) with S = 1/2 and S z = 1/2 for three fermions. These general formulas incorporate in a compact form both the orthogonal basis of spin functions that spans the spin space for a given S, as well as any linear superposition of them. In addition, Appendix A describes the process of mapping the many-body CI wave functions onto these spin eigenfunctions.
Appendix B discusses the mathematics of the Heisenberg model for 4 localized fermions in a ring-like rectangular configuration (DWPA case), while Appendix C discusses the corresponding case for an open chain arrangement (DWLI case).
Appendix D discusses the mathematics of the Heisenberg model for 3 localized fermions in a triangular (DWPA case) and linear (DWLI case) configuration, both associated with tilted wells.
Finally, Appendix E discusses the mathematics of the more general t − J model for 3 localized fermions in the case of a double trap with symmetric wells.

Many-Body Hamiltonian
A DWLI trap can be treated as a strictly-1D problem along a single direction. The DWPA trap which consists of two parallel needle-like wells, however, cannot be treated solely along one direction (e.g., the x-direction). Instead it requires consideration of the y coordinate as well. To treat both cases in a unified way, we consider a many-body Hamiltonian for N fermions of the form where r i − r j denotes the relative vector distance between the i and j fermions (e.g., 6 Li atoms). This Hamiltonian is the sum of a single-particle part H(i), which guarantees the needle-like shape of the individual wells, and the two-particle contact interaction.
The external confining potential [in H(i)] that models a double well (DW) is based on a two-center-oscillator (TCO) model [12,26] exhibiting a variable smooth neck along the x-direction. Along the x direction, this TCO model allows for an independent variation of both the separation d and of the barrier height V b between the two wells; see Fig. 3. Along the y-direction, the confinement consists of that of a single harmonic oscillator. The values of the frequencies ω x1 (left well), ω x2 (right well) and ω y that confine the two wells along the x and y directions, respectively, are also allowed to  [26]. The CPD in (d) displays the distribution of the two down spins when the fixed spin-up fermion (see the black arrow) is placed at r 0 = (0,0.8 µm). The zero of energy in the spectra corresponds to the ground-state total energy of the corresponding non-interacting system, that is, to 202.86 ω x in (a) and 203.52 ω x in (b). The difference in the non-interacting energies between (a) and (b) is due to the different interwell barrier. For two wells at infinite separation, the total energy for 4 non-interacting fermions is equal to 2 ω x +2 ω y = 202 ω x . vary independently; here we choose ω x1 = ω x2 = ω x . The needle-like shape of each individual trap is enforced by assuming that ω x << ω y (DWLI case) or ω x >> ω y (DWPA case). The TCO further allows consideration of a tilt ∆ between the left and right wells. Fig. 3 illustrates the TCO confining potentials in the x direction used in Fig. 4 below for the study of N = 3 fermions in tilted and symmetric double wells.
As we mentioned in the Introduction, we use the CI method for determining the solutions of the many-body problem specified by the Hamiltonian (1). The CI method expresses the fermionic many-body wave function as a supperposition of Slater determinants, and it is well known in quantum chemistry and in few-body physics of electrons; for a basic description of the CI method, see Ref. [49]. Thus a detailed description of the CI method will not be repeated here. Specific adaptations by us of this method to a few electrons in 2D semiconductor quantum dots and rotating bosons in the lowest Landau level have been reported in Refs. [12,32] and [50], respectively. An earlier application by us of this method to the case of N = 2 trapped ultracold fermions was reported in Ref. [26]. The reader can find an expanded description of the CI method in the literature mentioned above.
Convergence in the CI calculations is reached through the use of a basis of up to eighty TCO single-particle states as needed. Note that the TCO single-particle states automatically adjust to the separation d as it varies from the limit of the unified atom d = 0 to that of the two fully separated traps (for sufficiently large d). We verified that for ω y /ω x = 100 (strictly-1D single trap), our CI calculations agree with the results of Table 2 of Ref. [30].
The matrix elements of H MB between the CI determinants are calculated using the Slater rules. An important ingredient in this respect are the two-body matrix elements of the contact interaction in the basis (of dimension K) formed out of the single-particle (space) eigenstates ϕ i (r) of the TCO Hamiltonian. Because the individual wells remain needle-like in all of our calculations here, the s-wave scattering between two ultracold fermions takes place primarily along a single dimension, either the y-dimension or the x-dimension. As a result, the parameter g xy in front of the δ(r i − r j ) function in Eq. (1) does not reflect the physical process of two-dimensional s-scattering. Rather it is an auxiliary theoretical parameter that allows us to treat the DWPA and DWLI traps on an equal footing. In particular, the actual 1D interparticle interaction strengths, g, are related to g xy as follows where u is a dummy variable and W is the lowest-in-energy single-particle state in the y (x) direction for the LI (PA) trap configurations, respectively. Note that the 1D strength g relates to the 3D s-scattering length a 3D via the relation [51], precisely as is done in the experimental studies of Ref. [15]; µ = M/2 is the relative mass and l ⊥ is the harmonic oscillator strength in the direction perpendicular to the needle.
For the CI calculations in this paper, we assume that the total-spin projection S z = 0 for 4 fermions or S z = 1/2 for 3 fermions. This suffices to provide the full energy spectrum, as long as the many-body Hamiltonian does not depend on S z . Naturally, the many-body wave functions characterized by a given total spin S are different for the different projections S z > 0. For lack of space, we will not consider here many-body wave functions with S z = 0 for four fermions or with S z = 1/2 for three fermions. For an earlier study of such wave functions in the case of four electrons in a double quantum dot, see Ref. [12].

Four fermionic ultracold atoms: CI results
We treat here three different types of traps: a SW quasi-1D trap [see We start with the four-atom double-well systems, followed by a comparison with the single-well trap (end of Secs. 3.1 and 3.2), which is used as a reference point to allow for a deeper appreciation of the richness of magnetic behaviors introduced by the double-well geometries.
Figs. 1 and 2 illustrate the evolution of the spectra of N = 4 6 Li atoms for the DWPA and DWLI cases, respectively, as a function of the separateness of (or alternatively the strength of tunneling between) the two 1D wells (resulting from both the effect of separation in distance, d, and the height of the interwell barrier V b ). The limiting case of the single-well ("united atom") quasi-1D trap (at d = 0) is displayed in Fig. 1(a). The opposite limit of strongly separated wells is displayed in Fig. 1(e) and Fig.  2(b), respectively. All spectra are shown in the range −1/( √ 2l 0y(x) ω y(x) ) ≤ −1/g ≤ 0, which covers the regime of strong interparticle contact repulsion. A salient common feature of all five energy spectra in Figs. 1 and 2 is the emergence of a separate band formed by six low-energy states as the interaction strength approaches infinity (i.e., as −1/g → −0); all six states become degenerate at −1/g = 0. Qualitative differences between these spectra amount only in the extent of the spreading out of the six curves; the most spread out case (with six clearly distinct lines) arises for the single well [ Fig.  1(a)], whereas the strongly separated cases display a characteristic 1-2-3 degeneracy in the whole range −1 ≤ −1/g ≤ 0 [ Fig. 1(e) and Fig. 2(b)]. The tendency towards the regrouping of the energy curves according to the 1-2-3 degeneracy pattern is also visible in the intermediate cases [ Fig. 1(d) and Fig. 2(a)].
In all instances, i.e., for both the DWPA and DWLI cases, as well as the SW case, there are two states with total spin S = 0, three states with S = 1, and one state with S = 2. These total-spin multiplicities, denoted here by G(N = 4, S), arise from the group-symmetry properties of the spin eigenfunctions [12,31] of N = 4 fermions with spin 1/2; for the multiplicities G(N, S) of total-spin degeneracies for any N fermions, see the branching diagram [31] in Appendix A. (The theory of spin-1/2 eigenfunctions is well known in quantum chemistry (see Ref. [31]) and has been used [12] previously in the field of quantum dots, ant it will not be repeated here. However, for a brief outline and a description of the general spin eigenfunctions for N = 3 and N = 4, see Appendix A. Importantly, in the cases studied in this paper [that is, for N = 4 (Figs. 1 and 2), and for N = 3 in a SW and in the DWPA trap ( Fig. 4), as well as in a DWLI trap ( Fig.  5)] the AFM lowest spin-state is the ground state and the energy level spacings decrease with increasing interatomic repulsion.
The similar behavior of the sixfold-multiplet bands irrespective of the different geometries of the double-well traps, i.e., DWPA versus DWLI, indicates an underlying physical process independent of dimensionality (2D versus 1D). This underlying physics involves spatial localization of the 6 Li atoms at extemporaneously created sites within each well and the ensuing formation of quantum UCWMs, as can be seen by an inspection of corresponding single-particle densities (SPDs, green surfaces in Figs. 1 and 2) and spin-resolved conditional probability distributions (SR-CPDs, angle-resolved pair correlations, red surfaces in Figs. 1 and 2).
The SPD is the expectation value of the one-body operator where |Φ CI N denotes the many-body (multi-determinantal) CI wave function. We note that the SPD is the sum of the spin-up and spin-down single-particle densities, defined as where σ and σ i denote up or down spins. In all cases the SPDs display four humps corresponding to the four localized fermions at the self-generated localization sites. The detailed arrangement of these sites varies in order to accomodate the geometry of the traps. For the DWPA case [ Fig.  1(f)] with two fermions in the left well and the other two in the right well (n L = 2, n R = 2), a 2D rectangle is formed. For the DWLI (2,2) case [ Fig. 2(c)], including the limiting case of the single well [ Fig. 1(b)], the four sites fall onto a straight line. Note the opening in the middle of the DWLI density [ Fig. 2(c)], in contrast to the case of the single well in Fig. 1 Although several distinct spin structures can correspond to the same SPD of a UCWM, the spin eigenfunction associated with a specific CI wave function can be determined with the help of the many-body SR-CPDs, P σσ 0 , which yield the conditional probability distribution of finding another fermion with up (or down) spin σ at a position r, given that a specific fermion with up (or down) spin σ 0 is fixed at r 0 . In detail, the spin-resolved two-point anisotropic correlation function is defined as Using a normalization constant we further define a related spin-resolved conditional probability distribution (SR-CPD) as In particular, by calculating the ratios of the volumes under the CPD humps and equating them to the corresponding ratios of the squares of the angle-dependent coefficients of the general expressions for the spin eigenfunctions, one can determine the numerical values of the coefficients that map the spin eigenfunction to a specific SR-CPD (for details, see Appendix A and Ref. [12]). As an example, the spin eigenfunction associated with the 4-fermion S = 0, S z = 0 CI ground state at −1/g = −0.1( √ 2l 0y ω y ) −1 in the case of well-separated DWPA parallel wells [see star in Fig. 1(e); for the corresponding SR-CPD, see Fig. 1 where the α's (β's) denote up (down) spin-1/2 fermions situated at the self-generated sites (the maxima of the humps in the SPDs or CPDs); the methodology and detailed calculations used in determining the angle θ in the general spin eigenfunction in Eq. (A.2) are described in Appendix A.
The equal in absolute value |C i | = 1/2, i = 1, . . . , 4 coefficients in front of the four primitives in Eq. (10) agree with the probability of 0.5 [i.e., 0.5 = 2 × C 2 ] for the socalled "antiferromagnetic" component (ααββ and ββαα) found in Ref. [19] [see Fig. 1(d) therein] for the case of a two-parabola DWLI double well in the high-barrier regime. They also agree with the probability for the "mixed" component (αβαβ and βαβα) reported in the same paper. We note that in our treatment, we can vary the height of the barrier independently from the separation of the wells, unlike the case in Ref. [19]. The use of the terms "antiferromagnetic", "mixed", and "ferromagnetic" to characterize the spin primitives of the Ising basis is borrowed here and in a paragraph below from Ref. [19] in order to facilitate the comparisons. This use is not repeated anywhere else in the paper; instead, as aforementioned, we employ the term "antiferromagnetic" to describe finite systems that have ground states with the lowest possible total spin.
The mapping to the spin eigenfunction in Eq. (10) reflects the fact that at the highbarrier (or large-separation) regime the 4-fermion problem can be viewed as that of two pairs of strongly interacting fermions within each well, each pair interacting weakly with the other one through the high barrier. In this case, as discussed below, the energetics of the 4-fermion system can be understood simply by adding the singlet and triplet energy levels of the left and right fermionic pairs. However, the CI wave functions exhibit strong entanglement between the left-and right-well fermionic pairs in addition to the entanglement between the two fermions within each well. This across-the-barrier entanglement is not weakening as a result of a higher barrier, and it is manifested in the mapping of the CI ground-state wave function onto the spin eigenfunction in Eq. (10).
Furthermore, the discussion above applies also to the excited states. For example, the SPD and SR-CPD of the first excited state with S = 0, S z = 0 in the DWPA trap of Fig. 1 [having an energy E = 2 ω y in Fig. 1(e)] is displayed in Figs. 1(h) and 1(i), repectively. For this case, following an analysis as described above (and in Appendix A), we find an angle θ = π/6, which is associated with a spin function of the form We note that the spin eigenfunctions in Eqs. (10) and (11) are orthogonal.
The two coefficients C 1 = C 2 = 1/(2 √ 3) in front of the first two primitives in Eq. (11) agree with the probability of 0.166 (i.e., 0.166 = 2 × C 1 2 ) found in Ref. [19] for the "antiferromagnetic" (ααββ and ββαα), as well as for the "mixed" (αβαβ and βαβα) primitives in the case of a two-parabola DWLI double well at the high-barrier regime.
The spin eigenfunction associated with the 4-fermion S = 0, S z = 0 CI ground  10) and (11). Specifically, the analysis of the SR-CPD described in detail in Appendix A yields an angle θ = −π/5.12 in Eq. (A.2), which is associated with the following spin eigenfunction where C 1 = 0.332411, C 2 = −0.575017, and C 3 = 0.242606. In the next section, we utilize the trends uncovered by the CI solutions for the spectra and wave functions of the many-body Hamiltonian, in order to develop a Heisenberg-model phenomenology. This development aims at providing tools for analyzing quantum magnetism in double (and multi-well) ultracold-atom traps.

Four fermionic ultracold atoms: The Heisenberg model
We have verified that the CI energy spectra presented in Figs. 1 and 2, as well as the SR-CPD-derived spin eigenfunctions [see, e.g., the functions in Eqs. (10), (11) and (12) where the symbol < ij > denotes that the summation is restricted to the nearestneighbor sites. The second term is a scalar, leading simply to an overall energy shift; for a detailed description of H H , see Appendix B and Appendix C. The DWPA case is associated with a rectangular 2D Heisenberg ring [see schematic (III) in Fig. 1 ] are given in Appendix B and Appendix C. They can reproduce all the trends in the energy spectra of the sixfold energy band, as well as the total-spin multiplicities G(N = 4, S) and spin eigenfunctions calculated via the CI method. In particular, in the limit of well-separated wells (i.e., for r = 0), one gets E 2 = E 4 = E 6 = 0, E 1 = E 3 = −s, and E 5 = −2s, which coincides with the aforementioned 1-2-3 spin-group-theoretical degeneracy pattern and relative gaps within the sixfold lowest-energy CI band. Note further that the Heisenberg modeling reproduces the two different SR-CPD-derived spin eigenfunctions in Eqs. (10) and (11), associated with the fully-separated-wells (r = 0, for both the DWPA and DWLI cases); compare with the eigenvectors in Eqs. (B.22) and (C.12). It is notable that both the CI spectra [see Figs. 1(e) and 2(b)] and the Heisenberg energies for fully separated wells exhibit two energy gaps, one twice as large as the other (e.g., −s and −2s in the Heisenberg model). This behavior can be understood from the spectrum of two unrelated single wells each containing a pair of two strongly interacting fermions. Indeed, the two lowest levels of two interacting fermions consist of a singlet state with energy E s and a triplet state with energy E t . The low-energy spectrum of the double well has then three levels, E 1 = 2E t , E 2 = E t + E s , and E 3 = 2E s , corresponding to whether both fermion pairs are in a triplet state, one pair is in a triplet with the other in a singlet state, or both pairs are in a singlet state; this results in the two energy gaps ∆E 12 = E t − E s and ∆E 13 The topology of the spin chain in Fig. 1(III) (DWPA) is indeed a closed ring, whereas the one in Fig. 2 For the single-well case, all six CI energies have distinct values; see the spectrum in Fig. 1(a). By using the open-Heisenberg-chain eigenvalues E i , i = 1, . . . , 6 in Eqs. (C.2-C.7) and fitting the ratios (E 4 − E i )/(E 4 − E j ) to the CI spectrum, we can determine the parameter f = r/s that describes the single well. For example, using the fully polarized, E f p = E 4 = 0, the ground-state, E gs = E 5 , and the 1st-excited, E 1st = E 3 , energies, we obtain the ratio which is independent of s and allows for the determination of f . Fitting to the CI spectrum, we get f ∼ 1.35. This value agrees with that resulting from the nearestneighbor exchange constants of harmonically trapped particles listed in Table I of Ref. [18]. Another study [20] gave a value of ≈ 1.4 for this ratio.
With the value of r/s = 1. 35, the open linear Heisenberg chain yields E 1 = −s (S = 1), E 2 = −0.334985s (S = 1), E 3 = −2.01501s (S = 1), E 4 = 0 (S = 2), E 5 = −2.55853s (S = 0), and E 6 = −0.79147s (S = 0), i.e., six distinct values, in agreement with the CI spectrum in Fig. 1(a). The corresponding angle in Eq. (A.2) is θ = −π/4.58. This value is slightly different from the value of −π/5.12 (corresponding to an r/s ≈ 1.62) that was determined above in Sec. 3.1 from an analysis of the CI CPD in Fig. 1(a). This slight discrepancy is due to the elimination of the space degrees of freedom when considering the mapping of the CI wave function onto the spin eigenfunctions. Naturally, the spin eigenfunctions have constant coefficients in front of the Ising-expansion primitives and by themselves are unable to reflect the influence of the extent of space distribution of the localized fermions. Indeed the localization of the four fermions is sharper in a double well with a high barrier compared to that in a single well; compare the SPD's in Figs. 1(b) (4 fermions inside the same well) and 1(f) (2 fermions in each well). The CI SPDs and SR-CPDs incorporate in their definition the space degrees of freedom and they account for the actual extent of partial or full particle localization (which varies with g). A detailed investigation of this matter is beyond the scope of this paper, but it will be examined in a future publication [52]. The von Neumann entanglement entropies, calculated from the single-particle density matrix [11,26]. Note the increased or constant entanglement with increasing repulsion, and the larger values for the symmetric DW congurations (e,i) compared to the nonsymmetric (tilted) ones (m,q). The zero of energy in the spectra corresponds to the ground-state total energy of the corresponding non-interacting system, that is, to 11.14 ω y in (a), 12.62 ω y in (f), 13.11 ω y in (j), and 13.62 ω y in (n). The presence of the interwell barrier in (n) accounts for the difference from the single-well value of (5 ω y +3 ω x )/2 = 12.40 ω y .

Three fermionic ultracold atoms: CI results and the Heisenberg model for tilted double wells
CI results for N = 3 ultracold 6 Li atoms in a DWPA trap are displayed in Fig. 4 for both symmetric [zero tilt, ∆ = 0, see schematic in Fig. 4(II) and spectra in Fig. 4(a) and Fig. 4(f)] and asymmetric wells with a moderate tilt ∆ = 0.5 ω y [see Fig. 4(IV) and Fig. 4(j)] and a strong tilt ∆ = 2.5 ω y [see Fig. 4(VI) and Fig. 4(n)].
The cases of asymmetric wells are amenable to straightforward interpretations based on pure Heisenberg models. The moderate tilt [ Fig. 4(IV), ∆ = 0.5 ω y ] generates a ground state with a (2,1) distribution of the atoms (two in the left well and one in the right, tilted upward, one), which are localized in the shape of a isosceles triangular UCWM [see the SPD in Fig. 4(k)]. The corresponding CI energy spectrum [ Fig. 4(j)] exhibits a three-fold lowest-energy band with a characteristic 1-2 degeneracy pattern, converging to the same energy for −1/g → 0. The total-spin multiplicities in this band are G(N = 3, S = 1/2) = 2 and G(N = 3, S = 3/2) = 1, in agreement with the branching diagram for three fermions (see Appendix A). This CI energy spectrum and the correponding SR-CPDs [see, e.g., the SR-CPD in Fig. 4(l)] are reproduced by a 3-site Heisenberg-ring Hamiltonian with J 12 = s and J 13 = J 23 = r; for the numbering of the three sites, see the schematics in Fig. 4 6). We note that, while the spin spatial distribution is analyzed here with the use of the SR-CPDs [see Eq. (7)], it is also reflected in the spatial spin-densities [see Eq. (6)] shown in Figs. 6(k-l); the latter agree with those displayed in Fig. 6 of Ref. [18]. Note that the sum of the up-and down-spin densities in Figs. 6(k-l) agrees with the total SPD in Fig. 6(j).
A qualitatively different behavior, bringing extra intricacies and opening igress to novel complex physical systems, is exhibited by the symmetric DWPA cases (∆ = 0) for N = 3 shown in Fig. 4. Indeed, the CI energy spectra in Fig. 4(a) and Fig. 4(f) show a sixfold lowest-energy band, comprising four S = 1/2 states, and two S = 3/2 states, i.e., twice as many as in the case of tilted wells [Figs. 4(j) and 4(n)]. In particular, for the higher barrier [ Fig. 4(f)] a characteristic 2-4 degeneracy appears, which is a doubling of are calculated from the single-particle density matrix [11,26]. This figure complements Fig. 4 in that it displays examples of all possible SR-CPDs for the ground state of N = 3 cold fermions in a single 1D well. It is straightforward to check that the SR-CPDs agree with a mapping of the CI ground state onto the spin eigenfunction of the schematic (II). (j-l) Ground-state results for N = 3 6 Li atoms in a strictly-1D single trap with ω x = 100 kHz and ω y = 1 kHz. (j) SPD (green surface). (k) spin-down density (red surface). (l) spin-up density (blue surface). The zero of energy in the spectrum corresponds to the ground-state total energy of the corresponding non-interacting system, that is, to 13.62 ω y in (a). the 1-2 degeneracy pattern in Fig. 4(j). This doubling of the number of energies is due to the conservation of parity, which requires consideration of a second triangle (246), which is the mirror of the original (135) one; see the schematic in Fig. 4(I) and in Fig.  7(a). In each of these mirror reflected configurations, two atoms localize in one well and one atom localizes in the other well; see the two sets of different colored spheres in Fig. 4(I). The formation of these triangular atomic configurations is reflected in the SR-CPDs shown on Fig. 4(c,d) for the lower-barrier symmetric DW case and Fig. 4(g,h) for the higher-barrier case. One may view this situation as having six available sites altogether (three in each well), with the 3 fermionic atoms localizing in either of the aforementioned triangular configurations, (135) and (246) [see Fig. 7(a)], with 2 atoms in one well and 1 atom in the other; in each case we may term the unoccupied (empty) sites as "holes". This mapping leads to the picture of a 3-atom UCWM that resonates between the two interlocking triangles.
We mention that the resonating behavior and the symmetrization of the many-body wave function in two-center/three-electron bonded systems is well known [33,34,35] in theoretical chemistry and in particular in the valence-bond treatment of the threeelectron bond which controls the formation of molecules like He + 2 and F − 2 . Furthermore we mention that the symmetry properties of the strictly-1D few-fermion problem with contact interactions have been also investigated in Refs. [53,54]

Three fermionic ultracold atoms: The t-J model for symmetric double wells
To model the exact-diagonalization results shown above, one must go beyond the aforementioned simple Heisenberg Hamiltonian model [see Eq. (15) and Appendix D]. Indeed, we find that a generalization of the so called t-J model allows us to capture all the salient characteristics uncovered by the CI calculations. The t-J model [5,44] modifies (away from the half filling) the antiferromagnetic Heisenberg Hamiltonian associated with the Mott insulator at half-filling (one electron per crystal site); it has attracted much attention, because it has been proposed for explaining the high-T c superconductivity arising in the case of underdoped insulators (away from the half filling when holes are present). A finite t-J-type Hamiltonian may be expressed as where H c is the coupling between the two simple Heisenberg rings defined over the sites (135) and (246); see the two 3 × 3 blocks on the diagonal (upper left and lower right) in Eq. (18). H c , represented by the two off-diagonal blocks in Eq. (18), is defined by the matrix t = (α0α0β0|H c |0α0α0β) = (α0α0β0|H c |0α0β0α), and t 2 = (α0α0β0|H c |0β0α0α), where the "0" indicates an empty site; e.g., α0α0β0 corresponds to a state where sites 1,3, and 5 are occupied and 2,4, and 6 are empty (for site designation see Fig. 4(I) and Fig. 7). The Hamiltonian H tJ is equivalent to a six by six matrix, we will focus on the case of symmetric double wells, i.e., we will set ∆ = 0, J 15 = J 24 = r, and J 13 = J 46 = s. For lower values of the interwell barrier [V b = 11.14 kHz, Fig. 4(a)], the 2-4 [(1-2) ×2] doubling degeneracy is lifted, with two lowest S = 1/2 curves and four higher in energy (and parallel) curves (two with S = 1/2 and two with S = 3/2) forming distinct subbands; then one distinguishes all 6 lines as separate lines [see the spectrum in Fig.  4(a) and also in Fig. 5(a)]. It is remarkable that the nontrivial spectrum in Fig. 4(a) can be reproduced by setting t 2 ∼ −4t/10 > −1/2, with t < 0 and s > |t|. Then one has for the energy gap between the two lowest states, ∆E 12 = E 1 − E 2 = 14t/5; these energies are centered around −s. The remaining energies group together forming a fourfold band, centered around zero. The energy gap between the two outer (both S = 3/2) members of the fourfold band is ∆E 56 = E 5 − E 6 = 16t/5, i.e., similar to the ∆E 12 gap, in agreement again with the pattern in Fig. 4(a). Furthermore, the gap between the two higher energies in the fourfold band, as well as that between the two lower energies of this band, is ∆E 35 = ∆E 46 = E 4 − E 6 = t/5, which is much smaller than the width, ∆E 56 , of the same band, again in agreement with the pattern in Fig. 4(a). Note that for t 2 = −t/2, ∆E 35 = ∆E 46 = 0 and a degeneracy pattern 1-1-2-2 develops in disagreement with the CI spectrum. Also, when t 2 = 0, the width ∆E 56 (= 4t) of the fourfold band is twice as large as the energy gap, ∆E 12 (= 2t), between the two lowest states, again in disageement with the CI spectrum in Fig. 4(a).
Similar trends pertaining to the doubling of the spectrum (from three to six states) in conjunction with the emergence of two resonating UCWMs apply also in the case of N = 3 ultracold fermions in a symmetric (∆ = 0) DWLI (linear arrangement) trap, as is illustrated in Fig. 5. These results can be interpreted again through the use of a corresponding t-J model with a similar parametrization.

Quantifying entanglement using a CI-based von Neumann entropy
The entanglement entropy S vN for three 6 Li atoms in a DWPA trap in the configurations, whose spectra are shown in Fig. 4(a,f,j,n), are displayed in Fig. 4(e,i,m,q), respectively.
For the CI many-body wave functions, we adopt as a measure of entanglement the von Neumann entropy [11,26], where ρ is the single-particle density matrix and C = − log 2 N , yielding S vN = 0 for an uncorrelated single-determinant state. The single-particle density matrix ρ is given by and it is normalized to unity, i.e., Tr ρ = 1. The Greek indices µ (or ν) count the spin orbitals (of dimension 2K) and where α(β) denote up (down) spins. Since the allowed maximum value for S vN in our CI calculations is log 2 (2K) − log 2 (3) = 5.70 (we use a typical basis of K = 78 single-particle space orbitals), it is notable that the calculated values in Fig. 4 remain smaller than ∼ 1, and in particular in the regime of strong correlations, i.e., for −1/g → −0. This reflects formation of a Wigner molecule. Additionally, we find increased or constant entanglement with increasing repulsion, and larger values for the symmetric DW configurations [ Fig. 4(e,i)] compared to the nonsymmetric (tilted) ones [ Fig. 4(m,q)]. For S vN entropies for three 6 Li atoms in a DWLI trap, see Fig. 6(a).

Summary and Outlook
In this paper, we have presented timely advances in the growing field of few-body ultracold atoms with the aim of enhancing understanding of experimental endeavors and lodging new directions of research in this area. We progressed in two main courses: (i) uncovering universal non-itinerant and fermionization-like aspects of the physics of ultracold few fermions trapped in double-well confinements, with various 1D and 2D trapping geometries, as a conduit for emulating quantum magnetism and related phenomena beyond the strictly 1D single-well (SW) case, and (ii) making headways in the development and implementation of benchmark numerical simulations (exact diagonalization of the full microscopic Hamiltonian with configuration interaction, CI, techniques) as tools for modeling theoretical and experimental results with effective spin-Hamiltonians (Heisenberg and t-J models). Our calculations for N = 3 and N = 4 ultracold fermionic 6 Li atoms in SW and double well (DW) traps with linear (DWLI) or parallel (DWPA) geometries, reveal formation of antiferromagnetic ordering for the lowest-energy bands over the entire range of interparticle contact repulsion studied here.
For N = 4 ultracold atoms in a symmetric DWPA trap with very strong interatomic repulsion, we find (via miscroscopic, CI, calculations) formation of a two-dimensional ultracold Wigner molecule (UCWM) of non-itinerant character. For the symmetric parallel DW trap the formation of the 2D UCWM leads to mapping of the interacting 4-atom trapped system onto a 2D rectangular Heisenberg ring cluster, whereas for a symmetric DWLI trap (as well as for a SW trap) we find a four-atom linear (1D) UCWM in juxtaposition with mapping onto a linear Heisenberg spin-chain. These mappings enable employment of the corresponding Heisenberg model Hamiltonian, whose solutions reproduce well the results of the microscopic, numerically-exact, calculations.
For N = 3 ultracold atoms in DWLI or DWPA traps with a finite tilt (detuning) between the two wells, the numerically calculated (CI) spectrum for strong interatomic repulsion is described well with the use of the aforementioned Heisenberg Hamiltonian. As noted already in the Introduction, the high measure of entanglement predicted for the set of lowest energy states of the three strongly repelling fermionic atoms, together with the controllable tilt between the two wells, motivate consideration of this double well system as a cold-atom quantum computing qubit.
In contrast to the asymmetric DW case, description of the N = 3 ultracold-atom CI spectra for symmetric (vanishing tilt) DWLI or DWPA traps, that manifest doubling of the number of states in the lowest band, as well as modeling the corresponding SR-CPDs, are not attainable with the simple Heisenberg model, requiring instead the more intricate t-J-type model [5,6,44], consisting of two coupled resonating triangular 2D UCWM Heisenberg clusters. The emergence of the t-J model for the description of quantum magnetism (in particular AFM ordering) in a trapped few-body ultracold atom system, strongly suggests its future role as a useful laboratory for exploration of the elementary building blocks of high-T c superconducting behavior [44,55,56]. CI CPDs We outline in this Appendix several properties of the many-body spin eigenfunctions which are useful for analyzing the trends and behavior of the spin multiplicities exhibited by the CI wave functions for N = 4 and N = 3 ultracold fermions. The spin multiplicities of the CI wave functions lead naturally to analogies with finite Heisenberg clusters [8,10] and to t-J-type models.
A basic property of spin eigenfunctions is that they exhibit degeneracies for N > 2, i.e., there may be more than one linearly independent (and orthogonal) spin functions that are simultaneous eigenstates of bothŜ 2 and S z . These degeneracies are usually visualized by means of the branching diagram [31] displayed in Fig. A1. The axes in this plot describe the number N of fermions (horizontal axis) and the quantum number S of the total spin (vertical axis). At each point (N, S), a circle is drawn containing the number G(N, S) which gives the degeneracy of spin states. It is found [31] that Specifically for N = 4 particles, there is one spin eigenfunction with S = 2, three with S = 1, and two with S = 0. In general the spin part of the CI wave functions involves a linear superposition over all the degenerate spin eigenfunctions for a given S.
For a small number of particles, one can find compact expressions that encompass all possible superpositions. For example, for N = 4 and S = 0, S z = 0 one has: [12,57]   where the parameter θ satisfies −π/2 ≤ θ ≤ π/2 and is chosen such that θ = 0 corresponds to the spin function with intermediate two-fermion spin S 12 = 0 and threefermion spin S 123 = 1/2; whereas θ = ±π/2 corresponds to the one with intermediate spins S 12 = 1 and S 123 = 1/2. For N = 3 and S = 1/2, S z = 1/2 one has [57]: For the general expressions for the remaining spin combinations, S and S z , for N = 4 and N = 3 fermions, see Refs. [12,57,58]. For each SPD corresponding to a given CI state of the system, one can plot four different spin-resolved CPDs, i.e., P ↑↑ , P ↑↓ , P ↓↑ , and P ↓↓ . This can potentially lead to a very large number of time consuming computations and an excessive number of plots. For studying the spin structure of the S = 0, S z = 0 states for N = 4 fermions and the S = 1/2, S z = 1/2 states for N = 3, however, we found that knowledge of a single CPD, is sufficient in the regime of Wigner-molecule formation. Indeed, the specific angle θ specifying the spin function X  Fig. 1(g). Integrating numerically under the humps of the CI CPD in Fig. 1(g), we specify the ratio x =Vol(4)/[Vol (2) Fig. 1(a,b,c). Note the labeling of the four sites in space, which is "4123" and not "1234". This results from our taking J 34 = 0, when opening the four-site ring, and it is consistent with our treatment of the four-site linear Heisenberg chain in Appendix C below.
Noting that hump No. 4 in Fig. 1(c) is again well isolated from the rest, and focussing on the numbering of the remaining humps of this SR-CPD, it is apparent that we need to use the same set of the quantities Π ↓↑ (i, 1), i = 2, 3, 4 as was the case with the previous example. Integrating under the humps of the CI CPD in Fig. 1(c), we find the numerical values for the volumes Vol(i), i = 2, 3, 4. In particular, we determine that x = 0.789. With this value of the ratio x, condition (A.7) yields an angle of θ = −π/5.12.
Example 3; case of the CPD in Fig. 1(i): As a third example, we choose an excited state (the one with S = 0 and S z = 0) in the double well. Illustrative calculations for the spectrum, densities, and CPDs for this case are displayed in Fig. 1(e,h,i). Note again the labeling of the four sites in space. Noting that hump No. 4 in Fig. 1(i) is again well isolated from the rest, and focussing on the numbering of the remaining humps of this SR-CPD, it is apparent that we need to use the same set of the quantities Π ↓↑ (i, 1), i = 2, 3, 4 as was the case with the previous examples. Integrating under the humps of the CI CPD in Fig. 1(i), we find the numerical values for the volumes Vol(i), i = 2, 3, 4. In particular, we determine that x = 0.20 in this case. With this value of the ratio x, condition (A.7) yields an angle of θ = π/6.
For further detailed applications of this procedure, see Refs. [12,57]. where the indices k in S k denote the locations of the four sites, which are associated with the four humps in the s.p. density of Fig. 1 (in a clockwise direction); see also schematic in Fig. B1(a).
For the case of all four fermions being trapped in a single well, one has an open linear 4-site Heisenberg chain, which is obtained from Eq. (B.1) by setting J 34 = 0, i.e., The general eigenvalues E i and corresponding eigenvectors V i of the matrix (B.5) are calculated easily using MATHEMATICA [59]. The eigenvalues are: where The corresponding (unnormalized) eigenvectors and their total spins are given by: and f = r/s. To understand how the Heisenberg Hamiltonian in Eq. (B.5) captures the behavior seen in the CI spectra of Fig. 1 (DWPA case), we start with the limiting case r → 0, which is applicable (see below) to the larger interwell barrier V b = 11.14 kHz. In this limit, one can neglect r compared with s, which results in a characteristic 1-2-3 degeneracy pattern within the band; namely one has E 2 = E 4 = E 6 = 0, E 1 = E 3 = −s, and E 5 = −2s.
Furthermore, the fact that all six curves in the CI lowest-energy band cross at the same point 1/g = 0 suggests that s ∼ F (−1/g) and r ∼ F (−1/g) with (x = −1/g) Of interest is the fact that the ability of the Heisenberg Hamiltonian in Eq. (B.5) to reproduce the CI trends is not restricted solely to energy spectra, but extends to the CI wave functions as well. Indeed when r → 0, the last two eigenvectors of the Heisenberg matrix (having S = 0) become The CI spectra and spin functions for the smaller barrier V b = 6.08 kHz can be analyzed within the framework of the 4-site Heisenberg Hamiltonian (B.5) when small (compared with J 14 = s), but nonnegligible, values of the second exchange integral J 12 = r are considered. In this case, the partial three-fold and two-fold degeneracies are lifted. Indeed in Figs. 1(d) (V b = 6.08 kHz), the CI lowest-energy band consists of six distinct levels. The corresponding (unnormalized) eigenvectors and their total spins are given by: and f = r/s as previously defined.
In the limit of r → 0 (high interwell barrier V b ), the energies in Eqs. (C.2)-(C.7) reproduce the characteristic 1-2-3 degeneracy pattern, which appears also in the case of the rectangular arrangement of the four sites; namely one has E 2 = E 4 = E 6 = 0, E 1 = E 3 = −s, and E 5 = −2s. Furthermore, for r → 0, the corresponding (unnormalized) eigenvectors and their total spins are given by:

Appendix D. Heisenberg model for 3 localized fermions in tilted wells
In the case of N = 3 strongly-interacting fermions in a single 1D well or a tilted doublewell with a parallel arrangement of the two 1D wells, the simple Heisenberg model is applicable. For a (2, 1) fermion configuration, the three sites form an isosceles triangle (see Figs. 4 and C1), and the associated Heisenberg-ring Hamiltonian H trg H is given by Eq. (15). To proceed, we use the three-dimensional Ising Hilbert subspace for total-spin projection S z = 1/2, which is spanned by the following set of basis states: |1 >→ ααβ, |2 >→ αβα, and |3 >→ βαα. In this subspace, the complete Heisenberg Hamiltonian in Eq. Note that the eigenvectors are independent of s and r, however which one is the ground state depends on these exchange constants through the expressions for the eigenvalues E i given in Eqs. (D.2)-(D.4). In particular, when the interwell barrier is high (r → 0) [see Fig. 4(j)] a characteristic 1-2 degeneracy develops with E 1 = E 2 = 0 and E 3 = −s. When s > 0, the ground-state vector is given by V 3 in Eq. (D.7).
The case of 3 fermions in a single well [forming a linear Wigner molecule, see Fig.  4(VI)] is described by the matrix Hamiltonian (D.1) when s = 0 (open Heisenberg chain). Then all three eigenvalues are different with E 1 = 0, E 2 = −3r/2, and E 3 = −r/2 [see Fig. 4(n)]. Thus, with r > 0, the ground-state vector for the (3,0) fermion arrangement is given by V 2 in Eq. (D.6) and is different from that of the (2,1) fermion arrangement, although the total spin remains the same, i.e., S = 1/2 [compare SR-CPDs in Figs. 4(p,l)]. The t-J Hamiltonian matrix in Eq. (E.1) exhibits a very rich behavior. In the following, we will limit our analysis to the case with r = 0, i.e., for large interwell barrier V b which is also the case of both spectra in Fig. 4(a) and 4(f). In this limit, the eigenvalues and eigenvectors (unnormalized) of the Hamiltonian matrix (E.1) are: