Dipolar fermions in a multilayer geometry M. Callegari,1,2 M. M. Parish,3,4,∗ and F. M. Marchetti2,† 1Department of Physics and Astronomy “Galileo Galilei”, University of Padova, Via Marzolo 8, 35131 Padova, Italy 2Departamento de F´ısica Teo´rica de la Materia Condensada & Condensed Matter Physics Center (IFIMAC), Universidad Auto´noma de Madrid, Madrid 28049, Spain 3London Centre for Nanotechnology, Gordon Street, London, WC1H 0AH, United Kingdom 4School of Physics & Astronomy, Monash University, Victoria 3800, Australia (Dated: January 26, 2016) 6 We investigate the behavior of identical dipolar fermions with aligned dipole moments in two- 1 dimensional multilayers at zero temperature. We consider density instabilities that are driven by 0 theattractivepartofthedipolarinteractionand,forthecaseofbilayers,weelucidatetheproperties 2 ofthestripephaserecentlypredictedtoexistinthisinteractionregime. Whenthenumberoflayers n is increased, we find that this “attractive” stripe phase exists for an increasingly larger range of a dipole angles, and if the interlayer distance is sufficiently small, the stripe phase eventually spans J thefullrangeofangles,includingthesituationwherethedipolemomentsarealignedperpendicular 6 to the planes. In the limit of an infinite number of layers, we derive an analytic expression for the 2 interlayer effects in the density-density response function and, using this result, we find that the stripe phase is replaced by a collapse of the dipolar system. ] s PACSnumbers: 67.85.Lm,68.65.Ac,71.45.Gm a g - t I. INTRODUCTION n a u Motivated by the prospect of novel many-body phases q generated by anisotropic long-range dipolar interactions, . t much attention has recently been devoted to ultracold a polar molecules and magnetic atoms [1–3]. Dipolar m fermions, in particular, can be used to simulate strongly - correlated phenomena in electron systems, including d n chargedensitymodulations(stripes)andunconventional o superconductivity [4–6]. c Quantum degeneracy has already been achieved for [ dipolarFermigasesofatomswithapermanentmagnetic 1 dipolemomentsuchaschromium[7],dysprosium[8]and FIG.1. (Coloronline)Schematicrepresentationofthesystem v erbium [9]. This has enabled the observation of dipole- geometry. A gas of dipolar fermions is confined to N 2D 8 driven Fermi surface deformations in the Fermi liquid layerslabelledbyanindexj =1,2,...,N andseparatedbya 8 phase [10]. However, to investigate many-body phenom- distance d. We assume that each layer has the same density 9 enaatstrongerdipole-dipoleinteractions,itappearsnec- of dipoles n. The dipole moment of each fermion is aligned 6 essary to use polar molecules, which generally possess by an electric field E = (E ,0,E ) in the x-z plane (ϕ = 0 0 x z largerdipolemoments—theelectricdipolemomentcan direction), at an angle θ with respect to the zˆ direction. . 1 be as large as 5.5 Debye in the case of 133Cs6Li [11]. 0 Thusfar,therehasbeenmajorprogresstowardsproduc- 6 ing quantum degenerate clouds of long-lived fermionic 1 strong electric field, the nature of the effective 2D dipo- : dipolar molecules using 40K87Rb [12–14], 23Na6Li [15], lar interaction within the plane may be externally ma- v 133Cs6Li [16], and 23Na40K [17, 18]. i nipulated: The dipole-dipole repulsion is maximized by X The dipole-dipole interaction can be further tuned aligning the dipole moments perpendicular to the plane, and enhanced by confining the polar molecules to two- r whileanisotropyandattractionaregraduallyintroduced a dimensional (2D)layers. Suchageometryhasbeenused byvaryingthedipoletilt(seeFig.1). Thispossibilityhas to suppress chemical reactions [19] and to stabilize the stimulated much theoretical work on dipolar fermionic gas against mechanical collapse, which arises in three gases in single and multilayer 2D geometries [3]. dimensions for a sufficiently strong interaction [20–26]. For the single-layer geometry and for small (but non- Furthermore, by aligning all the dipole moments with a zero) tilting angles, the weakly interacting system corre- sponds to a Landau Fermi liquid with deformed Fermi surface [27], similarly to the case in 3D. With increasing ∗ [email protected] dipolar interaction (or cloud density), the system is then † [email protected] predictedtoundergoatransitiontoauni-directionalden- 2 sitymodulatedphase[28–31],wherethemodulationsare II. MULTILAYER SYSTEM AND MODEL perpendicular to the direction of the dipole tilt. Such a “stripe” phase has also been shown to exist in the We consider a gas of polar fermionic molecules in a isotropic case where the dipoles are aligned perpendic- multilayer geometry, as shown in the schematic picture ular to the layer, thus requiring the system to sponta- in Fig. 1. The molecules have a dipole moment D and neously break the rotational symmetry [32]. This result are confined to N two-dimensional layers, each labelled hasrecentlybeensupportedbydensityfunctionaltheory by an index j = 1,2,...,N, and equally separated by calculations [33], which predict a transition to a stripe a distance d. We assume that the dipoles are aligned phase followed by a transition to a triangular Wigner by an external electric field E = (E ,0,E ) in the x-z x z crystal at higher coupling. Quantum Monte Carlo cal- plane, which is tilted at an angle θ with respect to the zˆ culations also find a Wigner crystal phase for the case direction. Withineachlayer,weparameterizethex-y in- of perpendicularly aligned dipoles, although at a much plane wavevector by polar coordinates q=(q,ϕ), where higher dipolar interaction than that obtained from den- ϕ=0 corresponds to the direction of the dipole tilt. sity functional theory [34]. In the limit qW (cid:28)1, where W is the layer width, the effective 2D intralayer interaction between dipoles takes For tilting angles greater than a critical angle, the at- the following form [43]: tractive part of the dipolar interaction can lead to p- wave superfluidity in the single-layer system [35], and vjj(q)≡v(q)=V0−2πD2qξ(θ,ϕ), (1) this phase may even coexist with stripe order [36]. A sufficientlystrongattractioneventuallydrivesamechan- whereξ(θ,ϕ)=cos2θ−sin2θcos2ϕ. Here,qcorresponds ical instability of the cloud towards collapse [28, 31, 32, totherelativewavevectorbetweentwodipoles. Thecon- 35,37]. However, interestingly, ifoneinsteadconsidersa stantV0 isashort-rangecontactinteractiontermthatin bilayer geometry, the additional layer stabilizes the col- general depends on the width W [43], yielding a natural lapse at large tilt angles — as long as the dipoles are UV cut-off. Since we are considering identical fermions, aligned out of the plane (θ < π/2) — to form a new the system properties will not depend on V0. stripephase, wherethedensitymodulationsareoriented Inthelimitwherethelayerwidthismuchsmallerthan along the direction of the dipole tilt [38]. thelayerseparation,W (cid:28)d,theinteractionbetweentwo dipoles in different layers j >l is given by [44]: Inthisarticle,weinvestigatesuchastripephase,which v (q)=−2πD2qe−(j−l)qd[ξ(θ,ϕ)+isin2θcosϕ] . (2) isgeneratedbytheattractivepartofthedipolarinterac- jl tionforlargeenoughdipoletiltangle. Westartbyeluci- The remaining interlayer interactions can be obtained dating its properties in the case of the bilayer, where we from the condition v (q) = v (−q) = v∗(q), which is provide a classical argument for how the stripes in each lj jl jl derived from the fact the the dipolar interaction is al- layer are shifted with respect to each other. Then we ways real in real space. Likewise, the momentum-space extend our results for the density-density response func- interactioniscomplexforθ (cid:54)=0sincethereal-spaceinter- tiontothemultilayergeometry. Weemployanapproach actionisnotinvariantunderthetransformationr(cid:55)→−r. based on the STLS scheme [32, 38, 39] that incorporates Assuming that each layer has the s√ame density n, exchange interactions only, which should be reasonable we define the Fermi wavevector k = 4πn. This al- for the “attractive” stripe phase [38], although note that F lows us to define the dimensionless interaction strength we neglect the possibility of pairing [40, 41] and other U = mD2k /(cid:126)2, where m is the fermion mass. The stripe phases driven by the strong repulsion [38, 42]. As F other parameters that can be independently varied are the number of layers N is increased, we find that the the dipole tilt angle θ and the dimensionless layer sepa- attractive stripe phase spans an increasingly larger re- ration k d. gion of the phase diagram. However, this stripe phase F eventually gives way to collapse in the N →∞ limit. A. Response function and STLS equations The paper is organized as follows: In Sec. II, we de- scribe the system geometry and introduce the STLS Similarly to Ref. [38], we make use of linear response scheme which allows us to evaluate the density-density theory to analyze density wave instabilities. In the mul- response function matrix in the multilayer geometry; in tilayersystem,thelineardensityresponseδntoanexter- Sec.III,wedescribethepropertiesofthedensityinstabil- nal perturbing potential Vext defines the density-density itiesdrivenbytheattractivepartoftheinterlayerdipolar response function matrix [45], interaction and, in Sec. IIIA, we explain via a classical model how the stripes in each layer are shifted with re- δn (q,ω)=(cid:88)χ (q,ω)Vext(q,ω), (3) specttoeachother. InSec.IVweextendtheresultstoa j jl l l generic number of layers N, while, in Sec. V we consider the N →∞ limit. The concluding remarks are gathered wherej,larethelayerindices. Adivergenceinthestatic in Sec. VI. density-density response function matrix χ (q,ω = 0) jl 3 signals an instability of the system. Specifically, the sys- andisnotjustlimitedtoexchangecorrelations. Assuch, tem is unstable towards forming a stripe phase when theSTLSschemehasproventobeapowerfulmethodfor the smallest eigenvalue χ˜ (q) of the static response treating strongly correlated electron systems such as the min function matrix first diverges at a critical value of the 2Delectrongas[46]. Wepreviouslyadaptedanimproved wavevector q =(q ,ϕ ). version of this scheme to the dipolar system, both in c c c While the response function is known exactly for the the single- [32] and double-layer [38] geometries. The non-interacting gas, typically one can only incorporate scheme is improved by imposing, at each iteration step, the effect of interactions approximately. A standard ap- theconditionthattheintralayerpaircorrelationfunction proachistherandomphaseapproximation(RPA),where is zero at zero distance, g (0)=0, where, jj one replaces the external potential with one that con- tains an effective potential due to the perturbed den- 1 (cid:90) dq sity: Vjext (cid:55)→ Vjext + (cid:80)lvjlδnl. However, RPA ne- gjl(r)=1+ n (2π)2eiq·r[Sjl(q)−δjl] . (8) glects exchange correlations, which are always impor- tant in the dipolar system, even in the long-wavelength This ensures that the intralayer static structure factor limit q → 0 [32]. This issue may be remedied using S (q) is dominated by Pauli exclusion in the long wave- jj a conserving Hartree-Fock approximation [30, 31, 42], length limit, q (cid:29) 2k , and that the system response is F but we choose a simpler and physically motivated ap- independent of the short-range contact interaction term proachwherecorrelationsareincludedvialocalfieldfac- V and the cut-off W. 0 tors G (q) [46]. This yields the inverse density-density jl In the following, we first review the bilayer case and response function matrix describe the instability to a stripe phase occurring for δ large tilt angles θ, where the modulations are oriented χˆ−1jl(q,ω)= Π(qj,lω) −vjl(q)[1−Gjl(q)] , (4) along the dipole tilt, i.e., along ϕ = 0. We then show how the instability to the ϕ=0 stripe phase can be well where Π(q,ω) is the non-interacting response function, described using exchange correlations only, and we use which, for equal density layers, reads as [45, 47] this to investigate its existence in the multilayer geome- try. m (cid:40)√ (cid:20) (cid:113) (cid:21)1/2 (cid:41) Π(q,iω)= 2 a+ a2+(ωb)2 −b , 2π(cid:126)2b III. THE ϕ=0 STRIPE PHASE IN BILAYERS with a = b2 −bkF2 −ω2 and b = q2. RPA corresponds 4 m m to taking the limit where the layer-resolved local field Thecaseoftwolayers(N =2)waspreviouslyanalysed factors G (q) in Eq. (4) are all zero. withintheSTLSself-consistentapproximationschemein jl The response function (4) can be related to the layer- Ref. [38]. Here, at sufficiently small tilt angles θ < θ , c resolvedstaticstructurefactorS (q)bythefluctuation- and by increasing the value of the dimensionless cou- jl dissipation theorem: pling strength U, there is an instability from the uni- form phase to a stripe phase with modulations along the (cid:126) (cid:90) ∞ S (q)=− dωχ (q,iω). (5) y-axis (ϕ = π/2 stripe phase). The instability to this jl jl πn 0 stripe phase is driven by intralayer correlations beyond exchange, which are induced by the repulsive part of the Inthenon-interactinglimit, thestaticstructurefactoris intralayer interaction potential v(q). By contrast, for diagonal,i.e.,Sj(0l)(q)=δjl S(0)(q),andcanbeevaluated θc < θ < π/2, the system develops an instability to a exactly (see App. A), where stripephasealongthex-axis(ϕ=0stripephase). While forasinglelayer,theattractivesliveroftheintralayerin- (cid:115) 2 (cid:18) q (cid:19) q (cid:18) q (cid:19)2 teraction produces a collapse of the dipolar Fermi gas at S(0)(q)= arcsin + 1− , (6) large tilt angles, the bilayer geometry stabilizes the col- π 2k πk 2k F F F lapseinfavourofaϕ=0stripephase,whichthusderives for q ≤2k , while S(0)(q)=1 for q >2k [45, 46]. from a competition between the intralayer attraction in F F the ϕ=0 direction and the interlayer interaction. To determine the local field factors, we consider the STLS approximation scheme, where we have the expres- Interestingly,thislatterstripephasecanbeaccurately sion [39, 45]: described using intralayer exchange correlations only. In fact, it was found for this phase that the intralayer pair 1 (cid:90) dk q·kv (k) correlation function g (r) deviated only slightly from G (q)= jl [δ −S (q−k)] . (7) jj jl n (2π)2 q2 v (q) jl jl the non-interacting case, while the interlayer correla- jj tion function g (r) ∼ 1. In terms of local correlations, 12 In principle, G (q) can be determined self-consistently the interlayer local field factor can thus be neglected, jl by solving Eqs. (4), (5) and (7). Note that this self- G (r) = 0, while the intralayer one G (q) ≡ G(q) is 12 jj consistentapproachincludesallcorrelationsbeyondRPA determined from the non-interacting intralayer structure 4 FIG. 3. (Color online) Schematic representation of a single dipole in the top layer interacting classically with an infinite bottomlayerofdipolesarrangedinastripemodulatedphase. As in Fig. 1, all dipoles are aligned by an electric field E. In theleft(right)panel,thesingledipoleisshiftedbyx¯(y¯)with respecttooneofthestripecrestsfortheϕ=0(ϕ= π)stripe 2 phase. FIG. 2. (Color online) Main panel: Phase boundary for the tion (4) determine the phase-shift η between layers, and ϕ = 0 stripe phase in the two-layer geometry as a function for the bilayer geometry we have found that [38]: of the tilt angle θ and the dimensionless interaction strength U at fixed interlayer distance kFd=2. Density modulations eiη =− v12(q)[1−G12(q)] . (10) areinthedirectionϕ=0ofthedipoletilt(schematicfigure). |v (q)[1−G (q)]| 12 12 The instability to this phase is to the right of the plotted boundaries: The full STLS results [38] ([green] open trian- For both ϕ = π/2 and ϕ = 0 stripe phases, we find gles)arecomparedtotheresultsobtainedwithexchange-only that the interlayer phase shift between the modulations correlations([blue]solidline),alsoreferredtoastheX-STLS is independent of the dipole interaction strength U and approximation scheme. Within X-STLS, the ϕ = 0 stripe thelayerdistanced. Inparticular,fortheϕ=π/2stripe phaseappearsforθ>θ (cid:39)0.79. Atθ=π/2,thegasisunsta- c phase, both the interaction and the local field factor are ble towards mechanical collapse for U (cid:38)1.57 ([red] diamond realandthusbothlayermodulationsarealwaysinphase, symbol and thick solid line), where the gas compressibility is i.e.,η =0. Ontheotherhand,fortheϕ=0stripephase, infinite. The density modulations in the two layers have a ifweconsideranexchange-onlyapproximation(X-STLS) phase shift η (lower inset) equal to 2θ, i.e., the shift between the modulations 2θ/qc (schematic figure). for which G12(q)=0, we obtain a phase shift of η =2θ. The phase shift for the ϕ = 0 stripe phase is plotted in the inset of Fig. 2. We see a very good agreement factor (6): between the full STLS results ([green] open triangles) andthesimplifiedX-STLSscheme([blue]solidline). The 1 (cid:90) dk q·kv(k)(cid:104) (cid:105) ϕ=0 result may at first appear counter-intuitive, but it G(q)= 1−S(0)(q) . (9) n (2π)2 q2 v(q) canbereproducedbyevaluatingtheclassicalinteraction energy between an infinite layer of dipoles in one layer We refer to this approximation scheme as the exchange- and a single dipole in the second layer, as we discuss onlySTLSapproximation(X-STLS).Forthesingle-layer next. case [32], this approach yielded an instability towards collapse at large θ that agreed with the predictions from Hartree-Fock calculations [28, 30, 31, 35]. A. Classical model The phase boundary between the normal phase and the ϕ = 0 stripe phase obtained from the full STLS For the bilayer geometry, a simple classical model can schemeisdisplayedinFig.2([green]opentriangles)and easily explain the phase shifts found in both ϕ = 0 and is compared with the results of the X-STLS approxima- ϕ=π/2stripephases. Letusconsiderthesimplifiedcase tion([blue]solidline). Thevalueofthecriticaltiltangle ofan infinitelayerof dipoleswhose densityis modulated for such a phase is θc (cid:39) 0.75 if evaluated within the full sinusoidallywithawavevectorq =2π/λ andamplitude c c STLS scheme, while it is slightly higher, θc (cid:39) 0.79, if ρ0. Wefurtherassumethatthisinteractsclassicallywith evaluated within the X-STLS approximation. It turns a single dipole positioned at r in the other layer. The 0 out that, for the phase boundary, the X-STLS approx- two layer densities are thus respectively given by: imation works particularly well at large angles, all the way up to θ = π/2, where, for U (cid:38) 1.57 ([red] diamond ρ(1)(r)=n+ρ0cos(qc·r) (11) symbol and thick solid line), the gas collapses because ρ(2)(r(cid:48))=δ(r(cid:48)−r ). (12) 0 the Fermi pressure is not high enough to counteract the strong dipolar attraction. For the ϕ = 0 (ϕ = π/2) stripe phase we have qˆ = c The eigenvectors of the density-density response func- xˆ (qˆ = yˆ) and r = (x¯,0) (r = (0,y¯)) — see the c 0 0 5 schematic representation of both geometries in Fig. 3. The classical interaction energy is given by (cid:90) E = drdr(cid:48)ρ(1)(r)v (r−r(cid:48))ρ(2)(r(cid:48)), (13) cl 12 where x2+y2+d2−3(xsinθ+dcosθ)2 v (r)=D2 12 (x2+y2+d2)5/2 istheinterlayerpotential(assumingthatthedistancedis much larger than the layer thickness W). For a uniform density distribution ρ(1)(r)=n, the classical interaction FIG.4. (Coloronline)Evolutionoftheϕ=0stripephasefor energy would be zero; therefore only deviations from the the trilayer N =3 geometry when varying the dimensionless average density in ρ(1)(r) contribute (either positively or layer separation dkF. Left panel: phase boundary in the θ versus U plane; the instability to the stripe phase is on the negatively) to E . cl right side of the plotted boundary. Right panel: Rescaled Considering the Fourier transforms ρ(1,2)(r) = stripewavevectorq /k atthephaseboundariesasafunction (cid:82) (2dπq)2ρ(1,2)(q)eiq·r and v12(r) = (cid:82) (2dπq)2v12(q)eiq·r, we of the tilt angle θ.c F can rewrite (13) to obtain (cid:90) dq the X-STLS approximation, the critical tilt angle for the E = ρ(1)(−q)v (q)ρ(2)(q), (14) cl (2π)2 12 ϕ = 0 stripe phase in bilayers is θ (N = 2) ∼ 0.79, c which is lower than that for collapse in the single layer, where θ (N =1)∼0.89. It is therefore natural to ask whether c (cid:20) (cid:21) the ϕ = 0 stripe phase will tend to dominate the phase δ(q+q )+δ(q−q ) ρ(1)(q)=(2π)2 nδ(q)+ρ0 c 2 c diagram as the number of layers N is increased. ρ(2)(q)=e−iq·r0 . IV. N LAYERS Thus, as v (0)=0, we get in general 12 Ecl = ρ20 (cid:2)v12(−qc)eiqc·r0 +v12(qc)e−iqc·r0(cid:3) , (15) temM,owtievantoewdabpyptlyhethreesXu-lStsToLbStaaipnperdoxfoimr athtieonbislcahyeemr seytso- the general case of finite N > 2 layers and evaluate the and specifically for the two stripe phases: occurrence of the ϕ = 0 stripe phase when varying the system parameters. In particular, by neglecting all cor- 4π2D2ρ (cid:18) x¯ (cid:19) Eϕ=0 =− 0e−2πd/λccos 2π −2θ relations except for the exchange ones, we assume that cl λc λc all off-diagonal local field factors are zero, Gj(cid:54)=l(q) = 0, 4π2D2ρ (cid:18) y¯ (cid:19) while the intralayer ones G(q) are evaluated according Ecϕl=π/2 =− λ 0e−2πd/λccos2θcos 2πλ . to Eq. (9). To locate the stripe instabilities, we extract c c the smallest eigenvalue of the static density-density re- Therefore, we can conclude that, for both stripe configu- sponsefunctionmatrix, χ˜ (q), anddeterminethecrit- min rations, the distance, x¯ or y¯, that minimizes the interac- ical wavevector q =(q ,ϕ ) at which it first diverges. If c c c tion energy E does not depend on the dipole strength the instability is for a specific angle ϕ , then it signals cl c D or the layer separation d. Further, for the ϕ = π/2 theformationofadensitywavewithmodulationsinthat stripe phase, the best configuration is the one where the direction and with a period set by q . Here, we always c single dipole in layer 2 aligns with the maximum den- find that ϕ =0, as in the bilayer case. c sity of layer 1, i.e., y¯ = 0. By contrast, for the At the stripe transition, the phase shifts η between min jl ϕ = 0 stripe phase, the optimal configuration is for a the stripes in different layers j,l are extracted from phase shift equal to twice the dipole tilt angle θ, i.e., the eigenvector associated with the smallest eigenvalue, 2πx¯ /λ = 2θ. An analogous calculation was carried χ˜ (q ). Wefindthatthebehaviorofthemultilayersys- min c min c out for the ϕ = π/2 stripe phase in the simplified limit tem is a natural extension of the bilayer case: the phase wherethedensitymodulationswereapproximatedasdis- shiftbetweennearestneighbourlayersη growsmono- jj+1 crete dipolar wires [42]. tonically with the tilt angle θ, although linearly only for We now wish to extend these results to multiple layer smallvaluesofθ. Moreover,thephaseshiftbetweenmore N > 2 configurations. We have seen that the presence distant layers is always proportional to η . jj+1 of a second layer stabilises the region of collapse at large To gain insight into the phase diagram of the multi- tilt angles θ < θ < π/2, replacing it with a novel stripe layer system, we first focus on the trilayer N = 3. In c phase oriented along ϕ = 0 [38]. Furthermore, within Fig. 4, we plot the phase boundaries for the instability 6 when the layer distance decreases, we find that θ even- c tually reaches zero at a critical distance d ([violet] star c symbols): Here, the ϕ = 0 stripe phase spans the en- tire range of tilt angles θ. For smaller distances, d<d , c stripe formation is always possible for sufficiently large but finite values of the interaction strength U, even for dipoles aligned perpendicular to the planes. Notethatwhendecreasingthevalueofdk ,eventually F ourexchange-onlyformalismbecomesquestionable,since it neglects interlayer correlations. In fact, for k d (cid:46) 1, F interlayer pairing (e.g., dimers in the two-layer configu- ration [48] and bound chains in multilayers [49]) is ex- pectedtodominateoverstripeformation,andthisisnot included in our approximation scheme. We then investigate whether, by fixing the layer dis- tance to a value k d > 1, the ϕ = 0 stripe phase can F dominatethephasediagramasthenumberoflayersN is increased. Tothisend,weplotinFig.5thephasebound- FIG. 5. (Color online) Phase boundaries (left panels) and aries (left panels) for the ϕ=0 stripe phase for different rescaledstripewavevectorattheboundariesq /k (rightpan- c F valuesofN. Weobservequalitativelydifferentbehaviour els)fortheϕ=0stripephasefordifferentvaluesofthenum- depending on whether the distance d is above or below ber of layers N and two fixed values of the layer distance: ∼ 1.47/k , corresponding to the critical distance d for dk =1.2 (top panels) and dk =2 (bottom panels). F c F F N → ∞, as derived in Sec. V. When d > d (N → ∞) c (lower panels of Fig. 5), the ϕ = 0 stripe phase exists for an increasingly larger range of dipole tilt angles as N increases, but the critical angle θ finally saturates to a c finitepositivevalue. Wheninsteadd<d (N →∞)(up- c per panels of Fig. 5), the stripe phase eventually spans thefullrangeofangles,includingthesituationwherethe dipole moments are aligned perpendicular to the planes. TheseresultsaresummarizedinFig.6,whereweplot, as a function of 1/N, the critical value of the interlayer distanced k atwhichtheϕ=0stripephasefirstspans c F the full range of dipole tilt angles. The data for N → ∞ in the figures are evaluated following the procedure explained in the next section. V. THE N →∞ LIMIT FIG.6. (Coloronline)Criticalvalueoftheinterlayerdistance We now show how the calculation for the ϕ=0 stripe d k atwhichθ =0(atU →∞)asafunctionoftheinverse instabilitycanbeextendedtothelimitofaninfinitenum- c F c number of layers, 1/N. In the limit N → ∞, dckF = 1.47. beroflayers. Thekeypointisthattheinterlayerinterac- Numerical data are represented as (red) circles, while the tionpotentialv (q)onlydependsonthelayerindexdif- jl dashed(blue)lineisalinearfittothedata. Inset: maximum ference|j−l|, sothat, forasystemwithperiodicbound- value of the stripe wavevector q /k for d=d as a function c F c aryconditionsandN (cid:29)1,wecanmakeatransformation of 1/N. Data are (orange) triangles, while the (turquoise) fromthelayerindexspacej =1,2,...,N totherecipro- dashedlineisanon-linearfitgivingf(1/N)=2.83(1/N)0.49. cal space p=2πm/N, where m=−N/2,...,N/2 [50]: N (cid:88) to a ϕ = 0 stripe phase (left panel) and the associated u˜ = e−i(j−l)pv . (16) p jl critical wavevectors qc/kF (right panel) for different val- j−l=−N ues of the layer distance dk . The qualitative behaviour F Of course, in the actual system, we do not have peri- we observe here is common to any number of layers N, odic boundary conditions, but this should not change including the case of a bilayer N = 2. For large enough the physics in the limit N →∞. Inserting Eq. (2) yields layer distance d, the ϕ = 0 stripe phase can only occur the analytical solution for values of the tilt angle greater than a critical an- (cid:20) (cid:21) gle θ , i.e., at the stripe instability boundary, we have ζ(θ,ϕ) ζ(θ,ϕ+π) c u˜ (q)=−2πD2q + , (17) U → ∞ for θ → θc ([green] circle symbols). However, p eip+qd−1 e−ip+qd−1 7 the phase boundary of such a collapsed phase are sum- marized in Fig. 7 and are qualitatively similar to the boundaries found for the ϕ = 0 stripe phase in the case of a finite number of layers. Below the critical layer dis- tance d k (cid:39) 1.47, the collapsed phase exists for any c F value of the dipole tilt angle, including θ =0. The behaviour of the infinite N multilayer system is reminiscentofthatexpectedforthe3DdipolarFermigas, where one also has collapse for sufficiently large dipolar interactions [21]. In particular, the θ = 0 case possesses thesamerotationalsymmetryasthe3Dgaswithaligned dipole moments. Therefore, it is instructive to compare FIG. 7. (Color online) Left panel: Phase boundary of the the onset of collapse for θ = 0 in the multilayer system collapsed region for an infinite number of layers N →∞ and tothatinthe3Dcase. First,onecandefinea3Ddensity different values of the rescaled layer density dkF. The insta- for the multilayer geometry: bility to collapse is signalled by an infinite compressibility of the dipolar gas, i.e., by a divergence of the static response n k3 function for q = 0. While for d > d , collapse occurs only n = = F . (21) c c 3D d 4π(k d) for a tilt angle larger than a critical value θ , for d<d , the F c c collapsedphasespanstheentirerangeofangles. Rightpanel: This then yields the corresponding 3D dimensionless in- We plot the asymptotic values of the tilt angle at the phase boundary for U → ∞, θ , as a function of the rescaled layer teraction parameter: c distance dk ; we find that the critical distance for infinite F layers is given by dc(N →∞)kF =1.47. U ≡ mD2n1/3 = U . (22) 3D (cid:126)2 3D (4πk d)1/3 F whereζ(θ,ϕ)=ξ(θ,ϕ)+isin2θcosϕ,andwherewehave In the 3D dipolar gas, the collapse instability occurs at taken the limit N → ∞ after evaluating the geometric U (cid:39) 2.4 [21]. Thus, to obtain a comparable U for 3D 3D series. collapse in the multilayer system, we require dk (cid:39)1.31 F We thus obtain the following expression for the eigen- and U (cid:39)6.12. In our 2D infinite layer configuration, we values of the inverse density-density response function: caninterprettheparameterdk asaneffectiveFermisur- F face“deformation”parameterifwetreat2π/dasaFermi 1 χ˜−1(q,ω)= −v(q)[1−G(q)]−u˜ (q). (18) momentuminthez direction. Thisassumesthatitisthe p Π(q,ω) p inter-particlespacingratherthefermionexchangethatis the key feature of the Fermi surface deformation in 3D To investigate the instabilities, we take the static limit, when considering the collapse instability. For the mul- ω = 0, and determine the values of p and q for which tilayer geometry, the ratio between the Fermi momenta χ˜−1 first hits zero. In practice, this means we must find p in the z and radial directions is then given by 2π/(k d), thevalueofpthatmaximizes−u˜ (q)foreachq. Solving F p thus yielding k /k ∼ 4.8 at the collapse instability. for the stationary points gives us two solutions: F,z F,r This is not so dissimilar from that obtained for the 3D (cid:34)(cid:0) (cid:1) (cid:35) Fermi gas within Hartree-Fock mean-field theory, where eqd∓1 sinθcosϕ p =2arctan +πδ , (19) k /k ∼ 2 [21]. In the layered system, we can effec- i=1,2 (eqd±1)cosθ i,2 F,z F,r tively tune the deformation such that the critical U 3D forcollapseisraised(lowered)byincreasing(decreasing) where the argument is positive if we assume 0 ≤ ϕ ≤ dk . Eventually, when d > d , the Fermi surface is not F c π/2. The first solution p corresponds to the maximum 1 sufficiently elongated along the dipole direction to pro- of −u˜ (q), where we have p duce collapse. (cid:18) cos2θ cos2ϕsin2θ(cid:19) u˜ (q)=−4πD2q − . (20) p1 eqd−1 eqd+1 VI. CONCLUDING REMARKS Thus we have now considerably simplified the problem, as we only have to find the zero of the maximum inverse Inthiswork,wehaveanalysedthedensityinstabilities eigenvalueofthestaticresponse,χ˜−1(q,0),asafunction of dipolar Fermi multilayer systems that are driven by p1 of q. the attractive part of the dipolar interaction. We have In this limit, we always find that (q ,ϕ ) = 0, i.e., by arguedthatsuchinstabilitiesaredominatedbyexchange c c increasing the number of layers to infinity, the gas be- correlations and can thus be described using a simpli- comes unstable towards mechanical collapse. Here, the fied exchange-only STLS approach. We find that the gas compressibility, which is proportional to the static attraction-driven ϕ = 0 stripe phase expands to fill the response function at q = 0, is infinite. The results for phase diagram with an increasing number of layers N. 8 However, at the same time, the stripe wavevector de- where n is the 2D density and n =Θ(k −k) the zero- k F creases so that the stripe phase is eventually replaced by temperature Fermi distribution function. Here, the inte- collapse as N → ∞. For the case θ = 0, the infinite gral simply corresponds to calculating the area A of the N limit resembles a 3D dipolar gas with a Fermi surface overlapregionbetweentwoidenticalcirclesofradiusk , F deformation that can be tuned by varying the interlayer as shown in Fig. 8. Due to the symmetry of the prob- distance dk . lem, we only need to consider half of the overlap region F Our predicted stripe phases should be assessable in as follows. experiments on polar molecules with sufficiently large Assuming q < 2k , we first determine the segment F dipole moments. One also needs to consider the is- area spanned by the angle 2θ in the left circle of Fig. 8: sue of losses in experiments when strong interactions are involved. The restricted motion in the 2D geome- A =k2θ try reduces the possibility of head-to-tail collisions be- seg F tweendipoles, whichunderlies thedominant lossprocess in dipolar gases, but such collisions are not necessar- ily suppressed when we have a large dipole tilt. How- k F ever, we expect that chemically stable molecules such as ✓ 23Na40K [18] will make it possible to probe this regime q of parameter space. In the future, it would be interesting to extend our results to finite temperature, where the proliferation of topological defects can melt the stripe phase [51]. Fur- thermore, one could investigate the effects of pairing using more sophisticated approaches to the multilayer FIG. 8. Two identical circles of radius k , where the cen- F system such as the Euler-Lagrange Fermi-hypernetted- ters are separated by a distance q in k-space. The area of chain approximation [52]. Finally, there is the intriguing the shaded overlapping region corresponds to the integral in question of how our predicted phase diagram connects Eq.(A1),definingthethenon-interactingstaticstructurefac- with other instabilities such as nematic phases or col- tor in 2D. lapse within the stripe phase. where cosθ = q/2k . Next, we determine the area of F ACKNOWLEDGMENTS the left triangle obtained by drawing a line through the points where the circles intersect: We are grateful to G. M. Bruun and P. A. Marchetti (cid:115) for useful discussions. FMM acknowledges financial 1 1 (cid:18) q (cid:19)2 support from the Ministerio de Econom´ıa y Competi- A∆ = 2kFqsinθ = 2kFq 1− 2k tividad (MINECO), projects No. MAT2011-22997 and F No. MAT2014-53119-C2-1-R. MMP acknowledges sup- Then we obtain: port from the EPSRC under Grant No. EP/H00369X/2. A=2(A −A ) seg ∆ Appendix A: Non-interacting static structure factor (cid:115) (cid:18) (cid:19) (cid:18) (cid:19)2 q q =2k2 arccos −k q 1− (A2) We start with the general expression for the non- F 2k F 2k F F interactingstaticstructurefactorintwodimensions[46]: Inserting (A2) into (A1), and using the fact that 1 (cid:90) dk S(0)(q)=1− n n , (A1) arcsinx = π/2−arccosx, we finally recover Eq. (6) in n (2π)2 k k+q the main text. [1] M. A. Baranov, Phys. Rep. 464, 71 (2008). 393, 550 (1998). [2] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and [6] S.A.Kivelson,I.P.Bindloss,E.Fradkin,V.Oganesyan, T. Pfau, Rep. Prog. Phys. 72, 126401 (2009). J. M. Tranquada, A. Kapitulnik, and C. Howald, Rev. [3] M.A.Baranov,M.Dalmonte,G.Pupillo, andP.Zoller, Mod. Phys. 75, 1201 (2003). Chemical Reviews 112, 5012 (2012). [7] B. Naylor, A. Reigue, E. Mar´echal, O. Gorceix, [4] G. Gru¨ner, Rev. Mod. Phys. 60, 1129 (1988). B. Laburthe-Tolra, and L. Vernac, Phys. Rev. A 91, [5] S. A. Kivelson, E. Fradkin, and V. J. Emery, Nature 011603 (2015). 9 [8] M. Lu, N. Q. Burdick, and B. L. Lev, Phys. Rev. Lett. [30] M. Babadi and E. Demler, Phys. Rev. B 84, 235124 108, 215301 (2012). (2011). [9] K.Aikawa,A.Frisch,M.Mark,S.Baier,R.Grimm, and [31] L. M. Sieberer and M. A. Baranov, Phys. Rev. A 84, F. Ferlaino, Phys. Rev. Lett. 112, 010404 (2014). 063633 (2011). [10] K. Aikawa, S. Baier, A. Frisch, M. Mark, C. Ravensber- [32] M.M.ParishandF.M.Marchetti,Phys.Rev.Lett.108, gen, and F. Ferlaino, Science 345, 1484 (2014). 145304 (2012). [11] L.D.Carr,D.DeMille,R.V.Krems, andJ.Ye,NewJ. [33] B.P.vanZyl,W.Kirkby, andW.Ferguson,Phys.Rev. Phys. 11, 055049 (2009). A 92, 023614 (2015). [12] K. K. Ni, S. Ospelkaus, M. H. G. De Miranda, A. Pe’er, [34] N. Matveeva and S. Giorgini, Phys. Rev. Lett. 109, B.Neyenhuis,J.J.Zirbel,S.Kotochigova,P.S.Julienne, 200401 (2012). D. S. Jin, and J. Ye, Science 322, 231 (2008). [35] G.M.BruunandE.Taylor,Phys.Rev.Lett.101,245301 [13] S. Ospelkaus, K.-K. Ni, D. Wang, M. H. G. de Miranda, (2008). B. Neyenhuis, G. Qu´em´ener, P. S. Julienne, J. L. Bohn, [36] Z.Wu,J.K.Block, andG.M.Bruun,Phys.Rev.B91, D. S. Jin, and J. Ye, Science 327, 853 (2010). 224504 (2015). [14] K. K. Ni, S. Ospelkaus, M. H. G. De Miranda, A. Pe’er, [37] J. K. Block and G. M. Bruun, Phys. Rev. B 90, 155102 B.Neyenhuis,J.J.Zirbel,S.Kotochigova,P.S.Julienne, (2014). D. S. Jin, and J. Ye, Nature 464, 1324 (2010). [38] F. M. Marchetti and M. M. Parish, Phys. Rev. B 87, [15] M.-S. Heo, T. T. Wang, C. A. Christensen, T. M. Rva- 045110 (2013). chov, D. A. Cotta, J.-H. Choi, Y.-R. Lee, and W. Ket- [39] K.S.Singwi,M.P.Tosi,R.H.Land, andA.Sjo¨lander, terle, Phys. Rev. A 86, 021602 (2012). Phys. Rev. 176, 589 (1968). [16] M. Repp, R. Pires, J. Ulmanis, R. Heck, E. D. Kuhnle, [40] A. Pikovski, M. Klawunn, G. V. Shlyapnikov, and M. Weidemu¨ller, and E. Tiemann, Phys. Rev. A 87, L. Santos, Phys. Rev. Lett. 105, 215302 (2010). 010701 (2013). [41] N. Matveeva and S. Giorgini, Phys. Rev. A 90, 053620 [17] C.-H. Wu, J. W. Park, P. Ahmadi, S. Will, and M. W. (2014). Zwierlein, Phys. Rev. Lett. 109, 085301 (2012). [42] J. K. Block, N. T. Zinner, and G. M. Bruun, New J. [18] J.W.Park,S.A.Will, andM.W.Zwierlein,Phys.Rev. Phys. 14, 105006 (2012). Lett. 114, 205302 (2015). [43] U. R. Fischer, Phys. Rev. A 73, 031602 (2006). [19] M.H.G.deMiranda,A.Chotia,B.Neyenhuis,D.Wang, [44] Q. Li, E. H. Hwang, and S. Das Sarma, Phys. Rev. B G.Qu´em´ener,S.Ospelkaus,J.L.Bohn,J.Ye, andD.S. 82, 235126 (2010). Jin, Nature Phys. 7, 502 (2011). [45] L. Zheng and A. H. MacDonald, Phys. Rev. B 49, 5522 [20] T. Miyakawa, T. Sogo, and H. Pu, Phys. Rev. A 77, (1994). 061603 (2008). [46] G. F. Giuliani and G. Vignale, Quantum Theory of the [21] T. Sogo, L. He, T. Miyakawa, S. Yi, H. Lu, and H. Pu, Electron Liquid (Cambridge University Press, 2005). New Journal of Physics 11, 055017 (2009). [47] F. Stern, Phys. Rev. Lett. 18, 546 (1967). [22] Y.Endo,T.Miyakawa, andT.Nikuni,Phys.Rev.A81, [48] A. G. Volosniev, N. T. Zinner, D. V. Fedorov, A. S. 063624 (2010). Jensen, andB.Wunsch,J.Phys.B:At.Mol.Opt.Phys. [23] J.-N. Zhang and S. Yi, Phys. Rev. A 80, 053614 (2009). 44, 125301 (2011). [24] J.-N. Zhang and S. Yi, Phys. Rev. A 81, 033617 (2010). [49] A.Volosniev,J.Armstrong,D.Fedorov,A.Jensen, and [25] J. Krieg, P. Lange, L. Bartosch, and P. Kopietz, Phys. N. Zinner, Few-Body Systems 54, 707 (2013). Rev. A 91, 023612 (2015). [50] L. S´wierkowski, D. Neilson, and J. Szyman´ski, Phys. [26] D. Baillie, R. N. Bisset, and P. B. Blakie, Phys. Rev. A Rev. Lett. 67, 240 (1991). 91, 013613 (2015). [51] Z. Wu, J. K. Block, and G. M. Bruun, “Liquid crystal [27] Z.-K.LuandG.V.Shlyapnikov,Phys.Rev.A85,023614 phases of two-dimensional dipolar gases and berezinskii- (2012). kosterlitz-thouless melting,” ArXiv:1509.02679. [28] Y.Yamaguchi,T.Sogo,T.Ito, andT.Miyakawa,Phys. [52] S.H.Abedinpour,R.Asgari,B.Tanatar, andM.Polini, Rev. A 82, 013643 (2010). Annals of Physics 340, 25 (2014). [29] K. Sun, C. Wu, and S. Das Sarma, Phys. Rev. B 82, 075105 (2010).