Virial relations for ultracold trapped Fermi gases with finite range interactions through the Bardeen–Cooper–Schrieffer Bose–Einstein-condensate crossover

We study the virial relations for ultracold trapped two-component Fermi gases in the case of short finite range interactions. Numerical verifications for such relations are reported through the Bardeen–Cooper–Schrieffer (BCS) Bose–Einstein-condensate (BEC) crossover. As an intermediate step, it is necessary to evaluate the partial derivatives of the many-body energy with respect to the inverse of the scattering length and with respect to the interaction range. Once the binding energy of the formed molecules in the BEC side is subtracted, the corresponding energy derivatives are found to have extreme values at the unitary limit. The value of the derivative with respect to the potential range in that limit is large enough to yield measurable differences between the total energy and twice the trapping energy unless the interacting system is described by extremely short potential ranges. The virial results are used to check the quality of the variational wavefunction involved in the calculations.


Introduction
In the absence of interaction, the virial theorem relates the energy per particle of a confined atomic gas with the trapping potential. If that potential is harmonic, the theorem states that the total energy per particle is twice the mean trapping potential energy For strongly interacting two-component Fermi gases, confined by a harmonic trap in the unitary limit, this relation was also shown to be valid experimentally and theoretically [1,2]. The first derivation of that theorem considered zero-range interactions and made use of the local density approximation. Further insight into the fundamental basis of this relation revealed several remarkable features of the unitary gas such as its scaling properties [3] or a mapping, using group theory, between the trapped and the free space problem [4]. Recently, the Hellmann-Feynman theorem was used to prove equation (1) at the unitary limit [5,6] and to generalize the virial relations for finite scattering lengths [7,8]. The physical origin of these relations lies in the variational stability imposed by the Schrödinger equation to the energy and states under small independent changes of the parameters that define the Hamiltonian that models the system. Although nonzero range interactions would arise in any first principles description of ultracold interacting particles when the range of the interparticle interaction is smaller than all the length scales in the system, the details of the interaction are believed to be unnecessary and the parameters defining contact pseudo-potentials, like the scattering length, are expected to be sufficient to characterize the system. The importance of this asseveration deserves a detailed study of the dependence of the predicted properties of an ultracold dilute gas on the parameters that define short finite range interactions. For general confinement potentials and finite range interactions, a general virial theorem has been stated as follows in [8].
Consider a Hamiltonian for a system of N particles with arbitrary statistics: where H and its domain depend on p parameters with length dimensions 1 , . . . , p , on h and on the mass of the particles. U ( r 1 , . . . , r N ) denotes a regular arbitrary function that allows the domains of H and H to coincide, with r i the position vector for the ith particle. Then, with E being the total energy. For N particles confined by a harmonic trap, where E tr = U is the trapping potential energy.

Virial relations for a short-range potential
In the present paper, we study 2N fermionic atoms of mass m in two equally populated hyperfine states (N = N ↑ = N ↓ = 165) confined by an isotropic three-dimensional harmonic trap of frequency ω, and interacting through an attractive finite range potential V = −|V 0 |e −r/r v . This potential is characterized by two parameters, its strength V 0 and its range r v . When the kinetic energy of the atoms is low enough, the scattering length is a proper parameter to describe the interacting system. For a given number of s-wave bound states and a given r v , there is a one-to-one relationship between the strength of the potential V 0 and the scattering length a. We consider the case where at most one bound state is admitted by the potential and find the ground state of the many-body Schrödinger equation approximately, via a variational Monte Carlo calculation, for several scattering lengths a and short potential ranges r v r ho ≡ √h /mω.
We then study the behavior of the total, internal and trapping energies as a function of both the length parameters a and r v to verify (4). The explicit expression of the Hamiltonian is and the corresponding virial relation becomes On the BEC side of the crossover, the total energy E can become extremely large compared to the total energy E on the BCS side due to the contribution of the binding energy of the formed molecules. This fact increases the numerical errors in the evaluation of the derivatives in (7). In order to isolate this two-body effect from many-body effects, we have found it convenient to take into account the behavior of the free space binding energy as follows. The two-body problem, is analytically solvable for s-states, so that the scattering length is explicitly given by with ζ = (2r v √ |V 0 |m/h), C = 0.577 215 664 901 . . . is the Euler constant and J ν and N ν represent the Bessel functions of the first and second kinds of order ν, respectively. This problem has the following bound states: where N is a normalization factor and y = ζ e −r/2r v . The boundary condition at the origin implies J x s (ζ ) = 0, so the corresponding eigenenergies ε (r v ) s fulfil the equation That is, x s is determined by ζ and We shall work with z 0 < ζ < z 1 with z 0 and z 1 being the first two zeros of the Bessel function J 0 . Under these conditions, just one bound state is admitted for each positive scattering length a. Given a and r v and using (9), we can write As a consequence, the ground-state binding energy ε (r v ) 0 of the two interacting particle system in otherwise free space satisfies the equation Thus, if we definẽ and the virial relation, equation (7), reads This expression is easier to verify numerically than (7). Note that, from a dimensional analysis, equations similar to (13) can be expected to be valid for other forms of the potential.

Approximate wavefunctions for the ground state of an ultracold Fermi gas
Approximate ground-state eigenfunctions for the Hamiltonian (5) were obtained variationally. The trial wavefunctions used have the Eagles-Leggett form through the BCS-BEC crossover regime. In this equation, A denotes the antisymmetrizing operator to be applied to all fermions of each species and, as a contrast to the unrestricted Eagles-Leggett formalism, we take the form of the paired functions as where ϕ(r i j ) is the s-wave ground-state solution of the trapped interacting two-body problem The variational parameter λ EL modulates the optimal shape of the atomic cloud. The evaluation of the mean value of the many-body Hamiltonian, equation (5), for the restricted Eagles-Leggett trial wavefunction was performed using Monte Carlo techniques that take advantage of the structure of the function [10].
For weak interactions, that is, for negative scattering lengths shorter than the mean separation between interacting atoms, lower variational energies are obtained using the lengthscaled ground-state solution of the non-interacting problem (which is a product of Slater determinants) multiplied by a Jastrow correlation function The scaling factors β and λ J were taken as variational parameters. For the many-body groundstate calculation, the inputs of the Slater determinants are the single-particle eigenstates of the non-interacting trapped system φ ho n , and the set of quantum numbers {n} are chosen to give the lowest energy compatible with Pauli's exclusion principle. The Fermi energy F for β = 1 defines an effective Fermi wavenumber k F = √ 2m F /h, whose inverse provides a natural unit for measuring the scattering length and potential range in the interacting many-body problem.

Numerical results for virial relations
In this paper, we concentrate on potential ranges in the interval k F r v ∈ [0.0087, 0.075]. 2 For each r v and several scattering lengths through the crossover, upper bounds of the energyẼ were obtained by optimizing the variational parameters λ EL or λ J and β according to the trial wavefunction. The optimization of the numerical subroutines allowed us to explore higher 6 Figure 1. Many-body energy per particle,Ẽ, equation (15), as a function of the inverse of the scattering length for a potential range k F r v = 0.0109. The energy units correspond to E IFG , the total energy per particle of the non-interacting atomic cloud, and k F denotes the Fermi wavenumber. Many-body energy per particle, E, as a function of the potential range r v at unitarity, 1/k F a = 0. The energy units correspond to E IFG , the total energy per particle of the non-interacting atomic cloud, and k F denotes the Fermi wavenumber. Note that at unitarity, E =Ẽ.
statistics with respect to previous calculations [9] and yield the evaluation of the energy with improved accuracy. In figures 1 and 2, we illustrate the obtained energies; the vertical size of the plotted points is comparable to or higher than the numerical errors. Note that the resulting curves for the energy dependence on the scattering length for a given range and on the range . Partial derivative of the many-body energy per particle,Ẽ, equation (15), with respect to the inverse of the scattering length as a function of the inverse of the scattering length for a potential range k F r v = 0.0109. The energy units correspond to E IFG , the total energy per particle of the noninteracting atomic cloud, and k F denotes the Fermi wavenumber. Partial derivative of the many-body energy per particle,Ẽ, with respect to k F r v as a function of the inverse of the scattering length for a potential range k F r v = 0.0109. The energy units correspond to E IFG , the total energy per particle of the non-interacting atomic cloud, and k F denotes the Fermi wavenumber.
for a given scattering length show a soft structure that allows a numerical interpolation or even an analytical local fitting curve. These interpolations were used to numerically compute the derivatives necessary to verify the virial relations. The results are illustrated in figures 3 and 4. Softened curves for the many-body energy per particle,Ẽ, and the mean value of twice the trapping energy mω 2 R 2 as predicted by the virial relation, equation (17). The dots represent the value of the trapping energy obtained directly from the variational wavefunctions. The energy units correspond to E IFG , the total energy per particle of the non-interacting atomic cloud; k F denotes the Fermi wavenumber.
For all the short interaction ranges explored, both derivatives, ∂Ẽ/∂(1/k F a) and ∂Ẽ/∂(k F r v ), get an extreme value at unitarity.
In figure 5, we show a comparison between the trapping energy curve as predicted by the virial relation, equation (17), and specific values of that energy evaluated directly from the variational functions for a potential range k F r v = 0.0109. Note that there is a good overlap between variational and virial results all over the crossover. The width of the continuous curve corresponds to the numerical errors in its derivation. The main source of error for the trapping energies evaluated directly from the variational wavefunctions results from the non-uniqueness of the variational parameters λ EL and β that yield similar variational energies. In fact, the agreement between the trapping energies evaluated directly and using the virial relation can be used as an additional criterion to select those parameters. As a reference, in figure 5, the many-body energy curveẼ(1/k F a) is also shown. The crossing of virial mω R 2 and total energyẼ curves does not occur at unitarity as a finite range effect. For a lower potential range, k F r v = 0.0087, we have verified that twice the trapping energy is closer to the total energy at the unitary limit.
Although, in this paper, we report the results for N = 165, we have verified that similar results are obtained for N = 120, 84 and 56 particles. That is, independently of the number of particles we considered, it was found necessary to take into account finite range effects to satisfy virial relations within numerical accuracy for potential ranges k F r v 0.0087.
The energy at unitarity as a function of the potential range r v can be used to estimate the parameter ξ = 1 + β Berscht , which is expected to be a universal number for r v → 0. For a harmonically trapped system, ξ is given in terms of the ratio between the total energy per particle E and the total energy per particle of the analogous ideal trapped system E IFG by ξ = (E/E IFG ) 2 . In the past years, there has been a great deal of effort in the community to calculate this parameter [12]. It was found that those calculations performed with lower ratios between the effective range and the interparticle spacing gave lower values of ξ . For k F r v = 0.0065, we obtained a strict upper bound to the energy, yielding ξ = 0.376 ± 0.008. If the interpolation curves for the energy at unitarity as a function of the potential range r v , figure 2, are used to extrapolate to the zero range limit the value of the ratio E/E IFG , then slight variations of the interpolation curves compatible with errors of the calculated energies will yield different limiting values of E/E IFG at r v = 0. Using this procedure, we have obtained that E/E IFG (r v = 0) is in the interval [0.470, 0.498], so that the homogeneous gas parameter ξ would be in the interval [0.221, 0.248]. However, numerical improvements should be made to support further this extrapolated result.

Conclusions
We have studied the virial relations for trapped particles interacting through a two-parameter (intensity and potential range) attractive potential. We have applied those relations to the ground state study of a balanced mixture of two species fermionic trapped atoms when the wavefunction is approximately obtained using simple form variational trial functions. In this way, we are able to quantify the quality of our wavefunctions and, even more important, we can compare the virial relations of short finite range interaction versus contact interactions. Along the calculations, the partial derivatives of the energyẼ as a function of the scattering length and of the potential range were numerically evaluated. It was found that in all the cases considered, those derivatives get large values at unitarity. As a consequence, an accurate determination of the coefficientẼ/E IFG at unitarity using finite range potentials requires an accurate extrapolation procedure. In a similar way, an accurate experimental determination of such a coefficient does require an extremely precise realization of the 1/k F a → 0 limit.