Collisionless drag for a one-dimensional two-component Bose-Hubbard model
Daniele Contessi, Donato Romito, Matteo Rizzi, Alessio Recati
CCollisionless drag for a one-dimensional two-component Bose-Hubbard model
Daniele Contessi, Donato Romito,
2, 3
Matteo Rizzi,
4, 5 and Alessio Recati
2, 1 Dipartimento di Fisica, Universit`a di Trento, 38123 Povo, Italy INO-CNR BEC Center, 38123 Povo, Italy Mathematical Sciences, University of Southampton,Highfield, Southampton, SO17 1BJ, United Kingdom Forschungszentrum J¨ulich, Institute of Quantum Control,Peter Gr¨unberg Institut (PGI-8), 52425 J¨ulich, Germany Institute for Theoretical Physics, University of Cologne, D-50937 K¨oln, Germany (Dated: September 2020)We theoretically investigate the elusive Andreev-Bashkin collisionless drag for a two-componentone-dimensional Bose-Hubbard model on a ring. By means of Tensor Network algorithms, wecalculate superfluid stiffness matrix as a function of intra- and inter-species interactions and of thelattice filling. We then focus on the most promising region close to the so-called pair-superfluidphase, where we observe that the drag can become comparable with the total superfluid density.We elucidate the importance of the drag in determining the long-range behavior of the correlationfunctions and the spin speed of sound. In this way we are able to provide an expression for the spinLuttinger parameter K S in terms of drag and the spin susceptibility.Our results are promising in view of implementing the system by using ultra-cold Bose mixturestrapped in deep optical lattices, where the size of the sample is of the same order of the numberof particle we simulate. Importantly the mesoscopicity of the system far from being detrimentalappears to favour a large drag, avoiding the Berezinskii-Kosterlitz-Thouless jump at the transitionto the pair superfluid phase which would reduce the region where a large drag can be observed. The dynamics of multi-component systems, rangingfrom neutron stars [1] to superconducting layers [2],is supposed to be crucially influenced by an inter-component dissipationless drag. Such an entrainmenthas been first discussed by Andreev and Bashkin in1975 to describe the three-fluid hydrodynamics of a two-component mixture of He and He superfluids [3]. Upto that point, the constitutive relations for the super-fluid momenta were assumed to involve the transport ofparticles of one kind only [4–7]. Andreev and Bashkin in-stead introduced a superfluid stiffness matrix n ( s ) αβ , whoseoff-diagonal elements make it possible that a velocity v β in one component generates a (super-)current j α in theother component even without collisions: j α = (cid:88) β n ( s ) αβ v β . (1)Although the Andreev-Bashkin (AB) effect should begeneric of multi-component systems, it has never beendirectly observed (e.g., in Helium due to the very lowmiscibility of the two isotopes) and its dependence onthe microscopic parameters has been obtained only forfew cases.The high degree of control reached in manipulatingultra-cold multi-component Bose gases has sparkled thehope of having a platform for a detailed experimentalstudy of the collisionless drag. However, beyond mean-field corrections play a crucial role to achieve a sizeableentrainment. Unfortunately, in the most standard config-urations, i.e. three-dimensional homogeneous (or quasi-homogeneous) systems [8, 9], it is difficult to enhance theeffect by increasing the interactions. Indeed strong inter- actions, although possible, lead to strong 3-body losses,which make the study of the AB effect hard. Fortunately,other routes are available to increase the role of quantumfluctuations while keeping the system stable.Very recently, configurations with reduced dimension-ality have been proposed: in particular, in [10] it hasbeen shown that approaching the molecular phase in adouble-layer dipolar gas system the drag can become in-creasingly large. In [11], one-dimensional mixtures closeto the so-called Tonks-Girardeau regime have been shownto exhibit a large entrainment.Another possibility is to consider Hubbard-like models,by putting cold-atoms in deep optical lattices. In thisway it is possible not only to realise strongly interactingsuperfluids with reduced 3-body losses, but also to studynew phases that do not appear in continuous systems.An analysis of the AB effect in a two-component single-band two-dimensional Bose-Hubbard model can be foundin [12], where the effect of the proximity to the Mottinsulating phase is discussed in detail.In the present Letter we give a detailed account ofthe AB drag in two-component Bose-Hubbard modelon a one-dimensional ring. The reason is many-fold.The ring geometry is very convenient to study super-current related phenomena, both theoretically and ex-perimentally. The discrete nature allows us to have afinite range of parameters in the attractive regime wherethe two-superfluid state is stable and the drag can bestrongly enhanced [10, 12]. Moreover, differently to allthe above mentioned works based on Quantum Monte-Carlo (QMC) techniques, here we employ a Tensor Net-work approach, thus paving the way to study time de- a r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p pendent phenomena as well.Aside from determining the strength of the AB effect,our microscopic results show that the Luttinger liquiddescription of the low-energy dynamics of such systemsmust include an inter-velocity interaction. We shed newlight upon the meaning of the spin Luttinger parameter K s , by providing an expression in terms of the superfluiddensity matrix and the susceptibility of the system. Thiseffect is rooted in the fact that the f -sum rule for thespin channel is not exhausted by single-phonon excita-tions (see [13] and the Supplementary Material of [11]).Furthermore, we show that finite-size effects might ac-tually increase the visibility of the collisionless drag, bycircumventing the sudden jump that characterises thephase transition to a pair-superfluid phase in the ther-modynamic limit. This is particularly relevant for typ-ical 1D cold atomic setups, where the particle numberis comparable to our numerical simulations. In this re-spect hyperfine states mixtures of K atoms, or K - Rb mixtures (e.g. [14–18]), whose inter-species interactioncan be tuned by exploiting Feshbach resonances, seemvery promising to achieve the regimes where the elusiveAB effect could be finally observed. FIG. 1. Sketch of the 2-component Bose-Hubbard ring withhopping parameters ˜ t α = t α e − i πφ α /L , on-site intra-speciesinteractions U α and inter-species interaction U AB . The twofluxes φ α pierce the ring and give rise to bosonic currents j α in the ground state. Model:
As sketched in Fig. 1, we consider a 2-speciesBose-Hubbard Hamiltonian, H = H A + H B + H AB , on aring with L sites: H α = L (cid:88) x =1 (cid:20) − (cid:16) ˜ t α b † x +1 ,α b x,α + h . c . (cid:17) + U α n x,α ( n x,α − (cid:21) ,H AB = U AB L (cid:88) x =1 n x,A n x,B , (2) where b † x,α ( b x,α ) is the bosonic creation (annihilation) op-erator and n x,α = b † x,α b x,α is the number operator at site x for the species α ∈ { A, B } . The two species can beinternal degrees of freedom, band indices or two differentatomic elements. The single species Hamiltonian H α ac-counts for the hopping between neighboring sites, with˜ t α = t α e − i πφ α /L and t α , φ α ∈ R + , and for the on-siterepulsion characterised by the parameter U α >
0. Thefluxes φ α piercing the ring are introduced in order to com-pute the superfluid currents and densities (see Eqs. (3)-(4)), and are equivalent to twisted periodic boundaryconditions (PBC) [19]. The Hamiltonian H AB describesthe inter-species on-site interaction – responsible for thecollisionless drag phenomenon – with strength U AB .We limit ourselves to a zero temperature, Z symmet-ric mixture: t α = t , U α = U and filling ν α ≡ N α /L = ν/
2, in terms of the number of atoms N α . The phasediagram of the 1D model is very rich and has not yetbeen determined with the same accuracy as in higher di-mensions [20, 21]: to our knowledge the most completeanalysis can be found in [22]. It is beyond the scope ofthe present Letter to bridge this gap and the full phasediagram will be reported elsewhere [23]. For the purposeof the present Letter we are interested in only two of thepossible phases: the two-superfluid (2SF) phase and thepair-superfluid (PSF) phase. The 2SF phase is charac-terised by both components A and B being superfluid.The low energy spectrum consists of two gapless linear(Goldstone) modes corresponding to a density (in-phase)and a spin (out-of-phase) mode. In the PSF phase, thetwo components are paired and the spin-channel acquiresa gap [24]. One of our goals here is to determine the su-perfluid density matrix in the 2SF phase while approach-ing the PSF, where the collisionless drag should saturateto its maximum value [10, 12].We resort to a Matrix Product States (MPS) ansatz todeal with the full many-body problem. The model (2) isindeed not exactly solvable and our numerical treatmentis an almost unbiased approach to it. We overcome thedifficulties related to PBC by employing a loop-free ge-ometry of the tensor network and shifting the topology ofthe lattice into the Matrix Product Operator (MPO) rep-resentation of the Hamiltonian. The idea consists essen-tially in introducing nearest-neighbour couplings alonga snake-like enumeration of the physical sites [25]. Inthis way, we can reliably compute relevant quantities forsystems up to L = 96 sites achieving expectation values(e.g., of densities and currents) homogeneous up to 0 . Computing the superfluid densities:
Thanks to thediscrete translational invariance of the model and makinguse of the Hellmann-Feynman theorem [25], the currentdensity of species α is given by: j α = 2 (cid:126) Im (cid:104) ˜ t α b † x +1 ,α b x,α (cid:105) = 12 π (cid:126) ∂E∂φ α . (3)By linearising the currents for small fluxes and usingthe relation v β = (cid:126) m ∗ πφ β L , with m ∗ = (cid:126) / t the “bandmass”, we can compute the AB [3] superfluid density ma-trix n ( s ) αβ of Eq. (1) as: n ( s ) αβ = lim φ α ,φ β → Lm ∗ π (cid:126) ∂j α ∂φ β = lim φ α ,φ β → Lm ∗ (2 π (cid:126) ) ∂ E∂φ α ∂φ β . (4)Within Tensor Network methods the total energy andcurrent densities, computed from short-range correla-tions, are among the most reliable quantities to be ex-tracted. n ( s ) A B / n U AB /U . . . . . − − . . C o ll a p s e tr a n s i t i o n P h a s e - S e p a r a t i o n n ( s ) A B / ( n ( s ) AA + n ( s ) A B ) . . . . . . − − . . FIG. 2. (color online). Superfluid drag in terms of the to-tal density (main panel) and of the total superfluid density(inset). The red dashed line indicates the theoretical predic-tion via Bogoliubov approximation [13]. The relevant pairingcorrelations are responsible for the non symmetric trend ex-hibited by the attractive regime ( U AB < Drag at half-filling:
In order to highlight the roleof pairing correlations – relevant in 1D for any value U AB < n ( s ) AB = n ( s ) BA as a func-tion of the inter-species interaction U AB /U for half to-tal filling ν = 0 . U/t = 2. In this regime, as ex-pected from previous analysis [22], the mixture is alwaysin the 2SF phase, until undergoes either phase separa-tion or collapse. The location of the phase transitions isstrongly dependent on the parameters of the configura-tion. In our case, the collapse occurs beyond the blackdashed line in the shaded region. This is in contrast withthe mean-field Bogoliubov prediction (red dashed line)of a symmetric behaviour for attractive and repulsive in-teractions [13] that sets the transitions at | U AB | = U .The data display an evident asymmetry between the tworegimes. In fact the attractive mixture experiences amuch steeper growth of the drag as a function of theinteraction strength. This substantial increase can beascribed to pairing correlations, which are not capturedby the Bogoliubov approach. Superfluid drag approaching the PSF phase:
The in-crease of the entrainment for U AB < n ( s ) AB = n ( s ) AA . Therefore, n ( s ) AB becomes a quarter ofthe total superfluid density n ( s ) = n ( s ) AA + n ( s ) BB +2 n ( s ) AB , butsimultaneously ceases to be interpreted as a drag coeffi-cient. The results for different system’s sizes are reportedin Fig. 3, where U/t = 10 and ν = 1 are chosen to makesure the system can undergo the transition [22]). In theinset we report the behaviour of the normalized drag asa function of L − for different values of the interaction.We estimate that the SF-to-PSF transition should occurin the thermodynamic limit for U AB /U ∈ [ − . , − . n ( s ) AA − n ( s ) AB , from a finite value to zero. We stress, how-ever, that mesoscopic samples like the ones accessible incold-atomic setups will display no jump, but rather asizeable value of n ( s ) AB , thus making the AB collisionlessdrag finally observable. Collisionless Drag and Luttinger liquid parameters:
In the 2SF phase, the low energy theory corresponds ( n ( s ) A B ) / ( n ( s ) AA + n ( s ) A B ) U AB /U . . . . . − . − . − . − . − . − . − . − .
05 0 P S F tr a n s i t i o n /L . . .
05 0 . FIG. 3. (color online). Normalized drag with respectto the total superfluid density for different system’s sizes,as a function of U AB for a system with U = 10 t and ν = 1. From bottom to top (from light to dark shades) L = 8 , , , ,
96 with corresponding bond dimensions χ =100 , , , , n ( s ) AB ) / ( n ( s ) AA + n ( s ) AB ) = 0 . L − for differentvalues of U AB /U = − (0 . , . , . , . , .
35) from bottomto top. to two coupled hydrodynamic Hamiltonians, usually re-ferred to as Luttinger liquids in the 1D context [26]. TheHamiltonian can be diagonalised by introducing the den-sity ( D ) and spin/polarisation ( S ) channels: H µ = 12 π (cid:90) (cid:20) c µ K µ ( ∂ x φ µ ) + c µ K µ ( ∂ x θ µ ) (cid:21) d x, (5)where φ D ( S ) = ( φ A ± φ B ) / √ θ D ( S ) = ( θ A ± θ B ) / √ c D ( S ) , and K D ( S ) are the so-called Luttinger parameters. There’san additional non-linear coupling between the densitiesof the two Luttinger liquids which can be perturbativelyaccounted for by a term proportional to U AB cos(2 √ θ S ).This term is irrelevant in the 2SF phase and relevant inthe PSF phase. As long as the Hamiltonian Eq. (5) holds,an algebraic decay characterises the correlation functions(a.k.a. quasi-long-range order) [26]: R D ( x ) = (cid:104) b † i + x,A b † i + x,B b i,B b i,A (cid:105) ∝ / | x | K D ,R S ( x ) = (cid:104) b † i + x,A b i + x,B b † i,B b i,A (cid:105) ∝ / | x | K S , (6)which is easily checked by using the leading termin the long-wavelength field representation b j,α ∝ exp ( iφ α ) [27]. The Luttinger parameters satisfy the re-lation [26]: K S = π (cid:126) χc S / , (7)where χ = (cid:2) ∂ e/∂ ( ν A − ν B ) (cid:3) − is the spin susceptibilitywith e the energy density (the same goes for the densitychannel with the compressibility rather than the suscep-tibility).On the other hand a hydrodynamic approach based onthe energy functional including the collisionless drag [3,10] provides the following relation between the spin-speedof sound and the superfluid densities [28]: c S = 2 n ( s ) AA − n ( s ) AB m ∗ χ . (8)By direct comparison, the following also holds: K S = (cid:118)(cid:117)(cid:117)(cid:116) π (cid:126) (cid:16) n ( s ) AA − n ( s ) AB (cid:17) χ m ∗ . (9)The drag thus appears in the constitutive relations ofthe Luttinger parameters. This fact, often overlookedin literature [22, 29–33], is crucial in obtaining consis-tent results in the perturbative approach of Luttingerliquids. Therefore in writing the two-species LuttingerHamiltonian a term proportional to n ( s ) AB ∂ x φ A ∂ x φ B mustbe included. The latter term should be indeed gener-ated under renormalisation group flow obtained by means . . . . . − . − . − . − . − . − . − . − . . . . . − . − . K S U AB /U P S F tr a n s i t i o n χ [ t / a ] U AB /t P S F tr a n s i t i o n FIG. 4. (color online). Luttinger parameter K S for thespin channel, as obtained from the hydrodynamic relation inEq. (9) (solid line with error shadow) and from the correla-tion functions, Eq. (6) (points with error bars). Points arereported until the algebraic fit makes sense: the shaded re-gion indicates where deviations become sizeable [25]. In thePSF the parameter K S must go to zero. The inset displaysthe behaviour of the susceptibility as function of the interac-tion, as estimated from various system sizes (color code as inFig. 3 except for L = 8 that is omitted). of the operator product expansion as it has been stud-ied for the Josephson tunneling between superfluids (see,e.g., [34] and reference therein). In Fig. 4 we report thevalues of K S obtained by hydrodynamics (Eq. (9)) andby the long-range behaviour of the correlation functions(Eq. (6)). The latter, however, stops to be algebraic andbecomes exponential once the system enters in the PSFphase [22, 25]: we quantify this change by measuringthe deviation from a pure algebraic behaviour in log-log scale (see Supp. Mat. [25] for a detailed discussion).The region where the deviation becomes appreciable is inagreement with the region predicted from the drag sat-uration (Fig. 3). Before the fluctuations region the twoestimates for K S give consistent results.For the sake of comparison with mesoscopic experi-ments accessible with ultra-cold gases in optical lattices,we report in the inset of Fig. 4 also the estimated spinsusceptibility for different system’s sizes. Conclusions:
In conclusion, we provide a numericalestimation of the AB superfluid drag via a Tensor Net-work approach for a 1D Bose mixture on a ring lattice.The drag is enhanced for attractive interactions due tothe relevance of paring correlations in 1D systems. Inparticular close to the 2SF-PSF transition, the entrain-ment can become of the same order of the full super-fluid density, which makes this regime more suitable foran experimental measure. Differently from the Mott-Superfluid case for a single species [34, 35], the ther-modynamic estimate of the PSF transition proves to bechallenging from a numerical point of view due to theslow convergence to the thermodynamic limit for thespin-channel [11]. We stress again that our results arerelevant for ultra-cold gases experiment where tens to afew hundreds of atoms are considered. For this reason werefrain here from attempting to extract a careful finite-size analysis which will be the subject of future work [23].Moreover our analysis suggest that the inclusion of thedrag has fundamental implications in the understandingof any hydrodynamic Luttinger liquid approach to thetwo species Bose-Hubbard model. In particular we showthat AB hydrodynamics provides a reliable expression forthe Luttinger parameters of the spin-channel in terms ofthe superfluid densities and the susceptibility of the sys-tem.
Note added: while completing the present work webecame aware of a DMRG-QMC comparison for thestudy of the pairing properties of one-dimensional hard-core bosons [36]. The results show the slow convergenceto the thermodynamic limit for the spin-channel as in oursoft-core case.
Acknowledgements:
We gratefully acknowledge in-sightful discussions with S. Giorgini, A. Haller, S. Man-mana, J. Nespolo and S. Stringari. M.R. acknowledgessupport from the Deutsche Forschungsgesellschaft (DFG)through the grant OSCAR 277810020 (RI 2345/2-1),the EU through the project PASQUANS, the Alexan-der von Humboldt foundation and the kind hospitalityof the BEC Center in Trento, where great part of thisstudy was performed. D.C. acknowledges hospitality atthe Johannes Gutenberg University Mainz, where pre-liminary stages of this work were conducted. A.R. ac-knowledges financial support from the Italian MIUR un-der the PRIN2017 project CEnTraL, the Provincia Au-tonoma di Trento and the Q@TN initiative. The MPSsimulations were run on the Mogon Cluster of the Jo-hannes Gutenberg-Universit¨at Mainz (made available bythe CSM and AHRP) and on the JURECA Cluster atthe Forschungszentrum J¨ulich, with a code based on aflexible Abelian Symmetric Tensor Networks Library, de-veloped in collaboration with the group of S. Montangero(Padua). [1] J. M. Lattimer and M. Prakash, Science , 536–542(2004).[2] Ji-Min Duan and Sungkit Yip, “Supercurrent drag viathe coulomb interaction,” Phys. Rev. Lett. , 3647–3650(1993).[3] A. Andreev and E. Bashkin, “Three-velocity hydrody-namics of superfluid solutions,” Sov. Phys.JETP , 164(1976).[4] I. M. Khalatnikov, “Hydrodynamics of Solutions of TwoSuperfluid Liquids,” Sov. Phys. JETP , 542–545 (1957).[5] I. M. Khalatnikov, “Hydrodynamics of Solutions of TwoSuperfluid Liquids,” JETP Lett. , 386 (1973). [6] Z. M. Galasiewicz, “Hydrodynamics of Solutions of TwoSuperfluid Liquids,” Phys. Lett. A , 149 (1973).[7] V. P. Mineev, “Hydrodynamics of Solutions of Two Su-perfluid Liquids,” Sov. Phys. JETP , 338 (1974).[8] D. V. Fil and S. I. Shevchenko, “Nondissipative drag ofsuperflow in a two-component bose gas,” Phys. Rev. A , 013616 (2005).[9] J. Linder and A. Sudbø, “Calculation of drag and super-fluid velocity from the microscopic parameters and ex-citation energies of a two-component bose-einstein con-densate in an optical lattice,” Phys. Rev. A , 063610(2009).[10] J. Nespolo, G.E. Astrakharchik, and A. Recati,“Andreev-bashkin effect in superfluid cold gases mix-tures,” New Journal of Physics , 125005 (2017).[11] L. Parisi, G. E. Astrakharchik, and S. Giorgini, “Spindynamics and andreev-bashkin effect in mixtures of one-dimensional bose gases,” Phys. Rev. Lett. , 025302(2018).[12] K. Sellin and E. Babaev, “Superfluid drag in the two-component bose-hubbard model,” Phys. Rev. B ,094517 (2018).[13] D. Romito, C. Lobo, and A. Recati, “Linear Re-sponse Study of Collisionless Spin Drag,” , 1–12 (2020),arXiv:2002.03955.[14] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wol-swijk, F. Minardi, M. Modugno, G. Modugno, M. In-guscio, and M. Fattori, “Self-bound quantum dropletsof atomic mixtures in free space,” Phys. Rev. Lett. ,235301 (2018).[15] G. Ferioli, G. Semeghini, L. Masi, G. Giusti, G. Mod-ugno, M. Inguscio, A. Gallem´ı, A. Recati, and M. Fat-tori, “Collisions of self-bound quantum droplets,” Phys.Rev. Lett. , 090401 (2019).[16] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas,P. Cheiney, and L. Tarruell, “Quantum liquid dropletsin a mixture of Bose-Einstein condensates,” Science ,301–304 (2018).[17] L. Tanzi, C. R. Cabrera, J. Sanz, P. Cheiney, M. Tomza,and L. Tarruell, “Feshbach resonances in potassium bose-bose mixtures,” Phys. Rev. A , 062712 (2018).[18] C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich,F. Ancilotto, M. Modugno, F. Minardi, and C. Fort,“Observation of quantum droplets in a heteronuclearbosonic mixture,” Phys. Rev. Research , 033155 (2019).[19] A. J. Leggett, “Topics in the theory of helium,” PhysicaFennica , 125–170 (1973).[20] A. B. Kuklov and B. V. Svistunov, “Counterflow super-fluidity of two-species ultracold atoms in a commensurateoptical lattice,” Phys. Rev. Lett. , 100401 (2003).[21] A. Kuklov, N. Prokof’ev, and B. Svistunov, “Com-mensurate two-component bosons in an optical lattice:Ground state phase diagram,” Phys. Rev. Lett. ,050402 (2004).[22] A. Hu, L. Mathey, I. Danshita, E. Tiesinga, Carl J.Williams, and C. W. Clark, “Counterflow and pairedsuperfluidity in one-dimensional bose mixtures in opticallattices,” Phys. Rev. A , 023619 (2009).[23] D. Contessi, A. Recati, M. Rizzi, in preparation.[24] This phase is sometimes considered the bosonic coun-terpart of the much more relevant phenomenon of paircondensation occurring in fermionic systems.[25] See Supplemental Material at . . . . [26] T. Giamarchi, Quantum Physics in One Dimension (Ox-ford University Press, Oxford, UK, 2004).[27] F. D. M. Haldane, “Effective harmonic-fluid approach tolow-energy properties of one-dimensional quantum flu-ids,” Phys. Rev. Lett. , 1840–1843 (1981).[28] Let us remind that the in a Galilean invariant systemsone has that at T = 0 the total superfluid density is equalto the total density, i.e., 2( n ( s ) AA + n ( s ) AB ) = n , and thereforethe numerator of Eq. 8 would be fully determined by thedrag [10].[29] A. Kleine, C. Kollath, I. P. McCulloch, T. Giamarchi,and U. Schollw¨ock, Phys. Rev. A , 013607 (2008).[30] A Kleine, C Kollath, I P McCulloch, T Giamarchi,and U Schollw¨ock, “Excitations in two-component bosegases,” New Journal of Physics , 045025 (2008).[31] L. Mathey, I. Danshita, and C.W. Clark, “Creating asupersolid in one-dimensional bose mixtures,” Phys. Rev.A , 011602 (2009).[32] E. Orignac, R. Citro, M. Di Dio, and S. De Palo, “Vortexlattice melting in a boson ladder in an artificial gaugefield,” Phys. Rev. B , 014518 (2017).[33] R. Citro, S. De Palo, M. Di Dio, and E. Orignac, “Quan-tum phase transitions of a two-leg bosonic ladder in anartificial gauge field,” Phys. Rev. B , 174523 (2018).[34] L. Benfatto, C. Castellani, and T. Giamarchi, “Berezin-skii–kosterlitz–thouless transition within the sine-gordonapproach: The role of the vortex-core energy,” in (2013)pp. 161–199.[35] M Gerster, M Rizzi, F Tschirsich, P Silvi, R Fazio, andS Montangero, “Superfluid density and quasi-long-rangeorder in the one-dimensional disordered Bose-Hubbardmodel,” New Journal of Physics , 15015 (2016).[36] B. Gr´emaud and G. G. Batrouni, “Pairing and pair su-perfluid density in one-dimensional hubbard models,”(2020), arXiv:2006.13950.[37] P. Silvi, F. Tschirsich, M. Gerster, J. J¨unemann,D. Jaschke, M. Rizzi, and S. Montangero, “The TensorNetworks Anthology: Simulation techniques for many-body quantum lattice systems,” SciPost Phys. Lect.Notes , 8 (2019).[38] M. Haller, A.and Rizzi and M. Filippone, “Drude weightincrease by orbital and repulsive interactions in fermionicladders,” Phys. Rev. Research , 023058 (2020).[39] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac,and M. Rigol, “One dimensional bosons: From condensedmatter systems to ultracold gases,” Rev. Mod. Phys. ,1405–1466 (2011).[40] U. Schollw¨ock, “The density-matrix renormalizationgroup in the age of matrix product states,” Annals ofPhysics , 96–192 (2011).[41] Rom´an Or´us, “A practical introduction to tensor net-works: Matrix product states and projected entangledpair states,” Annals of Physics , 117 – 158 (2014). Supplemental material
Currents on the ring
The number current density is related to the numberdensity through the continuity equation. In a stationaryone-dimensional system, it translates to: (cid:42) dˆ j α ( x )d x (cid:43) = − (cid:28) ∂ ˆ n α ( x, t ) ∂t (cid:29) = i (cid:126) (cid:68) [ˆ n x,α , ˆ H ] (cid:69) . (10)in which ˆ n x,α = b † x,α b x,α is the local number density op-erator. The computation of the commutator yields[ˆ n α , ˆ H ] = − tL L (cid:88) x =1 [( b † x +1 ,α b x,α e − π i φαL − h.c. )+ − ( b † x,α b x − ,α e − π i φαL − h.c. )] , (11)for a unitary lattice spacing. Considering the finite differ-ence between the currents flowing on neighbouring linksalong the ring we get: (cid:68) ˆ j x + ,α − ˆ j x − ,α (cid:69) = 2 t (cid:126) (cid:16) Im (cid:104) b † x +1 ,α b x,α e − π i φ α /L (cid:105) + − Im (cid:104) b † x,α b x − ,α e − π i φ α /L (cid:105) (cid:17) , (12)that leads to Eq. (3) by identifying, thanks to transla-tional invariance (cid:68) ˆ j x + ,α (cid:69) = j α .Fig. 5 illustrates the effect of the AB drag on thespecies’ (super)currents: in presence of φ A only, the cur-rent j B is constantly zero in absence of inter-species in-teraction, U AB = 0, while it increases monotonically as U AB increases. The drag density n ( s ) AB is proportional tothe slope of the j B curve in the limit of φ A →
0. Theplot highlights the smallness of the drag effect at a genericpoint in parameter space – here the same ν = 0 . U/t = 2 and L = 32. Since the total mo-mentum is anyway conserved, the total current j A + j B must be independent of U AB , as confirmed by the datain the inset. Realization of periodic boundary conditions
Periodic Boundary Conditions usually pose a challengeto Tensor Network methods due to two main reasons:the increased amount of entanglement to be accountedfor when bi-partitioning the system in segments (twicethe boundaries), and the presence of loops impairing theexistence of a canonical form if naively mirrored in thenetwork structure [37]. While the first is unavoidable, thelatter could be circumvented by either employing TreeTensor Networks [35, 37] or – as we did here – by mapping − . − . . . − . − . . . U AB U AB U AB U AB j A + j B j α [ t / ¯ h ] φ A FIG. 5. (color online). The superfluid currents j A and j B in presence of φ A only for a L = 32 system. The four setsof curves (from dark to light color) correspond to U AB /U =0 . , . , . , .
75, with
U/t = 2 at half-filling (same regime asFig. 2). The drag density n ( s ) AB is proportional to the slope ofthe j B curve in the limit of φ A →
0. The inset demonstratesthe global momentum conservation via the constant value ofthe total current for all values of U AB . the ring onto a chain with tailored next-nearest neighborcouplings and boundary nearest neighbor ones [38].The loop therefore disappears from the MPS form andgets adsorbed inside the MPO. The basic idea of such“snake” pattern is shown in Fig. 6 for a very little sys-tem with an even number of sites. Despite a slightlybigger MPO and a bit of extra book-keeping in the post-processing of observable measurements, the numericalcalculations demonstrated to be efficient and stable [38]. FIG. 6. (color online). Sketch of the scheme used to mimic pe-riodic boundary conditions with a computational open bound-ary conditions ansatz.
Numerical derivatives and errors
The first source of error in the MPS ansatz that hasbeen taken into account is the maximum local bosonicoccupation number. The maximum total (species A andB) number of particles at site x is chosen to be N maxx =6. Because of the low filling of the computed systems,we checked that the probabilities of the high-occupationstates beyond our truncation are negligible.In order to compute the superfluid densities as inEq. (4), we resorted to a four-point numerical approxima-tion of the derivative, while the currents have been mea-sured from the MPS wavefunction according to Eq. (3): ∂j α ∂φ β (cid:12)(cid:12)(cid:12)(cid:12) (cid:39) ¯ j α [ − δ β ] − j α [ − δ β ] + 8 ¯ j α [ δ β ] − ¯ j α [2 δ β ]12 δ β + O ( δ β )(13)where the mean current of species α over the links ofthe lattice is ¯ j α = 1 /L (cid:80) Lx =1 j x +1 / ,α . The O ( δ β ) errordescending from the finite-difference derivative is negli-gible with respect to the one associated with the inho-mogeneity of the currents j x +1 / ,α along the ring, shownin Fig. 7. We have estimated the latter via the standarddeviation σ of the measured values and then propagatedthe uncertainty in Eq. (13). Moreover, the resulting er-ror has been compared and eventually combined withthe error estimation due to the MPS ansatz. This hasbeen accounted for by comparing the values of the finalquantity at different bond dimensions’ values. The wholeprocedure leads to the error bars reported in Figs. 2-3.Concerning the numerical error on the energy, we madeuse of the standard deviation of the ground state energyestimations along the last sweep of the MPS optimization(e.g., Chapter 5 of [37]). This should be closely relatedto the accumulation of the truncated probabilities in therenormalization procedure.Lastly, some presented quantities are the thermody-namic estimations unless explicitly specified. The errorsresult from the fit of those quantities as function of theinverse of the system’s size L . More details on the ex-traction procedure of the Luttinger parameters from thecorrelation functions are provided in the following sec-tion. Extraction of the Luttinger parameters from thecorrelation functions
The Luttinger parameters K D,S are extracted by ex-ploiting the quasi-long-range order behaviour of the fol-lowing correlation functions: G α ( x ) = (cid:104) b † i + x,α b i,α (cid:105) ∝ | d | − KD − KS ,R D ( x ) = (cid:104) b † i + x,A b † i + x,B b i,B b i,A (cid:105) ∝ | d | − KD ,R S ( x ) = (cid:104) b † i + x,A b i + x,B b † i,B b i,A (cid:105) ∝ | d | − KS . (14) j x + / , A [ − t / ¯ h ] L = 64 ) . σ σ FIG. 7. (color online). The j x + a/ ,A current as function of x , indicated by the site labels. The system’s configuration is: t = 1, U = 10, U AB = − ν = 1, φ A = 0 .
002 and φ B = 0.The red solid line is the average of the values and the shadedareas represent 1 σ and 1 . σ error respectively. Notice thatalmost all the values are included in the 1 σ region that is wellinside the 0 .
1% error (dashed blue line).
Here we expressed the algebraic decay in terms of thenatural measure of the distances between sites on a ringgeometry, i.e., the chord function [39]: d ( x/L ) = Lπ sin (cid:16) πxL (cid:17) , (15)where L is the number of sites and x ∈ N the lineardistance between the sites. For very large rings the ex-pressions further simplifies as in Eq. (6) of the main text.The single-body correlations, G α , have a mixed den-sity/spin character, consistently with the fact that theimaginary part of their nearest-neighbor value gives backthe species current j α as in Eq. (3). The two contribu-tions can be instead isolated with the help of two-bodycorrelations: R D concerns the superfluid character ofpairs of A − B particles and therefore the density chan-nel, while R S relates to particle-hole pairs and thereforethe spin channel [12, 22].Away from commensurate effects, possibly leading to aMott insulator, the density channel is always superfluid,i.e., R D scales algebraically. A change in R S and G α fromalgebraic to exponential decay – or equivalently a drop of K S to 0 – happens instead when entering the PSF phase,due to the opening of a gap in the spin channel. This isillustrated in Fig. 8, where the correlations measured fordifferent system sizes ( L = 8 , , , ,
96) are reportedfor two sample parameter values deep in the 2SF (blue)and PSF (orange) phases.In order to compute the Luttinger parameter K S re-ported in Fig. 4, we employed the following procedure:First, we extracted at each finite size L the exponents of − − − − − h ( a † ) ( a ) x i d ( x/L ) Algebraic fitExponential fit U AB /t = − (2SF) U A B / t = − ( P S F ) (a) − − h ( a † b † ) ( a b ) x i d ( x/L ) Algebraic fitAlgebraic fit U A B / t = − ( S F ) U A B / t = − ( P S F ) (b) − − − − − h ( a † b ) ( a b † ) x i d ( x/L ) Algebraic fitExponential fit U A B / t = − ( S F ) U A B / t = − ( P S F ) (c) FIG. 8. (color online). The correlation functions of Eq. (14) (the (a) panel is G α , (b) is R D and (c) is R S ) as function of theconformal distance for a unitary total filling system ν = 1, t = 1 and U = 10. The orange points concern a regime in which thesystem is in PSF phase while the blue ones concern the 2SF phase and we represent with a color gradient from dark to lightdifferent system’s sizes from L = 8 to L = 64. The dashed lines are exponential and algebraic fits depending on the expectedbehaviour of the functions for U AB /U = − . U AB /U = − .
5. The − − − h ( a † b ) ( a b † ) x i d ( x/L ) FIG. 9. (color online). The effect of the filtering procedurefor the R S correlations for a L = 96 ring with t = 1, U =10 t , U AB = − . U . The dark-blue points include all thecorrelations extracted from the MPS ansatz and the light-blue are the ones left from discarding those which are too farfrom a computational point of view thanks to the procedureexplained in the main text. The algebraic fit is representedwith a solid black line and the area between the two dashedblack lines is the error on the fit. the algebraic fits in Eq. (14) and checked their mutualconsistency. Then, we extrapolated the thermodynamiclimit of K S ( L ) as function of 1 /L , as we also did for otherquantities across the manuscript.It is important to recall here that any MPS Ansatzleads to correlation functions expressed as a sum of expo-nential functions (as dictated by the underlying transfermatrix structure, see [40, 41]), which could mimic even al-gebraic decays in finite systems. Due to our specific MPSstructure of Fig. 6, however, the same physical distancemight be represented by pairs of points at a pretty dif- ferent computational distance. Therefore, some of thesepairs return more precise estimates than others: for thisreason, when performing the (algebraic) fits, we filteredout data coming from pairs further than the half of thering length ( x > L/
2) in the computational basis. Theresult is shown in Fig. 9, where it is evident that thefiltering allows for cleaner fits. − − . − − . − − . − − . − − . − − . . . . . . . . . l og h ( a † b ) ( a b † ) x i log( d ( x/L )) FIG. 10. (color online). Fit of the correlations log[ R S ( x )] asfunction of the logarithm of the conformal distance with afitting function as Eq. (16). The configuration of the systemis the same of Fig. 9 except for the interpecies interactions:red-like point are for U AB /U = − . U AB /U = − .
1. The correspondent fitting functions (solidand dashed lines respectively) consider only the dark pointswith correspondent error bars. The bending of the correla-tions from a pure algebraic trend can be captured by increas-ing the α amplitude. Finally, in order to better locate the 2SF to PSF tran-sition region, we resorted as in [35] to an estimate ofthe curvature of the fitting function in the log-log plane,0quantified by α :log[ R S ( x )] = αξ − K S ξ + const ., (16)where ξ = log[ d ( x/L )]. If α is zero the fit is purely alge-braic and the system is in the 2SF phase (see dashed blueline in Fig. 10). Otherwise, if α assumes a finite value, thefit bends (down) from the algebraic trend and hence itsignals an exponential kicking in. This allows us in princi-ple not to make any assumption on a critical value K ∗ S atwhich the phase transition happens. In Fig. 11 we showthe α values as a function of the inter-species interaction U AB , for various system sizes – with the same color codeas in Fig. 3 ( U/t = 10). We can clearly spot two regimes:for U AB /U ≥ − . α appears to converge towards zerofor larger system sizes, while for U AB /U ≤ − .
25 the(negative) curvature α increases in amplitude with grow-ing system sizes. This latter behaviour hints at the factthat an algebraic fit of R S ceases to make sense beyondthis region, which we therefore identify as the one wherethe 2SF-PSF transition takes place.Larger system sizes and a more accurate sampling ofthe putative transition region would be needed in orderto precisely locate the critical point – e.g., the thermo- dynamic limit might emerge for a very large number ofparticles [11]. We leave such a more complete study ofthe transition for future work. α U AB /U − . − . − . − − . − . − . − . − . − . − . − . − . − . − . − . − . PSF