Pairing in a system of a few attractive fermions in a harmonic trap

We study a strongly attractive system of a few spin-1/2 fermions confined in a one-dimensional harmonic trap, interacting via two-body contact potential. Performing exact diagonalization of the Hamiltonian we analyze the ground state and the thermal state of the system in terms of one-- and two--particle reduced density matrices. We show how for strong attraction the correlated pairs emerge in the system. We find that the fraction of correlated pairs depends on temperature and we show that this dependence has universal properties analogous to the gap function known from the theory of superconductivity. In contrast to the standard approach based on the variational ansatz and/or perturbation theory, our predictions are exact and are valid also in a strong attraction limit. Our findings contribute to the understanding of strongly correlated few-body systems and can be verified in current experiments on ultra-cold atoms.

Introduction. Quantum engineering, i.e. coherent controlling of atomic systems, is a rapidly developing field of modern physics. Its development is important because it gives us the possibility of a deep study of quantum mechanics. Present experimental methods in quantum engineering allow one to prepare and control the system in a state with a rigorously chosen number of interacting atoms. Spectacular experiments on ultra-cold fewbody fermionic systems confined in a one-dimensional harmonic trap [1][2][3][4] have opened alternative way to study strongly correlated quantum systems. It is believed that these kind of studies can shed a new light on our understanding of the crossover from the few-body world to the situation with a macroscopic number of quantum particles [4]. They can help to answer the fundamental question on emerging of collective behavior from the two-body interactions.
Recently, harmonically trapped effectively one dimensional few-body systems were intensively studied on theoretical [5][6][7][8][9] and experimental [1][2][3][4] level. For example, it is known that in a strictly one-dimensional geometry, in the limit of infinite repulsion (known as Tonks-Girardeau limit) the bosonic system can be mapped on the system of noninteracting fermions [10][11][12][13]. Moreover, in this limit, the ground-state of the one-dimensional system of a few spin-1/2 fermions is highly degenerated [6,9] and a combination which is fully antisymmetric with respect to the position of fermions is identical to the ground state of the system of spinless fermions [5,14]. This degeneracy of the ground-state and its consequences for the experimental results were also studied [9]. The mapping mentioned above, in the limit of infinite repulsion was probed experimentally for two-fermion system [2], for which the exact solution is known [15][16][17][18].
In the recent paper [3] the attractive system was probed experimentally and formation of fermionic pairs was observed. This experiment opened an alternative path to the experimental exploration of the few-body systems. Therefore the quantitative predictions based on the first-principle theory of a few-body system are necessary. We should mentioned here that pairing in a system of fermionic atoms in the mean-filed Hertree-Fock approach was studied recently in [19,20].
The Model. In this Letter we study a few fermions of mass m, confined in a one-dimensional harmonic trap of frequency Ω, mutually interacting via short-range deltalike attractive potential. In the second-quantization formalism the Hamiltonian of the model has a form: where Ψ σ (x) annihilates a fermion with spin σ at a point x, and g measures a strength of the contact interactions. For convenience we expressed all quantities in harmonic oscillator units, i.e. energies inhΩ, lengths in h/mΩ, temperatures inhΩ/k B etc. The dimensionless interaction constant g = (m/h 3 Ω) 1/2 g 1D is proportional to the effective one-dimensional interaction g 1D between fermions of opposite spins. The Hamiltonian (1) commutes separately with the number of fermions in each spin componentN σ = dx Ψ † σ (x)Ψ σ (x). Here, we consider a situation where both species are distinguishable. Therefore, the field operators of different spins do commute, Ψ ↑ (x),Ψ † ↓ (x ) = Ψ ↑ (x),Ψ ↓ (x ) = 0. There are no any symmetry conditions related to the exchange of particles with different spins. However, particles of the same spin satisfy the standard anticommutation relations Ψ σ (x),Ψ † σ (x ) = δ(x − x ) and Ψ σ (x),Ψ σ (x ) = 0. We focus our attention on correlations between fermions with opposite spins in a balanced system N ↑ = N ↓ = N/2.
The Method. We fix the number of fermions in each component, then we decompose the field operator Ψ σ (x) in the single-particle basis φ n (x) of harmonic oscillator eigenstates with sufficiently large cutoff, n ≤ n max , and we diagonalize the full many-body Hamiltonian (1) in this basis via Implicitly Restarted Arnoldi Method [21]. In this way we find N eigenstates |G i and the corresponding lowest eigenenergies E i . In our case N = 60.
Properties of the ground state. To get some intuition, first we discuss properties of the system at zero temperature, i.e. we analyze the ground-state |G 0 . In the limit of vanishing interaction g, the ground state of the system is simply formed by filling up the lowest single-particle states of the harmonic potentials by pairs of opposite spin fermions. This fact is directly manifested in the spectral decomposition of the reduced two-fermion density matrix (normalized to Tr ρ (2) = 1). Let us notice that coordinates x 1 and x 1 (x 2 and x 2 ) correspond to the spin-down (spin-up) fermions. From now on we will use this convention throughout our Letter. Distinction between the two species is possible only because in our Hamiltonian there are no terms responsible for spin flipping mixing different species.
All nonzero eigenvalues of ρ (2) are equal to λ i = 1/N ↑ N ↓ . They determine the relative populations of the corresponding eigenvectors, i.e. the two-particle orbitals ψ i (x 1 , x 2 ) describing possible pairings of fermions with opposite spins. Situation changes dramatically when attractive interaction are switched on. Then, we observe that one eigenvalue starts to dominate in the system (Fig.  1). According to the Penrose-Onsager criterion, an occurrence of the dominant orbital in the spectral decomposition of the reduced density matrix is an indication of the condensation of fermionic pairs in the selected orbital. For convenience we denote the dominant orbital and its eigenvalue by ψ 0 and λ 0 respectively. To get better understanding of the structure of the ground state of the system it is convenient to relate the orbitals ψ i (x 1 , x 2 ) to the orbitals ϕ i (x) -the eigenvectors of the reduced single-fermion density matrices (normalized to one) where σ ∈ {↑, ↓}. Since we assumed exact symmetry between 'up' and 'down' fermions, ρ (1) is exactly the same for both spin components. By projecting two-fermion orbitals of ρ (2) on the product of the single-fermion orbitals of ρ (1) we find that the set of all orbitals ψ i (x 1 , x 2 ) can be divided into two classes. For a given two-fermion orbital ψ i belonging to the first class one can find exactly two single-fermion orbitals ϕ j and ϕ k (j = k) such that Note that the two-fermion orbital has no any symmetry with respect of interchanging spin up and spin down fermions, i.e. x 1 with x 2 . In the second class each twofermion orbital ψ i (x 1 , x 2 ) has a structure of the correlated pair -the superposition of states with two particles of opposite spin occupying the same singe-particle state: where j |α i (j)| 2 = 1. From our numerical analysis it follows that the dominant orbital, ψ 0 (x 1 , x 2 ), always belongs to the second class. Moreover it is symmetric with respect to exchange of spin-up and spin-down fermions. Moreover, the structure of the dominant orbital is unique in the sense that only in this orbital we find non vanishing expectation value of the operator where the operatorb σi annihilates fermion in the orbital ϕ i of the reduced single-fermion density matrix ρ (1) . Expectation value S 0 = G 0 |Â|G 0 is a direct counterpart of the superconducting order parameter S = dr Ψ ↑ (r)Ψ ↓ (r) known from the BCS theory [22]. Op-eratorÂ can be understood as describing the process of destruction of a correlated fermionic pair occupying the same one-particle state and putting the two fermions anywhere. This operator can be viewed as the number conserving destruction operator of the correlated pair. The strong correlation in the dominant orbital can be observed in the probability density of finding fermions with opposite spins (Fig. 2). In the position representation we find strong correlation between x 1 and x 2 . In contrast, in the momentum space we find strong anti-correlation between p 1 and p 2 . Moreover, for strong enough attraction we find that probability of finding a pair with opposite momenta does not depend on the value of the momentum (Fig. 3). The dominant orbital supports such correlated pairs that both particles are very likely to be found in the same position with opposite momenta. This observation is one of the fundamental features of correlated pairs known from the theory of superconductivity. Therefore, the dominant orbital |ψ 0 = j α 0 (j)b † ↓jb † ↑j |Ω is an analog of a Cooper pair [23].
Non zero temperatures. Now, we discuss a situation when the system studied is in contact with a thermostat of given temperature T . Therefore, instead of the manybody ground state |G 0 , we assume that the state of the system is well described by the density matrix operator ρ T = Z −1 exp(−Ĥ/T ), where Z is a canonical partition function assuring a proper normalization Trρ T = 1. Numerically, the density matrix operatorρ T is constructed from a finite number of N = 60 eigenstates |G i of the full many-body Hamiltonian as followsρ In analogy to the zero temperature case, we define a reduced twofermion density matrix and we diagonalize it. We perform the whole spectral analysis of ρ T in analogy to the zero temperature limit. First, we find temperature dependent single-particle orbitals annihilated by temperature dependent operators b σi . These states serve as a basis for two-fermion orbitals of ρ T . We find that structural properties of the orbitals are not affected by finite temperatures, i.e. the division of all orbitals into two classes with respect to their representation by finite-temperature single-fermion orbitals holds. In particular coefficients α 0 (j) practically do not depend on the temperature. We also checked that the dominant orbital of ρ (2) T remains the only orbital with non-vanishing expectation value of the particle conserving destruction operatorÂ (defined with temperature dependentb σi ). It means that the superconducting order parameter S(T ) = Tr ρ (2) TÂ is well defined and can be directly related to the superconducting fraction represented by the dominant eigenvalue at given temperature λ 0 (T ): With increasing temperature we observe a rapid decay of the occupation of the dominant orbital, and in consequence, the superconducting fraction tends to zero. It is worth to notice that the dominant eigenvalue λ 0 (T ) of the reduced two-fermion density matrix for the system studied is closely related to the gap function ∆ known from the theory of superconductivity for isotropic system of interacting fermions in the thermodynamic limit. along the lines x1 = x2 and p1 = −p2 for position and momentum domain respectively. We observe that for strong attraction the density distribution in momentum space becomes flat. This fact shows that probability of finding correlated fermions with opposite momentum does not depend on the value of the momentum. This is one of the features of superconducting fraction well known in the many-body theory of superconductivity.
Therefore it is quite obvious that the behavior of λ 0 (T ) is in one-to-one correspondence with the behavior of the gap function ∆ known from the BCS theory. In Fig. 4a we plot a temperature dependence of the superconducting fraction λ 0 (T ) for N ↓ = N ↑ = 4 and different interactions g. As it is seen, this quantity crucially depends on the interaction constant g. Nevertheless, it has some fundamental universality which indicates correspondence to the BCS theory. To show this universality we define two natural quantities which depend on g, i.e. superconducting fraction at zero temperature, λ 0 (0), and the characteristic temperature T 0 at which, by definition, the superconducting fraction drops by some chosen fac- tor η = λ 0 (T 0 )/λ 0 (0). We checked that our predictions are insensitive to any reasonable choice of the value of η and in the following we set η = 85%. Then, by plotting the normalized data, i.e. λ 0 (T )/λ 0 (0) vs. T /T 0 , we find that all numerical points collapse to the universal curve which does not depend on model details (see Fig.4b). We checked that this universal behavior of λ 0 (T ) is independent of the definition of the characteristic temperature T 0 . Moreover, we find that the characteristic temperature T 0 is a linear function of interaction g (Fig. 4c). On the other hand, the dependence of the zero temperature superconducting fraction λ 0 (0) on the interaction constant g is exactly the same as the dependence of the gap-function in the BCS theory (Fig. 4d) λ 0 (0) = Λe −1/M|g| .
Amazingly, the value of the parameter M obtained from the numerical fit to the data is equal to 2.001 ± 0.013, which is in perfect agreement with the intuitive picture of the theory of superconductivity. In the case of homogenous and infinite system, the theory predicts that for few first orbitals ψi at T = 0 and g = −10 for N = 8 (the gray scale corresponds to the orbital index i as indicated in the legend). The hight of the peak is equal to the probability of finding simultaneously two particles (described by ψi) above chosen harmonic oscillator level n. In the inset we show the probability that the pair detected above a level n comes from the dominant strongly correlated orbital ψ0. Note, that for large n the pair, if detected, is the Cooper-like pair with probability approaching one.
M should be equal to the density of states at the surface of the Fermi sea. In the one-dimensional case studied here, the size of the Fermi surface does not depend on the number of particles and is simply equal to 2obviously, the Fermi level can be occupied only by two fermions with opposite spins. Finally, let us notice that dependence of the characteristic temperature T 0 on g is linear and different than in the BCS theory for homogeneous and 3D system. This comes from the fact that the number of particles here is small and the system is one-dimmensional. Final remarks. We believe that our predictions can be verified in current experiments with ultra-cold fermions in optical traps. By tilting the trap potential (like in experiments of Jochim group [1]) the particles of energies above selected harmonic oscillator level are allowed to escape and can be detected. The probability of simultaneous detection of two particles above selected energy n is related to the two-particle cumulative distribution function C(n) (CDF). In Fig. 5 the CDF for a few first orbitals ψ i of the two-fermion density matrix ρ (2) (i = 0, . . . , 3) at T = 0 and g = −10 for N = 8 is shown. The gray scale corresponds to the orbital index i as indicated in the legend. The hight of the peak is equal to the probability of finding simultaneously two particles (described by ψ i ) above chosen harmonic oscillator level n. In the inset we show the probability that the pair detected above a level n comes from the dominant strongly correlated orbital ψ 0 . Note, that for large n the pair, if detected, is the Cooper-like pair with probability approaching one. The figure clearly shows that if the pair is detected above the harmonic oscillator level n = 9 (both particles must have energies larger than n) then with probability about 85% it has to be the correlated Cooper-like pair occupying initially the dominant orbital. This way the correlated pair can be distilled from the trap. Measuring the two particle correlation function of the pair in the momentum domain (Fig. 2) is the ultimate proof of detection of the Cooper pair. We checked that this is true also at finite temperatures.