Ultracold few fermionic atoms in needleshaped double wells: spin chains and resonating spin clusters from microscopic Hamiltonians emulated via antiferromagnetic Heisenberg and  models
Abstract
Advances with trapped ultracold atoms intensified interest in simulating complex physical phenomena, including quantum magnetism and transitions from itinerant to nonitinerant 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 onedimensional, 1D, wells, and (ii) two coupled parallel wells, forming a trap of twodimensional, 2D, nature. The spectra and spinresolved conditional probabilities reveal for both cases, under strong repulsion, atomic spatial localization at extemporaneously created sites, forming quantum molecular magnetic structures with nonitinerant character. These findings usher future theoretical and experimental explorations into the highlycorrelated behavior of ultracold stronglyrepelling 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 spinchain 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 like model, akin to that used in investigations of high T superconductivity. The highly entangled states discovered in the microscopic and model calculations of controllably detuned, asymmetric, double wells suggest threecoldatom DW quantum computing qubits.

19 May 2016
Published: New J. Phys. 18, 073018 (2016)
1 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 finitesize [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 and ultracold Li atoms confined in a singlewell onedimensional optical trap [13]. Moreover, progress aiming at bottomup approaches to fermionic manybody systems, addressing entanglement, quantum information, and quantum magnetism in particular, is predicated on experimental developments of which the recently created doublewell (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 onedimensional systems trapped within a single well (see, however, Refs. [19, 20] for doublewell 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 manybody wave functions. However, the recently demonstrated ability to create needleshaped double well traps [15], and the anticipated nearfuture further development of small arrays of such needlelike quasi1D 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 LippmanSchwinger equation approach, see Refs. [27, 28] and [29], respectively. Similarly, states of ultracold fermions in a single strictly1D trap, away from the fermionization limit, using an exact diagonalization of the manybody Hamiltonian have been reported [30].
In this paper, using largescale configurationinteraction (CI) calculations as means for exact diagonalization of the microscopic Hamiltonian, we report that for (even number) strongly interacting ultracold fermions in a doublewell trap with parallel arrangement [DWPA [26]; see Figs. 1(II,III)] the manybody problem can be reduced to that of a 2D rectangular AFM Heisenberg ring. The associated mapping between the manybody wave function and the spin eigenfunctions [31] for and 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 longrange Coulombic repulsion. Such molecular structures are usually referred to as Wigner molecules (WMs) [32]. For (odd number) ultracold fermions, fewbody quantum magnetism requires introduction of a more complex type [5, 6] model; here the  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 manybody wave function in twocenter/threeelectron bonded systems is well known [33, 34, 35] in theoretical chemistry and in particular in the valencebond treatment of the threeelectron bond which controls the formation of molecules like He, NO, and F. The concept of the threeelecton 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 stronglyinteracting, highly correlated fermionic atoms and the formation [26] of quantum 1D and 2D moleculelike 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 fewfermion molecular structures irrespective of the nature of the repelling twobody interaction, i.e., contact versus longrange 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 highdegree 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 wellknown relevant fields are: (i) the socalled direct exchange mechanism [3, 4] (related to ferromagnetism in the mixedvalency manganites of perovskite structure), and (ii) the  model [5, 44] which modifies (away from the half filling) the antiferromagnetic Heisenberg Hamiltonian associated with the Mott insulator at halffilling. The  model has attracted much attention, because it has been proposed for explaining the highT superconductivity arising in the case of underdoped insulators [44]. Due to the antiferromagnetic aspect, our resonating model Hamiltonian for fermions (see Section 4.2 below) represents a finite variant of the  model. The emergence of the  model in this work suggests future investigations of fundamental aspects associated with the physics of highT 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 doublewell trap with linear arrangement of the two needlelike 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 Wignermolecule (UCWM) formation and of the associated simpleHeisenbergchain and  resonatingspinchains 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 meanfield approaches (e.g. HartreeFock), 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 threeelectron system in semiconductor double quantum dots has recently attracted a major attention in conjunction with the fabrication and implementation of pulsegated fast hybrid qubits for solidstatebased quantum computing [47, 48]. These advances and the fascinating physics of doublewelltrapped 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 3atom 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., for fermions and for fermions.
The plan of the paper is as follows:
A statement of the manybody Hamiltonian, including a description of the doublewell employed by us, is given in Section 2.
In Section 3, we describe investigations concerning four ultracold Li atoms in doublewell traps with both parallel and linear arrangement of the individual needlelike wells. A comparison with the case of four fermions in a quasi1D 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 4fermion phenomenological Hamiltonian.
Section 4 presents our studies concerning the case of three ultracold 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 3fermion Heisenberg model. Going beyond the Heisenberg model, Sec. 4.2 introduces the 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 manybody wave functions.
The appendices provide detailed information concerning the mathematical formalism associated with the spin eigenfunctions and the finite Heisenberg and models. In particular,
A provides a brief description of the branching diagram that describes the multiplicities (number of degenerate spin states) of a given total spin . This Appendix also presents in the Ising basis the general formulas that describe a spin eigenfunction (i) with and for four fermions and (ii) with and 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 , as well as any linear superposition of them. In addition, A describes the process of mapping the manybody CI wave functions onto these spin eigenfunctions.
B discusses the mathematics of the Heisenberg model for 4 localized fermions in a ringlike rectangular configuration (DWPA case), while C discusses the corresponding case for an open chain arrangement (DWLI case).
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, E discusses the mathematics of the more general model for 3 localized fermions in the case of a double trap with symmetric wells.
2 ManyBody Hamiltonian
A DWLI trap can be treated as a strictly1D problem along a single direction. The DWPA trap which consists of two parallel needlelike wells, however, cannot be treated solely along one direction (e.g., the direction). Instead it requires consideration of the coordinate as well. To treat both cases in a unified way, we consider a manybody Hamiltonian for fermions of the form
(1) 
where denotes the relative vector distance between the and fermions (e.g., Li atoms). This Hamiltonian is the sum of a singleparticle part , which guarantees the needlelike shape of the individual wells, and the twoparticle contact interaction.
The external confining potential [in ] that models a double well (DW) is based on a twocenteroscillator (TCO) model [12, 26] exhibiting a variable smooth neck along the direction. Along the direction, this TCO model allows for an independent variation of both the separation and of the barrier height between the two wells; see Fig. 3. Along the direction, the confinement consists of that of a single harmonic oscillator. The values of the frequencies (left well), (right well) and that confine the two wells along the and directions, respectively, are also allowed to vary independently; here we choose . The needlelike shape of each individual trap is enforced by assuming that (DWLI case) or (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 direction used in Fig. 4 below for the study of fermions in tilted and symmetric double wells.
As we mentioned in the Introduction, we use the CI method for determining the solutions of the manybody problem specified by the Hamiltonian (1). The CI method expresses the fermionic manybody wave function as a supperposition of Slater determinants, and it is well known in quantum chemistry and in fewbody 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 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 singleparticle states as needed. Note that the TCO singleparticle states automatically adjust to the separation as it varies from the limit of the unified atom to that of the two fully separated traps (for sufficiently large ). We verified that for (strictly1D single trap), our CI calculations agree with the results of Table 2 of Ref. [30].
The matrix elements of between the CI determinants are calculated using the Slater rules. An important ingredient in this respect are the twobody matrix elements of the contact interaction
(2) 
in the basis (of dimension ) formed out of the singleparticle (space) eigenstates of the TCO Hamiltonian.
Because the individual wells remain needlelike in all of our calculations here, the wave scattering between two ultracold fermions takes place primarily along a single dimension, either the dimension or the dimension. As a result, the parameter in front of the function in Eq. (1) does not reflect the physical process of twodimensional 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, , are related to as follows
(3) 
where is a dummy variable and is the lowestinenergy singleparticle state in the direction for the LI (PA) trap configurations, respectively.
Note that the 1D strength relates to the 3D scattering length via the relation [51],
(4) 
precisely as is done in the experimental studies of Ref. [15]; is the relative mass and is the harmonic oscillator strength in the direction perpendicular to the needle.
For the CI calculations in this paper, we assume that the totalspin projection for 4 fermions or for 3 fermions. This suffices to provide the full energy spectrum, as long as the manybody Hamiltonian does not depend on . Naturally, the manybody wave functions characterized by a given total spin are different for the different projections . For lack of space, we will not consider here manybody wave functions with for four fermions or with 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].
3 Four fermionic ultracold atoms in a doublewell trap
3.1 Four fermionic ultracold atoms: CI results
We treat here three different types of traps: a SW quasi1D trap [see Fig. 1(I)], a DWPA trap [parallel arrangement, Fig. 1(II)], and a DWLI trap [linear arrangement, Fig. 2(I)]. We start with the fouratom doublewell systems, followed by a comparison with the singlewell 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 doublewell geometries.
Figs. 1 and 2 illustrate the evolution of the spectra of 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, , and the height of the interwell barrier ). The limiting case of the singlewell (“united atom”) quasi1D trap (at ) 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 , 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 lowenergy states as the interaction strength approaches infinity (i.e., as ); all six states become degenerate at . 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 123 degeneracy in the whole range [Fig. 1(e) and Fig. 2(b)]. The tendency towards the regrouping of the energy curves according to the 123 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 , three states with , and one state with . These totalspin multiplicities, denoted here by , arise from the groupsymmetry properties of the spin eigenfunctions [12, 31] of fermions with spin ; for the multiplicities of totalspin degeneracies for any fermions, see the branching diagram [31] in A. (The theory of spin1/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 and , see A. Importantly, in the cases studied in this paper [that is, for (Figs. 1 and 2), and for in a SW and in the DWPA trap (Fig. 4), as well as in a DWLI trap (Fig. 5)] the AFM lowest spinstate is the ground state and the energy level spacings decrease with increasing interatomic repulsion.
The similar behavior of the sixfoldmultiplet bands irrespective of the different geometries of the doublewell 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 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 singleparticle densities (SPDs, green surfaces in Figs. 1 and 2) and spinresolved conditional probability distributions (SRCPDs, angleresolved pair correlations, red surfaces in Figs. 1 and 2).
The SPD is the expectation value of the onebody operator
(5) 
where denotes the manybody (multideterminantal) CI wave function.
We note that the SPD is the sum of the spinup and spindown singleparticle densities, defined as
(6) 
where and denote up or down spins.
In all cases the SPDs display four humps corresponding to the four localized fermions at the selfgenerated 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 (, ), 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(b).
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 manybody SRCPDs, , which yield the conditional probability distribution of finding another fermion with up (or down) spin at a position , given that a specific fermion with up (or down) spin is fixed at . In detail, the spinresolved twopoint anisotropic correlation function is defined as
(7) 
Using a normalization constant
(8) 
we further define a related spinresolved conditional probability distribution (SRCPD) as
(9) 
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 angledependent 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 SRCPD (for details, see A and Ref. [12]). As an example, the spin eigenfunction associated with the 4fermion , CI ground state at in the case of wellseparated DWPA parallel wells [see star in Fig. 1(e); for the corresponding SRCPD, see Fig. 1(g)] is given by [ in Eq. (LABEL:x00)]
(10) 
where the ’s (’s) denote up (down) spin fermions situated at the selfgenerated 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. (LABEL:x00) are described in A.
The equal in absolute value , coefficients in front of the four primitives in Eq. (10) agree with the probability of 0.5 [i.e., ] for the socalled “antiferromagnetic” component ( and ) found in Ref. [19] [see Fig. 1(d) therein] for the case of a twoparabola DWLI double well in the highbarrier 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 largeseparation) regime the 4fermion 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 4fermion 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 rightwell fermionic pairs in addition to the entanglement between the two fermions within each well. This acrossthebarrier entanglement is not weakening as a result of a higher barrier, and it is manifested in the mapping of the CI groundstate wave function onto the spin eigenfunction in Eq. (10).
Furthermore, the discussion above applies also to the excited states. For example, the SPD and SRCPD of the first excited state with , in the DWPA trap of Fig. 1 [having an energy 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 A), we find an angle , which is associated with a spin function of the form
(11) 
We note that the spin eigenfunctions in Eqs. (10) and (11) are orthogonal.
The two coefficients in front of the first two primitives in Eq. (11) agree with the probability of 0.166 (i.e., ) found in Ref. [19] for the “antiferromagnetic” ( and ), as well as for the “mixed” ( and ) primitives in the case of a twoparabola DWLI double well at the highbarrier regime. The third coefficient for the “ferromagnetic” primitive in Eq. (11) yields a probability of 0.666 (), again in agreement with Ref. [19].
The spin eigenfunction associated with the 4fermion , CI ground state at in the case of wellseparated wells in the DWLI linear configuration [see star in Fig. 2(b); for the corresponding SRCPD, see Fig. 2(d)] is given by the same spin eigenfunction as in Eq. (10). This is due to the fact that the left and right pairs of fermions are isolated from each other in their respective wells.
Returning to the case of four fermions in a single quasi1D trap [Figs. 1(a,b,c)], the spin eigenfunction associated with the , CI ground state at [see star in Fig. 1(a); for the corresponding SRCPD, see Fig. 1(c)] is found to have a different form from those in Eqs. (10) and (11). Specifically, the analysis of the SRCPD described in detail in A yields an angle in Eq. (LABEL:x00), which is associated with the following spin eigenfunction
(12) 
where , , and .
In the next section, we utilize the trends uncovered by the CI solutions for the spectra and wave functions of the manybody Hamiltonian, in order to develop a Heisenbergmodel phenomenology. This development aims at providing tools for analyzing quantum magnetism in double (and multiwell) ultracoldatom traps.
3.2 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 SRCPDderived spin eigenfunctions [see, e.g., the functions in Eqs. (10), (11) and (12)] are related to those of a 4site Heisenberg Hamiltonian , with the four fermions being located at the humps of the SPDs and SRCPDs, namely to [see, e.g., Eqs. (30) and (31)]
(13) 
where the symbol 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 , see B and C. The DWPA case is associated with a rectangular 2D Heisenberg ring [see schematic (III) in Fig. 1], while the DWLI and SW cases represent open linear spin chains [see schematic (III) in Fig. 2]. Due to the and reflection symmetries, has only two different exchange constants. In particular, in general, for the rectangular Heisenberg ring in the DWPA case, the interwell exchange constants and the intrawell ones . For the open 1D linear configuration of the DWLI and SW traps, , , and . The energy eigenvalues and eigenvectors of [Eq. (13)] are given in B and C. They can reproduce all the trends in the energy spectra of the sixfold energy band, as well as the totalspin multiplicities and spin eigenfunctions calculated via the CI method. In particular, in the limit of wellseparated wells (i.e., for ), one gets , , and , which coincides with the aforementioned 123 spingrouptheoretical degeneracy pattern and relative gaps within the sixfold lowestenergy CI band. Note further that the Heisenberg modeling reproduces the two different SRCPDderived spin eigenfunctions in Eqs. (10) and (11), associated with the fullyseparatedwells (, for both the DWPA and DWLI cases); compare with the eigenvectors in Eqs. (51) and (64).
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., and 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 and a triplet state with energy . The lowenergy spectrum of the double well has then three levels, , , and , 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 and .
The topology of the spin chain in Fig. 1(III) (DWPA) is indeed a closed ring, whereas the one in Fig. 2(III) (DWLI) is that of an open ring. The corresponding Heisenberg Hamiltonians are given in Eqs. (LABEL:hh1matgen) and (53), respectively; note that they have different matrix elements. The similarities between these two cases arise from the fact that the spin eigenfunctions onto which the CI wavefunctions map (as we show in both the DWPA and DWLI cases) have the same group structure, differing only in the coefficients of their components [see, e.g., Eq. (LABEL:x00) in A]; the multiplicity of the four fermions spin eigenfunctions onto which the CI spectrum maps (in both the DWPA and DWLI cases) is six (for all arrangements of 4 fermions, see A and Fig. 8).
For the singlewell case, all six CI energies have distinct values; see the spectrum in Fig. 1(a). By using the openHeisenbergchain eigenvalues , in Eqs. (5459) and fitting the ratios to the CI spectrum, we can determine the parameter that describes the single well. For example, using the fully polarized, , the groundstate, , and the 1stexcited, , energies, we obtain the ratio
(14) 
which is independent of and allows for the determination of . Fitting to the CI spectrum, we get . 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 for this ratio.
With the value of , the open linear Heisenberg chain yields (), (), (), (), (), and (), i.e., six distinct values, in agreement with the CI spectrum in Fig. 1(a). The corresponding angle in Eq. (LABEL:x00) is . This value is slightly different from the value of (corresponding to an ) 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 Isingexpansion 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 SRCPDs 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 ). A detailed investigation of this matter is beyond the scope of this paper, but it will be examined in a future publication [52].
4 Three fermionic ultracold atoms in a doublewell trap
4.1 Three fermionic ultracold atoms: CI results and the Heisenberg model for tilted double wells
CI results for ultracold Li atoms in a DWPA trap are displayed in Fig. 4 for both symmetric [zero tilt, , see schematic in Fig. 4(II) and spectra in Fig. 4(a) and Fig. 4(f)] and asymmetric wells with a moderate tilt [see Fig. 4(IV) and Fig. 4(j)] and a strong tilt [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), ] 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 threefold lowestenergy band with a characteristic 12 degeneracy pattern, converging to the same energy for . The totalspin multiplicities in this band are and , in agreement with the branching diagram for three fermions (see A). This CI energy spectrum and the correponding SRCPDs [see, e.g., the SRCPD in Fig. 4(l)] are reproduced by a 3site Heisenbergring Hamiltonian
(15) 
with and ; for the numbering of the three sites, see the schematics in Fig. 4(V) and in D. For (case of a high barrier ), the eigenenergies of are and , reproducing the abovementioned 12 CI degeneracy pattern. The CIcalulated CPDs are also in full agreement with the eigenvectors of the Hamiltonian. For example, the CI groundstate SRCPD in Fig. 4(l) [at the point ] is found to map onto the 3fermion general spin eigenfunction [see Eq. (25)] for , i.e., that is to the function
(16) 
This CIderived spin function is schematically portrayed in Fig. 4(V) and agrees with the Heisenberg eigenvector in Eq. (86).
A larger tilt generates a (3,0) CI ground state, associated with a linear UCWM [see the SPD in Fig. 4(o)]. The CI energy spectrum [Fig. 4(n)] and the correponding CPDs [see, e.g., Fig. 4(p)] are related to a 3site openlinearchain Heisenberg Hamiltonian, obtained from Eq. (15) by setting . This Hamiltonian has three different eigenenergies , , and , in agreement with the threefold CI band. The groundstate CIderived spin function is schematically portrayed in Fig. 4(VII) [ in Eq. (25)] and agrees with the Heisenberg eigenvector in Eq. (85).
Fig. 6 complements Fig. 4 in that it displays examples of all possible SRCPDs for the ground state of cold fermions in a single 1D well. It is straightforward to check in detail that the all SRCPDs agree with a mapping of the CI ground state onto the spin eigenfunction of the schematic in Fig. 6(I), i.e., with the Heisenberg vector in Eq. (85). We note that, while the spin spatial distribution is analyzed here with the use of the SRCPDs [see Eq. (7)], it is also reflected in the spatial spindensities [see Eq. (6)] shown in Figs. 6(kl); the latter agree with those displayed in Fig. 6 of Ref. [18]. Note that the sum of the up and downspin densities in Figs. 6(kl) 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 () for shown in Fig. 4. Indeed, the CI energy spectra in Fig. 4(a) and Fig. 4(f) show a sixfold lowestenergy band, comprising four states, and two 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 24 degeneracy appears, which is a doubling of the 12 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 SRCPDs shown on Fig. 4(c,d) for the lowerbarrier symmetric DW case and Fig. 4(g,h) for the higherbarrier 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 3atom UCWM that resonates between the two interlocking triangles.
We mention that the resonating behavior and the symmetrization of the manybody wave function in twocenter/threeelectron bonded systems is well known [33, 34, 35] in theoretical chemistry and in particular in the valencebond treatment of the threeelectron bond which controls the formation of molecules like He and F. Furthermore we mention that the symmetry properties of the strictly1D fewfermion problem with contact interactions have been also investigated in Refs. [53, 54]
4.2 Three fermionic ultracold atoms: The  model for symmetric double wells
To model the exactdiagonalization results shown above, one must go beyond the aforementioned simple Heisenberg Hamiltonian model [see Eq. (15) and D]. Indeed, we find that a generalization of the so called  model allows us to capture all the salient characteristics uncovered by the CI calculations. The  model [5, 44] modifies (away from the half filling) the antiferromagnetic Heisenberg Hamiltonian associated with the Mott insulator at halffilling (one electron per crystal site); it has attracted much attention, because it has been proposed for explaining the highT superconductivity arising in the case of underdoped insulators (away from the half filling when holes are present). A finite type Hamiltonian may be expressed as
(17) 
where is the coupling between the two simple Heisenberg rings defined over the sites (135) and (246); see the two blocks on the diagonal (upper left and lower right) in Eq. (18). , represented by the two offdiagonal blocks in Eq. (18), is defined by the matrix , and , where the “0” indicates an empty site; e.g., 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 is equivalent to a six by six matrix,
(18) 
where the upper index denotes explicitly the dependence on the tilt. When and , one recovers the isolatedtriangle Hamiltonian, , in Eq. (15). Below we will focus on the case of symmetric double wells, i.e., we will set , , and .
For and [case of the very high interwell barrier, kHz, in Fig. 4(f)], reproduces () the characteristic CI 12 degeneracy pattern found earlier using the simple 3site Heisenberg model [compare Fig. 4(f) and Fig. 4(j)]; see the six eigenvalues , in Eqs. (88)(93).
For lower values of the interwell barrier [ kHz, Fig. 4(a)], the 24 [(12) ] doubling degeneracy is lifted, with two lowest curves and four higher in energy (and parallel) curves (two with and two with ) 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 , with and . Then one has for the energy gap between the two lowest states, ; these energies are centered around . The remaining energies group together forming a fourfold band, centered around zero. The energy gap between the two outer (both ) members of the fourfold band is , i.e., similar to the 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 , which is much smaller than the width, , of the same band, again in agreement with the pattern in Fig. 4(a). Note that for , and a degeneracy pattern 1122 develops in disagreement with the CI spectrum. Also, when , the width of the fourfold band is twice as large as the energy gap, , 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 ultracold fermions in a symmetric () DWLI (linear arrangement) trap, as is illustrated in Fig. 5. These results can be interpreted again through the use of a corresponding  model with a similar parametrization.
5 Quantifying entanglement using a CIbased von Neumann entropy
The entanglement entropy for three 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 manybody wave functions, we adopt as a measure of entanglement the von Neumann entropy [11, 26],
(19) 
where is the singleparticle density matrix and , yielding for an uncorrelated singledeterminant state.
The singleparticle density matrix is given by
(20) 
and it is normalized to unity, i.e., . The Greek indices (or ) count the spin orbitals (of dimension )
(21) 
and
(22) 
where denote up (down) spins.
Since the allowed maximum value for in our CI calculations is (we use a typical basis of singleparticle space orbitals), it is notable that the calculated values in Fig. 4 remain smaller than , and in particular in the regime of strong correlations, i.e., for . 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 entropies for three Li atoms in a DWLI trap, see Fig. 6(a).
6 Summary and Outlook
In this paper, we have presented timely advances in the growing field of fewbody 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 nonitinerant and fermionizationlike aspects of the physics of ultracold few fermions trapped in doublewell confinements, with various 1D and 2D trapping geometries, as a conduit for emulating quantum magnetism and related phenomena beyond the strictly 1D singlewell (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 spinHamiltonians (Heisenberg and  models). Our calculations for and ultracold fermionic Li atoms in SW and double well (DW) traps with linear (DWLI) or parallel (DWPA) geometries, reveal formation of antiferromagnetic ordering for the lowestenergy bands over the entire range of interparticle contact repulsion studied here.
For ultracold atoms in a symmetric DWPA trap with very strong interatomic repulsion, we find (via miscroscopic, CI, calculations) formation of a twodimensional ultracold Wigner molecule (UCWM) of nonitinerant character. For the symmetric parallel DW trap the formation of the 2D UCWM leads to mapping of the interacting 4atom 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 fouratom linear (1D) UCWM in juxtaposition with mapping onto a linear Heisenberg spinchain. These mappings enable employment of the corresponding Heisenberg model Hamiltonian, whose solutions reproduce well the results of the microscopic, numericallyexact, calculations.
For 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 coldatom quantum computing qubit.
In contrast to the asymmetric DW case, description of the ultracoldatom 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 SRCPDs, are not attainable with the simple Heisenberg model, requiring instead the more intricate type model [5, 6, 44], consisting of two coupled resonating triangular 2D UCWM Heisenberg clusters. The emergence of the  model for the description of quantum magnetism (in particular AFM ordering) in a trapped fewbody ultracold atom system, strongly suggests its future role as a useful laboratory for exploration of the elementary building blocks of highT superconducting behavior [44, 55, 56].
Appendix A Spin eigenfunctions for 4 and 3 fermions. Comparison with CI CPDs
We outline in this Appendix several properties of the manybody spin eigenfunctions which are useful for analyzing the trends and behavior of the spin multiplicities exhibited by the CI wave functions for and ultracold fermions. The spin multiplicities of the CI wave functions lead naturally to analogies with finite Heisenberg clusters [8, 10] and to type models.
A basic property of spin eigenfunctions is that they exhibit degeneracies for , i.e., there may be more than one linearly independent (and orthogonal) spin functions that are simultaneous eigenstates of both and . These degeneracies are usually visualized by means of the branching diagram [31] displayed in Fig. 8. The axes in this plot describe the number of fermions (horizontal axis) and the quantum number of the total spin (vertical axis). At each point , a circle is drawn containing the number which gives the degeneracy of spin states. It is found [31] that
(23) 
Specifically for particles, there is one spin eigenfunction with , three with , and two with . In general the spin part of the CI wave functions involves a linear superposition over all the degenerate spin eigenfunctions for a given .
For a small number of particles, one can find compact expressions that encompass all possible superpositions. For example, for and , one has: [12, 57]
where the parameter satisfies and is chosen such that corresponds to the spin function with intermediate twofermion spin and threefermion spin ; whereas corresponds to the one with intermediate spins and .
For and , one has [57]:
(25)  
For the general expressions for the remaining spin combinations, and , for and fermions, see Refs. [12, 57, 58].
For each SPD corresponding to a given CI state of the system, one can plot four different spinresolved CPDs, i.e., , , , and . 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 states for fermions and the states for , however, we found that knowledge of a single CPD, is sufficient in the regime of Wignermolecule formation. Indeed, the specific angle specifying the spin function [Eq. (LABEL:x00)] for or the spin function [Eq. (25)] for fermions can be determined through a procedure exemplified in the following through two examples for fermions:
Example 1; case of the CPD in Fig. 1(g): The same labeling that numbers the sites determines also the lefttoright ordering of the localized electrons in each of the six primitive spin functions , , etc., that span the eigenfunction in Eq. (LABEL:x00). Namely, the fermion localized at the hump No. 1 corresponds to the far left position in the primitive, the fermion localized at the hump No. 2 corresponds to the second from the left position in the primitive, the fermion localized at the hump No. 3 corresponds to the third from the left position in the primitive, and the fermion localized at the hump No. 4 corresponds to the far right position in the primitive. The numbering of the humps does not necessarily follow the cardinal ordering 1,2,3,4, as will become evident below from the second example concerning a linear Heisenberg chain. An inspection of Eq. (LABEL:x00) shows that only the first three primitive spin functions in can be associated with [compare the CPD in Fig. 1(g)], namely , , and ; these are the only primitives in Eq. (LABEL:x00) with a down spin in the site labeled as 1 [see diagram in Fig. 1(III)]. From these three primitives, only the first and the second contribute to the partial conditional probability of finding another fermion with spindown in site No. 4, while the first fermion is fixed at site No. 1. Taking the squares of the coefficients of and in Eq. (LABEL:x00), one gets
(26) 
Similarly, one finds
(27)  
and
(28) 