Drag effect and topological complexes in strongly interacting two-component lattice superfluids

The mutual drag in strongly interacting two-component superfluids in optical lattices is discussed. Two competing drag mechanisms are the vacancy-assisted motion and proximity to the quasi-molecular state, in which an integer number $q$ of atoms (or holes) of one component might be bound to one atom (or hole) of the other component. Then the lowest energy topological excitation (vortex or persistent current) becomes a composite object consisting of $q$ circulation quanta of one component and one circulation of the other. In the SQUID-type geometry, the value of $q$ can become fractional. These topological complexes can be detected by absorptive imaging. We present both the mean field and Monte Carlo results. The drag effects in optical lattices are drastically different from the Galilean invariant Andreev-Bashkin effect in liquid helium.

Multi-component quantum mixtures in optical lattice (OL) are a source of new and rich many-body physics. The one-component bosons in OL have been exhaustively studied theoretically [1] and experimentally [2], especially in the context of the quantum phase transition between superfluid (SF) and Mott insulator. The study of multi-component boson systems in OL has just began. Theoretical investigations predict variety of new quantum phases with unusual properties [3,4,5,6,7,8]. Two interesting recent examples include topological excitations -vortices and persistent currents with nonstandard winding properties in two-component superfluids (2SF) [5,6].
Crucial, but largely unaddressed effect is the impact of strong interaction on properties of superfluid phases where each component ψ a has its finite expectation value ψ a . Interesting manifestation of the strong interaction is the inter-component drag similar to the Andreev-Bashkin effect [9] in 2SF helium mixtures. In general, the drag effect between non-convertible species at zero temperature is represented by the cross-terms in the expansion of the ground state energy in terms of small gradients of the superfluid phases ∇ϕ a , a = 1, 2, with ρ ab standing for the superfluid stiffnesses. The cross-term ρ 12 is responsible for the drag. It is due to interaction effects and is not confined to some particular term in the full microscopic many-body Hamiltonian. Depending on its sign, this term describes either a mutual unidirectional flow or a counterflow of the components. The Galilean invariance argument, often attributed to Landau, imposes two constraints on ρ ab . These constraints are responsible for the Andreev-Bashkin effect in superfluid mixtures of liquid helium isotopes in which ρ 12 is uniquely related to the ratio of bare m 1 and effective m * 1 masses of minority atoms in the host superfluid of the majority component. The Galilean transformation to a frame moving with velocity V requires that the phase of each component changes as ϕ i → ϕ i −(m i / ) V · r, where m i are the bare masses. The energy density (1) transforms as δE → δE − PV, where P/ = N 1 ∇ϕ 1 + N 2 ∇ϕ 2 is the momentum density expressed in terms of particle densities N 1,2 of each component. This yields In other words, conservation of the total momentum requires that the difference of the bare and effective masses is compensated by the flux of the other component. Note that ρ 12 > 0 since m * 1,2 > m 1,2 . In the case of strong mass renormalization, (m * 1 /m 1 ) − 1 ≥ 1, quite spectacular effects should be expected [10] from the topological excitations -vortices. Specifically, the lowest energy single-circulation vortex of the majority component (ρ 22 ≫ ρ 11 ) should carry several circulation quanta q = 1, 2, .. of minority component. The equilibrium value of q is obtained by minimizing the factor m 2 q 2 + 2 (m * 1 − m 1 ) q in the energy of the vortex complex (or persistent current). These q+1 vortex complexes exhibit transformations with respect to the value of q depending on external conditions that determine the value of m * 1 . If the interaction is weak, ρ 12 can be calculated as an expansion in the gas parameter [11].
In this paper, we address the drag effect in a lattice 2SF in strongly interacting limit, and show that it is radically different from the Galilean-invariant case. The lattice plays a central role in violating the relation [9] between ρ 12 and m 1 /m * 1 (and the constraints (2)). We also argue that the value of q is affected by proximity of the 2SF to the quasi-molecular phase.
In OL, in contrast to the Galilean-invariant system, the lattice provides a preferred reference frame, so that the (hydrodynamic) properties of the two-component mixture are determined not by the relative velocity of components but by their individual velocities with respect to the lattice. Furthermore, the effective mass in OL is formed largely by the width and depth of laser-generated potential wells rather than by a trailing cloud of the second component. Another crucial difference is that in OL number of vacancies is a conserved quantity. Below we perform the mean field and Monte Carlo analysis of the mutual drag in 2SF in three different physical situations: a soft-core system close to molecular condensation, a hard core system with finite intercomponent exchanges, and a hard core system with vacancy-assisted motion without the intercomponent exchanges.
Drag due to proximity to the quasi-molecular state. Here we discuss a generic mechanism leading to the q + 1-topological complexes in the 2SF. Strong drag effect occurs if a two-component boson system is close to a transition into the quasi-molecular state in which the only broken symmetry has the order parameter Φ q ∼ exp(iϕ (q) ) ∼ ψ 1 ψ q 2 = 0 (or ψ 1 ψ †q 2 = 0). In pure molecular state with undefined individual phases ϕ 1,2 (that is, ψ 1,2 = 0), the phase-gradient energy is given by the molecular superfluid phase ϕ (q) as δE = dxρ q (∇ϕ (q) ) 2 /2, with ρ q being molecular superfluid stiffness. The molecular order parameter persists in the 2SF phase so that the additional broken U (1) symmetry emerges continuously [6]. The two phases ϕ 1,2 become well defined in the 2SF state with the molecular phase being locked as This locking can be understood as a consequence of virtual processes of transformation of a (q +1)-molecule into q B-atoms and one A-atom. The corresponding contribution to the energy functional is ∆E ∼ dxΦ q ψ * 1 ψ * q 2 + H.c.. This term (cf. the diatomic molecules with q = 1 [12]) ensures the relation (3) in the longwave limit. Then the energy (1) becomes with ρ ′ ij continuously changing from zero in the molecular phase to some finite values in the 2SF phase. It is important that the molecular stiffness ρ q is not a critical property of the system -it does not change while crossing the phase boundary. Thus, at least close to the phase boundary, minimization of the vortex energy gives ϕ 1 = −qϕ 2 , that is, the q + 1 vortex. In reality, the relations |ρ ′ ab | ≪ |ρ 12 | ≈ ρ 1,2 can hold quite far from the phase boundary. This implies that the q + 1 topological excitation exist deep in the 2SF phase. We demonstrate this numerically for q = 1 (see Fig.1

below).
It is convenient to introduce the drag coefficient k as a ratio k = ρ 12 /ρ 11 of the cross-stiffness to the smallest diagonal stiffness, ρ 11 ≤ ρ 22 . Then, as the minimization of the energy (1) shows, when |k| > 0.5, a vortex of the dominant component can lower its energy if it carries the circulation of the other component q = ±1. In symmetric case (ρ 11 = ρ 22 ), the integer q closest to k determines the q + 1 vortex (or persistent current) as the minimal topological excitation. It is important to note that even small |k| causes attraction between either vortices of equal circulations (k < 0) or between vortex and anti-vortex (k > 0) in different components, so that if both exist they will form a complex. Crossing the boundary |k| > 0.5 has strong impact on mechanisms of vortex creation and stability. For example, stirring the component with the largest stiffness (ρ 22 ) above the threshold will cause creation of the complex instead of a single vortex of the stirred component. Also, a single vortex of the component 2 becomes unstable with respect to inducing creation of vortex of the other component.
The Hubbard lattice model with molecular phases, has been extensively studied analytically [4,8] and numerically [6,7]. Here U α,α ′ is the interactions matrix, t α describes the nearest-neighbor jumps of component α; a † α,i , a α,j are the construction bosonic operators, and n αi = a † α,i a α,i are the on-site occupancies. As discussed in Refs. [6,7], the quasi-molecular phase (U 12 < 0), namely, the paired superfluid, is in many respects isomorphic to the super-counterfluid state (U 12 > 0) [4] . Both states can undergo second order phase transition into the 2SF phase so that the order parameter Φ q=1 = ψ 1 ψ 2 (or Φ q=1 = ψ 1 ψ † 2 ) remains finite and robust. Obviously, in the 2SF phase, the q = 1 composite vortices are the lowest topological excitations. As pointed out in ref. [13], the Hamiltonian (5) also allows molecular phases with arbitrary integer value of q. This issue, though, requires separate analysis.
Hard core limit U ab → ∞ of the Hamiltonian (5). This limit can exhibit quite interesting physics of strong quantum fluctuations even far from any phase transition [5]. Obviously, when N A + N B = 1 (N A , N B are the average on-site occupancies of the species A, B), the system in the hard core limit (HC) is the Mott insulator. Its ground state is degenerate with respect to possible permutations of bosons A and B. This degeneracy, which is a consequence of the HC approximation, is lifted by any infinitely small inter-component exchanges. Accordingly, the two-component HC model should be considered as a limit of the model in which the inter-component interaction V int = U 12 is finite and increasing. In contrast to free space, increasing V int leads to decrease of all superfluid stiffnesses because all transport is suppressed as ∼ t 1 t 2 /V int . Furthermore, in the limit V int → ∞ all stiffnesses are equal in magnitude. This is clearly at variance with the free space constraints (2) which prohibit uniform decrease of all stiffnesses at fixed densities.
The two-component Hamiltonian with residual softcore inter-component repulsion is represented in terms of the HC construction operators a † i , a i and b † i , b i with Pauli commutation relations for the A and B components with summation < ij > over the nearest-neighbor sites. At total filling 1, this Hamiltonian has two phases -2SF, where both SF order parameters are defined, and supercounterfluid (SCF), where the only SF order is observed in a i b † j . Transition between these two phases is con- tinuous in the universality class U(1) [6] and occurs in the symmetric case t 1 = t 2 = t at some value V int = V c , V c /t ∼ 1. As discussed above, the drag effect is strong in the 2SF phase even far from the transition. We proved this by performing the Worm algorithm [14] Monte Carlo simulations of the two-color J-current model [6,7,15] at zero temperature on a 2D square lattice. This model is a discrete-time grand-canonical analog of the Hamiltonian (6) with the hard-core constraints. The stiffnesses were determined from the statistics of the winding numbers similarly to Refs. [6,7,16]. The SCF phase was identified by observing ρ 11 = ρ 22 = −ρ 12 . The negative value of ρ 12 is due to counterflow of the componentseach winding of A-worldline is accompanied by opposite winding of B-worldline. In Fig.1, the drag coefficient k is plotted as a function of the relative interaction strength.
As can be seen, the domain 1/2 < |k| < 1 in the 2SF (between the dashed lines), where the composite 1+q-vortex with q = 1 has lower energy than any single vortex, is not restricted to the vicinity of the critical point V c but occupies about half of the phase diagram. Here ρ 12 < 0, indicating that both components participate in the counterflow even in the 2SF state.
Vacancy assisted drag. If the total filling is different from 1, system is always in 2SF phase at T = 0. In this case, another mechanism contributes to the drag -the vacancy assisted transport. Atoms tunnel to the unoccupied sites (vacancies) much faster than the rate of the A-B exchange with large V int . The vacancies stimulate mass flow in one direction and move in the opposite one. As a result, both components A and B move in one direction, which means that ρ 12 > 0. This situation implies crossover when ρ 12 changes sign at some special point with no drag, ρ 12 = 0. Since no symmetry change takes place, this is not a phase transition. The crossover from k < 0 to k > 0 takes place as V int increases at fixed number of vacancies.
Note that the drag coefficient k must increase when the number of vacancies This counterintuitive result stems from the nature of vacancies. In one component case, conservation of the number of vacancies N V makes them similar to particles. The HC limit links the flow of vacancies with the opposite flow of atoms. In the two-component case, the situation is similar with one crucial difference -a vacancy is not uniquely associated with a particular sort of atoms. Thus, motion of a single vacancy through a lattice in one direction leads to flows of both components in the opposite direction. This implies strong drag with positive k. When x v increases, system becomes more like a low density and, thus, weakly interacting mixture of two sorts of atoms with correspondingly small k.
To analyze the mutual drag and the possibility of complex vortices in the vacancy dominated regime we modified the HC model by imposing the additional constraint a i b i = 0 on (6) and introducing the chemical potentials term −µ 1 N A − µ 2 N B for each component to keep control of the filling factors. As discussed in Ref. [5], this limit can exhibit long range phase separation as well as short scale fluctuative phase separation corresponding to minority particles acquiring large cloud of vacancies.
If 1 − N B > N A , it is convenient to introduce a description in which the vacuum corresponds to all sites filled by B particles. Then, the number n = 1 − N B of B holes is shared between N A atoms and remaining x v = n − N A > 0 vacancies. In the limit N A ≪ N B ∼ 1 transport of vacancies can be considered as transport of B holes with the effective Hamiltonian where v † i , v i are the Pauli operators for B holes. In order to describe the mutual drag within the mean field approximation, one should replace the field operators a, v by the functions a = √ x v exp (−iϕ 2 ) with the slowly varying phases and perform the gradient expansion. [The minus in front of iϕ 2 indicates that flow of holes and actual flow of mass are opposite]. This automatically generates the term ∼ t 1 x 1 x v (∇(ϕ 1 + ϕ 2 )) 2 in effective energy from the first term in Hamiltonian (7). Obviously, the ratio of the stiffnesses becomes k = ρ 12 /ρ 11 = 1 which corresponds to positive crossterm typical for the vacancy assisted transport meaning that the mean field captures well the physics of the vacancy assisted transport. However, the prediction k = 1 and, therefore, q = 1, is not supported numerically.
We have performed Worm algorithm [14] Monte Carlo simulations of the two-color J-current model [6,7,15] at zero temperature in 2D square lattice in the HC limit with partial filling. This model is similar to described above (6) with an additional requirement of no double occupancy. We have found that (Fig.2), for x v 0.15, k < 0.42 ± 0.02, and, thus, no topological complexes can exist as the lowest energy topological excitations in this regime. At this point we do not have a simple explanation for this variance between numerical and mean field results. Most likely, the mean field result is not applicable for large x v in the symmetric mixture when the vacancies cannot be uniquely identified with the holes in the majority component.
As the number of vacancies is tuned to become x v ≤ 0.15, all stiffnesses exhibit large error bars which can be attributed to the regime of strong quantum fluctuations [5] associated with the degeneracy of the ground state in the HC limit. The precise nature of this effect requires separate analysis. For finite V int , depending on N A , N B , the ground state can exhibit various types of ordering including the checkerboard insulator [8]. Then, decreasing x v at N A = N B → 0.5 will result in the first order phase transition with strong fluctuations, similar to those in Fig.2, due to the domain formation.
Fractional q. In the case of finite drag with |k| < 0.5 fractional phase circulation q = k can be observed when persistent current is interrupted by a Josephson junction which lifts the requirement of the integer of 2π windings by creating the phase jump across the junction. Then, phase winding is determined solely by the minimization of energy.
Detection. Finally, the (q + 1) -vortex complexes can be observed by absorptive imaging technique similar to imaging of vortices in one-component Bose-Einstein condensates [17]. Typical pattern should include extra q fringes in one component.
In summary, we explored generic mechanisms of drag effect in quantum bosonic mixtures in optical lattice with hard and soft core interaction. Strong mutual drag can result in composite topological structures. The drag in lattice is not controlled by particle effective masses. The simplest mean field approximation does not adequately describe the strong drag.