Critical velocity, vortex shedding and drag in a unitary Fermi superfluid
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] J a n Critical velocity, vortex shedding and drag in a unitary Fermi superfluid
F. Ancilotto , , L. Salasnich , and F. Toigo , Dipartimento di Fisica e Astronomia ”Galileo Galilei” and CNISM,Universit`a di Padova, Via Marzolo 8, 35122 Padova, Italy CNR-IOM Democritos, via Bonomea, 265 - 34136 Trieste, Italy
We study the real-time motion of a microscopic object in a cold Fermi gas at unitary conditions byusing an extended Thomas-Fermi density functional approach. We find that spontaneous creationof singly quantized vortex-antivortex pairs occurs as a critical velocity is exceeded, which leads to adrag between the moving object and the Fermi gas. The resulting force is linear in the velocity forsubsonic motion and becomes quadratic for supersonic motion.
PACS numbers: 05.30.Fk, 03.75.Ss, 67.85.-d
I. INTRODUCTION
A Fermi gas of atoms at unitary conditions is charac-terized by a divergent s-wave scattering length [1], whichmakes it a unique system, being at the same time di-lute and strongly interacting. In this regime all scalesassociated with interactions disappear from the problemand the energy of the system is expected to be propor-tional to that of a non interacting fermions system. Theunitary Fermi gas (UFG) at low energies is known tobe superfluid, and the most striking manifestation of itsphase coherence has been the experimental observationof vortices in a rotating UFG of Li atoms [2].Vortex dipoles in a Fermi gas were studied theoreti-cally, within an extended Thomas-Fermi approach, in theBCS limit [3]. Vortex solutions within the same approachhave been studied in Ref. [4], where the rotation of theFermi system was predicted to generate a giant vortexin the presence of strong anharmonicity in the confiningpotential. Recently [5] a time-dependent Bogoliubov de-Gennes (BdG) approach has been used to study the 3-D,real-time formation of vortices in a UFG. Surprisingly,one of the main conclusion in Ref. [5] is that the sys-tem remains superfluid even when stirred at supercriticalspeed. Single vortex solutions in the Ginzburg-Landauregime of a trapped superfluid Fermi gase have been stud-ied in Ref. [6] and Ref. [7]. In a dilute Fermionic super-fluid a vortex state is characterized by a strong densitydepletion along the vortex core. The depletion is howevernot complete, according to the BdG calculations of Refs.[8, 9].Recently it has been remarked [10] that the superfluidunitary Fermi gas can be efficiently described at zero tem-perature by phenomenological density functional theory.Density functionals of different flavours have been pro-posed by different theoretical groups. Bulgac and Yuhave introduced a superfluid density functional based ona Bogoliubov-de Gennes approach to superfluid fermions[11, 12]. Papenbrock and Bhattacharyya [13] have in-stead proposed a Kohn-Sham density functional with aneffective mass to take into account nonlocality effects.Here we adopt instead the extended Thomas-Fermi func-tional of the UFG that we have proposed few years ago [14] and which has been used recently to successfully ad-dress a number of properties of such system [14–18].Our extended Thomas-Fermi (ETF) density functionalapproach for the UFG [14, 15] is based on the followingextended hydrodynamics equations: ∂n∂t + ∇ · ( n v ) = 0 , (1) m ∂ v ∂t + ∇ [ m v + U ( r , t ) + ξ ~ m (3 π ) / n / − λ ~ m ( ∇ √ n ) √ n ] = (2)where n ( r , t ) is the time-dependent scalar density fieldand v ( r , t ) the time-dependent velocity field. Here U ( r , t ) is the external potential and ξ = 0 .
40 and λ = 1 / λ is found to be crucial to accuratelydescribe surface effects, in particular in systems with asmall number of atoms, where the Thomas-Fermi (localdensity) approximation fails [14]. We have shown thatwhen fast dynamical processes occur and/or when shockwaves come into play such term is necessary also in thelarge N limit [18], where results in quantitative agreementwith experiments of real-time collision of strongly inter-acting Fermi gas clouds at unitarity [19] can be obtained[18].We use this method here to study the motion of a mi-croscopic 2-dimensional object with circular shape in aUFG and the associated process of vortex shedding, oncethe object velocity exceeds a critical value. Experimen-tal realization of this geometry would imply, for instance,moving a far-detuned laser beam through a trapped con-densate. II. CRITICAL VELOCITY AND DRAG IN THEUNITARY FERMI GAS
For a dilute Fermi gas the long wavelength elementaryexcitations are sound waves, and the Landau criterionfor the critical velocity v c = ( ǫ/p ) min for breakdown ofsuperfluidity gives v c = c . However, for fluid flowingpast an object (or an atomic impurity as well), the localvelocity near the object surface can become supersoniceven when the fluid velocity far from the object, v ∞ , issubsonic (the maximum velocity of, e.g., a 2-dimensionalflow of an incompressible fluid past a circular obstacleis reached on the perimeter of the obstacle, where it is ∼ v ∞ ).An estimate for the critical velocity for a UFG can beobtained as follows. At stationary conditions Eq. (2)provides the Bernoulli equation for the UFG: (in the fol-lowing α ≡ ξ ~ m (3 π ) / ) − λ ~ m ( ∇ √ n ) √ n + αn / + U ( r ) + m v = const (3)By assuming that the quantum term in the previousequation is negligible in the spirit of the long-wavelengthlimit approximation, and calling n the (uniform) densityfar from an impenetrable object, where v ∼ v ∞ , onefinds: n ( r ) = [ n / + m α ( v ∞ − v )] / (4)outside the object and n = 0 within it (the interactionU between the object and the fluid only provides theexcluded volume condition).The above equation may be written in terms of thelocal c loc ≡ r ξ v F ( r ) = q αn ( r ) / / (3 m ) (5)and bulk c ≡ r ξ v F, ∞ = q αn / / (3 m ) (6)speeds of sound as: c loc = c + ( v ∞ − v ) / v ∼ c loc local instabilities develop, leading to therelease of vortices [20]. The maximum local velocity ofa fluid flowing past an impenetrable cylindrical object,normally to its axis, is v ∼ v ∞ (on the surface of thecylinder and tangent to it). Using this value in (7) onecan solve the equation for v ∞ , thus providing an approx-imate value for the critical velocity of the fluid flowingpast a stationary object (or, equivalently, of a movingobject in the fluid at rest): v c = c / √ v c ∼ p / c obtained for the BEC case, using a similar argument, inRef. [20]. For a spherical object the equatorial velocityis instead v ∼ v ∞ /
2, giving for the critical velocity of asphere moving in a UFG v c = p / c . By combining the two equations (1) and (2) one can de-rive the equation for the momentum conservation (sum-mation over repeated indices is implied): ∂ t J k + ∂ i T ik + ρ∂ k ( U/m ) = 0 (9)where J k = ρv k is the supercurrent density and the stresstensor T ik is defined as: T ik ≡ ρv i v k + δ ik ( 2 α m ρ / ) −
14 ( ~ m ) ρ∂ i ∂ k ln ( ρ ) (10)(to derive the above equation the following identity hasbeen used: ρ − ∂ i ( ρ∂i∂ k ln ( ρ )) = 2 ∂ k ( ∂ i ∂ i √ ρ/ √ ρ )).Note that although in a superfluid there is no frictionalviscosity, nonetheless shear stress may arise from density(pressure) gradients (third term in Eq. (10). As a conse-quence, vortex formation and drag are possible even withno viscosity.The force exerted on a condensate by an obstacle mov-ing through it can be calculated from the rate of momen-tum transfer to the fluid. By integration, one finds thedrag force (per unit mass): F k = ∂ t Z Ω d Ω J k = − Z Σ d Σ n i T ik − Z Ω d Ω ρ∂ k U (11)Here Σ is the object external surface and n is a unitvector directed along the outward normal. In the case ofan homogeneous flow past an impenetrable object onlythe first term contributes. In the present case, however,where a partially penetrable object is used (see the fol-lowing) both contributions are present.In our calculations we used the full expression (10) forthe stress tensor to compute the drag (11). Notice thatthe quantum term in (10) is expected to be negligibleonly when v ∞ ≫ ( ~ /mR ), where R is the diameter ofthe moving object: this is not the case here since thesetwo terms are of comparable magnitude. III. METHODS AND CALCULATIONS
By using a Madelung transformation, equations (1)and (2) can equivalently be written in the form of atime-dependent nonlinear Schr¨odinger equation (NLSE)[14] involving the complex order parameter Ψ( r , t ) = p n ( r , t ) e iθ ( r ,t ) : i ~ ∂∂t Ψ = ˆ H Ψ (12)where ˆ H ≡ − ~ m ∇ + 2 U ( r − V t ) + 2 α | Ψ | / (13)The link between the two descriptions is provided bythe definition v ( r , t ) ≡ ~ m ∇ θ ( r , t ) (14)for the velocity field associated with the phase of theorder parameter θ ( r , t ). V in Eq. (13) is the (constant)velocity of the moving object.By means of a Galilean transformation to the referenceframe moving with the object (which will thus appear asstationary in our simulations), the NLSE to be solvedbecomes: i ~ ∂∂t Ψ = h ˆ H + i ~ V · ∇ i Ψ (15)For simplicity, we model the microscopic object insidethe moving Fermi gas by means of a repulsive cylindri-cal ”wall” and consider the flow motion perpendicularto it. Due to translational invariance along the axis ofthe cylinder (which we define as the z axis), the problemthus reduces to finding the density and the velocity ofthe fluid in the ( x, y ) plane.We have numerically solved this 2-D NLSE equation toobtain the long-time dynamics of the fluid due to the mo-tion of the object moving along the x direction. From thetime-dependent calculated density n ( x, y, t ) and velocity v ( x, y, t ) we computed the drag acting on the object, ac-cording to (11).Our simulation region is a square cell of side a =3 . µm . We used a uniform mesh to represent the wave-function Ψ, on 512 ×
512 points in the ( x, y ) plane. Wehave used the Runge-Kutta-Gill fourth-order method [21]to propagate in time the solutions of the NLSE. To accu-rately compute the spatial derivatives appearing in ourNLSE, we used a 13-point finite-difference formula [22].To avoid the outgoing sound waves and/or emitted vor-tices to interfere with the fluid dynamics after being re-flected on on the grid boundaries, we use an exponentialabsorbing buffer located in the periphery of the cell, asdescribed in Ref. [23]. The waves can travel freely in theundamped region (which occupy most of the simulationcell) but are quickly attenuated as they enter the exter-nal region, thus preventing unwanted interferences whichmight spoil the results.In the following we take m equal to the mass of a Liatom. We consider two different density values (”highdensity” and ”low density” in the following) for the UFG,namely n = 8600 µm − and n = 880 µm − , corre-sponding to interparticle distances d ∼ × − µm and d ∼ . µm , respectively.At equilibrium, the repulsive cylindrical ”wall” resultsin the formation of a circular cavity void of atoms. Theobject potential U and the associated initial density pro-file are shown in Fig. (1). IV. RESULTS AND DISCUSSION
We made a series of time-dependent calculations, solv-ing Eq. (15) for different values of the speed V , andfollowing the dynamics of the systems for several µs . Wefind that there is a critical value v c (which will be quanti-fied in the following) for the object speed V which sepa-rates two distinct regimes. When V < v c the fluid profile FIG. 1: Equilibrium density profile at t=0. The dotted linerepresents (in arbitrary units) the repulsive potential due tothe object. rapidly evolves with time into a stationary configuration.Both the density and the velocity field for such final con-figuration are fore-to-aft symmetric. This implies thatthe drag force on the object is zero, which is another ver-sion of the well-known D’Alembert paradox in classicalfluids.In Fig. (2) (lower curve) we show the calculated dragforce F x (computed using Eq. (11)) as a function of timewhen V < v c . After a transient where the fluid acco-modates under the sudden acceleration of the object, thedrag goes eventually to zero when the symmetric (sta-tionary) configuration is reached.Above v c , however, the Fermi gas dynamics changesdramatically: linear (antilinear) vortices are sponta-neously created in pairs close to the surface of the object,together with the emission of sound waves.We show in Fig. (3) a sequence of images taken dur-ing the real-time evolution of the system when V > v c .The object is moving from left to right. First a localizedbow wave moving with supersonic velocity is emitted infront of the cylinder and rapidly moves ahead. This isthe result of a ”shock” wave produced by the sudden ac-celeration of the object (we recall that the initial stateis with fluid at rest and the object moving at constantvelocity).Then a vortex pair is emitted in the rear, trailing be-hind (meaning that the pair velocity is lower than that ofthe moving object). The initial pair is soon overtaken bya pair emitted successively: the vortex lines apparentlyform a temporarily bound state (but they will eventuallysplit at later times, not shown).As shown in the upper curve of Fig. (2), in this case FIG. 2: Drag force during the real-time evolution of the sys-tem. Upper curve:
V > v c , lower curve: V < v c .FIG. 3: From left to right, top to bottom: configurationsat increasing times for V /c = 0 .
76. Only a portion of thesimulation cell is displayed. the instantaneous calculated drag force relaxes, after atransient, towards a finite nonzero value (with oscilla-tions that reflect the quasi-periodic emission of vortexpairs).In the simulations of Fig. (3) the velocity V is greaterthan v c , but still below the speed of sound c .The sequence shown in Fig. (4) is instead obtained for FIG. 4: From left to right, top to bottom: configurations atincreasing times for velocity
V /c = 1 .
44. Only a portion ofthe simulation cell is displayed. a supersonic motion of the object,
V /c = 1 .
44. It ap-pears that the vortex shedding frequency is considerablyincreased with respect to the previous case (similarly towhat happens for a BEC, where the shedding frequency is[24] ∝ V . Vortex-antivortex pairs are emitted in a semi-continuous way on different part of the rear section of theobject, leading to parallel rows of connected vortical linesthat eventually decay into separated vortices. Note alsothe appearance of a rather structured bow sound wavepattern moving along with the object. As time proceeds,the wake region behind the object becomes turbulent dueto the superposition of more and more vortices and soundwaves.In order to better characterize the vortex structure, wehave selected (from another simulation where V is justabove v c , so that vortex pairs are emitted with low fre-quency and thus are well separated from one another)a configuration after a pair of vortex lines have beenemitted and moved away from the object, and closelyanalyzed the vortex structure. We find that the vor-tex is singly quantized, and that the two vortices of thepair have opposite polarization. The velocity field (notshown) follows very closely the ideal vortex velocity pro-file, v ( r ) = ~ / (2 mr ) (here r = p x + y ).The vortex line structure has an empty core, as shownin Fig. (5), while the core size scale is set by k − F . Thisis in contrast with calculations based on Bogoliubov-DeGennes calculations [8, 9] where a partially filled core(between 0.2 and 0.3 n ) is predicted due to the presenceof some normal liquid coming from pair breaking in thecore region, where the velocity is higher. We will discussin the following the reason for such discrepancy. FIG. 5: Density profile in the vortex core region.FIG. 6: Current density distribution near the vortex core.Squares: our results. Solid line: BdG calculations from Ref.[9].
We computed the superfluid current density j = ρ v cir-culating around the vortex core (see Fig. (6)), and founda peak value at r max , in a remarkable overall agreementwith the T=0 Bogoliubov-De Gennes calculations of Ref.[9] (shown in Fig. (6) with a solid line).We find the agreement with the results of Ref. [9] par-ticularly rewarding since it shows that an important su-perfluid observable related to vorticity can indeed be de- scribed accurately by our Density Functional approach.The scale of the maximum circulating current is set bythe critical velocity which is determined by pair-breakingon the BCS side and by collective excitations on the BECside [9]. This is why our approach, which cannot obvi-ously describe single-particle processes like pair-breaking,is nonetheless able to reproduce the current pattern. For r < r max the kinetic energy cost associated with the cur-rent flow becomes larger than the condensation energy.Thus r max gives an estimate of the distance from the corecenter below which superfluidity is partially suppressed.We note that the peak position r max in Fig. (6) coin-cides with the value of r where the vortex velocity fieldbecomes equal to the local sound velocity, i.e. when ~ mr ∼ q αn ( r ) / / m (16)where n ( r ) is the vortex core density profile shown in Fig.(5).In Ref. [16] we have found that, at unitarity, our ap-proach leads to a maximum Josephson current across abarrier which is practically the same as the one obtainedfrom BdG. This suggests that, at unitarity, the currentis limited by Landau’s criterion for the creation of col-lective excitations, and not by single particle excitations.It is then not surprising that also in the present case ourmaximum current is close to that obtained by a muchmore demanding microscopic calculation based on BdGequations ([8, 9])It must be said that in Ref. [8, 9] the vortex struc-ture is imposed on the condensate order parameter ∆.Therefore near the core the current goes to zero because∆ vanishes and the superfluid density vanishes with it,since ρ s ∝ ∆ . However in the core region, where thesuperfluid velocity exceeds the critical value for the cre-ation of single particle excitations, ρ s < ρ , so that thetotal density is non zero even at the core where ρ s = 0.On the contrary, in our approach which does not takeinto account single particle excitations, all the fluid issuperfluid, and therefore the vanishing at the core of thesuperfluid current due to vorticity implies that the totaldensity is zero there.Results similar to those reported above were obtainedfrom the simulations of the lower density system, themain difference being the shallower vortex cores (whichscale as k − F ). The calculated superfluid current densityaround the core of an isolated vortex, plotted as a func-tion of k F r , is indistinguishable from the one shown inFig. (6).We show in Fig. (7) the calculated drag force onthe object, as a function of the velocity V . Each valueis obtained as a time-average [25] of curves like theone shown in Fig. (2) (the average is taken over atime interval where the drag force has already reacheda plateau). From these results, a value of the criticalvelocity v c ∼ . c is obtained, in agreement with thesimple estimate (8). FIG. 7: Average drag force for the high and low density sys-tems, plotted as a function of V . We also find that the drag force first increases linearlywith the velocity (”Stoke’s law”, usually associated withlaminar drag) and then turns to a quadratic behavior(”Newton’s law”, usually associated with turbulent flow),since at supersonic velocities there is also a contributionto the drag associated with sound radiation. A similarbehavior is observed in the case of an impenetrable cylin-der moving inside a dilute BEC [25].The phenomenology of the dissipative motion of anobject displaced through a UFG, as it appears from ourcalculations, is qualitatively similar in many aspects, inspite of the different non-linear interactions, to the be-havior observed in the case of an object moving in a BEC[20, 25]: the occurrence of vortex emissions in pairs andthe associated density patterns are similar in the twosystems, and also the behavior of the drag as a functionof the object velocity. There are important differences though: the emitted vortices in the UFG are doublyquantized, as expected from fermion pairing; the pre-dicted critical velocity is different; the vortex core struc-ture is also different, scaling in the present case as theinverse Fermi momentum.
V. CONCLUSIONS
In conclusion, we have numerically studied the mo-tion of an object in the ultracold unitary Fermi gas. Wedescribed the system by using an extended density func-tional approach, which has been used recently to success-fully describe a number of static and dynamical proper-ties of cold Fermi gases.We find that quantized vortices are spontaneously gen-erated in pairs during the time evolution at supercriticalvelocities. Moreover, the profile of the current density asa function of the distance from the core is quantitativelyclose to the one found by the much more demanding so-lution of the BdG solutions. We explain this agreementby observing that at distances larger than the one corre-sponding to the maximum current density, in both treat-ments the superfluid density coincides with the total den-sity, since the speed of the fluid is below its critical valuefor the creation of single particle excitations. At shorterdistances, in the core region, the current density is againsimilar in the two treatments: in both it vanishes imply-ing that the superfluid density is zero at the center ofthe vortex. The main difference is that while in our casethe superfluid density coincides with the total density, inthe BdG treatment the density may have a contributionfrom a ”normal” component related to single particle ex-citations and this component provides a nonzero densityat the center.We acknowledge Fondazione CARIPARO (project ofexcellence ”Macroscopic Quantum Properties of Ultra-cold Atoms under Optical Confinement”) and Universityof Padova (progetto di ateneo ”Quantum Informationwith Ultracold Atoms in Optical Lattices”) for partialsupport. [1] “The Many-Body X Challenge Problem”, formulated byG.F. Bertsch, see R. A. Bishop, Int. J. Mod. Phys. B ,issue: 10-11, iii (2001).[2] M.W. Zwierlein et al., Nature (London) , 1047(2005).[3] S. Gautam, arXiv:1205.5670v1.[4] E. Lundh and A. Cetoli, Phys. Rev. A , 023610 (2009).[5] A. Bulgac, Y.L. Luo, P. Magierski, K.J. Roche, and Y.Yu, Science , 1288 (2011).[6] M. Rodriguez, G-S. Paraoanu, and P.Torma, Phys. Rev.Lett. , 100402 (2001).[7] G.M. Bruun and L. Viverit, Phys. Rev. A , 063606(2001).[8] A. Bulgac and Y. Yu, Phys. Rev. Lett. , 190404 (2003). [9] R. Sensarma, M.Randeria and T-L. Ho, Phys. Rev. Lett. , 090403 (2006).[10] A. Bulgac, M. McNeil Forbes, and P. Magierski, in BCS-BEC Crossover and the Unitary Fermi Gas, LectureNotes in Physics, vol. , Ed. by W. Zwerger (Springer,Berlin, 2011).[11] A. Bulgac and Y. Yu, Phys. Rev. Lett. , 190404 (2003).[12] A. Bulgac, Phys. Rev. A , 040502(R) (2007).[13] T. Papenbrock, Phys. Rev. A , 041603(R) (2005); A.Bhattacharyya and T. Papenbrock, Phys. Rev. A ,041602(R) (2006).[14] L. Salasnich and F. Toigo, Phys. Rev. A , 053626(2008); L. Salasnich and F. Toigo, Phys. Rev. A ,059902(E) (2010). [15] S.K. Adhikari and L. Salasnich, Phys. Rev. A , 043616(2008); L. Salasnich, Laser Phys. , 642 (2009); S.K.Adhikari and L. Salasnich, New J. Phys. , 023011(2009).[16] F. Ancilotto, L. Salasnich and F. Toigo, Phys. Rev. A , 033627 (2009); L. Salasnich, F.Ancilotto, N. Maniniand F.Toigo, Laser Physics , 636 (2009).[17] L. Salasnich and F. Toigo, J. Low Temp. Phys. , 239(2011).[18] F. Ancilotto, L. Salasnich, and F. Toigo, Phys. Rev. A , 063612 (2012).[19] J. Joseph, J. Thomas, M. Kulkarni, and A. Abanov,Phys. Rev. Lett. , 150401 (2011). [20] T. Frisch, Y. Pomeau, and S. Rica, Phys. Rev. Lett. ,1644 (1992).[21] Mathematical Methods for Digital Computers , Ed. by A.Ralston and H. S. Wilf, vol. 1, p. 117 (Wiley, New York,1960).[22] M. Pi, F. Ancilotto, E. Lipparini, and R. Mayol, PhysicaE, , 297 (2004).[23] D. Jin and W. Guo, Phys. Rev. B , 094524 (2010).[24] T. Winiecki, J.F.McCann, and C.S.Adams, Phys. Rev.Lett. , 8186 (1999).[25] T. Winiecki, B. Jackson, J.F. McCann, and C.S.Adams,J. Phys. B: At. Mol. Opt. Phys.33