Curvature Induced Topological Defects of p-wave Superfluid on a Sphere
CCurvature Induced Topological Defects of p -wave Superfluid on a Sphere Ruihua Fan,
1, 2
Pengfei Zhang, and Zhe-Yu Shi
1, 3, ∗ Institute for Advanced Study, Tsinghua University, Beijing, 100084, China Department of Physics, Peking University, Beijing, 100871, China School of Physics and Astronomy, Monash University, Victoria 3800, Australia (Dated: February 26, 2018)We study the ground state of spinless fermions living on a sphere across p -wave Feschbach reso-nances. By construsting a microscopic model of fermions on a general curved surface, we show thatthe Guassian curvature induces an emergent magnetic field coupled to the p ± ip order parameters.In the case of a sphere, the magnetic field corresponds to a Dirac monopole field, which causestopological defects in the superfluid ground state. Using the BCS mean field theory, we calculateits many-body ground state self consistently and give the phase diagram. The ground state mayexhibit two types of topological defects, two voritces on the south and north pole or a domain wallwhich separates p θ + ip φ and p θ − ip φ superfluids. Introduction.
Recent development of the p -waveFeschbach resonance technique in ultracold fermonicatoms[1–5] has opened up possibilities to many novelquantum states. Among these states, the chiral p -wavesuperfluid[6–10] has attracted significant amount of inter-est, not only because it is a nontrivial topological phaseof matter, but also because the system may provide aplatform for fault-tolerant quantum computation[11–13].On the other hand, more and more newly developedtrapping techniques allow the experimentalists to con-fine cold atomic gases on various extraordinary geome-tries such as helices[14, 15] and spheres[16, 17]. Thisprovides a unique opportunity to explore quantum many-body systems in non-Euclidean spaces and study the ef-fects of the spatial curvatures on these systems.In this manuscript, we consider a system of spinlessfermions living on a sphere interacting through a p -wavecontact potential. Previous studies of such system on aflat plane shows that the ground state is a p x ± ip y pairedsuperfluid with uniform order parameter in space[6–10].If the system lives on a sphere, it has been known thatthe system must support topological defects[18, 19]. Inref.[19], the authors argued based on the nature of thechiral superfluid and used the Poincar´e-Hopf index the-orem to reach this conclusion. In this work, we demon-strate it from a few-body point of view. We use a two-body calculation to show that the two-body bound statewith nonzero angular momentum would experience anemergent monopole field when moving on the sphere[20].And the magnet flux of the monopole necessarily causetopological defects such as vortices and domain walls onthe sphere.Our main results are summarized in the ground statephase diagram fig.(1). The ground state may presentthree different phases depending on the radius of thesphere R and the interaction strength indicated by scat-tering area s . (a) phase corresponds to the configurationof a pair of vortices located on the south and north poleof the sphere separately. While both (b) and (c) phasecorrespond to a solution with a domain wall between the FIG. 1: A schematic plot of the phase diagram. Warm colorsdenotes p + ip order parameter. Cold colors denotes p − ip order parameter. R is the radius of the sphere, s is the lowenergy p -wave scattering area. Three phases correspondingto different configurations of topological defects appear, (a):a vortex pair on the south and north poles. (b): a domainwall on the southern or northern hemisphere. (c): a domainwall on the equator. p θ + ip φ and p θ − ip φ superfluids. The difference is thatin (b) phase the domain wall is close to one pole and lo-cates on the southern or northern hemisphere, and in (a)phase it is always on the equator. Model.
While most of the previous studies treatedthis subject phenomenologically[21, 22], in this work, wetarget this problem from a microscopic model. We be-gin with a Hamiltonian describing spinless fermions with p -wave interaction on a general two-dimensional surface M : a r X i v : . [ c ond - m a t . qu a n t - g a s ] D ec H = (cid:90) M d s H , (1) H = ψ † (cid:18) − D ν D ν m f − µ (cid:19) ψ + d † α (cid:18) − D ν D ν m f + (cid:15) − µ (cid:19) d α − iλ ( d ν ψ † D ν ψ † + d † ν ψD ν ψ ) . (2)Here d s = √ gdx dx is the area element of the surfacewith g = det g µν . ψ ( d ν ) is the annihilation operatorfor the open channel fermions(closed channel molecules). D ν ( D ν ) stands for the covariant(contravariant) deriva-tive with respect to the ν th coordinate. We use theEinstein summation convention and all the Greek lettersrange over { , } .Note that if we set g µν = δ µν (planar case), the Hamil-tonian becomes the well studied two-channel model forthe 2D p -wave superfluid[6, 8, 9]. While for generalcurved surfaces, our model shows several remarkable fea-tures because of the nontrivial metric g µν .First of all, to preserve the covariance of the Hamil-tonian we introduce both the covariant and contravari-ant molecular operators d ν and d ν , which are related by d µ = g µν d ν .Secondly, to guarantee the particle number conserva-tion, the commutation relations for both fields become { ψ ( x ) , ψ † ( y ) } = 1 √ g δ (2) ( x − y ) , (3)[ d µ ( x ) , d ν ( y )] = 1 √ g δ µν δ (2) ( x − y ) . (4)Given these relations, it is straightforward to define theconserved particle number operator, N = N F + 2 N B = (cid:90) M d s (cid:18) ψ † ψ + 2 d † ν d ν (cid:19) . (5)Moreover, the covariant and contravariant derivative inthe Hamiltonian lead to an extra term due to the connec-tion of the surface. Take the unit sphere (radius R = 1)as an example. We shall use the usual spherical coordi-nates x = θ, x = φ . The fermion kinetic energy termcan be computed, ψ † D ν D ν ψ = − ψ † L ψ, (6)where L = − ( θ ∂∂θ (sin θ ∂∂θ ) + θ ∂ ∂φ ) is the totalangular momentum operator.While for the molecular field, since the operator D ν D ν is acting on a contravariant vector field d α , the deriva-tive will give an extra connection term. Consequently, itcannot be simply replaced by L operator. Inspired bythe two-body calculation[20], we recombine d α operatorsand define d †± = √ ( d † ± i sin θd † ), which represents the close channel molecules with different chiralities(angularmomentums). Using these two operators, the kinetic en-ergy can be expressed in a simple form, d † α D ν D ν d α = − d †± ( L ± ˆ r × A ) d ± (7)with A = − cot θ ˆ e φ .Besides the common L term, there emerges an extragauge field A for the molecular field, and the effectivecharges for molecules with different chiralities are oppo-site. The corresponding magnetic field ∇ × A = ˆ e r /R is exactly the magnetic field of a Dirac monopole withcharge 1. Consequently, the superfluid ground state nec-essarily supports topological defects such as vortices anddomain walls. Two-body problem.
We first consider a two-fermionproblem on the sphere, which helps to relate the bareparameters (cid:15) and λ to physical parameters. A two-bodystate can be generally written as | ψ (cid:105) = (cid:20)(cid:90) ϕ ( r , r ) ψ † ( r ) ψ † ( r ) + (cid:90) χ ± ( r ) d †± ( r ) (cid:21) | (cid:105) , (8)where ϕ ( r , r ) is the two-body wave function and χ ± ( r )is the dimer wave function.Substituting it into ( E − H ) | ψ (cid:105) = 0, we obtain theSchr¨odinger equation for ϕ ( r , r ) and χ ± ( r ), (cid:18) E − L m f − L m f (cid:19) ϕ = − iλχ ± ( r c ) D ± δ (2) ( r ) , (9) (cid:18) E − ( L ± ˆ r × A ) m f − (cid:15) (cid:19) χ ± = i λD ± ϕ (cid:12)(cid:12) r =0 . (10)where r c = r + r and r = r − r are the center-of-massand relative coordinates, and D ± is defined as D ± = √ ( D ± i sin θD ).The term with δ -function on the R.H.S. of Eq. (9) maybe regarded as a contact pseudopotential[23, 24], whichcan be replaced by a short-range Bethe-Peierls boundarycondition. This helps to relate (cid:15) and λ to physical pa-rameters such as scattering area s and effective range r .Both s and r describe the low energy scattering proper-ties of the p -wave interaction[25–27], which lead to uni-versal physics for many-body system[28–32]. It can beshown that,1 λ = − m f π log Λ r , (cid:15)λ = − m f s + m f π Λ , (11)where Λ is the momentum cutoff.Note that the physical parameters and the renormal-ization relation (11) take the exactly same form as thetwo-body problem on a 2D plane[6, 20], which reflectsthe short-range nature of the interaction.Detailed results of the two-body problem can be foundin ref.[20]. Here we make a few remarks on the interestingmany-body effects from the two-body consideration.i). As mentioned previously, the emergence of mag-netic field A will significantly change the ground stateconfigurations of order parameter. Especially in the s → + ∞ (BEC) limit, where two fermions can form adeeply bound state, the many-body system can be viewedas a cloud of weakly interacting bosonic molecules. Andthe bosons should condense to the lowest Landau levelon the sphere. In the monopole field A , the lowest Lan-dau level is three-fold degenerate and the eigen func-tions(monopole harmonics) are given by Wigner-D ma-trices D ∗ m, ± [33]. This also serves as a benchmark for ourmany-body calculation.ii). From Eq.(10), one can see that two dimer wavefunctions χ + and χ − are coupled to each other throughtwo-body wave function ϕ . The coupling correspondsto the physical process that a d + molecule breaks intotwo fermions and recombines into a d − molecule on thesphere. This transition is a unique feature of curved sur-faces because the angular momentum conservation wouldforbid the transition process on the plane. Thus, on themany-body level, we may expect the p θ + ip φ Cooperpairs can also be transformed into the p θ − ip φ Cooperpairs due to the curvature effect, which can enrich possi-ble phases of the ground state.
Many-body ground state.
To study the ground state ofthe many-body system, we apply the BCS theory. Byapproximating the molecule field d α by its expectationvalue ∆ α ( x ) ≡ (cid:104) d α ( x ) (cid:105) (correspondingly, ∆ ± ( x ) is theorder parameter of p θ ± ip φ superfluid), we obtain themean-field Hamiltonian, H BCS = ∆ ∗ α (cid:18) − D ν D ν m f + (cid:15) − µ (cid:19) ∆ α + 12 (cid:0) ψ † ψ (cid:1) (cid:18) ˆ h − iλ ∆ ν D ν − iλ ∆ ∗ ν D ν − ˆ h (cid:19) (cid:18) ψψ † (cid:19) , where ˆ h = L m f − µ .The mean-field Hamiltonian is diagonalized througha Bogoliubov transformation α i = (cid:82) d s u ∗ i ( x ) ψ ( x ) + v ∗ i ( x ) ψ ( x ) † . ( u i , v i ) T is the normalized Nambu spinorfield on the sphere and satisfies the BdG equations, (cid:18) ˆ h OO † − ˆ h (cid:19) (cid:18) u i v i (cid:19) = E i (cid:18) u i v i (cid:19) , (12)where O = − iλ ∆ ν D ν − iλ ( D ν ∆ ν ). It is straightforwardto show that O satisfies O † = − O ∗ which preserves theparticle-hole symmetry.The gap equation is obtained through the variation ofthe ground state energy δ (cid:104) H (cid:105) g /δ ∆ ∗ α = 0, which gives (cid:18) ( L ± ˆ r × A ) m f + (cid:15) − µ (cid:19) ∆ ± = iλ (cid:88) n u i D ± v ∗ i . (13)And the number equation is given by (cid:104) N (cid:105) = (cid:90) d s (cid:18) (cid:104) ψ † ψ (cid:105) + 2∆ ∗ α ∆ α (cid:19) , (14) where (cid:104) ψ † ψ (cid:105) = (cid:80) i | v i | . The summation (cid:80) i should berestricted to 0 ≤ E i ≤ Λ / m f [34].To solve the BdG equation together with the gap equa-tion and the number equation self consistently, we expandthe order parameter as∆ ± ( θ, φ ) = (cid:88) m ∆ m ± ( θ ) e imφ . (15)Note that both the BdG equation and the gap equationpreserve the SO(2) rotation symmetry along z axis. Thisallows us to simplify the numerics by assuming the or-der parameter only contains one m component. And thesymmetry will assure the other m components remainzero during the iterative process.Thus, for different quantum number m , we may cal-culate the order parameter ∆ m ± ( θ ) separately. And theground state is obtained by comparing their free energies. Numerical results.
Since the effective range for p -wavescattering is approximately a constant across the Fesh-bach resonance[1–5], we fix the ratio of effective rangeand the radius of the sphere r /R = 0 .
01 in all cal-culations. Then the ground state properties are fullycontrolled by two dimensionless parameters 1 /k F R and − /k F s , where k F is related to the fermion density n by k F = √ πn .Here we only present the results of m = 0 and m = ± m = 0 chan-nel ∆ ± ( θ ) for different − /k F s . Firstly, we could seethat the order parameters always vanish at the northand south pole. Although the order parameters are real,the emergent gauge field gives nonzero velocity[19], whichleads to a 2 π phase winding. This is why the m = 0 so-lutions correspond to a phase with two vortices locate onthe north and south poles respectively. Secondly, as wetune the interaction, there are two significant phenomenademonstrating the effects of sphere from two aspects.One is the change of vortex size. When we tune theinteraction towards the BEC( − /k F s → −∞ ) or theBCS( − /k F s → ∞ ) side, the vortices grow larger, whichis consistent with the trend of healing length in a pla-nar system. In the extreme BEC or BCS limit, the vor-tex size saturates when it is comparable with the systemsize R . And the order parameters are well approximatedby monopole harmonics D ∗ , ± ∝ sin θ . This is becausein both limits the interaction terms in the Ginzburg-Landau theory are negligible and the minimization of ki-netic energy forces the order parameters to keep in thelowest Landau level, whose wave function is given by themonopole harmonics.The other is the change of the ratio between ∆ + and∆ − . In subplot (a) and (b), we can see that the order FIG. 2: The order parameter ∆ ± ( θ ) in m = 0 channel for different interaction strength 1 /k F s . For figure (a) and (b)1 /k F R = 0 .
14, for figure (c) 1 /k F R = 1 . ± as functions of θ . 1 /k F R =0 .
07 (a): − /k F s = 0 .
76 (b): − /k F s = 3 . parameter ∆ + is much larger than ∆ − , which means theground state is a p θ + ip φ superfluid. This is also consis-tent with the previous study on p -wave superfluid on theplane[6, 9]. While this conclusion changes completely inthe deep BCS regime. Since the Cooper pair size growsexponentially in this limit, it will soon become compara-ble with the radius R . This makes the Cooper pair feelsthe curvature of the sphere. As a result, the two-bodyprocess that couples ∆ + pairs to ∆ − pairs becomes verysignificant in BCS side. Indeed, our numerics shows thatthe system is in p φ pairing instead of p θ + ip φ in BCSlimit.In fig.(3) we plot the order parameters in m = 1 chan-nel ∆ ± ( θ ) for different interaction strength − /k F s (the m = − m = 1 phase only exists in BCSside, we only present the result in s < m = ± p θ ± ip φ phases.Similar to the m = 0 phase, the size of the domain wallalso increases towards the BCS limit. In the deep BCSregime, the domain wall moves to the equator and its sizesaturates when it is comparable with the radius of thesphere. And the order parameters ∆ ± are well approxi- mated by monopole harmonics D ∗ , ± ∝ (1 ∓ cos θ ) e iφ .After comparing the energies of different m channels,we obtain a phase diagram of the system. As shown infig.(1), the ground state may exhibit topological defectssuch as vortex pairs( m = 0) or domain walls( m = ± − /k F s and radius1 /k F R .In the BEC limit, the system is described by repulsiveinteracting molecules condensing at the lowest Landaulevel. Then it is straightforward to show that the in-teraction energy of m = 0 channel is lower. Numericssuggest this remains ture for the whole BEC side andthe phase diagram on this side can be understood.On the BCS side, the system may undergo phase tran-sitions from the m = ± m = 0 phase,when we change the radius R . Recall that the domainwall energy is proportional to its length. This makes thedomain wall phase energetically unfavored for large R ,which gives the phase transition in 1 /k F R → /k F R → ∞ limit, this argument breaksdown since the healing length is now comparable to thesystem size R . To see what happens in this limit, we didan calculation based on Ginzburg-Landau theory (see ap-pendix). Indeed, it shows that the vortex pair phase haslower free energy in small R limit.Furthermore, the m = ± re-flection symmetry is broken, and the domain wall movesto the northern or southern hemisphere. The Z phasetransition can be explained by a competition betweenthe domain wall energy and the kinetic energy of thesuperfluid[19]. Conclusions.
We consider a microscopic model de-scribing fermions with p-wave interaction on a sphereand study its ground state properties. We show that thecurvature of the sphere induces an emergent magneticmonopole field coupled to the closed channel molecules,which necessarily leads to topological defects in super-fluid ground state. By solving the BdG together withgap and number equations self consistently, we identifythree phases with different types of topological defectsand obtain the phase diagram.
Acknowledgment
We are grateful to Hui Zhai, Ran Qiand Zhenhua Yu for helpful discussion. Especially wewould like to thank Hui Zhai for helpful suggestions andcareful reading of this manuscript. ∗ Electronic address: [email protected][1] C. A. Regal, C. Ticknor, J. L. Bohn, and D. S. Jin, Phys.Rev. Lett. , 053201 (2003).[2] J. Zhang, E. G. M. van Kempen, T. Bourdel, L.Khaykovich, J. Cubizolles, F. Chevy, M. Teichmann, L.Tarruell, S. J. J. M. F. Kokkelmans, and C. Salomon,Phys. Rev. A , 030702(R) (2004).[3] K. G¨unter, T St¨oferle, H. Moritz, Michael K¨ohl, and T.Esslinger, Phys. Rev. Lett. , 230401 (2005).[4] Y. Inada, M. Horikoshi, S. Nakajima, M. Kuwata-Gonokami, M. Ueda, and T. Mukaiyama, Phys. Rev.Lett. , 100401 (2008).[5] C. Luciuk, S. Trotzky, S. Smale, Z. Yu, S. Zhang, and J.H. Thywissen, Nat. Phys. , 599-605 (2016).[6] V. Gurarie, L. Radzihovsky, and A. V. Andreev, Phys.Rev. Lett. , 230403 (2005).[7] C.-H. Cheng and S.-K. Yip, Phys. Rev. Lett. , 070404(2005)[8] V. Gurarie, and L. Radzihovsky, Phys. Rev. B , 212509(2007).[9] J. Levinsen, N. R. Cooper, and V. Gurarie, Phys. Rev.A , 063616 (2008).[10] G. Cao, L. He, and P. Zhuang, Phys. Rev. A , 013613(2013).[11] A.Yu. Kitaev, Ann. Phys. N.Y. , 010506 (2007).[13] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S.Das Sarma, Rev. Mod. Phys. , 1083 (2008).[14] E. Vetsch, D. Reitz, G. Sagu´e, R. Schmidt, S. T.Dawkins, and A. Rauschenbeutel, Phys. Rev. Lett. ,203603 (2010).[15] D. Reitz and, A. Rauschenbeutel, Opt. Commun. ,4705 (2012).[16] V. Bendkowsky, B. Butscher, J. Nipper, J. P. Shaffer, R.L¨ow, and T. Pfau, Nature , 1005 (2009).[17] A. Gaj, A. T. Krupp, J. B. Balewski, R. L¨ow, S. Hoffer-berth and T. Pfau, Nat. Commun. , 4546 (2014).[18] N. Read, and D. Green, Phys. Rev. B , 10267 (2000).[19] S. Moroz, C. Hoyos, and L. Radzihovsky, Phys. Rev. B , 024521 (2016).[20] Z.-Y. Shi, and H. Zhai, arXiv:1510.05815 (2015).[21] C. Hoyos, S. Moroz, and D. T. Son, Phys. Rev. B ,174507 (2014).[22] S. Moroz and C. Hoyos, Phys. Rev. B , 064508 (2015).[23] K. Huang, and C. N. Yang, Phys. Rev. , 767 (1957).[24] T. D. Lee, K. Huang, and C. N. Yang, Phys. Rev. ,1135 (1957).[25] B. Gao, Phys. Rev. A , 012702 (2009).[26] P. Zhang, P. Naidon, and M. Ueda, Phys. Rev. A , 062712 (2010).[27] E. Braaten, P. Hagen, H.-W. Hammer, and L. Platter,Phys. Rev. A , 012711 (2012).[28] Z. Yu, J. H. Thywissen, and S. Zhang, Phys. Rev. Lett. , 135304 (2015).[29] S. M. Yoshida and M. Ueda, Phys. Rev. Lett. , 135303(2015).[30] M.-Y. He, S.-L. Zhang, H. M. Chan, Q. Zhou, Phys. Rev.Lett. , 045301 (2016).[31] Y.-C. Zhang, S. Zhang, arXiv:1611.09499.[32] P. Zhang, Z. Yu, arXiv:1611.09454.[33] D ∗ , ± ∝ sin θ, D ∗ , ± ∝ (1 ∓ cos θ ) e iφ [34] In atomic gases, the momentum cutoff is set by the in-verse of the van der Waals radius, which remains a con-stant across Feshbach resononses. In our calulations, wealso fix the cutoff Λ = 80 /R (cid:29) k F .[35] Due to the limit of the computation power, it is difficultto keep a high density to approach the BEC limit withhigh precision. APPENDIX: GINZBERG-LANDAU IN BCSLIMIT
The Hamiltonian H BCS is largely simplified if we as-sume ∆ ± to be some function with undetermined numer-ical factors. In 1 /k F R → ∞ limit, inspired by numerics,we assume ∆ + = − ∆ − = ∆ sin θ/λ for m = 0 channel,where ∆ is a constant. Set the mass of the atom to be 1,We have:∆ + ψ † iD − ψ † + ∆ − ψ † iD + ψ † = √ ψ † ∂ φ ψ † /λ (16)Hence the Hamiltonian can be written as: (cid:90) d Ω 12 Ψ † H Ψ + 2∆ λ sin θ ( 14 − µ + (cid:15) ) (17)Ψ is the two component spinor ( ψ, ψ † ) T and H =( L − µ ) σ z − i √ ∂ φ σ y . Expand the result in sphericalharmonics and diagonalize the matrix give us the energyfor ground state: E = − (cid:88) lm (cid:114) ( l ( l + 1)2 − µ ) + 8∆ m + 16 π ∆ λ ( 14 − µ + (cid:15) ) (18)Similarly, for m = 1, the numerics suggest an assump-tion: ∆ ± = ∆ e iφ (1 ∓ cos θ ) / √ λ . Then we have:∆ + ψ † iD − ψ † + ∆ − ψ † iD + ψ † = − ∆ ψ † iL + ψ † /λ (19)We have used the fact that L ± = e ± iφ ( ± ∂ θ + i cot θ∂ φ ).Similar procedure give the energy: E = − (cid:88) lm (cid:114) ( l ( l + 1)2 − µ ) + 4∆ ( l − m )( l + m + 1)+ 16 π ∆ λ ( 14 − µ + (cid:15) ) (20)In BCS limit the ∆ is small. We could expand energyfor ∆ and get a Ginsberg-Landau free energy: E = − r ∆ + b (21)The energy have its minimum at ∆ = rb and lead to E = − r / b .Using n (cid:88) i = n ( n + 1)(2 n + 1)6 , it is easy to show that the second order (∆ ) contributionfrom the summation is the same for any l after sum over m . Then r is the same for m = 0 and m = 1. For b , b m =0 = (cid:88) l l ( l + 1)(2 l + 1) (cid:0) l + 3 l − (cid:1)
15 ( l + l − µ ) / (22) b m =1 = (cid:88) l l ( l + 1) (cid:0) l + 6 l + 4 l + 1 (cid:1) (cid:16) ( l + l − µ ) (cid:17) / (23) Then a subtraction tell us that: δb = b m =0 − b m =1 = (cid:88) l l ( l + 1) (cid:0) l + 3 l − l − (cid:1) (cid:16) ( l + l − µ ) (cid:17) / (24)Recall that a larger b gives a higher energy. For fixeddensity, the contribution comes mainly near l ( l +1) ∼ µ ,then for large density, the m = 1 wins while for extremelysmall density(large 1 /k F R ), mm