Ordering in one-dimensional few-fermion clusters with repulsive interactions

By using a diffusion Monte Carlo (DMC) technique, we studied the behaviour of a mixture of spin-up and spin-down cold fermionic atoms confined in one-dimensional harmonic potentials. We considered only small balanced clusters (up to eight atoms) or arrangements in which the difference between the population of fermions with different spins is one (with a maximum total number of nine). The atom–atom interactions were modeled by a contact repulsive potential. We focused on the ordering of the confined fermions, finding that the probability of having a pure antiferromagnetic arrangement of atoms decreases strongly with the number of fermions, as to be almost negligible for our largest clusters. Exactly the same can be said of a pure ferromagnetic state.


Introduction
Arguably, one of the most interesting recent experimental developments in the field of quantum physics is the possibility of creating effectively one-dimensional few-body clusters of cold fermions with a perfectly controlled number of atoms [1][2][3][4]. Not only their number, but the nature of their interactions can be fixed by manipulating the strength of an applied magnetic field in the vicinity of a Feshbach resonance [5]. In the case of a quasi-one dimensional traps, we can also change those interactions by changing the transverse confining frequency to generate a connement-induced resonance (CIR) [6].
That richness of potential behaviour has triggered the appearance of many theoretical studies devoted to those systems [7][8][9][10][11][12][13][14][15][16][17], all of them having in common the use of numerical approximations to the solution of the following strictly one-dimensional (1D) Hamiltonian: Here, N is the total number of particles of mass m, i.e. the sum of the number of atoms with spin up (N ↑ ) and of those with spin down (N ↓ ). The coupling constant, g 1D , is given by g ma 2 , with a 1D the 1D scattering length of the pair of atoms with unequal spins. The value of g 1D can take any value ranging from −¥ to 0 to ¥, even tough this work we only concerned ourselves with repulsive interactions (g D 1 > 0). The system is confined via a harmonic potential of frequency ω. The high aspect ratio of the traps used in the experiments makes equation (1) a reasonable approximation to the description of the few-atom clusters.
All previous theoretical works have been devoted mainly to obtain the energy of sets of spins confined in the 1D traps. This is reasonable since this is an experimental observable, at least in the case of attractive interactions [3]. However, another common output is the density profiles of spin-up and spin-down fermions [7,9,13]. For balanced (N N =   ) clusters, both spin densities are exactly the same [9]. On the other hand, for imbalanced (N N ¹   ) ones, the density of the minority spin component at the center of the cluster is increased with respect to the spin density of the majority spin at the same spot. In the impurity case of (N >  1, N ↓ =1), this has been interpreted as a incipient phase separation that would lead to the onset of a 1D form of Stoner ferromagnetism [13]. In this work, we will further explore that possibility by studying the ordering of atoms of different spins when confined in the 1D harmonic trap. Instead of considering the impurity case, we will study both balanced Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. and slightly imbalanced (N ↑ =N ↓ + 1) clusters with N ranging from 3 to 9. We calculated also other magnitudes in order to compare our results with those of the previous literature to check the validity of the numerical method used.

Method
Our goal was to study the ground state of a series of 1D clusters described by the Hamiltonian in (1). To solve the associated Schrödinger equation we used the diffusion Monte Carlo (DMC) method [18], a stochastic technique that allows us to obtain that ground state within certain statistical uncertainties. The only downside is that DMC needs an initial approximation to the wavefunction of the system to work, the so-called trial function. In this work, we used: where z i is the position of the atom i and D  and D  are Slater determinants that depend on the coordinates of the spin-up and spin-down fermions, respectively. Those determinants would include the eigenfunctions of the one-body harmonic Schrödinger equation with the N ↑ (or N ↓ ) lowest energies. This means that we studied only the ground state configurations of the spin arrangements. However, instead of calculating those determinants, we used the standard approximation and evaluated instead those of the equivalent Van der Monde matrices [8].
This speeds up the calculation considerably. The nodes of the total wavefunction are those of the Slater determinants, that might or might not be the ones of the exact solution. Since we employed a fixed-node diffusion Monte Carlo (FN-DMC) scheme, the position of those nodes did not change throughout our simulations. FN-DMC only provides exact estimators of the energies if the nodal surfaces coincide with those of the eigenfunctions of the many-body Hamiltonian considered. Unfortunately, we do not know a priori if that is the case, but we can be sure that those energies are at least upper bounds to their ground state values [21]. To completely define z z , , N 1 F  ( ) , we chose z ij y ( )as the exact solution of the two-particle Hamiltonian in the unconfined (ω=0) case [19]. For g D 1 > 0, this means [20]: with R m a variational parameter and z ij the distance between particles i and j. k was derived for each value of g 1D from the transcendental equation: For convenience and to follow the standard practice in all previous literature, all the lengths will be measured in units of the oscillator length, σ= m  w , and the energies will be scaled by w (ÿ, Planck's constant). Taking that into account, we can define a new coupling constant, g g D 1 ws = , parameter that will be used in the remaining of the paper. Also, to simplify the notation, a set of N ↑ =n and N ↓ =m atoms will be referred to as a n+m cluster.

Results
Our first concern will be to compare the energies given by DMC to those obtained by other methods. Our benchmark was the energy per particle of the 2+1 system, displayed in figure 1. To obtain those energies, in 2 we considered a D  that included the first two harmonic oscillator wavefunctions, and D  =1 was set as a single Gaussian, as corresponds to the ground state of a cluster that includes two spin-up and one spin-down fermions. We found our results comparable [9,11,13,15] or even a little bit lower [7] than those in the literature for the same g range. This supports the good quality of the nodes of (2) and validates the use of FN-DMC to treat this problem. The line in figure 1 corresponds to the McGuire analytical expression of [10], that is an upper bound to the DMC results.
Also, as one would expect, our energies approach asymptotically the value for a Tonks-Girardeau (TG) gas (g = ¥) for large values of g, even though there are noticeable differences even for gʼs as large as 50. The TG value was obtained considering that the energy of a 2+1 cluster in the TG regime is degenerate to that of a 3+0 arrangement for g=0 [8]. On the other hand, in figure 2, we show the energies for larger clusters, energies that are in good agreement with those of [15] for 2+2, 3+3 and 4+4 clusters in the range of gʼs covered there. In that figure, we also displayed the energies for g=0 and g = ¥ for comparison. We found that the bigger the cluster, the larger the difference between the g=50 and the TG limit. This is due to the fact that, for strictly 1D systems, the ratio between the kinetic and potential energies increases, thus there is more difference with respect to the Tonks regime.
Another way to characterise the clusters is to display their spin-up and spin-down density profiles. In figure (3), we show them both for the 4+4 and 5+4 cases, using a mixed estimator for the density. Given the quality of the ground state energies discussed above, we believe that the use of pure estimators would not change significatively the results displayed. This is supported by the fact that the results we obtained from variational Monte Carlo calculations using the trial functions defined above and not shown for simplicity are very close to those derived from FN-DMC. In addition, the profiles for our smallest clusters are similar to those found previously in literature [9]: the distributions of both components are alike for balanced clusters, and differ from each other for imbalanced ones. We can see also that, at least for small and intermediate values of g, the 'excess' of the majority component in imbalanced clusters is accumulated at the wings, that being similar to what happen to very large clusters with attractive interactions [22,23]. The results for g=0 are also very similar to those derived in [24,25].
In that regime of small and intermediate values of g, that increasing in the repulsive interactions produces also a broadening of all profiles that progressively erases the distinct peaks (one for each particle) that appeared in the non-interacting case. This flattening is more pronounced that for smaller clusters [9]. In that g range, we observe also a certain alternation in the position of the peaks of the density profiles between the majority and minority components. That could hint at the creation of an antiferromagnetic state, in which the two neighbours of a given spin-down atom will be two spin-up fermions. Unfortunately, even though this pure antiferromagnetic phase does exists, the probability of finding it is relatively small (see below). In the regime of larger gʼs, the smoothing-out of the profiles is reversed since the behaviour of the system would trend to that of a n m + ( ) + 0 cluster for g=0. This reappearance of the peaks can be seen for balanced clusters in figure (3)  when we go from g=10 to g=50. We have then a structure similar to that of that of [24] for non-interacting fully polarised clusters.
Apart from the previous results, the DMC algorithm allows the easy study of the spatial distribution of the atoms for relatively large clusters. The studying of this ordering is relatively difficult with any other technique and, to our knowledge, have been done only for clusters up to 4 atoms [14] when g  ¥. For instance, for a 2+1 cluster, the spin ordering inside the cluster can be , , ). We considered the first two to be equivalent to each other for symmetry reasons. So in this case, we have only three possible spin configurations, of which one of them could be considered a pure 'antiferromagnetic' (AF) state, and the other two can be labelled 'ferromagnetic' (F), since all the spin-up fermions are close together. For larger clusters, we will have also other spin distributions, nor purely AF nor F. Following [14], those configurations will be labelled 'mixed' (M). In general, for arrangements of N spins, the total number of possible configurations is C N =N!/ (N N   ! !), of which only two can be considered ferromagnetic, and one (or two, equivalent to each other, for balanced arrangements) antiferromagnetic. So if the spin orderings were purely random, the probability of having a spin-separated F configuration will be 2/C N , and of having an AF one either 1/C N (imbalanced clusters) or 2/C N (balanced case). However, one would expect those probabilities to change due to the existence of repulsive interactions between atoms with different spins and to their fermionic nature.
To illustrate that point, we show in figure 4 the probabilities of the different configurations for a 4+4 cluster. The main features of this figure are common to all the other sets of fermions considered in this work. First of all, we can see that for all values of g, the spin-separated ferromagnetic state is almost non existent, since its probability is always very close to zero. At the same time, the probability of finding a pure AF ordering is between 0.2 for g=0 and ∼ 0.15 for g=100. The remaining configurations are mixed, this state being by far the most common. This is very easily understood if we take into account the fact that the label 'mixed' lumps up many different spin orderings. We can see also that the changes in the different probabilities with the strength of the interaction are relatively minor, a 25% at most. On top of that, all the probability shifts are produced between g=0 and g~20, being unchanged for larger gʼs.
Our simulation results for other clusters together with the values derived for purely random distributions are given in figure 5. A couple of things are immediately apparent. First, that the general trend is captured by the 'entropic' (purely combinatorial) term: the probability of having AF and F orderings decreases strongly with the number of particles in the cluster, something more obvious for the spin-separated F phase than for the AF one. Second, that AF distributions are more common than a random ordering would suggests, and more so for smaller values of g, something presumably due to the particle-particle interactions. We can see also that an AF phase is more probable than a mixed one only for a 2+2 cluster, in line with the results obtained in [14] for a TG gas. Our simulations indicate as well that the AF distribution is more favoured in balanced clusters than in imbalanced ones, as could be expected by combinatorial reasons. However, this effect is smoothed out for larger clusters, and probably will be absent in very big ones.
Finally, with the help of figure 6, we will try to explain the short range spin ordering inside the clusters. There, we display a radial distribution function (that we called g(z)), defined as the probability of finding another spin of a given type at a fixed distance from the one serving as reference. To simplify matters, the value of the    integrals of all distribution functions were set to one. What we see there is that for g=0 distributions (upper panel in figure 6), the probability of having a different spin as a first neighbour is larger than that of having a same-spin closest pair, even tough this last one is not negligible (except at z ij =0, where the probability is exactly zero as a consequence of Pauli's exclusion principle). This favours but not implies the existence of AF clusters in large proportions since we do not see an alternating structure of maxima and minima at distances further away from zero. This is true both of balanced and imbalanced clusters. A similar picture can be deduced from the lower panel of figure 6 for a larger value of g: the probability of a ,   ( ) first neighbour pair is still larger than that of the other possibilities, and this translates into a probability of having an AF cluster higher than for a random spin distribution. We observe also that for g=50, the pair distribution functions go to zero at z=0 in all cases, something that is expected in the TG limit. This is due to the fact that for g  ¥, the interparticle part of the wavefunction for different spins should trend to zero, being equal to zero only in the TG limit [8,16].

Discussion
In this work, we have used a DMC calculation to study both the energetic and the spin ordering of small 1D clusters with repulsive interactions. As indicated above, the energies obtained were similar to those in the previous literature, what validates this numerical tool to study these systems, and testifies to the quality of the nodes of the trial function. As a bonus, DMC can readily give us the probability of the different spin configurations compatible with each ground state, information that is more difficult to get from other techniques. From those probabilities, we can infer that the possibility of having pure AF or F phases at T=0 decreases very fast with the size of the cluster. This means that the hopes of having a ground-state 1D form of Stoner ferromagnetism are negligible, except for the rather peculiar 2+1 case, in which we can have only AF and F phases (see above). We concluded also that the spin ordering is mostly determined by the configurational entropy of each cluster, related to the total number of configurations possible for a set of N ↑ and N ↓ fermions, and that the interaction energy plays a relatively minor role.