Dissipationless Vector Drag--Superfluid Spin Hall Effect

Dissipationless flows in single-component superfluids have a significant degree of universality. In He4, the dissipationless mass flow occurs with a superfluid velocity determined by the gradient of the superfluid phase. However, in interacting superfluid mixtures, principally new effects appear. In this Letter, we demonstrate a new kind of dissipationless phenomenon arising in mixtures of interacting bosons in optical lattices. We point out that for a particular class of optical lattices, bosons condense in a state where one of the components' superflow results in dissipationless mass flow of the other component, in a direction different from either of the components' superfluid velocities. The free-energy density of these systems contains a vector product-like interaction of superfluid velocities, producing the dissipationless noncollinear entrainment. The effect represents a superfluid counterpart of the Spin Hall effect.

Dissipationless flows in single-component superfluids have a significant degree of universality. In 4 He, the dissipationless mass flow occurs with a superfluid velocity determined by the gradient of the superfluid phase. However, in interacting superfluid mixtures, principally new effects appear. In this Letter, we demonstrate a new kind of dissipationless phenomenon arising in mixtures of interacting bosons in optical lattices. We point out that for a particular class of optical lattices, bosons condense in a state where one of the components' superflow results in dissipationless mass flow of the other component, in a direction different from either of the components' superfluid velocities. The free-energy density of these systems contains a vector product-like interaction of superfluid velocities, producing the dissipationless noncollinear entrainment. The effect represents a superfluid counterpart of the Spin Hall effect.
As shown in [1], in the presence of intercomponent interactions, the free-energy density describing a binary superfluid mixture should necessarily include a scalar product of the superfluid velocities v a and v b , i.e., f = ρ a v 2 a /2 + ρ b v 2 b /2 + ρ ab v a · v b . The resulting superflow j a = ∂f /∂v a = ρ a v a + ρ ab v b indicates that even if v a = 0, there will still be a nonzero superflow of component a with the superfluid velocity v b .
In the condensed matter context, the effect became of great interest with the advent of optical lattices, which allow for precise control of strongly correlated superfluids [14,15]. The strength of the Andreev-Bashkin drag is controlled by the optical lattice parameters in combination with on-site interactions [16]. It was shown that the effect, in relative terms, can be arbitrarily strong and that the Andreev-Bashkin drag coefficient ρ ab can also become negative. In the latter case, one deals with a counterflow, where the flow of one component generates a mass flow of the other component in the opposite direction [8][9][10][11]. It was pointed out that the effect should lead to the formation of new superfluid states where only dissipationless coflow (paired superfluids) or only counterflow (supercounterfluids) can exist [8][9][10][11][12][13][16][17][18]. At the same time, even relatively weak drag substantially changes rotational responses [19,20].
In this Letter, we demonstrate the existence of a new dissipationless phenomenon in superfluid mixtures, where the dissipationless superfluid entrainment is not collinear with superfluid velocities. Namely, we show that the superflow-superflow interaction can, in general, be described by a nontrivial tensor ρ ij αβ which enters the bilinear free-energy density: Here Greek subscripts and Roman superscripts, respectively, label components and Cartesian directions of the superfluid velocity vector v i α and the superfluid stiffness tensor ρ ij αβ = ρ ji βα . The latter describes both kinetic (α = β) and drag (α = β) phenomena. While physical properties of the system are encoded in the tensor ρ ij αβ , as we elaborate below, a direct interpretation of individual coefficients may be deceptive since they depend on the choice of the coordinate system.
In this Letter we consider a two-dimensional twocomponent system-i, j ∈ {x, y} and α, β ∈ {a, b}-and study the drag-related elements ρ ij ab . In such a case, the quantities are coordinate system independent [21]. All other pairwise combinations vary under rotation, and it turns out that it is always possible to find a Cartesian coordinate system in which ρ xy ab + ρ yx ab = 0. If, in addition, the difference ρ xx ab − ρ yy ab is negligible-which in principle can be guaranteed in certain situations-one finds where f 0 = α ij ρ ij αα v i α v j α /2 represents the standard kinetic contribution to the free-energy density. The corresponding superflows read where j 0α = ∂f 0 /∂v α and e i denotes the unit vector in the i th direction. The existence of the vector product-like arXiv:2103.16531v2 [cond-mat.quant-gas] 27 Aug 2021 contribution to f given by ρ ⊥ constitutes a novel superfluid effect where a nonzero superflow of one component induces a perpendicular superflow response of the other component. This new phenomenon, which we coin vector drag, may be viewed as an intercomponent-interactiondriven counterpart of the Spin Hall effect [22] originating in the spin-orbit coupling, discussed particularly in hybrid structures involving superconductors [23,24]. Let us now provide both analytical and numerical evidence for the existence of a nonzero ρ ⊥ , and thus the existence of the vector-drag phenomenon. Our starting point in deriving the effective model in Eq. (1) is a twocomponent Bose-Hubbard-type model on a rectangular lattice, given by the Hamiltonian Hereb iα (b † iα ) denotes the bosonic annihilation (creation) operator of component α at site i, andn iα =b † iαb iα is the corresponding particle number operator. The mass of the α-component boson is m α , and the rectangular lattice consists of N × N sites-with lattice vectors a x = l x e x , a y = l y e y and lattice constants l x , l y -on which we impose periodic boundary conditions. For simplicity, we assume constant particle number densities n α and restrict ourselves to on-site and nearest-neighbor interactions whilst allowing for nearest-neighbor and next-nearestneighbor hopping. The problem is amenable to analytical treatment only in the weakly interacting regime. Here, like in the case of the Andreev-Bashkin effect [25][26][27][28][29], the drag effects are expected to be inherently small. We will first demonstrate the existence of vector drag in the weakly interacting regime analytically. Then we investigate the effect in the strongly correlated regime by employing large-scale quantum Monte-Carlo calculations.
The standard Andreev-Bashkin effect can be analytically calculated in macroscopic weakly interacting systems. That was previously done for square and triangular lattices and in a continuum [25][26][27][28]. We begin by employing a similar analytic approach to establish the new phenomenon: the vector drag. Consider a weakly interacting regime where at low enough temperatures both components are condensed into the zero-momentum mode. In such a case, the Hamiltonian (5) can be approximated and subsequently diagonalized in momentum space. The corresponding zero-temperature free-energy density (1) is then obtained as the ground state energy.
As discussed above, the vector drag should be the most transparent when ρ xy ab = −ρ yx ab , which is guaranteed for systems being invariant under a reflection in either of the two lattice vectors combined with an exchange of components a ↔ b [21]. Our aim is to construct a microscopic model which exhibits ρ ⊥ = 0 in addition to satisfying the above-mentioned symmetry. A simple choice of parameters obeying these conditions is m a = m b = m, n a = n b = n, and where r i indicates the i th lattice site position. The resulting model is illustrated in Fig. 1. While in general ρ xx ab = ρ yy ab , it turns out that when both ρ xx ab and ρ yy ab are nonzero and of the same sign-which for the considered parameter region is the case-one can in principle completely eliminate ρ xx ab − ρ yy ab . This is achieved by rescaling the ratio l x /l y by ρ yy ab /ρ xx ab while at the same time keeping all other model parameters fixed, which leaves the l x , l y -independent ρ ⊥ unchanged [21]. In this way we can realize the model effectively described by Eq. (3). We find that for the lattice geometry illustrated in Fig. 1(a), this rescaling leads to l x < l y in the region of interest, which is consistent with keeping the nearest-neighbor interactions in the x direction only-typical interatomic interactions rapidly decay with increased separation distance. Nevertheless, omitting the nearest-neighbor interactions along the y direction is merely a simplification to reduce the number of system parameters and does not change the main result. Namely, we would like to stress that the presence of vector drag is not limited to systems having nearest-neighbor interaction in one direction only. In what follows we will restrict ourselves to m = 1, n = 1/2, U = 1, U = 0.9, and V = 0.9V .
Our findings indicate that in the weakly interacting , (e), (h) shows the lx/ly-dependent collinear-drag coefficient ρ , computed using the ratio lx/ly plotted in the right column (c), (f), (i), which is adjusted such that ρ xx ab = ρ yy ab required for the system to be described by Eq. (3). The upper (a)-(c) and middle (d)-(f) row display how the vector drag ρ ⊥ , Andreev-Bashkin drag ρ , and lx/ly depend on t/U and t /U when V = 0.2 and V = −0.2, respectively. Even though the presence of the nearest-neighbor interactions modifies the amplitude of ρ , it does not change its character determined by the sign. On the contrary, both the sign and magnitude of the vector drag strongly depend on the sign and magnitude of V . Note that the sign change of V affects the ratio ρ yy ab /ρ xx ab implying differences between lx/ly visible in (c) and (f). The same quantities, ρ ⊥ , ρ and lx/ly, but for t/U = t /U and V = ±0.1, ±0.2, are presented in the bottom row (g)-(i).
regime, the macroscopic system harbors a substantial vector drag when V = 0, at least within our approach. In Fig. 2 we present the analytically derived drag coefficients ρ ⊥ and ρ versus t/U and t /U , for different values of V , in panels (a), (d), and (g) and in (b), (e), and (h), respectively. The latter quantity is calculated for l x /l y adjusted such that ρ xx ab = ρ yy ab , which entails l x < l y apart from the region where t t, see panels (c), (f), and (i) of Fig. 2 . Since both the magnitude and sign of ρ ⊥ are determined by the magnitude and sign of V , the nearestneighbor interactions are important for the vector-drag phenomenon in the considered regime.
After demonstrating the effect analytically in the weak- 56. An alternative ratio lx/ly-for which (ρ xx ab − ρ yy ab )/2 is completely eliminated-and the corresponding ρ are plotted in panels (d) and (e), respectively. In the latter panels, the missing data points around t/U ≈ 0.08 are due to ρ xx ab and ρ yy ab having opposite signs. All results were obtained using V /U = 0, 0.1, 0.2, t/U = t /U , N = 10, and β = N/t.
coupling regime, we proceeded to study the model, Eqs. (5) and (6), deep inside the strongly correlated regime. For this purpose, we performed large-scale wormalgorithm Monte Carlo simulations [30,31]. Using a generalization of the Pollock-Ceperley formula [21], which in its original form allows computing the superfluid density through winding number statistics [32], we obtained the drag-related elements of the more general tensor: Here β is the inverse temperature and w i a w j b are winding number correlations-the winding number w i α is the net number of times α-type particles cross the periodic boundary along the i th direction. Similarly as in the weakly interacting case, we consider the system with constant particle densities equal for both components, i.e., n a = n b = n = 1/2.
The results obtained for the parameters N = 10, t = t , β = N/t, and V /U = 0, 0.1, 0.2 are presented in Fig. 3. Here ρ ⊥ , ρ , and (ρ xx ab − ρ yy ab )/2 are plotted against t/U in panels (a)-(c), respectively. At small values of t/U , the system is in an insulating phase which is evident from the coefficients being identically zero. Then, at t/U ≈ 0.053, the system undergoes a transition into a superfluid state with negative drag coefficients. By further increasing t/U , ρ ⊥ and ρ increase in value and eventually change sign, whereas (ρ xx ab − ρ yy ab )/2 reveals a more complex dependence on t/U . The behavior of the Andreev-Bashkin drag ρ is very similar to that of the conventional twodimensional two-component Bose-Hubbard model [16]. Note that the magnitudes of both ρ and ρ ⊥ are maximal in the strongly correlated regime very close to the phase transition point, which is similar to the fact that in the previously studied systems the Andreev-Bashkin drag is also maximal close to insulating phases [16].
For l x /l y = 0.56, the difference (ρ xx ab − ρ yy ab )/2 is comparable to ρ ⊥ in magnitude and should not be ignored. Nevertheless, by adjusting the l x /l y ratio as previously mentioned we can prompt ρ xx ab = ρ yy ab in situations where 0 < ρ yy ab /ρ xx ab < ∞. This l x /l y ratio, as a function of t/U , is shown in panel (d) of Fig. 3, and the resulting ρ can be found in panel (e). In the vicinity of t/U ≈ 0.08 the ratio ρ yy ab /ρ xx ab does not fulfill the above-mentioned criteria which is why there are missing data points.
The analytical result for the weakly interacting macroscopic limit, and the unbiased numerical data for the strongly correlated regime, both demonstrate the new type of dissipationless transport manifested in the nonzero value of the vector-drag coefficient ρ ⊥ . In the strong coupling regime, substantial nonzero vector drag ρ ⊥ exists in case of exclusively on-site interactions, i.e., for V = 0, for the system size considered. The vectordrag coefficient's dependence on the system size with purely on-site interactions is further discussed in the Supplemental Material [21], where it is demonstrated that the effect is still present also in larger systems.
In summary, we demonstrated a new dissipationless transport phenomenon in superfluid mixtures in optical lattices. Namely, we have shown that for a class of optical lattices, the free-energy density of an interacting superfluid mixture should contain a vector product-like interaction between the superfluid velocities. This implies that in such a mixture the superflow of each component is not collinear with either of the components' superfluid velocity, i.e., a superflow of one of the components induces a superflow of the other component in the orthogonal direction, in addition to the standard drag. It may be viewed as an intercomponent-interaction-induced superfluid counterpart of the Spin Hall effect.
The strength of the effect was investigated both ana-lytically in the weakly interacting macroscopic limit and numerically in the strongly correlated regime using largescale worm-algorithm quantum Monte-Carlo simulations. One way to realize this model is through a bilayer optical lattice loaded with dipolar bosons [33,34]. The effect should also be present in multicomponent superconductors. There, at the level of the Ginzburg-Landau model, it should manifest through the presence of terms in the form of vector-like product of components of supercurrents, i.e., mixed terms fourth order in fields and second order in gradients. Investigation of these aspects will be presented in future studies.

DIAGONALIZATION OF THE HAMILTONIAN
We consider a two-component Bose-Hubbard-type model on a two-dimensional rectangular lattice consisting of N × N sites with periodic boundary conditions. The Hamiltonian readŝ where α, β ∈ {a, b} whilst i, j run over all lattice sites. Hereb iα (b † iα ) is the bosonic annihilation (creation) operator of component α at site i, andn iα =b † iαb iα is the corresponding particle number operator. We have chosen to restrict the model parameters to where σ = ±1 and r i indicates the i th lattice site position, while l s and e s are the lattice constant and unit vector in the s ∈ {x, y} direction, respectively. We switch to the momentum representation by Fouriertransforming the operatorsb jα = kb kα e −ik·rj /N with k = k x e x + k y e y . In the macroscopic limit at a sufficiently low temperature, when each component macroscopically occupy the k = 0 mode, one can replacê b 0α ,b † 0α → b † 0αb 0α 1/2 = N √ n 0α and express the particle number density as n α = n 0α + k =0b † kαb kα /N 2 . By omitting terms containing more than two of the remainingb † k =0α andb k =0α , we finally obtain Here H 0 is a constant and introducing ξ s = k s l s , U αβ = U αβ √ n α n β , and V s αβ = V s αβ √ n α n β we defined Since a superfluid flow is related to a phase gradient along the lattice, we may incroporate the corresponding superfluid velocity v α = v x α e x + v y α e y in Peierls-like factors, i.e., t ijα → t ijα e −i∆φ ij α . Here where m α is the α-type particles mass, = 1, and the integration taken along a straight line between i th and j th lattice sites. Note that since ∆φ ii α = 0, the introduced phase nonuniformity modifies only the kinetic part of the Hamiltonian and can be incorporated by the substitution k → k − m α v α in Eq. (S6). The Hamiltonian (S4) can be diagonalized similarly as in Refs. [1][2][3][4] and cast into the form where H 0 is a constant and β k,± are bosonic annihilation operators of ±-type quasiparticles. To do so, it is convenient to rewrite Eq. (S4) into the matrix form where the following basis vectors and ε kα = ε kα + U kαα are introduced. The Hamilto-nianĤ can be diagonalized via a unitary transformation O k with the new basis vectors Ψ k =Ô † k ψ k . The requirement that the new operators β kσ , β † kσ stored in Ψ k satisfy the standard bosonic commutation relations leads toÔ † kσ 3Ôk =σ 3 , wherě with 1 n being the n-dimensional identity matrix. In result, O k diagonalizes M k = M kσ3 where the corresponding eigenvalues can be determined by solving the characteristic equation |M k − E1 8 | = 0. Consequently, one finds four eigenenergies ±E k,σ=± where the positive ones describe quasiparticle excitations (see also [1][2][3][4]).

FREE ENERGY EXPANSION
The free-energy density of the system at temperature T = 0K is given by an expectation value of the Hamiltonian (S8) in a state free of quasiparticle excitations [1,3,4]. Since the superfluid velocities are assumed to be small, the eigenenergies E k,σ can be determined by expansion in terms of small v a and v b , where in the zeroth order E kab , see also Refs. [1][2][3][4] The subsequent higher order terms in the expansion can be found by expanding ε kα with k → k − m α v α in terms of small v α and solving the characteristic equation order by order. Note, that in order to prevent spatial collapse or separation between atomic clouds, the interactions have to satisfy the relation U ka U kb > U 2 kab . For small superfluid velocities the free-energy density can-up to a constant-be cast into the following form where ρ αβ := ij ρ ij αβ e i e T j , and the drag-related densities are found to be given by with G(k) = 2 ε ka ε kb U 2 kab /E vanishes when U kab → 0. Note that both G(k) and ε kα can be viewed as functions of ξ s which is l s independent. Therefore, the values of Andreev-Bashkin terms ρ ii ab can be simply modified by adjusting the ratio l x /l y . The drag-related elements (S16) are calculated in the thermodynamic limit where we transit from summation to integration, V −1 k =0 I(k) → (2π) −2 1BZ I(k)d 2 k, noting that the integrand I(k) vanishes for k = 0.

GENERALIZATION OF POLLOCK-CEPERLEY EQUATION
We will consider a n-component superfluid inhabiting a d-dimensional lattice with periodic boundary conditions. The set of lattice vectors {a i } i=1,...,d need in general not be orthonormal, i.e., a i · a j = δ ij and |a i | = l i = 1, where l i are the corresponding lattice constants. The number of lattice sites along the direction of the i th lattice vector is given by N i , and the corresponding side length is therefore L i = N i l i . Introducing the matrices N := i N i e i e T i and M := i a i e T i , where {e i } i=1,...,d is a set of orthonormal coordinate vectors (here e 1 = e x , e 2 = e y , etc.), the lattice volume may be expressed as V = det(N )|det(M)|.
Following the derivation in [5], we set the superfluid components in motion via a component-dependent Galilean transformation. Compared to a motionless system, this will introduce Peierels phase factors in the hopping amplitudes, i.e., t ij α → t ij α e −i∆φ ij α , with ∆φ ij α given by (S7). For the purpose of this derivation, it is sufficient to consider a uniform velocity field such that ∆φ ij α = m α v α · (r i − r j ). Due to the periodic boundary conditions, each term in the partition function gains a factor exp[−im α v α · (±N i a i )] for each αcomponent particle crossing the boundary in the ±a i /|a i | direction. Introducing the winding number w i α ∈ Z, as the flux of type α particles through the boundary perpendicular to a i , the net phase factor may be ex- Next, we decompose the partition function Z(v α ) in terms of fixedwinding-number partition functions Z {wα} , i.e., Z(v α ) = {wα} Z {wα} e iθ(vα,wα) . In this notation, the partition function of the motionless system corresponds to Z(v α = 0), and the ratio Z(v α )/Z(0) is identified to be the expectation value of the phase factor, namely e iθ(vα,wα) . For small velocities v α , this quantity may be approximated as where the winding-number correlation tensor W αβ := w α w T β has been introduced. The linear terms w i α vanish since the microscopic model is assumed to be invariant under a parity inversion. In addition, the ratio Z(v α )/Z(0) = e −β∆F is also related to the free-energy difference ∆F = F (v α ) − F (0) of having set the system in motion. Here β is the inverse temperature, and due to the previous assumption of uniform velocity fields, the free-energy difference may be expressed as the free-energy density (S15) times the volume, i.e., ∆F = Vf . In the limit of small velocities v α , the free-energy exponential may be accurately approximated by the two first terms of its series expansion, Combining Eq. (S17) and Eq. (S18), we obtain which generalizes the relation which was first derived in [6].
We find that in the case of a N ×N rectangular lattice, the drag coefficients are given by which is an analog of Eq. (S16) but here valid also in the strongly interacting regime.

WORM-ALGORITHM MONTE CARLO
In this Letter, we use continuous-time worm-algorithm Monte Carlo [7] to extract completely unbiased windingnumber correlations for the two-component Bose-Hubbard type model, described in the main text. Wormalgorithm Monte Carlo efficiently samples diagonal elements of the density matrix at thermal equilibrium as world-line configurations in real space and imaginary time. This is achieved by combining the partition function sector with the Green's function sector, so that, in order to go from one partition function world-line configuration to another, one needs to pass through the Green's function sector.
The d + 1 dimensional system containing the worldlines is periodic in the spatial dimensions due to imposed periodic boundary conditions. It is also periodic in imaginary time to reflect the cyclic property of the trace, present in the density matrix being sampled. We use separate worms for the two components, which are able to change their winding numbers w i α by winding in spatial dimensions. In contrast, the winding number in imaginary time determines the total particle number, and since we here consider a fixed number of particles, the latter type of winding is forbidden.

SYMMETRY CONSIDERATIONS
The novel dissipationless transport phenomenon is especially transparent when the representation of the interspecies interaction of Eq. (S15) reads Here ρ := (ρ xx ab +ρ yy ab )/2 and ρ ⊥ := (ρ xy ab −ρ yx ab )/2 are coordinate system independent drag-coefficients. To achieve this, we first consider rotation and reflection transformations, which in two dimensions can be represented by the orthogonal matrices respectively. Here c φ = cos φ and s φ = sin φ, so that R 1 rotates the system by an angle θ, whilst R 2 reflects against a line making the angle ϕ with the x-axis. Since matrices transform according to M = R k M R T k (k ∈ {1, 2}), drag coefficients of the superfluid stiffness tensor will in the rotated system r = R k r be given by ρ x x ab + ρ y y ab = ρ xx ab + ρ yy ab ≡ 2 ρ , (S23) ρ x x ab − ρ y y ab = (ρ xx ab − ρ yy ab )c 2θ − (ρ xy ab + ρ yx ab )s 2θ , (S25) ρ x y ab + ρ y x ab = (ρ xy ab + ρ yx ab )c 2θ + (ρ xx ab − ρ yy ab )s 2θ . (S26) From this we note that ρ and ρ ⊥ indeed are invariant under rotations, but more importantly that we need ρ xx ab = ρ yy ab and ρ xy ab = −ρ yx ab to consistently cancel other contributions to f ab . Note, that since the model is homogeneous, the superfluid stiffness tensor and its components are not affected by translation.
The explicit dependence of ρ ii ab on l i , as revealed by Eq. (S16) and Eq. (S20), makes it possible to cancel ρ xx ab − ρ yy ab by scaling the ratio l x /l y by ρ yy ab /ρ xx ab , as long as ρ yy ab and ρ xx ab are nonzero and of the same sign. However, to fix ρ xy ab + ρ yx ab = 0, we will consider a reflection in the x-or y-axis (2ϕ = πn , n ∈ Z), combined with a component exchange operation, i.e., a ↔ b. The former of the two operations yield ρ xx ab ± ρ yy ab → ρ xx ab ± ρ yy ab and ρ xy ab ± ρ yx ab → −(ρ xy ab ± ρ yx ab ), but in combination with the latter we find ρ xx ab ± ρ yy ab → ρ xx ab ± ρ yy ab and ρ xy ab ±ρ yx ab → ∓(ρ xy ab ±ρ yx ab ). Hence ρ , ρ ⊥ , and ρ xx ab +ρ yy ab are invariant with respect to this transformation, whereas the sign of ρ xy ab + ρ yx ab is negated. From the plain reflection we further observe that it is crucial that the system breaks reflection symmetry in both the x-and y-direction as the vector-drag coefficient otherwise will vanish. Thus, by incorporating these symmetries into the model in addition to adjusting the l x /l y ratio, we obtain the desired form of the interacting part of the free-energy density (S21).

LATTICE ASYMMETRY
We can make some illustrative remarks here which should not be mistaken for a rigorous description of the phenomena. The ordinary drag phenomenon arises when elementary collinear hopping processes on a lattice, shown in Fig. S1(b), are not equally likely. However, in such a situation the perpendicular component of the drag vanishes on symmetry grounds: elementary perpendicular hopping processes, i.e., the first particle hops to the right "pushing" the second particle up or down, are equiprobable, Fig. S1(c).
Unlike the traditional symmetric setting for which the vector drag vanishes, our system is asymmetric in respect to interaction-stimulated perpendicular hopping, Fig. S1(d). This results in elementary non-collinear hopping processes being non-equiprobable, Fig. S1(e,f). The asymmetry of the elementary processes leaves the imprint on the superfluid hydrodynamics resulting in a nonzero vector-drag coefficient ρ ⊥ .

FINITE SIZE SCALING
Here we investigate how the vector-drag coefficient ρ ⊥ scale with the system size N , in the absence of nearestneighbor interactions, i.e., V = V = 0. From Fig. S2 we notice a shift in the peak position towards lower values of t/U when N is increased, which initially is accompanied by a decrease in the peak value. However, at N = 20, the trend changes and the peak value starts growing when further increasing N . This suggests that the vector drag in the model with exclusively on-site interactions is not a finite-size effect.  S2. The vector-drag coefficient ρ ⊥ in the absence of nearest-neighbor interaction, i.e., V = V = 0, for different system sizes N , as a function of t/U . The inset shows the peak value for each system size, which shifts towards lower t/U with increased N . Initially the peak values decrease with increased N , but at N = 20 the trend changes: the peak value starts increasing with increased N . The remaining parameters used are t = t , U /U = 0.9, β = N/t , and na = n b = 1/2.