Superfluidity and Density Order in a Bilayer Extended Hubbard Model
Tuomas I. Vanhala, Jildou E. Baarsma, Miikka O. J. Heikkinen, Matthias Troyer, Ari Harju, Päivi Törmä
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] A p r Superfluidity and Density Order in a Bilayer Extended Hubbard Model
Tuomas I. Vanhala, Jildou E. Baarsma, Miikka O. J. Heikkinen, Matthias Troyer, Ari Harju, and P¨aivi T¨orm¨a ∗ COMP Centre of Excellence, Department of Applied Physics, Aalto University, FI-00076 Aalto, Finland Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland (Dated: August 15, 2018)We use cellular dynamical mean field theory to study the phase diagram of the square latticebilayer Hubbard model with an interlayer interaction. The layers are populated by two-componentfermions, and the densities in both layers and the strength of the interactions are varied. Wefind that an attractive interlayer interaction can induce a checkerboard density ordered phase andsuperfluid phases, with either interlayer or intralayer pairing. Remarkably, the latter phase doesnot require an intralayer interaction to be present: it can be attributed to an induced attractiveinteraction caused by density fluctuations in the other layer.
PACS numbers: 71.10.Fd,67.85.Fg,67.85.Hj,03.75.Ss
I. INTRODUCTION
Layered lattice systems may exhibit interesting physicssince they allow effects typical to two and three dimen-sions to interplay. For instance, the study of superflu-idity in layered systems has been an active area of re-search since the discovery of high temperature super-conductivity in the cuprates . Another intriguing phe-nomenon in layered systems is the condensation of in-terlayer electron-hole pairs, or excitons, induced by re-pulsive interactions between electrons . Ultracold gasesprovide one more system where various types of layeredsystems can be realized. There the possibility to con-trol inter- and intralayer tunnelings (hoppings) has beenavailable for some time but the interactions have beenlimited to on-site intralayer type. The interesting newdevelopment are dipolar atoms and molecules thatprovide long-range dipole-dipole forces and thereby bringin the possibility of interlayer interactions. In this article,we predict within the bilayer extended Hubbard model acompetition of density order and two types of superflu-ids: an interlayer one analogous to exciton condensates,and an intralayer superfluid that is induced by interlayer interactions.Motivated by the variety of experimental systems, sev-eral theoretical and computational studies have exploredthe physics of bilayer or other two-band Hubbard models.Superfluidity and possible modes of pairing were inves-tigated using the dynamical cluster approximation withrepulsive intralayer interactions and by Monte Carlosimulations with attractive ones . The magnetic prop-erties of the repulsive model were explored in . Theinclusion of an interlayer or interband interaction opensthe possibility of exciton formation and condensation,which has been studied using determinantal Monte Carlosimulations , exact diagonalization , DMFT andvarious other theoretical approaches .Here we consider an extended bilayer Hubbard model,where both layers are populated by two (spin) speciesof fermionic particles. The interaction between the twolayers is attractive, whereas, unlike in most earlier stud- ies, interactions within the layers are vanishing or veryweak. Our results are relevant for experiments with dipo-lar atoms in a bilayer optical lattice, where the dipolemoments are aligned by an external field such that the in-terlayer interaction is attractive . Especially Dy isa promising candidate: Feshbach resonances are observedand using them the s -wave interaction can be tuned tocounteract the on-site dipolar interactions . More-over, by adjusting the interlayer separation the ratio be-tween intra- and interlayer interactions can be tuned sothat nearest-neighbor interactions within the layers arenegligible. Our results are equally relevant to the case ofa repulsive interlayer interaction, since there is an exactparticle-hole transformation between the repulsive andattractive case, and thus a connection to, for instance,exciton condensation can be made .The most striking feature we find this model to ex-hibit is an intralayer superfluid phase, where the Cooperpairs are bound by an interaction induced , via the in-terlayer attraction, by density fluctuations in the otherlayer . This kind of induced superfluidity has previouslybeen proposed by Kuroki and Aoki as a candidatemechanism for superconductivity in strongly correlatedelectron systems. Using Monte Carlo methods on thelattice and the bosonization technique in the continuum,they show that, in one dimensional systems, the inducedinteraction can indeed produce slowly decaying intrabandpairing correlations. In this article we consider, in twodimensions, the competition between the intralayer su-perfluid and other types of long range order. In additionto the induced intralayer superfluid, we find a densityordered phase and an interlayer superfluid. All of thesephases can be reached for a fixed interaction strength bytuning only the densities of the layers, making this modela very promising one to realize experimentally. II. MODEL AND METHOD
The hopping Hamiltonian of the two-component bi-layer square lattice Hubbard model is, H h = − t X σ,l, h i,j i c † σli c σlj − t ⊥ X σ,i ( c † σ i c σ i + h.c. ) − ǫ X σ,i ( n σ i − n σ i ) − µ X σ,i ( n σ i + n σ i ) , (1)where c † σli creates a particle of spin component σ at site i of layer l , the corresponding density operators are de-fined by n σli = c † σli c σli and h i, j i means nearest neighboursummation within the layers. The external field ǫ is re-lated to the interlayer polarization P i = h n i − n i i andthe chemical potential µ can be used to tune the totaldensity ρ i = h n i + n i i /
2, where n il = n il ↑ + n il ↓ . Notethat ǫ does not need to be a physical field but can beunderstood as an energy offset or a difference betweenchemical potentials of the layers. In this article we as-sume the interlayer hopping t ⊥ to be zero, which can berealized experimentally by using a very deep lattice in theperpendicular direction, and take the intralayer hoppingas our unit of energy, t = 1.The interaction Hamiltonian is written in the particle-hole symmetric form H I = V X i,σ,σ ′ (cid:18) n σ i − (cid:19) (cid:18) n σ ′ i − (cid:19) + U X i,l (cid:18) n ↑ li − (cid:19) (cid:18) n ↓ li − (cid:19) , (2)so that half filling corresponds to ǫ = µ = 0. We focuson the case where the onsite interaction U is tuned tozero (e.g. using a Feshbach resonance), and set V = − Dy atoms in a bilayer with intralayerlattice spacing d = 225 nm and the interlayer spacing d ⊥ = d/ V =7 E R and V , ⊥ = 20 V , respectively, with E R the recoilenergy (see appendix A).To study superfluidity in the model, we use cellulardynamical meand field theory (CDMFT) with thecontinuous time auxiliary field (CT-AUX) impuritysolver that is capable of treating general density-densityinteractions. This method maps the lattice problem ontoa finite size cluster that is self-consistently embedded inthe lattice, producing a self energy local to the cluster.Correlation effects within the cluster are treated exactly,and thus the quality of the approximation can be con-trolled by examining the dependence of the results oncluster size.The smallest possible cluster that includes both in-tralayer and interlayer pairing correlations includes onesite in both layers. By symmetry arguments (see ap-pendix B), the relevant superfluid order parameters FIG. 1. The interlayer order parameter ∆ ⊥ (yellow or lightgrey) and intralayer order parameter ∆ l for the more denselayer (red or dark grey) from translation invariant two siteDMFT, as functions of the interlayer density polarizing field ǫ and the chemical potential µ . Half filling is given by µ = 0.The intralayer order parameter is zero where the interlayerorder parameter is nonzero, and vice versa. The inverse tem-perature is 1 /k B T = 15 and the interaction strength V = − t = 1. The datapoints are marked byblack dots, and the surfaces interpolate between them. TheCooper pairs forming the interlayer superfluid are directlybound by the interlayer interaction, while the attraction nec-essary for the intralayer superfluid is indirectly induced bythe other layer. within this cluster are the intralayer pairing order param-eters ∆ l = h c ↑ li c ↓ li i for both layers l and the interlayerpairing order parameter ∆ ⊥ = h c ↑ i c ↓ i i = h c ↓ i c ↑ i i .Here we also consider clusters which consist of L × L × L = 2. We treat the system in the Nambu formalismand include the anomalous Green’s functions between allorbitals of opposite spins.It is interesting to note that doing a particle-hole trans-formation in layer 2 leaves the Hamiltonian invariantapart from ǫ and µ switching roles and a sign changeof V . In other words, there exists a transformation be-tween the attractive and the repulsive model. In thistransformation, the ∆ l remain invariant and the inter-layer order parameter becomes ∆ ⊥ = sgn( i ) D c ↑ i c †↓ i E measuring pairing correlations of excitonic particle-hole-pairs. III. PHASE DIAGRAM
The intralayer and interlayer superfluid order param-eters are shown as a function of the external fields µ and ǫ in Fig. 1. These results were obtained within two-site DMFT ( L = 1) that excludes translational symme-try breaking. We find that the system prefers interlayersuperfluidity for low values of the interlayer density po-larizing field ǫ , while higher polarizations drive it intothe intralayer superfluid state. The intralayer superfluidregion is larger for the more dense layer, i.e. the layer fur-ther away from half filling, and therefore the intralayerorder parameter is plotted for that layer.We have thus predicted that, when the symmetry ofthe layers is distorted by the interlayer density polariz-ing field, an intralayer superfluid emerges that does notrequire intralayer interactions . An explanation of this in-triguing finding could be that the superfluidity is inducedby second order effects mediated by the interlayer interac-tion. To test this hypothesis we calculate the interactionstrength in layer l ′ mediated by density fluctuations inlayer l using a simple mean-field theory , U ind l ′ ∝ V X k ,n G k nl G k nl ∝ V k B T X k n F ( ε k − µ l ) [ n F ( ε k − µ l ) − , (3)where in the first line G k nl is the non-interacting Green’sfunction for the particles in layer l and the summationruns over all energy and momentum states, while the fac-tor 2 results from summing over spin states. In the secondline n F are the Fermi distribution functions, where ε k arethe particle dipersions and µ l determines the density inlayer l . The strength of U ind in one layer thus dependson the density in the other layer and turns out to bestronger if mediated by the layer closer to half filling.This explains why the superfluid is more easily formedin the more dense layer, i.e. in the layer further awayfrom half filling, while intuition gained from e.g. the sin-gle band attractive Hubbard model would suggest theopposite.Now, the crucial question is whether the superfluidswe predict survive the competition with density order,highly typical for lattice systems. To investigate this, wehave obtained the phase diagram of the system also usingthe L = 2 cluster, allowing the possibility of a checker-board density ordered phase, which can be identified bya nonzero value of the order parameter D = 1 N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)X σli sgn( i ) h n σli i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (4)where sgn( i ) is an alternating sign that has the oppositevalues for any pair of neighbouring sites inside one layer,and N is the total number of sites on the lattice. Thedensity order is present close to half filling and zero po-larization, as can be seen in the phase diagrams in Fig. 2aand 2b. Compared to Fig. 1, it can be seen that in quitea considerable region where we first found a superfluidphase, actually the density ordered phase turns out to bethe ground state of the system. However, we still find aregion where the interlayer superfluid phase is the mostfavourable and, strikingly, the intralayer superfluid phase induced by density fluctuations in the opposite layer re-mains present in the phase diagram.Above the critical temperatures of the superfluids, thetransition from the density ordered phase to the normalphase is of the second order and gets sharper as the tem-perature is lowered. At the temperature used in Fig. 2we cannot distinguish this transition from a first orderone, which can also be seen in the large gap present inthe density - interlayer polarization ρ − P plot. In fact,the competition between the superfluid phases and thedensity order is expected to cause a true first order tran-sition for almost the whole phase boundary. This is alsomanifest in the fact that the density ordered region isslightly larger above the superfluid critical temperatures,although it never totally covers the superfluid regionspresent in Fig. 2. The rapid variation of the density andpolarization as a function of ǫ and µ suggests phase sep-aration, in particular between the density ordered phaseand the intralayer superfluid, for which DMFT solutionscan be obtained at the same external polarizing field ǫ but with large differences in polarization P .The translation invariant solution is in itself an inter-esting model for two band superfluidity where the inter-band coupling is the dominant one. The order parame-ters and densities near the intralayer-interlayer transitionare shown in Fig. 3. It can be seen that the transitionbetween these two superfluid states is a sharp, first or-der transition. The interlayer superfluid order parameter∆ ⊥ is zero when the intralayer superfluid order param-eter ∆ l is nonzero and vice versa. The discontinuity inthe polarization and the interlayer order parameter grad-ually decreases with increasing density, and the transitionturns into a second order one when the density is highenough so that the intralayer superfluid is not present.Fig. 3 can be compared to Fig. 2 of reference , whereanalogous behaviour is found in a single band attractiveHubbard model as the spin polarization within the singleband is increased by a magnetic field. In their case thetransition is always a second order one, since the com-petition between two superfluid phases as in our case isabsent.The intralayer-interlayer superfluid transition is causedby a competition between the different modes of pairing.For instance, if we force the intralayer pairing fields tobe zero, the interlayer superfluid region is extended andexhibits a second order transition to the normal statealong the whole phase boundary. Similarly, the intralayersuperfluid is also present at zero polarization if the in-terlayer pairing is excluded. Nevertheless, we did notfind any hysteretic behaviour in the simulations near theintralayer-interlayer transition, which supports the relia-bility of the numerics. ρ P a) L=2U=0 intralayer SFinterlayer SFnormal statedensity order µε b) L=2U=0 µ c) L=1U=0 µ d) L=1 U=0.1
FIG. 2. a) The phase diagram from eigth-site CDMFT at V = − /k B T = 15 as a function of total density ρ andpolarization P . b) and c) Comparison of the ǫ - µ phase diagrams including density order with cluster sizes L = 1 and L = 2.d) The effect of a small intralayer repulsion on the phase diagram. IV. EFFECT OF THE CLUSTER SIZE ANDCOMPARISON TO MEAN FIELD RESULTS
The critical temperatures of the superfluids and thedensity ordered phase at selected points in the ( ǫ, µ ) planeare listed in Table I. There is only little variation be-tween different cluster sizes, which is a sign that neglect-ing nonlocal quantum fluctuations beyond two sites hasonly a small quantitative effect on the results. This isin contrast to the case of exotic superfluidity in quasi-1D systems or d-wave superfluidity in 2D where non-local correlations play a crucial role.It is possible to allow translational symmetry breakingalso within the two-site CDMFT by treating the systemin a partial real space formalism by including twoinequivalent impurity problems representing the high-density and low-density sublattices of the checkerboarddensity order. A comparison of the phase diagrams ob-tained within this approximation and using the L = 2cluster is presented in Fig. 2b and 2c. We find that theresults are in good agreement, which further supports theconclusion that the effect of the cluster size is not veryimportant.The effects of different approximations on the transi-tion from the normal phase to the density order as a func-tion of the inverse temperature are illustrated in Fig. 4.As expected, the critical temperature is slightly lowerfor larger clusters, but the variation is small and the ef-fect of the cluster size is insignificant when the systemis not very close to the critical point. In contrast, themean field approximation gives considerably higher criti-cal temperatures, emphasizing the importance of moreaccurate methods. For a weaker interaction strength V = − L = 1) is T c = 0 . ± .
01. In this case the system ε ∆ ⊥ ∆ ∆ ε ρ −1.29P FIG. 3. The superfluid order parameters (left) and the den-sity ρ and interlayer polarization P (right) from translationinvariant two-site DMFT as a function of the polarizing field ǫ in the intralayer-interlayer transition region. The interactionstrength is V = −
3, the chemical potential µ = 0 . /k B T = 15. is well in the perturbative region, and while the meanfield treatment still gives considerable errors, includingsecond order corrections (see appendix C) already yieldsa good approximation. For V = −
3, however, even withthe second order corrections mean field theory is clearlyinadequate.
V. EFFECT OF ONSITE REPULSION ANDINTERLAYER HOPPING
The presence of a small repulsive onsite interaction U does not destroy the intralayer superfluid. For example,the superfluidity persists up to U ≈ . ǫ = 1 . , µ = 0 . TABLE I. Critical temperatures T c (in units of t/k B ) atselected points in the ( ǫ, µ ) plane for different cluster sizes L × L ×
2. The interaction strength is V = − ǫ, µ ) (1 . , .
45) (0 . , .
6) (0 . , . T c (L=1) 0 . ± .
003 0 . ± .
007 0 . ± . T c (L=2) 0 . ± .
003 0 . ± .
007 0 . ± . T c (L=4) 0 . ± . β V=−10 2 400.20.40.60.8 β D V=−3
MFMFCL=1L=2L=4
FIG. 4. The order parameter D of the density ordered phaseas a function of the inverse temperature β = 1 /k B T frommean field theory (MF), mean field including second ordercorrections (MFC), and CDMFT with cluster size L × L × µ = ǫ = 0. Fig. 2.The µ − ǫ phase diagram for U = 0 . t ⊥ is shown inFig. 5. The possibility to tunnel between the layers turnsthe sharp first order transition between the interlayer andintralayer superfluid states into a smooth crossover. Thisis expected, since in this case an intralayer Cooper paircan turn into an interlayer one and vice versa thus mixingthe two types of superfluids. Note that the smootheningis not very large for t ⊥ = 0 .
1, and for t ⊥ = 0 .
01 it isnegligible.
VI. CONCLUSIONS
In conclusion, we have studied the bilayer extendedHubbard model for attractive interlayer interactions us-ing a beyond mean field approach. We found density order starting from half filling and, when the density isincreased, an interlayer superfluid that appears beforereaching the normal state. This superfluid, where theCooper pairs are formed by particles from different lay-ers, is conceptually related to bilayer exciton condensateswhere pairing happens between particles in one layer andholes in the other. Notably, the freedom to polarize adensity difference between the layers revealed also a noveltype of superfluid with intralayer pairing that is inducedby density fluctuations in the opposite layer. The pre-dicted phases could be experimentally realized, for in-stance, with ultracold dipolar atoms or molecules. Aprospect for future research is to study if some othertype of density order than the checkerboard consideredhere, possibly incommensurate, would be present nearthe density order phase boundary. This would be partic-ularly interesting as there is a possibility of simultaneoussuperfluidity and breaking of translation symmetry, i.e.a supersolid phase.
ACKNOWLEDGMENTS
This work was supported by the Academy of Fin-land through its Centers of Excellence Programme (2012-2017) and under Projects No. 263347, No. 251748,No. 272490, and by the European Research Coun-cil (ERC-2013-AdG-340748-CODE and ERC-2011-AdG-290464-SIMCOFE). T.I.V. is grateful for the supportfrom the Vilho, Yrj¨o and Kalle V¨ais¨al¨a Foundation.Computing resources were provided by CSC – the FinnishIT Centre for Science and the Triton cluster at Aalto Uni-versity.
Appendix A: Hopping parameters and energies inthe extended Hubbard model
In the Hubbard model the hopping parameter t de-scribing hopping from site j to the nearest-neighbor site j ′ is given by t = − Z d rφ ∗ j ( r ) (cid:20) − ¯ h m ∇ + V latt ( r ) (cid:21) φ j ′ ( r ) , (A1)where V latt ( r ) is the lattice potential. To calculate theinterlayer hopping t ⊥ in our model, j and j ′ are takenon different layers, whereas for the intralayer hopping t the sites j and j ′ are in the same layer. The interactionparameters are given by U = 4 π ¯ h a S m Z d r | φ j ( r ) | + Z Z d rd r ′ | φ j ( r ) | | φ j ( r ′ ) | U ( r − r ′ ) (A2) V j,j ′ = Z Z d rd r ′ | φ j ( r ) | | φ j ′ ( r ′ ) | U ( r − r ′ ) , (A3) ρ P ∆ ⊥ ∆ t ⊥ = . εµ εµ εµ
012 00.5100.040.08 µε t ⊥ = . εµ εµ εµ
012 00.5100.040.08 µε t ⊥ = . εµ εµ εµ
012 00.5100.040.08 µε t ⊥ = . εµ εµ εµ
012 00.5100.040.08 µε FIG. 5. A comparison of the superfluid order parameters ∆ l and ∆ ⊥ , the density ρ and the polarization P for different valuesof the interlayer hopping t ⊥ . These results were obtained within translation invariant two-site DMFT, which excludes thedensity ordered phase. The parameters were t = 1, V = − U = 0 and the inverse temperature β = 15. Note that in the tworightmost columns the axes µ , ǫ are rotated to better display the sharpness of the transition. where a S is the s -wave scattering length and U ( r ) is thedipole-dipole interaction , which in the case of all dipolesbeing aligned reads U ( r ) = C dd π − θr , (A4)where θ is the angle between r and the direction of thepolarization. Here the C dd the strength of the dipolarinteraction which in case of dysprosium is C dd = µ µ with the magnetic dipole moment µ = 10 µ B . The on-site interaction U results from the contact interactionand the dipole-dipole interaction, whereas the intersiteinteractions V j,j ′ result from the dipole-dipole interac-tions only. To make an order of magnitude estimationfor the hoppings and interactions we use in this article,we calculated the above coefficients using the tight bind-ing approximation (The onsite wavefunctions φ i ( r ) weretaken to be harmonic oscillator ground state wave func- tions). We found that for Dy in a layer with latticespacing d = 225 nm and lattice height V = 3 . E R ,with E R the recoil energy, the onsite interaction can betuned to zero, U = 0, since there are several Feshbachresonances available to control a S . An interlayer interac-tion V = − t can be reached when the spacing betweenthe layers is d ⊥ = d/ V , ⊥ = 20 V . In that case theinterlayer hopping t ⊥ and intralayer nearest-neighbor in-teraction V k are small enough (less than 10%) in com-parison to t and V respectively, to be taken to zero inour model. A stronger interaction, compared to the hop-ping, can be reached by taking larger lattice heights inall directions. The stronger interaction we consider inthis article can be reached when the spacing between thelayers is d ⊥ = d/ V = 7 E R and V , ⊥ = 20 V in the parallel and perpendicular directions,respectively. Appendix B: Choice of the superfluid orderparameters
In this appendix we provide details for the choice ofthe superfluid order parameters taking into account thesymmetries of the system. We first discuss the case t ⊥ =0 and comment on the general case in the end.Within the L = 1 cluster, it is possible to have a sixcomponent pairing order parameter which consists of theintralayer order parameter ∆ l = h c ↑ li c ↓ li i in each layer l and of the four component matrix M defined by ψ il = c ↑ li c ↓ li ! , M = (cid:10) ψ i ψ Ti (cid:11) . (B1)As the interlayer hopping vanishes, and the interlayerinteraction only couples the total densities of the layers,the Hamiltonian is invariant under rotations ψ ′ il = U l ψ il in both layers l where U l is a unitary matrix. The SU (2)-part of this symmetry corresponds to conservation ofspin, and the U (1)-part corresponds to conservation ofparticle number, in both layers separately. Letting U l = a bc d ! , (B2)the transformation for ∆ l is∆ ′ l = D c ′↑ li c ′↓ li E = h ( ac ↑ li + bc ↓ li ) ( cc ↑ li + dc ↓ li ) i = h ( ad − bc ) c ↑ li c ↓ li i = det( U l )∆ l . (B3)Thus the ∆ l are invariant in pure SU (2) spin rotationsand only gain a phase det( U l ) in general unitary trans-formations.The transformation for the interlayer order parametersis given by M ′ = D ψ ′ i ( ψ ′ i ) T E = D U ψ i ( U ψ i ) T E = U D ψ i ( ψ i ) T E U T = U M U T . (B4)By the singular value decomposition it is always possibleto choose U and U so that M ′ is diagonal and only hasnonnegative real elements. Thus, without losing general-ity, the order parameter reduces to two nonnegative realvalued fields and the two, in general complex valued, ∆ l .By doing a further transformation with the matrices U = ! , U = I, (B5)∆ gets a minus sign and the interlayer order parame-ter matrix M ′ becomes purely off diagonal, the nonzerocomponents being h c ↑ i c ↓ i i and h c ↓ i c ↑ i i .This is the representation where we perform our nu-merics. Because of technical limitations we also considerthe intralayer order parameters to be real valued. As we have fixed the interlayer order parameters to be non-negative, this still leaves the signs of the intralayer orderparameters as nontrivial factors that cannot be rotatedaway by symmetry. However, from the simulations wefind that the interlayer order parameters (the matrix M ′ )are always zero, when the intralayer order parameters arenonzero, and vice versa. Under this assumption it is al-ways possible to find a transformation that flips the signof one of the nonzero components, leaving the other oneinvariant. Thus all sign configurations are equivalent,and we are free to choose the order parameters to bepositive, including the intralayer ∆ l . Furthermore, wefind that the two nonzero components of M ′ are alwaysequal, and thus we just call this value the interlayer orderparameter ∆ ⊥ ≡ h c ↑ i c ↓ i i = h c ↓ i c ↑ i i .Also within larger clusters, our CDMFT implementa-tion includes all anomalous Green’s functions betweenopposite spin components within the cluster. The t ⊥ = 0case can be thougth of as a limiting case of a modelwhere the particles have a small probability to tunnel be-tween the layers. This tunneling breaks the conservationof particle number in both layers separately, thus leav-ing only two exactly conserved spin components. Thusit is expected that the superfluidity can be treated ina picture where the Cooper pairs are only formed be-tween particles of opposite spin, and not between par-ticles of the same spin in different layers. In the caseof nonzero t ⊥ the intralayer and interlayer order pa-rameters take finite values simultaneously, and the rel-ative signs of the order parameters become important.We find that the preferred sign configuration is suchthat ∆ ⊥ ≡ h c ↑ i c ↓ i i = − h c ↓ i c ↑ i i = h c ↑ i c ↓ i i and∆ ∆ < Appendix C: Mean field treatment of the densityordered phase
To make comparisons with the CDMFT solution of themodel, we have also considered a mean field type approx-imation for the density ordered phase. The lowest ordermean field approximation can be obtained by decompos-ing the interlayer density-density interaction as (cid:0) n A − (cid:1) (cid:0) n B − (cid:1) = n A n B − ( n A + n B ) + = ( h n A i + δ A )( h n B i + δ B ) − ( n A + n B ) + ≈ h n A i δ B + h n B i δ A + h n A i h n B i − ( n A + n B )= h n A i n B + h n B i n A − h n A i h n B i − ( n A + n B ) , where A and B stand for spin, site and layer indices, δ A = n A −h n A i and δ B = n B −h n B i are fluctuations fromthe mean value and terms quadratic in the fluctuationshave been neglected as well as the constant 1 /
4. Thisleads to the mean field Hamiltonian H = H h + V P i,σ,σ ′ ,l (cid:0) h n σli i − (cid:1) n σ ′ l ′ i − V P i,σ,σ ′ h n σ i i h n σ ′ i i , (C1) + +1 2 FIG. 6. The two-particle irreducible diagrams contributing tothe self energy up to second order in the interaction strength V . The solid lines stand for the full interacting propagator G , and the dashed lines represent the interaction vertices.The second order diagram depends on the imaginary timedifference τ − τ thus bringing frequency dependence to theself energy. The layer indices must fulfill the condition l ′ = l ,since there are no interactions within the layers, while thespin index σ ′ can freely be summed over. where l ′ = 1 when l = 2 and vice versa. The resultingself energy is diagonal in spin, site and layer indices, andis given in Matsubara frequency space byΣ (1) σli ( iω n ) = − V X σ ′ (cid:18) h n σ ′ l ′ i i − (cid:19) . (C2)Subsequently, the propagator G σlij ( τ ) = D c σli ( τ ) c † σlj (0) E can be calculated as G σl ( τ ) = 1 β X ω n exp ( − iω n τ ) ( − iω n + T − Σ σl ( iω n )) − , (C3)where G , T and Σ are matrices in the site indices, T isthe hopping matrix of the square lattice including the chemical potential contributions, and β is the unitlessinverse temperature. In practice the matrix inversion isdone in Fourier space with a unit cell that allows thebreaking of translation invariance. The self-consistencycondition is that the density h n σli i = 1 − G σlii ( τ = 0+)agrees with the density in Eqn. C2.It is possible to derive the mean field theory from theBaym-Kadanoff (or Luttinger-Ward) functional formal-ism, which can also be used to systematically includehigher order corrections to the approximation . Inthis formulation, the self energy is expressed as a dia-grammatic series expansion in terms of the interactionvertices and interacting propagator lines. As only thetwo-particle irreducible diagrams enter this series, thecontributions up to second order in V consist of only twodiagrams, which are depicted in Fig. 6.Including only the first order diagram gives the meanfield self energy C2. The second order diagram includescontributions which are not diagonal in the site indices,see Fig. 6. However, we neglect these contributions andapply a local approximation where both vertices of thediagram are on the same site i = j . Evaluating the sec-ond order diagram gives the correction termΣ (2) σli ( τ − τ ) = − V G σlii ( τ − τ ) ·· P σ ′ G σ ′ l ′ ii ( τ − τ ) G σ ′ l ′ ii ( τ − τ ) , (C4)where again l ′ = l . We calculate this term directly inan imaginary time grid and then perform a numericalFourier transformation to Matsubara frequencies. Thisyields a self-energy Σ σli ( iω ) = Σ (1) σli ( iω ) + Σ (2) σli ( iω ) whichcan be used to calculate the propagator according to Eqn.C3. The then obtained propagator is used to calculatethe self energy again, and the process is iterated until aconverged, self-consistent solution is found. ∗ paivi.torma@aalto.fi Patrick A. Lee, Naoto Nagaosa, and Xiao-Gang Wen,“Doping a Mott insulator: Physics of high-temperaturesuperconductivity,” Rev. Mod. Phys. , 17–85 (2006). J. P. Eisenstein and A. H. MacDonald, “Bose–Einstein con-densation of excitons in bilayer electron systems,” Nature , 691–694 (2004). Immanuel Bloch, Jean Dalibard, and Wilhelm Zwerger,“Many-body physics with ultracold gases,” Rev. Mod.Phys. , 885–964 (2008). D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, andP. Zoller, “Cold Bosonic Atoms in Optical Lattices,” Phys.Rev. Lett. , 3108–3111 (1998). S. Ospelkaus, A. Pe’er, K. K. Ni, J. J. Zirbel, B. Neyen-huis, S. Kotochigova, P. S. Julienne, J. Ye, and D. S. Jin,“Efficient state transfer in an ultracold dense gas of het-eronuclear molecules,” Nat Phys , 622–626 (2008). M. A. Baranov, M. Dalmonte, G. Pupillo, and P. Zoller,“Condensed Matter Theory of Dipolar Quantum Gases,” Chemical Reviews , 5012–5061 (2012). C. Trefzger, C. Menotti, B. Capogrosso-Sansone, andM. Lewenstein, “Ultracold dipolar gases in optical lat-tices,” Journal of Physics B: Atomic, Molecular and Opti-cal Physics , 193001 (2011). Axel Griesmaier, J¨org Werner, Sven Hensler, J¨urgen Stuh-ler, and Tilman Pfau, “Bose-Einstein Condensation ofChromium,” Phys. Rev. Lett. , 160401 (2005). L. Santos, “Dipolar Gases – Theory,” in
Quantum GasExperiments – Exploring Many-Body States , edited byP. T¨orm¨a and K. Sengstock (Imperial College Press, Lon-don, 2015) pp. 293–309. E. A. L. Henn, J. Billy, and T. Pfau, “Dipolar Gases –Experiment,” in
Quantum Gas Experiments – ExploringMany-Body States , edited by P. T¨orm¨a and K. Sengstock(Imperial College Press, London, 2015) pp. 311–325. T. A. Maier and D. J. Scalapino, “Pair structure and thepairing interaction in a bilayer Hubbard model for un-conventional superconductivity,” Phys. Rev. B , 180513 (2011). Aleksander Zujev, Richard T. Scalettar, George G. Ba-trouni, and Pinaki Sengupta, “Pairing correlations inthe two-layer attractive Hubbard model,” New Journal ofPhysics , 013004 (2014). Robert R¨uger, Luca F. Tocchio, Roser Valent´ı, andClaudius Gros, “The phase diagram of the square latticebilayer Hubbard model: a variational Monte Carlo study,”New Journal of Physics , 033010 (2014). Louk Rademaker, Steve Johnston, Jan Zaanen, and Jeroenvan den Brink, “Determinant quantum Monte Carlo studyof exciton condensation in the bilayer Hubbard model,”Phys. Rev. B , 235115 (2013). T. Kaneko, S. Ejima, H. Fehske, and Y. Ohta, “Exact-diagonalization study of exciton condensation in electronbilayers,” Phys. Rev. B , 035312 (2013). J. Kunes, “Phase diagram of exciton condensate indoped two-band Hubbard model,” ArXiv e-prints (2014),arXiv:1410.5198 [cond-mat.str-el]. Jan Kunes and Pavel Augustinsk´y, “Excitonic instabil-ity at the spin-state transition in the two-band Hubbardmodel,” Phys. Rev. B , 115134 (2014). Yuh Tomio, Kotaro Honda, and Tetsuo Ogawa, “Exci-tonic BCS-BEC crossover at finite temperature: Effects ofrepulsion and electron-hole mass difference,” Phys. Rev. B , 235108 (2006). Louk Rademaker, Jeroen van den Brink, Jan Zaanen, andHans Hilgenkamp, “Exciton condensation in strongly cor-related electron bilayers,” Phys. Rev. B , 235127 (2013). Mingwu Lu, Nathaniel Q. Burdick, and Benjamin L. Lev,“Quantum Degenerate Dipolar Fermi Gas,” Phys. Rev.Lett. , 215301 (2012). K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler,R. Grimm, and F. Ferlaino, “Bose-Einstein Condensationof Erbium,” Phys. Rev. Lett. , 210401 (2012). A. Pikovski, M. Klawunn, G. V. Shlyapnikov, andL. Santos, “Interlayer Superfluidity in Bilayer Systems ofFermionic Polar Molecules,” Phys. Rev. Lett. , 215302(2010). Bo Yan, Steven A. Moses, Bryce Gadway, Jacob P. Covey,Kaden R. A. Hazzard, Ana M. Rey, Deborah S. Jin, andJun Ye, “Observation of dipolar spin-exchange interactionswith lattice-confined polar molecules,” Nature , 521–525 (2013). G. M. Bruun and E. Taylor, “Quantum Phases of a Two-Dimensional Dipolar Fermi Gas,” Phys. Rev. Lett. ,245301 (2008). S. M¨uller, J. Billy, E. A. L. Henn, H. Kadau, A. Griesmaier,M. Jona-Lasinio, L. Santos, and T. Pfau, “Stability ofa dipolar Bose-Einstein condensate in a one-dimensionallattice,” Phys. Rev. A , 053601 (2011). Kristian Baumann, Nathaniel Q. Burdick, Mingwu Lu,and Benjamin L. Lev, “Observation of low-field Fano-Feshbach resonances in ultracold gases of dysprosium,”Phys. Rev. A , 020701 (2014). Tian-Sheng Zeng and Lan Yin, “Supersolidity of a dipolarFermi gas in a cubic optical lattice,” Phys. Rev. B ,174511 (2014). Hideo Aoki and Kazuhiko Kuroki, “Superconductivity ina two-band Hubbard model,” Phys. Rev. B , 2125–2136(1990). Kazuhiko Kuroki and Hideo Aoki, “Realization of negative-U superconductivity in a class of purely repulsivesystems: Interacting carrier and insulating bands,” Phys.Rev. Lett. , 3820–3823 (1992). Kazuhiko Kuroki and Hideo Aoki, “Superconductivity ina repulsively interacting two-band Fermi gas,” Phys. Rev.Lett. , 2947–2950 (1994). Antoine Georges, Gabriel Kotliar, Werner Krauth, andMarcelo J. Rozenberg, “Dynamical mean-field theory ofstrongly correlated fermion systems and the limit of infinitedimensions,” Rev. Mod. Phys. , 13–125 (1996). Thomas Maier, Mark Jarrell, Thomas Pruschke, andMatthias H. Hettler, “Quantum cluster theories,” Rev.Mod. Phys. , 1027–1080 (2005). Gabriel Kotliar, Sergej Y. Savrasov, Gunnar P´alsson, andGiulio Biroli, “Cellular Dynamical Mean Field Approach toStrongly Correlated Systems,” Phys. Rev. Lett. , 186401(2001). E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov,M. Troyer, and P. Werner, “Continuous-time Monte Carlomethods for quantum impurity models,” Reviews of Mod-ern Physics , 349–404 (2011). Emanuel Gull, Peter Staar, Sebastian Fuchs, PhaniNukala, Michael S. Summers, Thomas Pruschke,Thomas C. Schulthess, and Thomas Maier, “Sub-matrix updates for the continuous-time auxiliary-fieldalgorithm,” Phys. Rev. B , 075122 (2011). E. Gull, P. Werner, O. Parcollet, and M. Troyer,“Continuous-time auxiliary-field Monte Carlo for quantumimpurity models,” EPL (Europhysics Letters) , 57003(2008). Akihisa Koga and Philipp Werner, “Polarized Superfluid-ity in the Imbalanced Attractive Hubbard Model,” Journalof the Physical Society of Japan , 064401 (2010). Gerald D. Mahan,
Many-Particle Physics , 2nd ed.(Plenum, New York, N.Y., 1993). M. O. J. Heikkinen, D. H. Kim, M. Troyer, and P. T¨orm¨a,“Nonlocal Quantum Fluctuations and Fermionic Superflu-idity in the Imbalanced Attractive Hubbard Model,” Phys.Rev. Lett. , 185301 (2014). T. A. Maier, M. Jarrell, T. C. Schulthess, P. R. C. Kent,and J. B. White, “Systematic Study of d-Wave Supercon-ductivity in the 2D Repulsive Hubbard Model,” Phys. Rev.Lett. , 237001 (2005). Minh-Tien Tran, “Inhomogeneous phases in the Falicov-Kimball model: Dynamical mean-field approximation,”Phys. Rev. B , 205110 (2006). M Snoek, I Titvinidze, C Tke, K Byczuk, and W Hof-stetter, “Antiferromagnetic order of strongly interactingfermions in a trap: real-space dynamical mean-field anal-ysis,” New Journal of Physics , 093008 (2008). G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko,O. Parcollet, and C. A. Marianetti, “Electronic struc-ture calculations with dynamical mean-field theory,” Re-views of Modern Physics , 865–951 (2006), arXiv:cond-mat/0511085. Gordon Baym and Leo P. Kadanoff, “Conservation Lawsand Correlation Functions,” Phys. Rev. , 287–299(1961). J. M. Luttinger and J. C. Ward, “Ground-State Energy ofa Many-Fermion System. II,” Phys. Rev.118