Comment on “Superfluid stability in the BEC-BCS crossover” arXiv:cond-mat/0608172v2 [cond-mat.supr-con] 28 May 2007 Daniel E. Sheehy and Leo Radzihovsky Department of Physics, University of Colorado, Boulder, CO, 80309 (Dated: August 7, 2006) We point out an error in recent work by Pao, Wu, and Yip [Phys. Rev. B 73, 132506 (2006)], that stems from their use of a necessary but not sufficient condition [positive compressibility (magnetic susceptibility) and superfluid stiffness] for the stability of the ground state of a polarized Fermi gas. As a result, for a range of detunings their proposed ground-state solution is a local maximum rather than a minimum of the ground state energy, which thereby invalidates their proposed phase diagram for resonantly interacting fermions under an imposed population difference. There has been considerable recent interest in paired superfluidity of fermionic atomic gases under an imposed spin polarization1,2 , i.e., when the numbers N↑ and N↓ of the two atomic species undergoing pairing are different. Along with the detuning δ of the Feshbach resonance (controlling the strength of the interatomic attraction), the population difference ∆N = N↑ − N↓ is an experimentally-adjustable “knob” that allows the study of novel regimes of strongly-interacting fermions. A crucial question concerns the phase diagram of resonantly interacting fermions as a function of δ and imposed ∆N . Two early theoretical studies that have addressed this issue are the work by Pao, Wu, and Yip3 on the one-channel model of interacting fermions and our work4 on the two-channel model of interacting fermions, with the details and extensions of the latter presented in our recent preprint.5 Apart from the Fulde-FerrellLarkin-Ovchinnikov (FFLO) phase appearing over a thin sliver of the phase diagram in Ref. 4 (explicitly not considered by Pao, Wu, and Yip), for the positive-detuning BCS and crossover regimes the phase diagrams in these manuscripts are in qualitative agreement6 . However, at negative detuning in the BEC regime, the phase diagrams are qualitatively different. In particular, in Fig. 4 of Ref. 3, stable superfluidity is claimed to exist above and to the left of a nearly vertical and positively sloped phase boundary that we re-plot here in Fig. 1 as a dashed line at negative detuning. In qualitative contrast, in our work, Ref. 4, we found that stable superfluidity exists only below and to the left of a negatively-sloped phase boundary (see Fig. 1). What is the source of this qualitative discrepancy? Although Ref. 3 and Ref. 4 use different models of resonantly-interacting fermions, the close relationship between the one- and two-channel models (particularly within the mean-field approximation) implies that they should yield qualitatively similar phase diagrams. Furthermore, we have extended our original two-channel model study4 to that of a one-channel model5 and, as expected, found results in qualitative agreement with those in our Letter4 , but in disagreement with that of Ref. 3. Indeed, as we explicitly show here, the origin of the discrepancy is that the criterion for stability of the superfluid phase used in Ref. 3 (based on positivity of magnetic susceptibilities) is a necessary but not sufficient condition for stability and does not ensure7,8,9 that the state
N N 1 Normal 0.8 Phase Separation 0.6  FFLO 0.4 0.2 SFM x 1 -2 -1 0 1 2 kF as FIG. 1: (Color online) Polarization ∆N/N vs. detuning − kF1as phase diagram of the one-channel model (appropriate for width γ ≫ 1) showing regimes of FFLO, superfluid (confined to the ∆N = 0 axis), magnetized superfluid (SFM ), and phase separation. The dashed lines are the purported phase boundaries reported in Ref. 3. As we show5 , for the entire shaded region single-component (uniform) solutions to the stationarity conditions Eqs. (7) are local maxima of EG (∆) (and therefore unstable, leading to phase separation), shown for the point marked with the “x” in Fig. 2. is even a local minimum of the energy. This thereby leads to incorrect phase boundaries in both the BEC and BCS regimes, although below we shall focus on the phase boundary inside the BEC regime, where the error is qualitative and most pronounced. Thus, much of what is claimed to be a “stable superfluid” in Ref. 3 (in the negative-detuning BEC regime) is actually unstable to phase separation. This is illustrated in the correct T = 0 mean-field phase diagram, Fig. 1, for the one-channel model (quantitatively consistent with other recent work10,11 and derived in detail in Ref. 5), plotted as a function of the dimensionless parameter −(kF as )−1 (proportional to the Feshbach resonance detuning δ, with kF the Fermi wavevector and as the s-wave scattering length) and the polarization ∆N/N = (N↑ − N↓ )/(N↑ + N↓ ). Our starting point is the single-channel model Hamiltonian (studied in Ref. 3) for two resonantly-interacting 2 species of fermion ĉkσ (σ =↑, ↓): Ĥ = X k,σ ǫk ĉ†kσ ĉkσ + λ X † † ĉk↑ ĉp↓ ĉk+q↓ ĉp−q↑ , V (1) kqp where ǫk = k 2 /2m, m is the fermion mass, and V is the system volume. A broad Feshbach resonance is modeled by an attractive interaction λ < 0, the magnitude of which increases with decreasing Feshbach resonance detuning. The equilibrium ground state of a many particle system in the grand-canonical ensemble (with chemical potentials µ↑ and µ↓ ) at T = 0 is characterized by the grand thermodynamic potential Ω(µ↑ , µ↓ ) defined by12 h i Ω(µ↑ , µ↓ ) = min hĤi , (2) with the grand-canonical Hamiltonian Ĥ ≡ Ĥ − µ↑ N̂↑ − µ↓ N̂↓ . (3) P † Here, N̂σ ≡ k ĉkσ ĉkσ is the number operator for fermion species σ and the minimization in Eq. (2) is over all possible ground states. The standard (BCS-type) mean-field approximation that we shall utilize (as done in Ref. 3) amounts to assuming a restricted class of possible many-body groundstates self-consistently parametrized by the ground-state expectation value ∆ = λhĉ↓ (r)ĉ↑ (r)i. (4) For pairing amplitude ∆ 6= 0, this class includes both weakly-paired BCS-type and strongly-paired molecular BEC-type pairing order. For ∆ = 0 it is an unpaired Fermi gas. Once this mean-field approximation has been made, it is straightforward4,5 to compute the variational groundstate energy, which is the expectation value EG (∆) = hĤi. Evaluating this expectation value, and converting momentum sums to integrals, we find Z ∆2 d3 k m 2 (ξ − E + ) ∆ + EG (∆)/V = − k k (2πh̄)3 2ǫk 4πas h̄2 Z d3 k + (Ek − h)Θ(h − Ek ), (5) (2πh̄)3 p where ξk ≡ ǫk − µ, Ek ≡ ξk2 + ∆2 , and we have defined the chemical potential µ = 21 (µ↑ + µ↓ ) and the chemical potential difference h = 12 (µ↑ − µ↓ ). We have also exchanged the bare interaction parameter λ for the vacuum s-wave scattering length as given by Z m 1 d3 k 1 = . (6) + 2 λ (2πh̄)3 2ǫk 4πas h̄ The determination of the mean-field phase diagram using Eq. (5) is conceptually quite simple. According to Eq. (2), the ground state at a particular µ↑ and µ↓ (or, equivalently, µ and h) is given by the minimization of EG (∆) with respect to the pairing amplitude ∆ that can be taken to be real. Any such minima of course satisfy the stationarity constraint or gap equation [equivalent to Eq. (4)] 0= ∂EG , ∂∆ (7a) where we emphasize that the derivative is taken at fixed µ and h. Since experiments are conducted at fixed atom number, we must augment Eq. (7a) with the number constraint equations Nσ = hN̂σ i. Examining Eq. (3), we see that the constraints can be rewritten as ∂EG , ∂µ ∂EG ∆N = − , ∂h N = − (7b) (7c) with the total particle number N = N↑ + N↓ and population difference ∆N = N↑ − N↓ . Our key point (apparently missed by the authors of Ref. 3) is that not every simultaneous solution of the gap and number equations, Eqs. (7), corresponds to a physical ground state of the system; the additional criterion is that the solution ∆ must also be a minimum of EG (∆) at fixed µσ . The verification that an extremum solution is indeed a minimum is particularly essential when there is the possibility of a first-order transition, with EG (∆) exhibiting a local maximum that separates local minima, as is the case for a polarized Fermi gas, studied here and in Refs. 3,4,5. Analyzing Eqs. (7), we find that for sufficiently large ∆N in the positive-detuning BCS regime and for ∆N = N in the BEC regime, a solution to Eqs. (7) may be found that minimizes EG (∆) at ∆ = 0, indicating a polarized normal phase. Also, for sufficiently low detuning in the BEC regime, a polarized molecular superfluid (SFM ) solution exists that minimizes EG (∆) at ∆ 6= 0. A more general analysis4,5,13 shows that a periodicallypaired FFLO solution is the ground state over a thin range of polarization values at sufficiently large positive detuning. However, over the large shaded portion of the phase diagram, at intermediate detuning and polarization, we find that it is not possible to satisfy Eqs. (7) with a (homogeneous, single component) minimum of EG (∆). For the corresponding range of parameters the system phase-separates14 into two coexisting ground states (that are degenerate minima of EG [∆]). The resulting phaseseparated state can be explicitly accounted for by generalizing the ground-state ansatz to include the possibility of such an inhomogeneous mixture5 . The contrasting strategy of Pao, Wu, and Yip3 is to find solutions of Eqs. (7) for all values of N , ∆N , and − kF1as , some of which do not correspond to ground states. The unphysical (unstable) solutions are then discarded 3 EG =F -0.2 -0.25 -0.3 -0.35 0.5 1 1.5 2
=F FIG. 2: Plot of EG (∆) at fixed µ = −0.013ǫF and h = 1.16ǫF , normalized to the Fermi energy ǫF , for coupling g = kF1as = 0.63, so that Eqs. (7) yield a stationary solution with ∆N/N = 0.2 and ∆ = 1.02. As seen in Fig. 2, however, this solution is a local maximum (saddle point) of EG (∆). At this point ( kF1as = 0.63 and ∆N/N = 0.2) in the phase diagram (illustrated by an “x” in Fig.1) the actual mean-field ground state is a phase-separated mixture of a superfluid and a normal state. (thereby determining the phase boundaries plotted in Ref. 3) based on criteria of the positivity of magnetic susceptibilities (atomic compressibilities) and the superfluid stiffness. However, as we now discuss, these stability criteria are necessary but not sufficient (i.e., not generally restrictive enough) to ensure that a solution to Eqs. (7) is indeed a minimum of EG (∆). In contrast, it can be shown that the converse is true, namely that the minimization of EG (∆) ensures that the local susceptibilities are positive definite.15 The stability criteria16 used by Pao, Wu, and Yip can be understood by examining Fig. 3 of Ref. 3. The solid lines in this figure correspond to solutions of Eqs. (7) at different values of the coupling g ≡ kF1as . In particular they plot h/ǫF (with ǫF the Fermi energy, related to the √ 3/2 density n = N/V by n = 34 cǫF with c = m3/2 / 2π 2 h̄3 ) as a function of the polarization nd /n = ∆N/N (at fixed density), where nd = ∆N/V is the magnetization. At positive and intermediate detunings (the bottom curves of Fig. 3 of Pao, Wu, and Yip), they find solutions satisfy∂h ing ∂n < 0 and correctly conclude that such solutions d n (having a negative magnetic susceptibility) are unstable. However, at sufficiently negative detuning (g > ∼ 0.5) Pao, Wu, and Yip. find solutions to Eqs. (7) with a positive ∂h > 0, and based on their criterion (ersusceptibility ∂n d n roneously) conclude that these solutions indicate a stable magnetic superfluid ground state. They then define a phase boundary in the BEC regime (the leftmost dashed curve of Fig. 1), to a stable magnetized superfluid, by ∂h where ∂n changes sign. d n However, our explicit calculation of EG (∆) (plotted in Fig. 2) for one such solution (with kF1as = g = 0.63 and nd /n = 0.2, indicated with an “x” in Fig. 1 and corresponding to a point on the uppermost solid curve of Fig. 3 of Ref. 3), purported by Pao, Wu, and Yip to be stable (to the left of their proposed stability boundary), shows that in fact this solution (indicated with a dot in Fig. 2) is a local maximum and therefore does not represent a ground state. This solution was obtained by numerically solving Eqs. (7) at g = 0.63 and ∆N/N = 0.2, yielding µ = −0.013ǫF, h = 1.16ǫF and ∆ = 1.02ǫF, the latter two values consistent with Figs. 2 and 3 of Pao, Wu, and Yip, showing that we are indeed reproducing a solution claimed to be stable by Pao, Wu, and Yip. Thus, the method used by Pao, Wu, and Yip has not correctly located the global minimum of the ground-state energy; indeed, it has not even found a local minimum. Although it might appear from the plot of EG (∆, µ, h) (Fig. 2) that the true ground state is an unpaired (∆ = 0) normal state, this state does not satisfy Eqs. (7b) and (7c); thus, it is also not the ground state at this coupling and polarization. Indeed, as noted above, we find that it is impossible to minimize EG (∆) while satisfying Eqs. (7b) and (7c) at this coupling and polarization (marked by an “x” in Fig. 1), nor anywhere inside the shaded region in the phase diagram Fig. 1, indicating the absence of a uniform solution. The true mean-field ground state everywhere in the shaded region is a phaseseparated mixture of two phases of different densities in chemical equilibrium such that the total number and polarization constraints are satisfied4,5 . We note that a ground state determined by minimizing EG (∆) at a particular µ↑ and µ↓ automatically satisfies the condition of having a positive magnetic susceptibility (compressibility). Indeed, it is straightforward to generally show15,17 that, since Ω(µ↑ , µ↓ ) is concave downwards12 , the eigenvalues of the susceptibility matrix   2 ! ∂2Ω ∂ Ω ∂N ∂N χ≡ ↑ ↑ ∂µ↑ ∂N↓ ∂µ↑ ∂µ↓ ∂N↓ ∂µ↓ = − ∂µ2↑ ∂2Ω ∂µ↑ ∂µ↓ ∂µ↑ ∂µ↓ , ∂2Ω ∂µ2↓ (8) are positive in the ground state. An equivalent stability criterion was derived in Ref. 18 by considering stability against local density variations. The procedure used by Pao, Wu, and Yip, however, did not amount to analyzing Eq. (8); as discussed above, the phase diagram of Pao, ∂h Wu, and Yip was obtained by computing ∂n which is d n not an equivalent condition. If the authors had instead d studied ∂n ∂h µ , they would have found that the solution plotted in Fig. 2 has a negative magnetic susceptibility in the grand-canonical ensemble (and therefore is unstable). However, we must emphasize the important point that any particular extremum solution may have a positive magnetic susceptibility and still not be the ground state. The simplest example of this is the normal Fermi gas state (∆ = 0), which satisfies the gap and numberconstraint equations everywhere in the phase diagram (including ∆N = 0) and has a positive magnetic susceptibility, but is only the actual ground state (a minimum of EG (∆)) at sufficiently large ∆N . Thus, if the authors of Ref. 3 had computed the eigenvalues of Eq. (8) instead ∂h of ∂n , they would have been able to discard some of d n 4 the erroneous solutions plotted in Fig. 3 of Ref. 3. However, in general , Eq. (8) is still not sufficient and the most correct scheme is to use Eq. (2), i.e., to find the global minimum, in the grand-canonical ensemble, of the mean-field ground-state energy. We conclude by noting that, although a mean-field analysis of the one-channel model is not expected to be quantitatively accurate near the resonance position where kF |as | → ∞, it is expected to yield a qualitatively correct description of a polarized resonantly-interacting Fermi gas. To summarize, we have shown that while for 1 2 3 4 5 6 7 8 9 10 11 12 13 M.W. Zwierlein, A. Schirotzek, C.H. Schunck, and W. Ketterle, Science 311, 492 (2006). G.B. Partridge, W. Li, R.I. Kamar, Y. Liao, and R.G. Hulet, Science 311, 503 (2006). C.-H. Pao, S.-T. Wu, and S.-K. Yip, Phys. Rev. B 73, 132506 (2006). D.E. Sheehy and L. Radzihovsky, Phys. Rev. Lett. 96, 060401 (2006). D.E. Sheehy and L. Radzihovsky, cond-mat/0607803, Ann. Phys. (N.Y.) (to be published). Although in this Comment we have focused on the BEC regime where the failure of the stability criteria used by Pao, Wu, and Yip3 is most drastic, the inadequacy is also obvious (although less drastic) on the BCS side, yielding a quantitatively incorrect phase boundary, as illustrated by dashed curves in Fig. 1. Unfortunately, following Ref. 3, this erroneous (insufficient for stability) criterion was also more recently used by other authors yielding a similarly incorrect phase diagram; see for example Ref. 8. As also discussed in Ref. 9, such authors incorrectly identify spinodals (where local stability of a metastable solution is lost) as being phase boundaries. Such an error yields, e.g., the rightmost dashed-line boundary in Fig. 1. M. Iskin and C.A.R. Sá de Melo, Phys. Rev. Lett. 97, 100404 (2006). A. Lamacraft and F.M. Marchetti, cond-mat/0701692 (unpublished). Z.-C. Gu, G. Warner, and F. Zhou, cond-mat/0603091 (unpublished). M.M. Parish, F.M. Marchetti, A. Lamacraft and B.D. Simons, Nat. Phys. 3, 124 (2007). M.M. Forbes, E. Gubankova, W.V. Liu, and F. Wilczek, Phys. Rev. Let. 94, 017001 (2005). In Refs. 4,5, we use a more general variational ansatz that is equivalent to assuming λhĉ↓ (r)ĉ↑ (r)i = ∆Q eiQ·r , allowing the possibility of a periodically-modulated FFLO-type ground state. Since the more general Q 6= 0 analysis, presented in Ref. 5 for the one-channel model, only yields stable FFLO states for a very thin window of parameters at positive detuning (indicated in red in Fig. 1), for simplicity (and because this assumption was also made by Pao, Wu, equal species number (∆N = h = 0) such analysis can simply proceed by solving the gap and number equations [Eqs. (7)], the existence of first-order transitions at h 6= 0 implies that the ground-state energy EG (∆) exhibits local maxima as a function of ∆, yielding solutions to Eq. (7a) that do not represent physical ground states. Acknowledgments — We gratefully acknowledge discussions with V. Gurarie and M. Veillette as well as support from NSF DMR-0321848 and the Packard Foundation. 14 15 and Yip3 ) here we assume Q = 0 at the outset. P.F. Bedaque, H. Caldas, and G. Rupak, Phys. Rev. Lett. 91, 247002 (2003). The concavity of Ω(µ↑ , µ↓ ) [Eq. (2)], implying positive eigenvalues of the susceptibility matrix Eq. (8) in the ground-state, is exactly true following general arguments12,17 and is also true within the mean-field approximation. To show the latter, we note that, within meanfield theory, Ω(µ↑ , µ↓ ) is a minimization with respect to ∆ of EG [∆, µ↑ , µ↓ ], which, in turn, is a concave function (at fixed ∆) of µ↑ and µ↓ since it is the exact groundstate energy for the mean-field Hamiltonian. The fact that the mean-field Ω(µ↑ , µ↓ ) is a minimization of the concave function EG [∆, µ↑ , µ↓ ] with respect to ∆ then implies it is also concave. To illustrate this, consider for simplicity setting h = 0 so that the thermodynamic potential Ω(µ) is a function only of µ. Then one can show that the number susceptibility dN/dµ = −d2 EG [∆0 (µ), µ]/dµ2 satisfies [with ∆0 (µ) the minimum of EG [∆, µ]] dN ∂ 2 EG ∂ 2 EG “ ∂∆ ”2 ˛˛ . =− + 2 ∆=∆0 dµ ∂µ ∂∆2 ∂µ 16 17 18 The second term of this equation is clearly positive at the minimum, and the first term is positive by virtue of the above concavity argument, guaranteeing that dN/dµ > 0 (positive compressibility) in the ground state. However, if ∆0 is only a stationary point (e.g., a local maximum), this equation still holds and shows that dN/dµ can be either positive or negative depending on the relative magnitude of the two terms. This change in sign of dN/dµ (as takes place, for dnd /dh, across the dashed curve in Fig.1) must not be interpreted as a transition to a stable superfluid phase as was done by Pao, Wu, and Yip3 . Pao, Wu, and Yip (Ref.3 ) also consider a criterion based on the positivity of the superfluid stiffness but find this criterion to be even less stringent than the positive magnetic susceptibility criterion that we focus on here. G.L. Sewell, Quantum Mechanics and its Emergent Macrophysics, Princeton University Press, 2002. L. Viverit, C.J. Pethick, and H. Smith, Phys. Rev. A 61, 053605 (2000).