Current-Phase Relation of a Bose-Einstein Condensate Flowing Through a Weak Link
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] D ec Current-Phase Relation of a Bose-Einstein CondensateFlowing Through a Weak Link
F. Piazza, L. A. Collins, and A. Smerzi CNR-INFM BEC center and Dipartimento di Fisica, Universit`a di Trento, I-38050 Povo, Italy Theoretical Division, Mail Stop B214, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 (Dated: March 2, 2018)We study the current-phase relation of a Bose-Einstein condensate flowing through a repulsivesquare barrier by solving analytically the one dimensional Gross-Pitaevskii equation. The barrierheight and width fix the current-phase relation j ( δφ ), which tends to j ∼ cos( δφ/
2) for weak barriersand to the Josephson sinusoidal relation j ∼ sin( δφ ) for strong barriers. Between these two limits,the current-phase relation depends on the barrier width. In particular, for wide enough barriers,we observe two families of multivalued current-phase relations. Diagrams belonging to the firstfamily, already known in the literature, can have two different positive values of the current atthe same phase difference. The second family, new to our knowledge, can instead allow for threedifferent positive currents still corresponding to the same phase difference. Finally, we show that themultivalued behavior arises from the competition between hydrodynamic and nonlinear-dispersivecomponents of the flow, the latter due to the presence of a soliton inside the barrier region. PACS numbers: 03.75.Lm, 74.50.+r, 47.37.+q
I. INTRODUCTION
The current-phase relation characterizes the flow of a su-perfluid/superconductor through a weak link [1, 2, 3].The latter is a constricted flow configuration that can berealized in different ways: i) apertures in impenetrablewalls mostly for helium, ii) sandwich or bridge structuresfor superconductors, and iii) penetrable barriers createdby laser beams for ultracold dilute gases. Much infor-mation about such systems can be extracted from thecurrent-phase relation, which, given a fluid, depends onlyon the link properties. For instance, with He, the tran-sition from the usual AC Josephson effect to a quantizedphase slippage regime [4] corresponds to the switchingfrom a sine-like current phase relation to a multivaluedone [5].A weak link configuration can be modelled very generallyupon taking a portion of a superfluid/superconductor tohave “different conduction properties” with respect tothe rest of the system. Two pieces of superconductorjoined by a third superconducting region with a smallercoherence length provide one example, whose current-phase relation in one dimension has been studied withthe Ginzburg-Landau equation [6].In the context of ultracold dilute gases, raising a repul-sive penetrable barrier across the flow yields an equiva-lent configuration. For instance, with Bose-Einstein con-densates (BEC), Josephson effect(s) have been theoreti-cally studied [7, 8, 9, 10, 11] and experimentally demon-strated using multiple well traps [12, 13, 14]. Theoret-ically, the current-phase relation has been studied for aflow through a repulsive square well with fermions acrossthe BCS-BEC crossover by means of one dimensionalBogoliubov-de Gennes equations [15], for weak barrierswith bosons in a local density approximation [16], andfor fermions on the BEC side of the crossover using anonlinear Schr¨odinger equation approach [17].In this manuscript, we study the current-phase relationfor a BEC flowing through a repulsive square well. The weak link configuration, and in turn the current-phaserelation, is then determined by the barrier height withrespect to the chemical potential and by the barrierwidth with respect to the healing length. Though wesolve a one-dimensional Gross-Pitaevskii equation, theresults presented in this manuscript are not just rele-vant for BECs, but also include the essential features ofcurrent-phase relations of superconducting or superfluidHe-based weak links when governed by the Ginzburg-Landau equation. For any barrier width, we find that inthe limit of zero barrier height, the current phase relationtends to j ( δφ ) = c ∞ cos( δφ/ c ∞ being the bulksound velocity, which corresponds to the phase across agrey soliton at rest with respect to the barrier. On theother hand, if the barrier height is above the bulk chem-ical potential at zero current, the limit of tunneling flowis reached either when the barrier height is much biggerthan the bulk chemical potential at zero current or whenthe barrier width is much larger than the bulk healinglength. In this regime, we recover the the usual Joseph-son sinusoidal current-phase relation and obtain an ana-lytical expression for the Josephson critical current as afunction of the weak link parameters. For barriers widerthan the healing lenght inside the barrier region, we ob-serve two families of multivalued (often called reentrant)current-phase relations. The first, already studied sincethe early works on superconductivity [2, 6], shows a posi-tive slope of the current when the phase difference is closeto π , thereby reaching a phase difference larger than π at least for small currents. The second family, appear-ing at a smaller barrier height, has instead a negativeslope of the current close to π , and in some cases canremain within the 0 − π interval across the whole rangeof currents. These two families can also be distinguishedby the maximum number of different positive currentscorresponding to the same phase difference: two for thefirst family, three for the second one. As the first kindof reentrant behavior was proven to be connected to theonset of phase-slippage in the AC Josephson effect [5],the second might then be connected to the appearanceof new features in the Josephson dynamics. We finallyobserve that the hysteresis characterizing both familiesof reentrant current-phase relations is always due to thecompetition between a hydrodynamic component of theflow and a nonlinear-dispersive component, the latter dueto the presence of a soliton inside the barrier region. Thetwo components can coexist only for barriers wide enoughto accomodate a soliton inside. In this spirit, we developa simple analytical model which describes very well reen-trant regimes of current-phase relations. II. THEORETICAL DEVELOPMENTSA. Stationary Solutions
We consider a dilute repulsive Bose-Einstein condensateat zero temperature flowing through a 1D rectangularpotential barrier. We look for stationary solutions of the1D GPE [18]: − ~ m ∂ xx Ψ + V ext ( x )Ψ + g | Ψ | Ψ = µ Ψ , (1)where Ψ( x ) = p n ( x ) exp[ iφ ( x )] is the complex orderparameter of the condensate, µ is the chemical poten-tial, and g = 4 π ~ a s /m with m the atom mass and a s > s -wave scattering length. The order param-eter phase φ ( x ) is related to the superfluid velocity via v ( x ) = ( ~ /m ) ∂ x φ ( x ). The piecewise constant externalpotential describes the rectangular barrier of width 2 d and height V : (cid:26) V ext ( x ) = V , | x | ≤ d ,V ext ( x ) = 0 , | x | > d . (2)We consider solutions of Eq. (1) which are symmetricwith respect to the point x = 0, therefore discarding casesin which a reflected wave is present [19]. Such symmetricsolutions in the presence of a barrier exist due to thenonlinearity in the GPE. We also restrict our analysis tosubsonic flows v ∞ ≤ c ∞ , with c ∞ = p gn ∞ /m being thesound velocity for a uniform condensate of density n ∞ .As boundary conditions, we fix the condensate density n ∞ and velocity v ∞ at x = ±∞ , thereby determining thechemical potential µ = gn ∞ + mv ∞ . Using the relationΨ = √ n exp[ iφ ], Eq. (1) can be split into a continuityequation, stating the constancy in space of the current j = n ( x ) v ( x ) = n ∞ v ∞ , and an equation for the densityonly µ = − ~ m ∂ xx √ n √ n + mj n + V ext ( x ) + gn , (3)where we have used the continuity equation v ( x ) = j/n ( x ). Its solution n ( x ) is expressed in terms of Jacobielliptic functions [6, 21]. Due to symmetry, we need onlyconsider half of the space x >
0. The solution outsidethe barrier x > d becomes n out ( x ) = n ∞ − A ∞ + A ∞ tanh (cid:20)r mg ~ A ∞ ( x − d ) + x (cid:21) , (4) where x = arctanh p ( n d − mv ∞ /g ) /A ∞ , A ∞ = n ∞ − mv ∞ /g ≥
0, and n d is the density at the barrier edge x = d . The solution inside the barrier x < d is: n ( x ) = n + A sn [ b x, k ]cn [ b x, k ] , ∆ ≥ A ≥ ,n ( x ) = n + A − cn[ b x, k ]1 + cn[ b x, k ] , else , (5)where n is the density at x = 0, ∆ = ( n − µ/g ) − mj /gn with ˜ µ = µ − V , A = 3 n / − ˜ µ/g − √ ∆ / A = p n − n ˜ µ/g + mj / gn ). Finally, the pa-rameters entering the Jacobi sines sn and cosines cn are b = q mg ( √ ∆ + A ) / ~ , k = ( √ ∆ / ( √ ∆+ A )) / , and b = p mgA / ~ , k = (( A − A − √ ∆ / / A ) / .Given n ∞ and v ∞ , we are left with two free parame-ters: n and n d , which are next determined by matchingthe density and its derivative at the barrier edge x = d .First, the derivative matching condition, enforced usingthe first integral of Eq. (3), allows us to write n d as afunction of n for any value of ∆: n d = g V (cid:20) n ∞ + n + mv ∞ n ∞ g (2 − n ∞ n ) − µg n (cid:21) . (6)The density matching equation n d = n i ( d ) ( i = 1 ,
2) isthen solved by a numerical root finding method, yelding n .Two bounded solutions are always found [22], an exampleis given in the upper panel of Fig. 1. In the following,the solution which tends to a plane wave for V → V and d , and at afixed density at infinity n ∞ , the solutions exist up to acritical injected velocity v ∞ = v c < c ∞ , at which theymerge and disappear. This behavior was found in thecase of a repulsive delta potential in [23]. Similarly, in[19], the same kind of merging was reported for a 1D BECflow through a repulsive square well when the width ofthe latter increases. B. Current-phase Relation
As pointed out in the Introduction, the current-phaserelation j ( δφ ) for a given superfluid only depends on theproperties of the weak link, in our case the barrier height V with respect to the chemical potential and the width2 d with respect to the healing length. For a fixed current j = v ( x ) n ( x ), the phase difference across the system iscalculated using the relation φ ( x ) = R x dy ( m/ ~ ) j/n ( y )and then renormalized by the phase accumulated by the aplane wave with the same boundary conditions in absenceof barrier (see lower panel in Fig. 1). Two different valuesof δφ are found, corresponding to the upper and lowersolutions.In this section, we will use dimensionless quantities,employing the chemical potential at zero current gn ∞ as the unit of energy, the bulk healing length ξ ∞ =FIG. 1: (Color online) Typical solutions for a barrierwith width 2 d = 4 ξ ∞ . Density (upper panel) and phase(lower panel) as a function of position are shown for boththe upper and the lower solution (see inset). Dashed linesin the lower panel correspond to the phase accumulatedby a plane wave in absence of barrier. ~ / √ mgn ∞ as the unit of length, and ~ /gn ∞ as theunit of time. Exploiting the symmetry of the systemabout x = 0, the phase difference can be written as δφ i =lim x →∞ [ R d jdy/n i ( y ) + R xd jdy/n out ( y ) − jx/n ∞ ] , i = 1 , δφ i = Z d jdyn i ( y ) − jd +2 " arcsin( j p n i ( d ) ) − arcsin( j √ . (7)The first two terms in Eq. (7) correspond to the phaseacquired inside the barrier while the third, which we cancall the pre-bulk term, gives the phase accumulated out-side the barrier, where the density has not yet reachedits bulk value n ∞ . We have taken the latter to be one.In Fig. 2, the current phase relation is shown for differentbarrier widths and heights. Each curve has a maximumat the point ( δφ c , j c ), with j c = n ∞ v c being the criticalcurrent at which the two stationary solutions merge anddisappear. The upper solutions constitute the part ofthe current-phase diagram which connects the maximumwith the point ( δφ = 0 , j = 0), while the lower onesbelong to the branch connecting the maximum to thepoint ( δφ = π, j = 0). Indeed, we will now show that,for any d and V →
0, the upper branch tends to a planewave, while the lower branch tends to a grey soliton. Inorder to have a finite n d in this limit, the term in squarebrackets in Eq. (6) must tend to zero, yielding a cubicequation with two coincident solutions n = 1 and a third n = j /
2, where we have set n ∞ = 1 for simplicity. For n = 1 we have ∆ = (1 − j / ≥
0, correspondingto the plane wave solution n ( x ) = n = 1. For n = j / A = j / − < n ( x ) = n +(1 − n ) tanh ( p / − n / x ). Therefore, in this limit FIG. 2: (Color online) Current-phase relation for a bar-rier with half width d = 0 . ξ ∞ (left panel) and d = 3 ξ ∞ (right panel), for different barrier heights V (see inset).The arrows sketch the behavior of the maximum of thecurves upon increasing V . δφ = 0 for any j , meaning that the upper branch isactually a vertical line, while for the lower branch wehave cos( δφ j √ . (8)This curve has a maximum at δφ = 0, corresponding to j = √
2, that is, the sound velocity c ∞ in dimensionlessunits.The arrows in Fig. 2 sketch the behavior of the maximumof the current-phase relation ( δφ c , j c ) as the height V isincreased at a fixed barrier width 2 d . For any width, thecurrent-phase diagram initially takes a cosine-like shape(Eq. (8)) when V ∼ δφ ) shape for sufficiently large V , characterizingthe Josephson regime of tunneling flow (blue trianglesin Fig. 2). Between these two limits, the behavior ofthe maximum is determined by the barrier width. Forthin barriers ( d . ξ ∞ ), as shown in the left panel ofFig. 2, the point ( δφ c , j c ) moves down-right, reachingthe Josephson regime keeping δφ c always smaller than π/
2. On the other hand, for sufficiently wider barriers, δφ c is able to reach values larger than π/ V during the down-right displacement ofthe maximum. The latter then moves down-left to finallyenter the Josephson regime, as shown in the right panelof Fig. 2. We note that in this way, while V is increased,the phase δφ c takes the value π/ V is much smaller than the chemicalpotential, indeed the current-phase relation is symmetri-cal with respect to π/
2, but not sinusoidal.
C. Josephson Regime
As described in the previous section, for strong enoughbarriers the flow enters the tunneling regime, and thecurrent-phase relation takes a sinusoidal form. In thefollowing, we will describe analytically this behavior, de-riving a relation between the Josephson critical currentand the barrier parameters V and d . Since we are nowinterested in tunneling flows, we will take V > gn ∞ andwill show that the Josephson regime is attained by eitherincreasing the barrier height V or its width 2 d .Since in this regime the injected velocity of a stationaryflow v ∞ must be much smaller than the sound veloc-ity c ∞ , the chemical potential can be written as ˜ µ ≃ gn ∞ − V <
0, upon neglecting the kinetic energy term mv ∞ /
2. Moreover, the density inside the barrier n beingexponentially small, we can write ∆ ≃ (2˜ µ/g ) − s > s = 4( n ˜ µ + mj /n ) /g , where we have neglected n . Thus, both the upper and lower solutions are of thekind n ( x ), with A ≃ n − mj / n ˜ µ , b ≃ p m | ˜ µ | / ~ ,and k ≃
1. The density in the Josephson regime has thusthe form n jos ( x ) = n + ( n + mj n | ˜ µ | ) sinh r m ~ | ˜ µ | x ! . (9)In order to write the density matching equation n d = n ( d ), we approximate n d by discarding both n and2 mv ∞ n ∞ /g in Eq. (6), obtaining a quadratic equationfor n . Further assuming that sinh q | ˜ µ | ǫ d ≫
1, with ǫ d = ~ / md being the kinetic energy associated withthe barrier length scale d , the solutions of the above equa-tion are of the form n + / − = ¯ n (1 ± √ − q ) with ¯ n =( gn ∞ / V ) n ∞ / sinh p | ˜ µ | /ǫ d , and q = mj / n | ˜ µ | .The Josephson critical current corresponds to the merg-ing of the two solutions at q = 1, and reads j jos = n ∞ r | ˜ µ | m gn ∞ V e − √ | ˜ µ | /ǫ d , (10)where we have used sinh( x ) ≃ exp( x ) /
2, for x ≫
1. Thecritical velocity for a bosonic and fermionic superfluidflowing through a repulsive square well has been calcu-lated in [16] within the local density approximation (forBEC case see also [24]). Analytical expressions for thecritical current of a BEC flow were found for both slowlyvarying and weak barriers [23]. Eq. (10) thus enrichesthe above set of analytical results by providing the criti-cal current for strong barriers.Finally, we calculate the current-phase relation usingEq. (7). Since j jos is exponentially small, only the firstterm in Eq. (7) contributes, and the integral can be per-formed analytically to obtain: δφ + = arcsin (cid:18) jj jos (cid:19) , δφ − = π − arcsin (cid:18) jj jos (cid:19) . (11)Thus, we recover the sinusoidal current-phase relationcharacterizing a Josephson regime of tunneling flow. FIG. 3: (Color online) Current-phase relation in thereentrant regime d = 20 ξ ∞ . Insets show in more detailthe shape of the diagram close to the critical point. Thearrows sketch the behavior of the maximum of the curvesupon increasing V .FIG. 4: (Color online) Typical density profiles in thereentrant regime d = 20 ξ ∞ for the upper and lower solu-tions. D. Reentrant Regimes
When the barrier width 2 d greatly exceeds the healinglength, the current-phase relation becomes multivaluedas shown in Fig. 3 for d = 20 ξ ∞ . These so-called reen-trant current-phase diagrams were first predicted for longsuperconducting weak links [2, 6]. Remarkably, they havebeen experimentally demonstrated with superfluid He[5].In our system, two kinds of multivalued diagrams arefound for a fixed barrier width 2 d . We designate the firstkind (green dots), appearing at larger barrier heights V ,as reentrant type I, and the second kind (orange dia-monds in Fig. 3), occurring for smaller V , as reentranttype II. They differ in the behavior of the lower branchat small currents j ≪
1. The phase decreases with in-creasing j in type II diagrams while it increases in typeI, reaching values larger than π for j ≪
1. These twofamilies of diagrams also differ in the number of positivevalues of the current j corresponding to the same phasedifference δφ . Indeed, diagrams of type I can have twovalues of j at the same δφ , while diagrams of type IIcan allow for three. The existence of these type II dia-grams, to our knowledge, has not been discussed in theliterature. In [5], only type I current-phase relations wereobserved and connected to the onset of phase-slip dissi-pation in the system. In this spirit, a type II diagrammight for instance lead to an instability of a differentkind.In the remainder of this section, we will develop an an-alytical model that captures the essential features un-derlying both kinds of reentrant current-phase diagrams.Examination of typical density profiles belonging to thereentrant regime (see Fig. 4) suggests to construct thelower solution by starting from the upper one at thesame current, then adding a grey soliton inside the bar-rier region. Since we are dealing with wide barriers, wedescribe the upper solution in the local density approxi-mation (LDA) (its current-phase relation in this approx-imation is also discussed in [16]). At a fixed current j ,the density inside the barrier | x | < d is constant andequal to n = ˜ µ/ g + 2˜ µ cos( ω/ / g , when j < j th , or n = ˜ µ/ g + 2˜ µ cos( π/ − ω/ / g , when j > j th with j th defined by j = 4˜ µ / mg , and ω = arccos | − mg j / µ | . The phase difference calculated withinLDA misses the pre-bulk term in Eq. (7), thus, for theupper branch, it is simply δφ = 2 mjd (1 /n − /n ∞ ) / ~ .Now, for the density profile of the lower branch n ( x ),we take a grey soliton (Eq. (4)) placed inside the barrierat x = 0, with a bulk density given by n and a bulkvelocity v = j/n while in the region | x | > d we keepthe density profile of the upper branch, that is, a con-stant density n ∞ and velocity j/n ∞ . Notice that in thissection n stands for the density of the upper solution at x = 0, as indicated in Fig. 4. The density at x = 0 forthe lower solution corresponds to the center of the dip inthe grey soliton density profile. Finally, using Eq. (7) weobtain the phase difference for the lower branch δφ = δφ + 2 arccos s mv gn ( d ) . (12)At a given current j , the overall phase difference cor-responding to the lower solution has two contribu-tions: 1) the “hydrodynamics phase” δφ coming fromLDA and 2) the “nonlinear-dispersive phase” δφ sol =2 arccos p mv /gn ( d ) accumulated across the grey soli-ton. While δφ is a monotonically increasing functionof j , δφ sol is instead monotonically decreasing, startingfrom π at zero current [25].Therefore, the hysteresis characterizing a reentrantcurrent-phase relation is due to the competition betweenthe hydrodynamic and the nonlinear-dispersive compo-nents of the flow, which can coexist only for barrierswide enough to accomodate a soliton inside. In par-ticular, we can derive a condition for the appearanceof type I reentrant behavior upon expanding Eq. (12)for small currents, and requiring δφ > π . Usingarccos( x ) ≃ π/ − x , for x ≪
1, we get δφ ≃ π + 2 jν ,where ν = md (1 /n − /n ∞ ) / ~ − p m/gn , and we FIG. 5: (Color online) Comparison between theLDA+soliton model (left panel) or the Deaver-Piercemodel (right panel) and the exact results.have taken n ( d ) ≃ n . The condition for type I reen-trance to appear is thus ν >
0. For j ≪
1, we have n ≃ ˜ µ/g ≃ n ∞ − V /g , and since within the LDA ap-proximation V ≪ gn ∞ , the condition ν > V gn ∞ dξ √ > , (13)with ξ = ~ / √ mgn being the healing length inside thebarrier region where the density is n . Equation (13),holding for V ≪ gn ∞ , has a clear physical meaning:in order to have a type I reentrant current phase dia-gram, the barrier width 2 d must be sufficiently largerthan 2 √ ξ , which is the characteristic size of a solitonplaced inside the barrier.In the left panel of Fig. 5, we compare the current-phaserelation calculated with the above model (solid lines) tothe exact results. Within the reentrant regime, for bothtype I and type II, the agreement is striking with onlyslight differences close to the the maximum ( δφ c , j c ). Onthe other hand, for thin/strong barriers, LDA, and inturn the above model, is in clear disagreement with theexact results. (See cases d = 3 ξ ∞ in the left panel ofFig. 5).The difference between type I and type II diagrams has aphysical interpretation within the above model, namelythat the hydrodynamic component of the flow dominatesfor all currents in type I reentrance (excluding the re-gion j ≃ j c [26]), while it is overcome by the nonlinear-dispersive part for sufficiently small current in type II.In the literature [2], multivalued current-phase relationsare typically modelled by describing the weak-link withan equivalent circuit containing a linear inductance in se-ries with a sinusoidal inductance, the latter correspond-ing to an ideal Josephson junction. When compared toour GPE exact results (see right panel in Fig. 5), thismodel in general fails to describe the curvature of bothbranches, and misses type II reentrance, as well as nearlyfree regimes (e.g. red squares in Fig. 2,3) since it does notallow the phase δφ c , corresponding to the maximum ofthe diagram, to be smaller than π/
2. It proves quite ac-curate only for sufficiently large barrier heights V , veryclose to the Josephson regime. In the He experimentof [5], this so called Deaver-Pierce model [27] agrees wellwith the measured current-phase relation. This mightbe due to the fact that these experiments are performedwith a fixed weak link configuration, moving the systemacross the transition between Josephson and type I reen-trant regime, upon changing the He healing length withtemperature, but always staying sufficiently close to theJosephson regime.
III. CONCLUSIONS
We studied the current-phase relation of a BEC flow-ing through a weak link created by a repulsive barrier.The link is thus modelled by two parameters, the barrierheight and width, which fix the current-phase relation.Even though we solved a simplified model, we believethat the results obtained will also be relevant for super-conducting and superfluid He-based weak links.We obtained analytical results for the weak barrier limit,for which the current-phase relation has a ∼ cos( δφ/ ∼ sin( δφ ) form characterizing the Josephson regime.In particular, we derived an expression for the Josephsoncritical current as a function of the link parameters.Finally, we found two kinds of multivalued current-phasediagrams which we show, by means of an analyticalmodel, to appear due to the competition between a hy-drodynamic and a nonlinear-dispersive part of the flow,which can coexist only for barriers wide enough to ac-comodate a soliton inside. The first kind, showing twodifferent positive currents at same phase difference, hasrecently been experimentally demonstrated with He andproven to be connected to the appearance of phase slip-page in the AC Josephson dynamics. The second one,new to our knowledge, can instead allow for three dif-ferent positive currents corresponding to the same phasedifference. We believe that this new kind of hysteresisin the current-phase relation can be associated with newfeatures emerging in the Josephson dynamics, which willbe studied elsewhere.
Acknowledgments
The Los Alamos National Laboratory is operated by LosAlamos National Security, LLC for the National NuclearSecurity Administration of the U.S. Department of En-ergy under Contract No. DE-AC52-06NA25396. [1] A. Barone, G. Paterno, P hysics and Applications of theJosephson Effect ( John Wiley and Sons, New York,1982).[2] K. K. Likharev, D ynamics of Josephson Junctions andCircuits (Gordon and Breach, New York, 1986).[3] R. E. Packard, Rev. Mod. Phys. , 641 (1998).[4] O. Avenel and E. Varoquaux, Phys. Rev. Lett. , 2704(1985).[5] E. Hoskinson, Y. Sato, I. Hahn, and R. E. Packard, Nat.Phys. , 23 (2006).[6] A. Baratoff, J. A. Blackburn, and B. B. Schwartz, Phys.Rev. Lett. , 1096 (1970).[7] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy,Phys. Rev. Lett. , 4950 (1997).[8] G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls,Phys. Rev. A , 4318 (1997).[9] I. Zapata, F. Sols, and A. J. Leggett, Phys. Rev. A ,R28 (1998).[10] D. Ananikian and T. Bergeman, Phys. Rev. A , 013604(2006).[11] M. L. Chiofalo, and M. P. Tosi, Europhys. Lett. , 326(2001).[12] F. S. Cataliotti, et al. , Science , 843 (2001).[13] Albiez M., et al. , Phys. Rev. Lett. , 010402 (2005).[14] S. Levy, et al. , Nature , 579 (2007).[15] A. Spuntarelli, P. Pieri, and G. C. Strinati, Phys. Rev.Lett. , 040401 (2007).[16] G. Watanabe, et al. , Phys. Rev. A , 053602 (2009).[17] F. Ancilotto, L. Salasnich, and F. Toigo, Phys. Rev. A , 033627 (2009).[18] F. Dalfovo, S. Giorgini, L. Pitaevskii, and S. Stringari, Rev. Mod. Phys. , 463 (1999).[19] P. Leboeuf, and N. Pavloff, Phys. Rev. A , 033602(2001).[20] The study of the stability of the two solutions is beyondthe scope of the present paper. This issue has been ad-dressed in [23].[21] Yu. G. Mamaladze and O. D. Che˘ishvili, Zh. Eksp. Teor.Fiz. , 169 (1966) [Sov. Phys. JETP , 112 (1966)];J. S. Langer, and V. Ambegaokar, Phys. Rev. , 498(1967); B. T. Seaman, L. D. Carr, and M. J. Holland,Phys. Rev. A , 033609 (2005); W. D. Li, Phys. Rev. A , 063612 (2006).[22] In general, the two solutions can both correspond to apositive ∆, both to a negative ∆, or one to a positiveand the other to a negative ∆.[23] V. Hakim, Phys. Rev. E , 2835 (1997).[24] A. M. Leszczyszyn, G. A. El, Yu. G. Gladush, and A. M.Kamchatnov, Phys. Rev. A , 063608 (2009).[25] By further approximating n ( d ) ≃ n we get δφ sol =2 arccos( j/ p gn /m ). Since n is a monotonically de-creasing function of j , the argument of arccos increasesfor increasing current. Therefore, starting from π at j = 0, δφ sol decreases monotonically with j .[26] Close to the critical point ( δφ c , j c ), and for every current-phase relation, the lower branch always has a phase de-creasing with increasing current, up to the critical pointitself, at which it meets the upper branch. In particular,for reentrant diagrams, this means that the dispersivepart of the flow always dominates over the hydrodynamicpart sufficiently close to the critical point.[27] B. S. Deaver, and J.M. Pierce, Phys. Lett. A , 81, 81