Core structure of two-dimensional Fermi gas vortices in the BEC-BCS crossover region
CCore structure of two-dimensional Fermi gas vortices in the BEC-BCS crossover region
Lucas Madeira, ∗ Stefano Gandolfi, and Kevin E. Schmidt Department of Physics, Arizona State University, Tempe, Arizona 85287, USA Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA (Dated: October 8, 2018)We report T = 0 diffusion Monte Carlo results for the ground-state and vortex excitation ofunpolarized spin-1/2 fermions in a two-dimensional disk. We investigate how vortex core structureproperties behave over the BEC-BCS crossover. We calculate the vortex excitation energy, densityprofiles, and vortex core properties related to the current. We find a density suppression at thevortex core on the BCS side of the crossover and a depleted core on the BEC limit. Size-effectdependencies in the disk geometry were carefully studied. I. INTRODUCTION
The study of cold Fermi gases has proven to be avery rich research field, and the investigation of low-dimensional systems has become an active area in thiscontext [1, 2]. Particularly, the two-dimensional (2D)Fermi gas has attracted much interest recently. It was theobject of several theoretical investigations [3–8], but itsexperimental realization, using a highly anisotropic po-tential, was a milestone in the study of these systems [9].Many other studies have been carried out since [10, 11].Quantum Monte Carlo (QMC) methods were successfullyemployed to compute several properties of the BEC-BCScrossover. These methods include diffusion Monte Carlo(DMC) [12, 13], auxiliary-field quantum Monte Carlo[14], and lattice Monte Carlo [15–17]. The fact that afully attractive potential in 2D always supports a boundstate, and the ability to vary the interaction strength overthe entire BEC-BCS crossover regime offers rich possibil-ities for the study of these systems.The presence of quantized vortices is an indication ofa superfluid state in both Bose and Fermi systems. Inthree-dimensional (3D) systems, much progress has beenmade [18–21], including the observation of vortex lat-tices in a strongly interacting rotating Fermi gas of Li[22]. With the recent progress on the 2D Fermi gases,it seems natural to also extend the theoretical study ofvortices to these systems. Interest is further augmentedin 2D, where a Berezinksii-Kosterlitz-Thouless transition[23, 24] could take place at finite temperatures, and pairsof vortices and antivortices would eventually condense toform a square lattice [25].We are interested in how the properties of a vortexchange over the BEC-BCS crossover. In this work wefocus on ultracold atomic Fermi gases, but it is note-worthy that a duality is expected between neutron mat-ter and superfluid atomic Fermi gases. In 3D, bothultracold atomic gases and low-density neutron matterexhibit pairing gaps of the order of the Fermi energy[26]. Neutron-matter properties depend on the interac-tion strength and, unlike the Fermi atom gases, the pos- ∗ [email protected] sibility of microscopically tuning interactions of neutron-matter is not available. However, we can study neu-tron pairing by looking at the BCS side of the crossover[27, 28]. Vortex properties are also of significant interestin neutron matter [29, 30] because a significant part ofthe matter in rotating neutron stars is superfluid, andvortices are expected to appear. Moreover, phases callednuclear pasta, where neutrons are restricted to 1D or 2Dconfigurations, are predicted in neutron stars [30, 31].We report properties of a single vortex in a 2D Fermigas. We considered the ground-state to be a disk withhard walls and total angular momentum zero, and thevortex excitation corresponds to each fermion pair hav-ing angular momentum (cid:126) . Hopefully, our results willmotivate experiments to increase our understanding ofvortices in 2D Fermi gases.This work is structured as it follows. In Sec. II we in-troduce the methodology employed. In Sec. II A we dis-cuss aspects of finite-size fermionic systems, we brieflyintroduce 2D scattering in Sec. II B, Sec. II C is devotedto the wave functions employed for the bulk, disk, andvortex systems, and we summarize the employed QMCmethods in Sec. II D. The results are presented in Sec. III.Sec. III A contains the ground-state energies in the diskgeometry and discussions on size-effects. In Sec. III B wepresent the vortex excitation energy. The determinationof the crossover region is done in Sec. III C. Density pro-files of the vortex and ground-state systems are shown inSec. III D. Properties of the vortex core are discussed inSec. III E. Finally, a summary of the work is presented inSec. IV. II. METHODS
Previous simulations of vortices in 3D bosonic systems,such as He, have often employed a periodic array ofcounter-rotating vortices, which enables the usage of pe-riodic boundary conditions. In the He calculations ofRef. 32, the simulation cell consisted of 300 particles infour counter-rotating vortices. If we had employed a sim-ilar methodology, we would need the same number offermion pairs, i.e., a system with 600 fermions. Thereare simulations of fermionic systems that have been per- a r X i v : . [ c ond - m a t . qu a n t - g a s ] M a y formed with this number of particles, but the variancerequired for a detailed optimization is beyond the scopeof this work. Instead, we considered a disk geometry sim-ilar to the one used in Ref. [33] for DMC simulations ofthe vortex core structure properties in He.
A. Finite-size systems
We are interested in the interacting many-body prob-lem, but it is useful to first consider the non interactingcase. In this section we compare the energy of finite-size2D systems to the results in the thermodynamic limit.First let us consider the case of N fermions in a squareof side L with periodic boundary conditions. The single-particle states are plane waves ψ k n ( r ) = e i k n · r /L , withwave vector k n = 2 πL ( n x ˆ x + n y ˆ y ) . (1)The eigenenergies are E n = (cid:126) k n / m , where m is themass of the fermion. At T = 0, all states with en-ergy up to the Fermi energy (cid:15) F = (cid:126) k F / m , where k F is the Fermi wave number, are occupied. A shellstructure arises from the fact that different combina-tions of n x and n y in Eq. (1) yield the same | k n | .The closed shells occur at total particle number N =(2 , , , , , , , · · · ). The free gas energy of a fi-nite system with N fermions, E bulk F G ( N ), is readily cal-culated by filling the lowest energy states described byEq. (1). In the thermodynamic limit, which correspondsto N, L → ∞ and n = N/L held constant, the en-ergy per particle of the free gas is E FG = (cid:15) F / k F = √ πn .Now let us consider the case of N fermions in a diskof radius R with a hard wall boundary condition, i.e.,the wave function must vanish at R . The single-particlestates are ψ νp ( ρ, ϕ ) = N νp J ν (cid:18) j νp R ρ (cid:19) e iνϕ , (2)where ( ρ, ϕ ) are the usual polar coordinates, N νp is a nor-malization constant, J ν are Bessel functions of the firstkind, and j νp is the p -th zero of J ν . The quantum num-ber ν can take the values 0 , ± , ± , · · · and p = 1 , , · · · .The corresponding eigenenergies are E νp = (cid:126) m (cid:18) j νp R (cid:19) . (3)This system also presents a shell structure, due to the en-ergy degeneracy of single-particles states with the same | ν | , with shell closures at total particle number N =(2 , , , , , , , , , , · · · ). Notice that the en-ergy levels of the bulk system are much more degeneratethan the ones of the disk. In practice this means thatmore shells are needed to describe a disk with a given N . The free gas energy for the disk, E disk F G ( N ), can be calculated analogously to the bulk case using the energylevels of Eq. (3). The thermodynamic limit for this casecorresponds to R → ∞ with n = N/ ( π R ) held constant,and E F G and k F go to the same expressions as the bulkones.The comparison between the free gas energy of finitesystems in the bulk case and in the disk geometry is notimmediate due to the presence of hard walls in the latter.In order to compare the free gas energy in both geome-tries, we define E disk0 ( N ) = E disk F G ( N ) − λ s (cid:114) nπN , (4)in which we separated the total energy E disk F G ( N ) into abulk component, E disk0 ( N ), and a surface term, the sec-ond term on the RHS. For further discussions on the func-tional form of the surface term, see Sec. III A. Figure 1shows E bulk F G ( N ) and E disk0 ( N ), with λ s = 17 . E F G k − F ,at the same density. The value of λ s , within a 0.2% error,was determined by fitting the data for 10 (cid:54) N (cid:54)
226 tothe functional form of Eq. (4).The disk presents a considerably higher free gas energy,if compared to the bulk system, due to the presence ofhard walls, but the difference between them is rapidlysuppressed as we increase the particle number. E F G ( A ) [ un i t s o f ε F / ] N E bulk E FIG. 1. (Color online) Free-gas energy for finite-size systemsas a function of the number of particles N , where the dottedlines are drawn to guide the eye. The (red) closed circles de-note the energy of the bulk system, E bulk FG ( N ), and the (green)open circles indicate the bulk energy component in the diskgeometry, E disk0 ( N ), as defined in Eq. (4). Local minima in E bulk FG ( N ) correspond to shell closures. B. Scattering in 2D
Two-body scattering by a finite-range potential V ( r )in 2D is described by the Schr¨odinger equation. We sep-arate the solutions into radial R ( r ) and angular P ( φ )parts, the latter being a constant for s -wave scattering.The two-body equation for an azimuthally symmetric ( s -wave) solution is (cid:20) − (cid:126) ∇ m r + V ( r ) (cid:21) u ( r ) = (cid:126) k m r u ( r ) , (5)where m r is the reduced mass of the system, and (cid:126) k / m r is the scattering energy. The scattering length a and effective range r eff can be easily determined fromthe k → u ( r ), and its asymptoticform y . We choose the solution y ( r ) = − ln (cid:16) ra (cid:17) , (6)and we match u and y , and their derivatives, outsidethe range of the potential.In 2D, the low-energy phase shifts δ ( k ), a , and effectiverange r eff , are related by [34]cot δ ( k ) ≈ π (cid:20) ln (cid:18) ka (cid:19) + γ (cid:21) + k r , (7)where γ = 0 . . . . is the Euler-Mascheroni constant,and the effective range is defined as [35] r = 4 (cid:90) ∞ ( y ( r ) − u ( r )) r dr. (8)Equation (7) is often called the shape-independent ap-proximation because it guarantees that a broad rangeof well-chosen potentials can be constructed to describelow-energy scattering. We consider the modified Poschl-Teller potential V ( r ) = − v (cid:126) m r µ cosh ( µr ) , (9)where v and µ can be tuned to reproduce the desired a and r eff .Bound-states occur for purely attractive potentials forany strength in 2D. If we continually increase the depthof V ( r ), a will eventually reach zero, and then it divergesto + ∞ when a new bound-state is created. The bindingenergy of the pair is given by (cid:15) b = − (cid:126) ma e γ . (10)We chose values of v and µ such that only one bound-state is present, and k F r eff is held constant at 0.006 [13].This choice guarantees that the systems studied in thiswork are in the dilute regime, since r (cid:29) r eff , where r = 1 / √ πn is of order of the interparticle spacing. C. Wave functions
The BCS wave function, which describes pairing ex-plicitly, has been successfully used in a variety of stronglyinteracting Fermi gases systems, such as: 3D [36] and 2D[13] bulk systems, vortices in the unitary regime [21], two-component mixtures [37, 38], and many other sys-tems. This wave function, projected to a fixed number ofparticles N (half with spin-up and half with spin-down),can be written as the antisymmetrized product [39] ψ BCS ( R , S ) = A [ φ ( r , s , r , s ) φ ( r , s , r , s ) . . .φ ( r N − , s N − , r N , s N )] , (11)where R is a vector containing the particle positions r i , S stands for the spins s i , and φ is the pairing function,which is given by φ ( r , s, r (cid:48) , s (cid:48) ) = ˜ φ ( r , r (cid:48) ) [ (cid:104) s s (cid:48) | ↑ ↓(cid:105) − (cid:104) s s (cid:48) | ↓ ↑(cid:105) ] , (12)where we have explicitly included the spin part to im-pose singlet pairing. The assumed expressions for ˜ φ de-pend on the system being studied (see Secs. II C 1, II C 2,and II C 3). Since neither the Hamiltonian or any op-erators in the quantities we calculate flip the spins, weadopt hereafter the convention of primed indexes to de-note spin-down particles and unprimed ones to refer tospin-up particles. Equation (11) reduces to ψ BCS ( R , S ) = A [ φ ( r , s , r (cid:48) , s (cid:48) ) φ ( r , s , r (cid:48) , s (cid:48) ) . . . φ ( r N/ , s N/ , r N/ (cid:48) , s N/ (cid:48) )] , (13)where the antisymmetrization is over spin-up and/orspin-down particles only. This wave function can be cal-culated efficiently as a determinant [40].In addition to fully paired systems, it is also possible tosimulate systems with unpaired particles [36], describedby single particle states Φ( r ). For q pairs, u spin-up,and d spin-down unpaired single particles states, N =2 q + u + d , we can rewrite Eq. (13) as ψ BCS ( R , S ) = A [ φ ( r , s , r (cid:48) , s (cid:48) ) · · · φ ( r q , s q , r q (cid:48) , s q (cid:48) )Φ ↑ ( r q +1 ) · · · Φ u ↑ ( r q + u )Φ ↓ ( r ( q +1) (cid:48) ) · · · Φ d ↓ ( r ( q + d ) (cid:48) )] . (14)We also included a two-body Jastrow factor f ( r ij (cid:48) ), r ij (cid:48) = | r i − r j (cid:48) | , which accounts for correlations betweenantiparallel spins. It is obtained from solutions of thetwo-body Schr¨odinger’s equation (cid:20) − (cid:126) ∇ m r + V ( r ) (cid:21) f ( r < d ) = λf ( r < d ) , (15)with the boundary conditions f ( r > d ) = 1 and f (cid:48) ( r = d ) = 0, where d is a variational parameter, and λ isadjusted so that f ( r ) is nodeless. The total trial wavefunction is written as ψ T ( R , S ) = (cid:89) i,j (cid:48) f ( r ij (cid:48) ) ψ BCS ( R , S ) . (16)
1. Bulk system
The assumed form of the pairing function for the bulkcase is the same as Ref. [36],˜ φ bulk ( r , r (cid:48) ) = n c (cid:88) n =1 α n e i k n · ( r − r (cid:48) ) + ˜ β ( | r − r (cid:48) | ) , (17)where α n are variational parameters, and contributionsfrom momentum states up to a level n c are included.Contributions with n > n c are included through the ˜ β function given by˜ β ( r ) = (cid:40) β ( r ) + β ( L − r ) − β ( L/
2) for r (cid:54) L/
20 for r > L/ β ( r ) = [1 + cbr ][1 − e − dbr ] e − br dbr , (19)where r = | r − r (cid:48) | and b , c , and d are variational param-eters. This functional form of β ( r ) describes the short-distance correlation of particles with antiparallel spins.We consider b = 0.5 k F , d = 5, and c is adjusted so that ∂ ˜ β/∂r = 0 at r = 0.
2. Disk
The pairing function for the disk geometry is con-structed using the single-particle orbitals of Eq. (2). Eachpair consists of one single-particle orbital coupled withits time-reversed state. This ansatz has been used beforein the 3D system [21], a cylinder with hard walls, andthe form presented here is analogous to that one if wedisregard the z components. We supposed the pairingfunction to be˜ φ disk ( r , r (cid:48) ) = n c (cid:88) n =1 ˜ α n N νp J ν (cid:18) j νp R ρ (cid:19) J ν (cid:18) j νp R ρ (cid:48) (cid:19) e iν ( ϕ − ϕ (cid:48) ) + ¯ β ( r , r (cid:48) ) , (20)where the ˜ α n are variational parameters, and n is a la-bel for the disk shells, such that different states with thesame energy are associated with the same variational pa-rameter. The ¯ β function is similar to ˜ β employed in thebulk system, but we modify it to ensure the hard wallboundary condition is met,¯ β ( r , r (cid:48) ) = N J (cid:0) j ρ R (cid:1) J (cid:16) j ρ (cid:48) R (cid:17) × [ β ( r ) + β (2 R − r ) − β ( R )] for r (cid:54) R r > R (21)and β has the same expression as the bulk case, Eq. (19).
3. Vortex
The vortex excitation is accomplished by consideringpairing orbitals which are eigenstates of L z with eigen-values ± (cid:126) . This is achieved by coupling single-particlestates with angular quantum numbers differing by one. In this case we used pairing orbitals of the form˜ φ vortex ( r , r (cid:48) ) = n c (cid:88) n =1 ¯ α n N νp N ν − p × (cid:26) J ν (cid:18) j νp R ρ (cid:19) J ν − (cid:18) j ν − p R ρ (cid:48) (cid:19) e i ( νϕ − ( ν − ϕ (cid:48) ) + J ν (cid:18) j νp R ρ (cid:48) (cid:19) J ν − (cid:18) j ν − p R ρ (cid:19) e i ( νϕ (cid:48) − ( ν − ϕ ) (cid:27) , (22)where n is a label for the vortex shells, and ¯ α are vari-ational parameters. The largest contribution is assumedto be from states with the same quantum number p forthe radial part [21]. Equation (22) is symmetric underinterchange of the prime and unprimed coordinates, asrequired for singlet pairing.The ¯ β function of Eq. (21) is not suited to describe thevortex state because it is an eigenstate of L z with angularmomentum zero. We tried different functional forms thathad the desired angular momentum eigenvalue, but noneof them resulted in a significant lower total energy. Thus,we chose to employ only the terms in Eq. (22). D. Quantum Monte Carlo
The Hamiltonian of the two-component Fermi gas isgiven by H = − (cid:126) m N ↑ (cid:88) i =1 ∇ i + N ↓ (cid:88) i = j (cid:48) ∇ j (cid:48) + (cid:88) i,j (cid:48) V ( r ij (cid:48) ) , (23)with N = N ↑ + N ↓ , and V ( r ij (cid:48) ) given by Eq. (9). TheDMC method projects the lowest energy state of H froman initial state ψ T , obtained from variational MonteCarlo (VMC) simulations. The propagation, which iscarried out in imaginary time τ , can be written as ψ ( τ ) = e − ( H − E T ) τ ψ T , (24)where E T is an energy offset. In the τ → ∞ limit, onlythe lowest energy component Φ surviveslim τ →∞ ψ ( τ ) = Φ . (25)The imaginary time evolution is given by ψ ( R , τ ) = (cid:90) d R (cid:48) G ( R , R (cid:48) , τ ) ψ T ( R (cid:48) ) , (26)where G ( R , R (cid:48) , τ ) is the Green’s function associated with H . The Green’s function contains two pieces, a diffusionterm related to the kinetic operator, and a branchingterm related to the potential. We solve an importancesampled version of Eq. (26) iteratively, using the Trotter-Suzuki approximation to evaluate G ( R , R (cid:48) , τ ), which re-quires the time steps ∆ τ to be small. We circumvent thefermion-sign problem by using the fixed-node approxi-mation, which restricts transitions across a chosen nodalsurface [41]. Hence our estimates of energy expectationvalues are upper bounds.We carefully optimized the trial wave function ψ T ,since it is used in three ways: an approximation of theground-state in the VMC calculations, as an importancefunction, and to give the nodal surface for the fixed-node approximation. The variational parameters [42]in Eqs. (17), (20), and (22) were determined using thestochastic reconfiguration method [43].Expectation values of operators that do not commutewith the Hamiltonian, for example the current and den-sity, were calculated using extrapolated estimators [44] (cid:104) Φ | ˆ S | Φ (cid:105) ≈ (cid:104) Φ | ˆ S | ψ T (cid:105) − (cid:104) ψ T | ˆ S | ψ T (cid:105) + O [(Ψ − ψ T ) ] , (27)where we combine the results of VMC and DMC runs. III. RESULTS
We define the interaction strength η ≡ ln( k F a ). Largevalues of η correspond to the BCS side of the crossover,while small η are on the BEC side. We probed 0 . (cid:54) η (cid:54) .
5, which encompasses the crossover region (seeSec. III C). For all systems the number density is n = N/ ( π R ), and k F = √ N / R . A. Ground-state energy and size-effects
We used the pairing function of Eq. (17), and N = 26,to calculate the ground-state energy per particle of thebulk systems. Our results (see Table I) are in agreementwith previous DMC calculations [13]. TABLE I. Comparison between the ground-state energy perparticle of the bulk ( E bulk ) and disk systems as a function ofthe interaction strength η . The parameters E and λ s , seeEq. (28), are related to our assumption of the functional formof the ground-state energy per particle in the disk geometry. η E bulk [ E FG ] E [ E FG ] λ s [ E FG k − F ]0.00 -2.3740(3) -2.32(3) 6(2)0.25 -1.3316(3) -1.31(3) 8(2)0.50 -0.6766(2) -0.65(2) 8(1)0.75 -0.2562(2) -0.25(2) 11(1)1.00 0.0233(2) 0.03(1) 11(1)1.25 0.2149(2) 0.22(2) 12(1)1.50 0.3523(2) 0.34(1) 13(1) Previous DMC simulations of 2D Fermi gases foundthat N = 26 is well suited to simulate bulk properties ofsystems in the region studied here [13]. However, the diskgeometry presents more intricate size-dependent effects.We investigated how the ground-state energy depends onthe disk radius R . In the thermodynamic-limit, R → ∞ ,the energy per particle should go to the bulk value. Sinceour system has hard walls, the energy has a dependence on the “surface” of the disk. Including this surface term,the energy per particle can be fit to E disk ( R ) = E + λ s π R , (28)where E and λ s are constants related to the bulk andsurface terms, and λ s / (2 π R ) can be viewed as a surfacetension.A few words about Eq. (28) are in order. The relationbetween the thermodynamic properties of a confined fluidand the shape of the container where it is confined hasbeen an active field of study. Our choice was inspired byfunctional forms (see for example Ref. [45]) where, asidefrom the constant term, thermodynamical properties areexpressed as functions of the various curvatures of thecontainer. The next correction to this functional form ofthe energy per particle would include a term proportionalto R − . We found that the inclusion of such a term doesnot significantly improve our description of the ground-state energy.In order to determine the number of particles nec-essary to simulate systems in the disk geometry, withcontrollable size effects, we performed simulations with26 (cid:54) N (cid:54)
70, and all particles paired, i.e., only even val-ues of N . The dependence of E with the system size wasinvestigated by fitting our data using Eq. (28) for differ-ent intervals of R or, equivalently, different intervals of N .We found that fitting the data for 58 (cid:54) N (cid:54)
70 re-sulted in a good agreement between E bulk and E , thatis, we were able to separate the bulk portion of the en-ergy from the hard wall contribution in the disk geome-try. The resulting parameters of the fitting procedure aresummarized in Table I, and Fig. 2 shows the energy perparticle as a function of R for all interaction strengthsstudied in this work.The E values agree with the bulk energies within theerror bars, except for η = 0 and η = 0 . E is of order 0.01 E F , independent of theinteraction strength. Thus the relative error can be quitelarge for systems where the absolute value of the bulkenergy is small, as it is observed for η = 1 .
0. This is animprovement if compared to a similar DMC calculationin 3D [21] which used the same procedure to calculatethe ground-state energy per particle of a unitary Fermigas, where the discrepancy between the result and theknown bulk value was ≈ N = 26, then we sim-ply could not rely on the results. In our simulationswith 26 (cid:54) N (cid:54)
38 the discrepancies between E and η =1.5 0.36 0.38 0.40 η =1.25 0.16 0.18 0.20 η =1.0-0.12-0.10-0.08 E d i sk [ un i t s o f E F G ] η =0.75-0.56-0.54-0.52 η =0.5-1.22-1.20-1.18 η =0.25-2.24-2.22-2.2010.6 10.8 11.0 11.2 11.4 11.6 11.8 12.0 R [units of k
F-1 ] η =0.0 FIG. 2. (Color online) Ground-state energy per particle E disk as a function of the disk radius R for several interactionstrengths. The curves correspond to the assumed functionalform of Eq. (28), with the parameters given in Table I. Errorbars are smaller than the symbols. E bulk were as large as 50%, and in some cases the uncer-tainty in λ s was bigger than the value itself. Results with58 (cid:54) N (cid:54)
70 are much more well-behaved, and they arewithin computational capabilities.It is also noteworthy to mention that the energy con-tribution of the surface term, due to the presence of hardwalls, is more significant for the BCS side than in theBEC limit (see the λ s values in Table I). This is expected,since the largest energy contribution in the BEC sideshould be from the binding energy of the pairs, Eq. (10),and they are smaller than the BCS pairs so that surfaceeffects are smaller. One of our goals is to obtain the vor-tex excitation energy, which is the difference between thevortex and the ground-state energies. Since both systemshave hard walls, we expect that the surface effects willtend to cancel. B. Vortex excitation energy
The energy per particle of the vortex system is ob-tained using the pairing functions of Eq. (22). The vor-tex excitation energy is given by the difference betweenthe energy of the vortex and ground-state systems, forthe same number of particles. We performed simulationswith 58 (cid:54) N (cid:54)
70 and averaged the results.In Fig. 3 we show the vortex excitation energy per par- ticle as a function of the interaction strength. The energynecessary to excite the system to a vortex state increasesas we move from the BCS to the BEC limit. The insetshows the vortex and ground-state energies per particlefor η = 1 .
5, although the other interaction strengths dis-play the same qualitative behavior. E e xc [ un i t s o f E F G ] η =ln(k F a) E [ un i t s o f E F G ] N FIG. 3. (Color online) Vortex excitation energy per particle E exc as a function of the interaction strength η . The insetshows the ground-state (squares) and vortex (triangles) ener-gies per particle as a function of the number of particles N for η = 1 . C. Crossover region
In 2D, the BCS limit corresponds to k F a (cid:29) k F a (cid:28)
1, however unlike 3D wherethe unitarity is signaled by the addition of a two-bodybound state, there is no equivalent effect with two-bodysector in 2D. Nevertheless, we can determine the interac-tion strength for which we can add a pair to the systemwith zero energy cost. The chemical potential µ can beestimated as µ = ∂E∂N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Even N , (29)for each interaction strength, where the even numbercondition implies that all particles are paired. For eachvalue of η we used a finite difference formula to evaluateEq. (29), for 58 (cid:54) N (cid:54)
70 (see Fig. 4).We found that µ = 0 at η ≈ .
90 for the ground-state ofthe disk. Previous DMC simulations of 2D bulk systems[13] found that the chemical potential changes sign at η ≈ .
65. Although the results differ, most probablydue to the different geometry employed in this work, itis safe to assume that the interaction strength interval0 (cid:54) η (cid:54) . µ = 0 is at a smallerinteraction strength, η ≈ . -2.5-2.0-1.5-1.0-0.50.00.5 0.0 0.5 1.0 1.5 µ [ un i t s o f E F G ] η =ln(k F a) VortexGround-state E [ un i t s o f E F G ] N FIG. 4. (Color online) Chemical potential of the ground-state(triangles) and vortex (circles) as a function of the interactionstrength. The chemical potential changes sign at η ≈ .
90 forthe ground-state, and η ≈ .
85 for the vortex state. In theinset we show the total energy as a function of the number ofparticles for the ground-state of η = 1.5. Other interactionstrengths with positive (negative) µ have positive (negative)slopes. D. Density profile
We calculated the density profile D ( ρ ) along the radialdirection ρ for both the vortex and ground-state systems.The normalization is such that (cid:90) D ( ρ ) d r = 1 , (30)where the integral is performed over the area of the disk.The results are obtained using the extrapolation proce-dure of Eq. (27), which combines both VMC and DMCruns. It is noteworthy to point out that, although thedensities observed in VMC and DMC simulations differ,they are much closer than previous results in 3D [21]. Inthat calculation it was needed to explicitly include a one-body term in the wave function to maximize the densityoverlap between DMC and VMC runs, whereas in thiswork no such term was employed.Figure 5 shows the density profile of both the vortexand ground-state systems for N = 70 and η = 1 .
5. Theoscillations in the density profiles are much more pro-nounced than in a similar DMC calculation of a unitaryFermi gas in 3D [21]. In this 3D calculation a cylindri-cal geometry was employed, with hard walls and periodicboundary conditions along the axis of the cylinder. Thedensity profiles were obtained by averaging the resultsover the z direction of the axis of the cylinder, we there-fore expect more fluctuations in 2D where the particlesare confined to a plane. For the ground-state, the densityoscillations are surface effects. They are present in boththe interacting and non-interacting systems, as it can beseen in Fig. 5.In Fig. 6 we show the density profiles of the other in- D en s i t y [ un i t s o f k F ] ρ [units of k F-1 ]VortexGround-stateFree-gas
FIG. 5. (Color online) Density profile along the radial di-rection ρ of the vortex (red squares) and ground-state (greencircles) for N = 70 and η = 1 .
5. Although there is a den-sity suppression at the vortex core of ≈ teraction strengths studied in this work, 0 (cid:54) η (cid:54) . ≈
30% at η = 1 . η (cid:54) D en s i t y [ un i t s o f k F ] η =1.25 η =1.0 η =0.75 ρ [units of k F-1 ] η =0.5 ρ [units of k F-1 ] η =0.25 ρ [units of k F-1 ] η =0.0 FIG. 6. (Color online) Density profile along the radial di-rection ρ of the vortex (red squares) and ground-state (greencircles) for N = 70 and 0 (cid:54) η (cid:54) .
25. It is interesting toobserve that the density at the vortex core diminishes as wego from the BCS to the BEC limit, and at η (cid:54) .
25 the coreis completely depleted.
The regions close to the walls exhibit a characteristicbehavior due to the hard wall condition we imposed, asit can be seen in Figs. 5 and 6. In order to estimate thenumber of particles outside this region, we can define theparticle number a distance R from the center of the diskas N ( R ) = N (cid:90) π dϕ (cid:90) R dρ ρ D ( ρ ) . (31)For the case of Figs. 5 and 6 where N = 70, if we set R ∼ k − F , N is approximately between 40 and 45 forthe ground-state, and between 35 and 40 for the vortexsystems. Hence the number of particles in this regime islarger than the usual value of N = 26 employed in bulksystems [13].Additionally, we performed simulations of the vor-tex systems with an odd number of particles, i.e., oneunpaired particle was added to a fully paired system,Eq. (14) with q = 34, u = 1, and d = 0. We set its an-gular momentum to zero, Eq. (2) with ν = 0 and p = 1.In the BEC limit we observed a non-vanishing density atthe center of the disk, which suggests that the unpairedparticle fills the empty vortex core region. On the otherhand, in the BCS limit the density close to the wall in-creased, while the density at the origin was unchanged.We chose a qualitative discussion of this phenomenonbecause the required variance for a detailed optimizationis beyond the scope of this work. Future calculationsshould include quantities such as the one-body densitymatrix, which may contribute to an accurate quantita-tive approach. E. Vortex core size
The probability current density operator can be writ-ten as J ( r ) = 12 N N (cid:88) j =1 (cid:2) v j δ ( r − r j ) + δ ( r − r j ) v j (cid:3) , (32)where the velocity operator is v j = p j /m → − i (cid:126) ∇ j /m .We are interested in the angular component as a functionof the radial coordinate, J ϕ ( ρ ), because the position ofits maximum can be used as an estimate of the vortexcore size, J max ≡ J ϕ ( ρ = ξ ).We followed the extrapolation procedure of Eq. (27).Figure 7 shows J ϕ ( ρ ) for N = 70 and 0 (cid:54) η (cid:54) .
5. Themaximum of the current increases as we go from the BCSto the BEC limit, its value at the BEC side, η = 0, beingmore than twice J max at the BCS side, η = 1 .
5. Theposition of the maximum is between ξ = 1.7 and 1.8 k − F at the BCS side of the crossover, i.e., 0 . (cid:54) η (cid:54) .
5; atthe BEC side, η = 0 .
25 and 0.5, ξ ∼ . k − F . The case η = 0 moves away from the trend of a smaller core as wego from the BCS to the BEC limit, with ξ = 2 . k − F . It isunclear if ξ or J max depend on the disk radius R , becausethe R values are closely spaced for 58 (cid:54) N (cid:54)
70, and nosignificant difference was observed in the maximum as wevaried N . Nevertheless, the relative results contribute tounderstanding how the vortex core evolves over the BEC-BCS crossover. J ϕ [ un i t s o f - h / m ] ρ [units of k F-1 ] η =1.50 η =1.25 η =1.00 η =0.75 η =0.50 η =0.25 η =0.00 FIG. 7. (Color online) Angular component of the probabilitycurrent J ϕ as a function of the radial coordinate ρ for severalinteraction strengths η . The position of its maximum providesan estimate of the vortex core size. The wave function that we employed for the vortexstate is an eigenstate of the total angular momentum op-erator. Since this operator commutes with the Hamilto-nian, the diffusion procedure does not change the eigen-value of the state. In addition, the calculation of theprobability current density operator allowed us to ver-ify that the vortex corresponds to a N (cid:126) / L = m (cid:90) ( r × J ) d r , (33)and the component of interest is L z = 2 πm (cid:90) ρ J ϕ ( ρ ) dρ. (34)In our definition of the probability current density opera-tor, we divide by the number of particles N , see Eq. (32).Thus, the evaluation of L z using Eq. (34) should yield (cid:126) /
2. We verified that, for all interaction strengths, thisis in agreement with our simulations.
IV. SUMMARY
We have investigated several properties of vortices in2D Fermi gases over the BEC-BCS crossover region. Wededicated a considerable portion of this work to carefullyunderstand and control size effects in the disk geometry,since it is very convenient to simulating a single vortex.Given that we were interested in the evolution of theproperties in the BEC-BCS crossover, determining thecrossover region was important to verify that the inter-action strengths studied in this work span the crossover.The vortex excitation energies and the density profilesare quantities that can be compared with experiments,once they become available. Interestingly, the observeddensity depletion of the vortex core goes from ≈ η = 1 .
5, to an empty core for η (cid:54) .
25, at the BEC limit. In 3D, Bogoliubov-de Gennestheory has been used to calculate the density suppressionat the vortex core throughout the BEC-BCS crossover[18–20]. Similar calculations in 2D could be comparedto our findings [46]. Also, determining the probabilitycurrent was essential to investigate the changes in thevortex core throughout the crossover region.In 3D the interplay between experiments, theory, andsimulations led to rapid advances in our comprehensionof cold Fermi gases. Hopefully, our results will motivateexperiments to increase our understanding of vortices in2D Fermi gases.
ACKNOWLEDGMENTS
We would like to thank G. C. Strinati for useful dis-cussions. This work was supported by the National Science Foundation under grant PHY-1404405. Thiswork used the Extreme Science and Engineering Discov-ery Environment (XSEDE), which is supported by Na-tional Science Foundation grant number ACI-1053575.The work of S.G. was supported by the NUCLEI Sci-DAC program, the U.S. DOE under Contract No. DE-AC52-06NA25396, and the LANL LDRD program. Wealso used resources provided by NERSC, which is sup-ported by the U.S. DOE under Contract No. DE-AC02-05CH11231. [1] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod.Phys. , 1215 (2008).[2] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[3] M. Randeria, J.-M. Duan, and L.-Y. Shieh, Phys. Rev.Lett. , 981 (1989).[4] M. Randeria, J.-M. Duan, and L.-Y. Shieh, Phys. Rev.B , 327 (1990).[5] D. S. Petrov, M. A. Baranov, and G. V. Shlyapnikov,Phys. Rev. A , 031601 (2003).[6] J.-P. Martikainen and P. T¨orm¨a, Phys. Rev. Lett. ,170407 (2005).[7] J. Tempere, M. Wouters, and J. T. Devreese, Phys. Rev.B , 184526 (2007).[8] W. Zhang, G.-D. Lin, and L.-M. Duan, Phys. Rev. A , 043617 (2008).[9] K. Martiyanov, V. Makhalov, and A. Turlapov, Phys.Rev. Lett. , 030404 (2010).[10] A. A. Orel, P. Dyke, M. Delehaye, C. J. Vale, and H. Hu,New Journal of Physics , 113032 (2011).[11] V. Makhalov, K. Martiyanov, and A. Turlapov, Phys.Rev. Lett. , 045301 (2014).[12] G. Bertaina and S. Giorgini, Phys. Rev. Lett. ,110403 (2011).[13] A. Galea, H. Dawkins, S. Gandolfi, and A. Gezerlis,Phys. Rev. A , 023602 (2016).[14] H. Shi, S. Chiesa, and S. Zhang, Phys. Rev. A , 033603(2015).[15] E. R. Anderson and J. E. Drut, Phys. Rev. Lett. ,115301 (2015).[16] L. Rammelm¨uller, W. J. Porter, and J. E. Drut, Phys.Rev. A , 033639 (2016).[17] Z. Luo, C. E. Berger, and J. E. Drut, Phys. Rev. A ,033604 (2016).[18] A. Bulgac and Y. Yu, Phys. Rev. Lett. , 190404 (2003).[19] R. Sensarma, M. Randeria, and T.-L. Ho, Phys. Rev.Lett. , 090403 (2006). [20] S. Simonucci, P. Pieri, and G. C. Strinati, Phys. Rev. B , 214507 (2013).[21] L. Madeira, S. A. Vitiello, S. Gandolfi, and K. E.Schmidt, Phys. Rev. A , 043604 (2016).[22] M. W. Zwierlein, J. R. Abo-Shaeer, A. Schirotzek, C. H.Schunck, and W. Ketterle, Nature , 1047 (2005).[23] V. L. Berezinsky, Sov. Phys. JETP , 493 (1971).[24] J. M. Kosterlitz and D. J. Thouless, J. Phys. C , L124(1972).[25] S. S. Botelho and C. A. R. S´a de Melo, Phys. Rev. Lett. , 040404 (2006).[26] J. Carlson, S. Gandolfi, and A. Gezerlis, Fifty Years ofNuclear BCS , edited by R. A. Broglia and V. Zelevinsky(World Scientific Publishing Company, 2013).[27] A. Gezerlis and J. Carlson, Phys. Rev. C , 032801(R)(2008).[28] A. Gezerlis and J. Carlson, Phys. Rev. C , 025803(2010).[29] F. V. De Blasio and O. Elgarøy, Phys. Rev. Lett. ,1815 (1999).[30] Y. Yu and A. Bulgac, Phys. Rev. Lett. , 161101 (2003).[31] D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys.Rev. Lett. , 2066 (1983).[32] M. Sadd, G. V. Chester, and L. Reatto, Phys. Rev. Lett. , 2490 (1997).[33] G. Ortiz and D. M. Ceperley, Phys. Rev. Lett. , 4642(1995).[34] N. N. Khuri, A. Martin, J.-M. Richard, and T. T.Wu, Journal of Mathematical Physics , 072105 (2009),http://dx.doi.org/10.1063/1.3167803.[35] S. K. Adhikari, W. G. Gibson, and T. K. Lim, TheJournal of Chemical Physics , 5580 (1986).[36] J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E.Schmidt, Phys. Rev. Lett. , 050401 (2003).[37] A. Gezerlis, S. Gandolfi, K. E. Schmidt, and J. Carlson,Phys. Rev. Lett. , 060403 (2009).[38] S. Gandolfi, Journal of Physics: Conference Series , , 553 (1988).[40] S. Gandolfi, A. Y. Illarionov, F. Pederiva, K. E. Schmidt,and S. Fantoni, Phys. Rev. C , 045802 (2009).[41] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra-jagopal, Rev. Mod. Phys. , 33 (2001).[42] See Supplemental Material at [ url to be inserted by pub-lisher ] for the variational parameters of Eqs. (17), (20),and (22).[43] M. Casula, C. Attaccalite, and S. Sorella, The Journalof Chemical Physics , 7110 (2004).[44] D. Ceperley and H. Kalos, Monte Carlo Methods inStatistics Physics Quantum Many-Body Problems , edited by K. E. Binder, Vol. 7 (Springer-Verlag, 1986).[45] P.-M. K¨onig, R. Roth, and K. R. Mecke, Phys. Rev.Lett. , 160601 (2004).[46] After our calculations were done, it was pointed out tous that pseudogap phenomena occurring in 2D and 3DFermi gases can be related in a universal way througha variable that spans the BEC-BCS crossover [47]. Fur-ther studies are necessary to determine if this universal-ity holds for other quantities, such as the density andthe probability current density per particle. This wouldprovide a very clean way of comparing 2D and 3D results.[47] F. Marsiglio, P. Pieri, A. Perali, F. Palestini, and G. C.Strinati, Phys. Rev. B91