Supersonic and subsonic shock waves in the unitary Fermi gas
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] O c t epl draft Supersonic and subsonic shock waves in the unitary Fermi gas
Luca Salasnich
Dipartimento di Fisica “Galileo Galilei” and CNISM, Universit`a di Padova, Via Marzolo 8, 35131 Padova, Italy
PACS – Degenerate Fermi gases
PACS – Shock wave interactions and shock effects
Abstract. - We investigate shock waves in the unitary Fermi gas by using the zero-temperatureequations of superfluid hydrodynamics. We obtain analytical solutions for the dynamics of alocalized perturbation of the uniform gas. These supersonic bright and subsonic dark solutionsproduce, after a transient time, an extremely large (divergent) density gradient: the shock. Wecalculate the time of formation of the shock and also simulate the space-time behavior of the wavesby solving generalized hydrodynamic equations, which include a reliable dispersive regularizationof the shock. We find that the shock spreads into wave ripples whose properties crucially dependon the chosen initial configuration.
One of the basic problems in physics is how density per-turbations propagate through a material [1,2]. In additionto the well-known sound waves, there are shock wavescharacterized by an abrupt change in the density of themedium [1, 2]. Shock waves are ubiquitous and have beenstudied in many different physical systems [1,2]. Ten yearsago shock waves have been experimentally observed alsoin atomic Bose-Einstein condensates (BECs) [3–6], andtheoretically investigated in various BEC configurations[7–14]. Very recently the observation of nonlinear hydro-dynamic waves has been reported in the collision betweentwo strongly interacting Fermi gas clouds of Li atoms [15].The experiment shows the formation of density gradients,which are nicely reproduced by hydrodynamic equationswith a phenomenological viscous term [15]. Nevertheless,the role of dissipation is questionable [16] since the ultra-cold unitary Fermi gas is noted as an example of an almostperfect fluid [17]. Indeed in Ref. [16] it has been shown,by solving zero-temperature time-dependent Bogolibov-deGennes equations, that the viscous term is not necessaryto reproduce the experimental results of Ref. [15].Here we investigate the formation and dynamics ofshock waves in the dilute and ultracold unitary (diver-gent inter-atomic scattering length) Fermi gas by usingthe zero-temperature equations of superfluid hydrodynam-ics [1]. At zero temperature fermionic superfluids in theBCS-BEC crossover can be modelled by hydrodynamicequations [18] and their generalizations with a gradientterm [19] that induces a dispersive regularization of theshock. In this paper we obtain analytical solutions for the dynamics of a localized perturbation of the uniformgas. We calculate the supersonic (or subsonic) velocity ofpropagation of these bright (or dark) perturbations. Weshow that bright perturbations evolve towards a shock-wave front, while dark perturbations produce the shock intheir back, and we calculate the period T s of formation ofthe shock. In addition, we study the space-time behaviorof the shock waves beyond this characteristic time T s byincluding a reliable quantum correction in the hydrody-namic equations [19]. In fact, according to the two-fluidmodel of Landau [1], the viscous term acts only on thenormal component of the fluid, and at zero temperaturethe normal component is zero. Morover, recent theoreticalmicroscopic calculations [20] suggest that the viscosity ofthe unitary Fermi gas is extremely small at very low tem-peratures because the transverse current does not coupleto collective modes. By solving numerically these general-ized equations of the unitary Fermi gas we show that theshock-wave front spreads into wave ripples whose proper-ties crucially depend on the “brightness” (bright or dark)of the chosen initial configuration. We expect that our re-sults are reliable when the normal density (with its viscousterm) is quite small, i.e. for a temperature T much smallerthan the critical temperature T c of the normal-superfluidtransition. For the unitary Fermi gas T c ≃ . T F , where T F = ǫ F /k B is the Fermi temperature with k B the Boltz-mann constant, ǫ F = ¯ h (3 π n ) / / (2 m ) the bulk Fermienergy of the ideal Fermi gas, n the total density, and m the mass of each atomic fermion. According to Ref. [21]for T < . T F the normal fraction is below 10% and,p-1uca Salasnichin practice, in this range of temperatures our approach isfully justified.At zero temperature the low-energy collective dynam-ics of a Fermi superfluid of neutral and dilute atoms atunitarity can be described by the equations of irrotationaland inviscid hydrodynamics ∂∂t n + ∇ · ( n v ) = 0 , (1) m ∂∂t v + ∇ h m v + µ ( n ) + U ( r ) i = , (2)where n ( r , t ) is the total density of the superfluid, v ( r , t )is its velocity field, and µ ( n ) = ξ ¯ h m (3 π ) / n / (3)is the bulk chemical potential of the system, with ξ ≃ . σ = ↑ , ↓ ). The term U ( r )models the external potential which traps the atoms.Let us consider the unitary Fermi gas with constant den-sity ¯ n with U ( r ) = 0. Experimentally this configurationcan be obtained with a very large square-well potential(or a similar external trapping), such that in the modelone can effectively impose periodic boundary conditionsinstead of the vanishing ones. A density variation alongthe z axis with respect to the uniform configuration ¯ n canbe experimentally created by using a blue-detuned (brightperturbation) or a red-detuned (dark perturbation) laserbeam [22]. In practice, we perform the following factor-ization n ( r ) = n ⊥ ( x, y ) n k ( z ) , (4)by imposing also n ⊥ ( x, y ) = ¯ n ⊥ (5) n k ( z ) = ¯ n k ρ ( z ) (6)such that n ( r ) = ¯ n ρ ( z ) (7)where ¯ n = ¯ n ⊥ ¯ n k , and ρ ( z, t ) is the relative density, i.e. thelocalized axial modification with respect to the uniformdensity ¯ n . We impose periodic boundary conditions alongthe z axis, namely ρ ( z = L z , t ) = ρ ( z = − L z , t ), with 2 L z the axial-domain length. We set v ( r , t ) = (0 , , v ( z, t ))with v ( z, t ) the velocity field such that v ( z = L z , t ) = v ( z = − L z , t ). Moreover we impose that the initiallocalized wave packet satisfies the boundary conditions ρ ( z = ± L z , t = 0) = 1 and v ( z = ± L z , t = 0) = 0.Because the dimensional reduction is done assuming theuniformess in x , y directions, we shall consider the propa-gation of a plane wave along the z axis.Inserting Eq. (7) into Eqs. (1) and (2) one finds the1D hydrodynamic equations for the axial dynamics of the superfluid, given by ˙ ρ + vρ ′ + v ′ ρ = 0 , (8)˙ v + vv ′ + c ls ( ρ ) ρ ρ ′ = 0 , (9)where dots denote time derivatives, primes space deriva-tives, and c ls ( ρ ) = c s ρ / (10)is the local sound velocity, with c s = c ls (1) = p ξ/ v F thebulk sound velocity, v F = p ǫ F /m is bulk Fermi velocityand ǫ F = ¯ h m (3 π ¯ n ) / the bulk Fermi energy.The bulk sound velocity c s is the speed of propagationof a small perturbation ˜ ρ ( z, t ) with respect to the uniformsuperfluid of density ¯ n . In fact, setting ρ ( z, t ) = 1+ ˜ ρ ( z, t ),with ˜ ρ ( z, t ) ≪ v ( z, t ) of the same order of ˜ ρ ( z, t ),from the linearization of Eqs. (8) and (9) we get the fa-miliar linear wave equation (cid:18) ∂ ∂t − c s ∂ ∂z (cid:19) ˜ ρ ( z, t ) = 0 , (11)for ˜ ρ ( z, t ) and a similar equation for v ( z, t ). Modelling theinitial perturbation with a Gaussian shape, i.e.˜ ρ ( z,
0) = 2 η e − z / (2 σ ) , (12)one finds [1, 2] from the linearized equations˜ ρ ( z, t ) = η e − ( z − c s t ) / (2 σ ) + η e − ( z + c s t ) / (2 σ ) , (13)with initial condition ˙˜ ρ ( z, t = 0) = 0. Thus, for the con-servation of the linear momentum, the initial wave packetsplits into two waves travelling in opposite directions withthe speed of sound c s . Obviously Eq. (13) is reliable onlyif | η | ≪
1. As expected, a small (infinitesimal) perturba-tion gives rise to sound waves.We can find wave solutions of Eqs. (8) and (9) witha generic initial density profile by following the approachdescribed by Landau and Lifshitz [1]. By supposing thatthe velocity v depends explicitly on the density ρ , i.e. v = v ( ρ ( z, t )), one has ˙ v = dvdρ ˙ ρ , v ′ = dvdρ ρ ′ . We now imposethat the two hydrodynamic equations reduce to the samehyperbolic equation ˙ ρ + c ( ρ ) ρ ′ = 0 , (14)where c ( ρ ) = v ( ρ ) + dvdρ ρ = v ( ρ ) + c ls ( ρ ) ρ (cid:18) dvdρ (cid:19) − . (15)It is quite easy to verify that, given a initial condition F ( z ) = ρ ( z, t = 0) for the density profile, the time-dependent solution ρ ( z, t ) of the hyperbolic equation (14)satisfies the following implicit, but algebraic, equation: ρ ( z, t ) = F ( z − c ( ρ ( z, t )) t ) . (16)p-2upersonic and subsonic shock waves in the unitary Fermi gas -0.5 -0.25 0 0.25 0.5 0.75 100.511.522.53 M -0.5 -0.25 0 0.25 0.5 0.75 1 η T s c s / σ supersonicsubsonic Fig. 1: Properties of the shock waves. Upper panel: Machnumber M as a function of the amplitude η of the pertur-bation (solid line). For M < − / ≤ η <
0) waves, while for
M > η >
0) waves. The dotted line separates the tworegimes. Lower panel: period T s of formation (breaking time)of the shock-wave front as a function of the amplitude η of theperturbation. T s is in units of σ/c s , where σ is the width ofthe perturbation and c s = p ξ/ v F is the bulk speed of sound,with v F the Fermi velocity. To determine the local velocity of propagation c ( ρ ),which is not equal to the sound velocity c ls ( ρ ), we observethat from Eqs. (15) we get dvdρ = ± c ls ( ρ ) ρ . (17)After separation of variables, and imposing that at infinitythe density is equal to one and the velocity field is zero,we finally get v ( ρ ) = ± c s ( ρ / − . (18)The velocity c ( ρ ) follows directly from the velocity v ( ρ )by using Eqs. (15). One finds c ( ρ ) = v ( ρ ) ± c ls ( ρ ), namely c ( ρ ) = ± c s (cid:16) ρ / − (cid:17) . (19)In conclusion, we have found that the density field ρ ( z, t )satisfies the implicit algebraic equation (16) with c ( ρ )given by Eq. (19). Note that a similar result, but with avery different local velocity, has been obtained by Damski[8] for the 1D BEC.Eqs. (16) and (19) contain the dynamics of the twowaves propagating to the left and to the right with ini-tial condition (12). Some properties characterizing thedynamics can be extracted from these equations. First ofall the two travelling waves have symmetric shapes duringthe time evolution. In addition, both amplitude and ve-locity of the extrema (maxima or minima, depending on the sign of η ) of the two waves are practically constantduring time evolution. In particular, the amplitude of theextrema is given by A ( η ) = 1 + η while the velocity of theextrema reads V ( η ) = c (1 + η ) = ± c s (cid:16) η ) / − (cid:17) . (20)Notice that taking η = 0 the velocity of the impulse ex-trema reduces to the sound velocity: V (0) = c (1) = c s = p ξ/ v F . Moreover, bright perturbations ( η >
0) movefaster than dark ones ( η < M = V ( η ) /V (0) of these perturbations in the unitaryFermi gas is simply M = 4(1 + η ) / − . (21)For M >
1, which means η > ≤ M <
1, which means η < M as function of the amplitude η of the perturbation. Notethat since 2 η is the amplitude of the initial condition, seeEq. (12), the region η ≤ − / η >
0) moving tothe right. The speed of impulse maximum V ( η ) is biggerthan the speed of its tails V (0). As a result the impulseself-steepens in the direction of propagation and a shockwave front takes place. The breaking-time T s requiredfor such a process can be estimated as follows: the shockwave front appears when the distance difference traveledby lower and upper impulse parts is equal to the impulsehalf-width 2 σ , namely [ V ( η ) − V (0)] T s = 2 σ . This, byusing Eqs. (19) and Eq. (10), gives T s = σ c s ((1 + η ) / − . (22)In the case of a dark perturbation ( η <
0) the tails ofthe wave packet move faster than the impulse minimum.The shock appears in the back of the travelling wave, andthe period of shock formation is simply T s = 2 σ/ ( V (0) − V ( η )). In the lower panel of Fig. 1 we plot the period T s as a function of the amplitude η of the perturbation. Thefigure shows that as η goes to zero the period T s goes toinfinity; in fact, in this limit the shock wave reduces to asonic wave (sound wave) which does not produce a shock.After the formation of the shock Eqs. (1) and (2) arenot reliable because their exact solutions given of Eqs.(16) and (19) are no more single-valued. To overcomethis difficulty we include a gradient quantum term in thehydrodynamic equations, which become ∂∂t n + ∇ · ( n v ) = 0 , (23) m ∂∂t v + ∇ h m v − λ ¯ h m ∇ √ n √ n + µ ( n ) + U ( r ) i = , (24)We stress that at zero temperature the simplest regular-ization process of the shock is a purely dispersive quantump-3uca Salasnich -300 -150 0 150 300 ρ ( z ) -300 -150 0 150 3000.811.21.41.6-300 -150 0 150 3000.811.21.41.6 -300 -150 0 150 3000.811.21.41.6 z/l F ω F t=0 ω F t=75 ω F t=150 ω F t=225 ω F t=300 ω F t=375 Fig. 2: Time evolution of supersonic shock waves. Initialcondition with σ/l F = 18 and η = 0 .
3. The curves givethe relative density profile ρ ( z ) at subsequent frames, where l F = p ¯ h / ( mǫ F ) is the Fermi length and ω F = ǫ F / ¯ h is theFermi frequency. gradient term. Historically, the gradient term with λ wasintroduced by von Weizs¨acker to treat surface effects innuclei [23]. This approach has been adopted for quantumhydrodynamics of electrons by March and Tosi [24], andalso by Zaremba and Tso [25]. In the study of the BCS-BEC crossover the gradient term has been considered byvarious authors [26]. The choice of the parameter λ in Eqs.(23) and (24) is still debated, here we choose λ = 1 /
4, avalue which gives good agreement with Monte Carlo cal-culations at zero and finite temperature (for details see[19, 21]). Moreover we set ξ = 0 . v c a normal componentwith a dissipative term appears in the fluid [1]. As dis-cussed in [27], for the unitary Fermi gas one has v c ≃ c s .Nevertheless, at ultrcold temperature the normal compo-nent is negligible and also the shear viscosity [17, 20]. Forthese reasons at zero temperature the shock waves are dis-persive and not dissipative [16].Eqs. (23) and (24) can be formally written (for anyvalue of λ , also λ = 0) in terms of a Galilei-invariantnonlinear Schr¨odinger equation [19]. Setting U ( r ) = 0and using Eqs. (3) and (7) we easily get from Eqs. (23)and (24) a 1D nonlinear Sch¨odinger equation. We solvethis equation by using a Crank-Nicolson finite-differencepredictor-corrector algorithm [28] with the initial condi-tion given by Eq. (12) and v ( z, t = 0) = 0. In fact, asalso shown by Damski [8], we have verified that the ini-tial velocity field v ( ρ ( z, t = 0)) and v ( z, t = 0) = 0 givepractically the same time evolution. In Fig. 2 we plot the time evolution of supersonicshock waves obtained with σ/l F = 18 and η = 0 .
3, with q ¯ h / ( mǫ F ) the Fermi length of the bulk system. The fig-ure displays the density profile ρ ( z ) at subsequent times.Note the splitting on the initial bright wave packet intotwo bright travelling waves moving in opposite directions.As previously discussed, there is a deformation of the twowaves with the formation of a quasi-horizontal shock-wavefront. Eventually, this front spreads into wave ripples.There is no qualitative difference with respect to a Bose-Einstein condensate [8] in the physical manifestation ofsupersonic shock waves in the zero-temperature unitaryFermi gas. Nevertheless, due to the very different equa-tion of state, there are large quantitative differences. Ournumerical analysis confirms that the breaking time T s de-creases by increasing the amplitude η , while the velocity V of the maxima of the travelling waves increases by increas-ing η . There is a good agreement between our analyticalformulas, Eqs. (20) and (22), and simulations: the rela-tive difference is within 5% for the velocity V and within20% for the breaking time T s .In Fig. 3 we plot the time evolution of subsonic shockwaves obtained again with σ/l F = 18 and η = − .
2. Alsoin this case the figure shows the splitting on the initialdark wave packet into two dark travelling waves movingin opposite directions. But here, as expected, the quasi-horizontal shock appears in the back side of the travellingwaves. Our simulations show that for η < T s are always dark,i.e. they never exceed the bulk density (compare wave rip-ples of Fig. 2 with those of Fig. 3). Note that dark shockwaves have been studied long ago [29] in a different phys-ical context: the discrete nonlinear Schr¨odinger equation.In that case the wave ripples can exceed the bulk density,probably due to the discrete nature of the Schr¨odingerequation. Also for dark shock waves our analytical pre-dictions on velocity V of the minima and breaking time T s are quite accurate with respect to numerical findings(similar relative differences of runs with η > Li atomic clouds. The shape of these waves changes dur-ing the time evolution giving rise to a shock-wave frontat a characteristic breaking time. We have determinedthe Mach number of these travelling waves as a functionof the perturbation amplitude, showing that supersonicbright and subsonic dark waves behave quite differently.The author thanks Sadahn Kumar Adhikari, TilmanEnss, Boris Malomed, Enzo Orlandini, Mario Salerno, andFlavio Toigo for useful suggestions and discussions.p-4upersonic and subsonic shock waves in the unitary Fermi gas -300 -150 0 150 300 ρ ( z ) -300 -150 0 150 3000.60.811.21.4-300 -150 0 150 3000.60.811.21.4 -300 -150 0 150 3000.60.811.21.4 z/l F ω F t=0 ω F t=75 ω F t=150 ω F t=225 ω F t=300 ω F t=375 Fig. 3: Time evolution of subsonic shock waves. Initial con-dition with σ/l F = 18 and η = − .
2. The curves givethe relative density profile ρ ( z ) at subsequent frames, where l F = p ¯ h / ( mǫ F ) is the Fermi length and ω F = ǫ F / ¯ h is theFermi frequency.REFERENCES[1] L.D. Landau and E.M. Lifshitz, Fluid Mechanics (Perga-mon Press, London, 1987).[2] G.G. Whitham,
Linear and Nonlinear Waves (Wiley, NewYork, 1974).[3] Z. Dutton, M. Budde, C. Slowe, and L.V. Hau, Science , 663 (2001).[4] M.A. Hoefer, M.J. Ablowitz, I. Coddington, E.A. Cornell,P. Engels, and V. Schweikhard, Phys. Rev. A , 023623(2006).[5] J.J. Chang, P. Engels, and M.A. Hoefer, Phys. Rev. Lett. , 1311 (2009).[6] R. Meppelink, S.B. Koller, J.M. Vogels, P. vander Straten, E.D. van Ooijen, N.R. Heckenberg, H.Rubinsztein-Dunlop, S.A. Haine, and M.J. Davis, Phys.Rev. A , 043606 (2009).[7] I. Kulikov and M. Zak, Phys. Rev. A , 063605 (2003).[8] B. Damski, Phys. Rev. A , 043610 (2004).[9] A.M. Kamchatnov, A. Gammal, and R.A. Kraenkel,Phys. Rev. A , 063605 (2004).[10] V.M. Perez-Garcia, V.V. Konotop, and V.A. Brazhnyi,Phys. Rev. Lett. , 220403 (2004).[11] B. Damski, Phys. Rev. A , 043601 (2006).[12] A. Ruschhaupt, A. del Campo, and J.G. Muga, Eur. Phys.J. D , 399 (2006).[13] L. Salasnich, N. Manini, F. Bonelli, M. Korbman, and A.Parola, Phys. Rev. A , 043616 (2007).[14] M.A. Hoefer, M.J. Ablowitz, and P. Engels, Phys. Rev.Lett. , 084504 (2008).[15] J. Joseph, J. Thomas, M. Kulkarni, and A. Abanov, Phys.Rev. Lett. , 150401 (2011).[16] A. Bulgac, Y.-L. Luo, and K.J. Roche, e-preprintarXiv:1108.1779.[17] C Cao, E Elliott, H Wu and J E Thomas, New J. Phys. , 075007 (2011).[18] S. Giorgini, L.P. Pitaevskii, and S. Stringari, Rev. Mov. Phys. , 1215 (2008).[19] L. Salasnich and F. Toigo, Phys. Rev. A , 642 (2009); L. Salasnich andF. Toigo, Phys. Rev. A , 059902(E) (2010).[20] H. Guo, D. Wulin, C. C. Chien, and K. Levin,Phys.Rev.Lett. , 020403 (2011).[21] L. Salasnich, Phys. Rev. A , 063619 (2010).[22] H.J. Metcalf and P. van der Straten, Laser Cooling andTrapping (Springer, New York, 1999); C.J. Foot,
AtomicPhysics (Oxford University Press, Oxford, 2005).[23] C.F. von Weizs¨acker, Zeit. Phys. , 431 (1935).[24] N.H. March and M. P. Tosi, Proc. R. Soc. A , 373(1972).[25] E. Zaremba and H.C. Tso, Phys. Rrev. B , 8147 (1994).[26] Y.E. Kim and A.L. Zubarev, Phys. Rev. A , 033612(2004); M.A. Escobedo, M. Mannarelli and C. Manuel,Phys. Rev. A , 063623 (2009); G. Rupak and T. Sch¨afer,Nucl. Phys. A , 52 (2009); S.K. Adhikari, Laser Phys.Lett. , 901 (2009); A. Csordas, O. Almasy, and P. Szep-falusy, Phys. Rev. A , 063609 (2010).[27] R. Combescot, M. Yu. Kagan, and S. Stringari, Phys.Rev. A , 042717 (2006); F. Ancilotto, L. Salasnich, F.Toigo, Phys. Rev. A , 033627 (2009).[28] E. Cerboneschi, R. Mannella, E. Arimondo, and L. Salas-nich, Phys. Lett. A , 495 (1998); G. Mazzarella andL. Salasnich, Phys. Lett. A , 3616(1997)., 3616(1997).