Correlated hopping of bosonic atoms induced by optical lattices

In this work we analyze a particular setup with ultracold atoms trapped in state-dependent lattices. We show that any asymmetry in the contact interaction translates into one of two classes of correlated hopping. After deriving the effective lattice Hamiltonian for the atoms, we obtain analytically and numerically the different phases and quantum phase transitions. We find for weak correlated hopping both Mott insulators and charge density waves, while for stronger correlated hopping the system transitions into a pair superfluid. We demonstrate that this phase exists for a wide range of interaction asymmetries and has interesting correlation properties that differentiate it from an ordinary atomic Bose-Einstein condensate.

. (a) Two atoms in state |+ collide and flip state in free space. (b) When these atoms are confined in a state-dependent optical lattice, their change of state must be accompanied by a tunneling event to a neighboring site. This leads to correlated hopping of atoms through the lattice.

Introduction
Ultracold neutral atoms are a wonderful tool to study many-body physics and strong correlation effects. Since the achievement of Bose-Einstein condensation in alkali atoms [1,2,3], we have been witnesses of two major breakthroughs. One is the cooling of fermionic atoms and the study of Cooper pairing and the BCS to BEC transition [4,5,6]. The other one is the implementation of lattice Hamiltonians using bosonic [7,8] and fermionic atoms in off-resonance optical lattices. Contemporary and supported by these experimental achievements, a plethora of theoretical papers has consolidated ultracold atoms as an ideal system for quantum simulations. The goal is two-fold: AMO is now capable of implementing known Hamiltonians which could describe real systems in Condensed Matter Physics, such as Hubbard models [8] and spin Hamiltonians [9,10]; but it is also possible to study new physical effects, such as the quantum Hall effect with bosons [11] and lattice gauge theories [12,13].
In this work we deepen and expand the ideas presented in [14], where we introduced a novel mechanism for pairing based on transport-inducing collisions. As illustrated in Fig. 1, when atoms collide they can mutate their internal state. If the atoms are placed in a state-dependent optical lattice, whenever such a collision happens the pair of atoms must tunnel to a different site associated to their new state. For deep enough lattices, as in the Mott insulator experiments [7], this coordinated jump of pairs of particles can be the dominant process and the ensemble may become a superfluid of pairs.
The focus of this work is to develop these ideas for all possible interaction asymmetries that can happen in experiments with bosonic atoms in two internal states, where the scattering lengths among different states can be different, (g ↑↑ = g ↓↓ ) [15]. We want to understand the different types of correlated hopping that appear when we move beyond the limited type of interactions considered in Ref. [14]. We will better study the resulting phases and phase transitions, using different analytical and numerical tools that describe the many body states, and evolving from few particles to more realistic simulations. The main results will be to show that when one considers more general asymmetries, a new type of correlated hopping appears, but that both the previous [14] and the newly found two-body terms cooperate in creating a superfluid state with pair correlations. Indeed, this new state will also be shown to be more than just a condensate of pairs, based on analytical and numerical studies.
Correlated hopping is not a new idea. It appears naturally in fermionic tight- binding models, where it has been used to describe mixed valence solids [16] and, given that they are able to mimic the attractive interactions between electrons, also high-T c superconductors [17,18,19,20,21,22]. In most of these works, the correlated hopping appears in the form n i a † j a k , indicating that the environment can influence the motion of a particle. This would seem substantially different from correlated motion of pairs of fermions. Nevertheless, even this more elaborate form of correlated hopping has been shown to lead to the formation of bound electron pairs [20,19] and it has been put forward as a possible explanation for high T c superconductivity [23,24].
This work is organized as follows. In Sec. 2 we introduce our model of correlated hopping (1), qualitatively discussing how it originates and what are the quantum phases that we expect it to develop. We present a possible implementation of this model which is based on optical superlattices and atoms with asymmetric interactions. Sec. 3 includes exact diagonalizations for a small number of atoms and sites. These calculations reveal the existence of insulating and coherent regimes, as well as pairing, and will be the basis for later analysis. In Sec. 4 we study the many-body physics of larger lattices with correlated hopping, using a variety of techniques: starting with insulating regime and following with the implementation of perturbation theory and the quantum rotor model. These methods suggest a number of possible phases, including a Mott insulator, a pair superfluid, a normal superfluid and a charge density wave state, and we estimate the parameters for which these phases appear. In Sec. 5 we develop two numerical methods to study our system, a Gutzwiller ansatz and an infinite Matrix Product State method. With these simulations we confirm the above mentioned phases and locate the quantum phase transitions, which are found to be of second order. Finally, in Sec. 6 we suggest some currently available experimental methods to detect and characterize these phases.

Correlated hopping model
We suggested in Ref. [14] that the combination of atomic collisions with optical superlattices can be used to induce correlated hopping. The basic idea is shown in Fig. 2b, where atoms are trapped in two orthogonal states called (+) and (−). The interaction terms change the state of the atoms, forcing them to hop to a different superlattice every time they collide. In this sense, interactions are responsible for transport.
In this section we introduce the most general model of correlated hopping that can be produced by such means with two-states bosons. This model is presented in the following subsection, where we explain qualitatively the role of each Hamiltonian term. Later on, in Sec. 2.2, we establish the connection between the parameters of this model and the underlying atomic model. This is the foundation for the subsequent analytical and numerical studies.

Lattice Hamiltonian
In this work we study the ground state properties of a very general Hamiltonian that contains different kinds of correlated hopping. More precisely, the model will be Here, c † i and c i are bosonic operators that create and anihilate atoms according to the site numbering from Fig. 2b-c, and the colons : A i B j : denote normal ordering of operators A i and B j .
Let us qualitatively explain the roles of the different terms in Eq. (1). The first and second terms, U and V, are related to on-site and next-neighbor interactions. When these terms are dominant, we expect the atoms in the lattice form an insulator. Such a phase is characterized by atoms being completely localized to lattice sites, having well-defined occupation numbers, the absence of macroscopic coherence and a gapped energy spectrum. Whether this insulating state is itself dominated by strong on-site interactions U or by nearest-neighbor repulsion/attraction V will decide whether it presents an uniform density, a Mott insulator (MI), or a periodic density pattern, a charge density wave (CDW), respectively.
The third term is the key feature of our model. It describes the tunneling of pairs between neighboring lattice sites, with amplitude t. Given U, V, j = 0, we expect the atoms to travel along the lattice in pairs forming what we call a pair superfluid (PSF). These pairs will be completely delocalized, establishing long range coherence along the lattice. The observable a 2 would be the figure of merit describing this kind of delocalization, while a vanishing a indicates the absence of the single-particle correlations appearing in a normal superfluid. Furthermore, we expect this phase to have a critical velocity, similar to that of an atomic condensate, and the energy spectrum should be gapless.
Unlike in Ref. [14], when one considers the most general kind of atomic interaction, a second kind of correlated hopping appears, described by the last term in Eq. (1). Here, individual atoms will hop only if there is already a particle in the site they go to (c † i c j (n i − 1)) or leave at least a particle behind ((n i − 1)c j c † i ). One might be induced to think that this term is equivalent to single-particle hopping with a strength that depends on the average density, thus giving rise to a single-particle superfluid (SF) phase. However, this does not seem to be the case. We will show that the correlated hopping j generates a mixed phase which contains features of both the ordinary BEC and the PSF created by t.

Relation to atomic parameters
We now establish the relation between the model in Eq. (1) and the dynamics of atoms in an optical superlattice. The actual setup we have in mind is shown in Fig. (2)ab and described in more detail in Appendix A.1. It consists on a three-dimensional lattice that is strongly confining along the Y and Z directions, creating isolated tubes. On top of this, we create an optical superlattice acting along the X direction [14]. This superlattice traps atoms in the dressed states |+ and |− , while the atomic interaction is diagonal in the basis of bare states |↑ and |↓ . The interaction will be described by a contact potential and parameterized by some real constants g αβ These interaction constants are functions of the s-wave one-dimensional scattering lengths between different species g αβ = 4π 2 a (1D) αβ /m. In general, the interaction strengths among different atomic components are different from each other, a situation that can be enhanced with Feshbach resonances. We will use a parameterization that makes the symmetries more explicit The total Hamiltonian combines the previous interaction with the kinetic energy and the trapping potential for one particle which is written in a different basis Since the superlattice potential V ± (x) is the dominant term, we may approximate the bosonic fields as linear combinations of the Wannier modes in this superlattice and in the dressed state basis, a process detailed in Appendix A. Note that out of all the terms in the interaction Hamiltonian (4), only the first one is insensitive to the state of the atoms. This is important because the asymmetries g 1 and g 2 , when expressed in the dressed basis, produce terms that change the state of the atoms during a collision. Once we introduce the effective interaction constants in the lattice where w(x) is the single site Wannier wavefunction, we arrive to the effective Hamiltonian in Eq. (1) with parameters U, V, t and j that relate to the microscopic model as follows Unlike in the specific case of Ref. [14], the most general situation contains not only two-body correlated hopping t, but also the terms proportional to j.

Preliminary analysis
In this section we study the eigenstates of Hamiltonian (1) for systems that we can diagonalize exactly. The goals are to characterize the effect of the different interaction and hopping terms, as well as to understand the structure of the ground state wavefunction. Although we are limited to a small number of particles, the following examples provide enough evidence of the roles of correlated hopping, nearest neighbor repulsion and the utility of different correlators to characterize the states.

A two-sites example
Let us take the simplest interesting case: four particles in two sites. We write the Hamiltonian in the basis {|40 , |22 , |04 , |31 , |13 }, where the notation |n 1 n 2 stands for n 1 particles in the first site and n 2 in the second and we restrict to n 1 + n 2 = 4, Notice that in this particular case, U gives rise to a global energy shift and does not affect the different eigenstates. This is consistent with later studies where we will see that on-site interactions just add a global, density dependent contribution to the energy. To better understand the role of the remaining terms, we will consider separately three limiting cases, two of superfluid nature and an insulating one.
Limit j = 0, t = V = 0, single-particle delocalization. In this case we take for simplicity V = 0 and diagonalize Eq. 9, finding the normalized ground state Note that this state is exactly a BEC of 4 particles spread over two sites This suggests that, at least in this small example, the correlated hopping proportional to j is equivalent to the single-particle hopping in the ordinary Bose-Hubbard model, giving rise to the delocalization of individual particles. However, as it will become evident later on, for larger systems and more particles this interpretation is wrong.
Limit j = 0, t ≫ |V |, pair delocalization. In the presence of two particle hopping, the lowest energy state has the form with coefficients In particular, for dominant pair hopping t ≫ |V | this is a state of delocalized pairs Observe that this wavefunction is not equivalent to what one would naïvely understand as a "pair condensate" from analogy with the single-particle case Instead the previous wavefunction is isomorphic to the BEC of two bosons under the replacement of each boson with two atoms. It is also interesting to remark that ψ BEC (2) has larger pair-correlations than the state in Eq. (14).
Limit |V | → ∞, an insulator. Reusing the previous wavefunction (12) and taking the limit of dominant nearest-neighbor interaction V, we obtain two possible states. For strong repulsion V → +∞, states |40 and |04 are favored, forming a charge density wave (CDW) with partial filling |ψ CDW ∝ |40 + |04 . On the other hand, for strong nearest-neighbor attractions V → −∞, the particles are evenly distributed forming a Mott insulator |22 .

Superfluidity of pairs
We have seen that four particles in a two-sites lattice recreate the exact wavefunction of an ordinary BEC under the replacement of single bosons with pairs. We can test this idea for slightly bigger lattices, diagonalizing numerically the Hamiltonian which only contains the pair hopping term (t = 0, V, U, j = 0). The resulting wavefunctions are compared side by side with the BEC-like ansatz we mentioned. In the case of two particles we get indeed the expected result |ψ g.s.
we find a disagreement between the ideal case of a BEC-like state with coefficients c 1 = 1/5, c 2 = c 3 = √ 2/5, and the exact diagonalization with c 1 ∼ 0.2735, c 2 ∼ 0.3073 and c 3 ∼ 0.1754. We observe that when compared to the ideal BEC, our paired state breaks the translational symmetry, revealing an effective attraction between different pairs, that favors their clustering.
In Fig. 3 we plot the projection between these states, namely the solution of Eq. 1 with only t = 0 and the ideal superfluid of pairs. In the nearby plot we also analyze two relevant correlators that will be also used later on in the manuscript, namely, a single-particle coherence and the pair correlator As it is evident from the wavefunction and from the plots, there is no single-particle coherence or delocalization because particles move in pairs. Hence, C 1 ∆ ∼ δ ∆0 . The other correlator, C 2 ∆ , which we identify with the delocalization of pairs is rather large and it only decreases with increasing the lattice size because the total pair density becomes smaller.

Analytical methods
We now study the many-body physics of our model for a much larger number of particles using exact analytical methods. We begin with the regime in which the interaction terms U and V dominate, obtaining the different insulator phases on the j, t = 0 plane. Then, using perturbation theory, we compute the phase boundaries of these insulating regions for growing j and t. Finally, we study the properties of the ground state and its excitations in the superfluid phase, with j = 0 and dominating t, proving indeed that this region describes a superfluid of pairs.

No hopping limit: insulating phases
To analyze the phase diagram it is convenient to work in the grand-canonical picture, in which the occupation is determined by the chemical potential µ. In this picture the ground state is determined by minimizing the free energy where N = k n k is the total number of particles, including both states |+ and |− . The free energy has a very simple form in the absence of tunneling This function is defined over positive occupation numbers n k ∈ {0, 1, 2, . . .}. A discrete minimization will determine the different insulating phases and the regions where the system is stable against collapse. For a translational invariant system with periodic boundary conditions, all solutions can be characterized as a function of two integers x t = (n, m), representing the occupations of even n 2k = n and odd sites n 2k+1 = m. The optimization begins by noticing that the bond energy of two sites has a quadratic form where physical solutions are in the sector with n, m ≥ 0. For these occupation numbers to remain bounded, the bond energy ε( x) has to increase as n, m or both grow. This gives us two conditions that need to be fulfilled to prevent collapse. If these conditions are not met, the ground state will be an accumulation of all atoms in the same site. In that case, the large interaction energies and the many-body losses induced by the large densities will cause the breakdown of our model and quite possibly of the experimental setup.
The first stability condition is found by studying ε( x) along the boundaries of our domain (n, m ≥ 0). Take for instance m = 0, this gives a total energy ε B = (U/4)n 2 − [(U + 2µ)/4]n. For this function to have a local minimum at finite n, we must impose The second condition comes from analyzing the interior of the domain.
Given that Eq. (21) and Eq. (22) are satisfied, the system is stable and we have two possibilities to attain the minimum energy: either at the boundaries, n = 0 or m = 0, or right on the eigenvector of A. Inspecting ε B and ε + we conclude that a positive value of V will lead to the formation of charge density waves (CDW) of filled sites alternating with empty sites If V ≤ 0 our energy functional will be convex and the minimum energy state will be a Mott insulator with n = m, when n + m is even, or a charge density wave with n = m ± 1, when n + m is odd. The actual choice between these two insulating phases is obtained by computing the energy of both states Having ǫ(2n + 1) − ǫ(2n) = 0 defines the value of µ at which the state with 2n particles every two sites, a Mott with n particles, stops being the ground state and becomes more favorable to acquire an extra particle to form a CDW. The boundaries of these insulating phases for t, j = 0 are given by Thus summing up, for µ(2n−1 → 2n) ≤ µ ≤ µ(2n → 2n+1) the optimal occupation is n particles per site, forming a Mott, while for µ(2n → 2n+1) ≤ µ ≤ µ(2n+1 → 2n+2) the occupation number is 2n + 1 particles spread over every two sites, having a CDW. The results of this section are summarized in Fig. 4.1.

Perturbation theory: insulator phase boundaries
The previous calculation can be improved using perturbation theory for t, j ≪ U, V around the insulating phases, obtaining the phase boundaries around the insulators as t and j are increased. This is done applying standard perturbation theory up to second order on both variables [25], using as unperturbed Hamiltonian the operator (19) and as perturbation the kinetic energy term We start calculating analytically the ground state energies of the first four insulating phases according to (19), considering the perturbation W up to second order in j, t. For the CDW with n i = 1 and n i+1 = 0 this energy is obviously zero For the MI with one particle per site we have virtual processes of the correlated hopping j, as environment-assisted hopping starts being allowed in an uniformly filled lattice For the CDW with n i = 2 and n i+1 = 1, we find some doubly occupied sites and contributions from the pair hopping t L. (31) In both cases the lowest region is a CDW with alternating 0 and 1 particle occupation, followed upwards by a Mott of one particle per site, a CDW with 1 and 2 particles and the highest area a Mott of two particles.
Finally, for the MI with two particles per site, a calculation detailed in [26], Here L is the total number of sites and all results presented in this section are for the case V < 0. At each value of j, t, the boundary of an insulating phase with average densityn is given by the degeneracy condition with a compressible state E(nL) = E(nL ± 1). Those points correspond to the chemical potential at which a hole, µ h (nL) = E(nL) − E(nL − 1), or a particle, µ p (nL) = E(nL + 1) − E(nL), can be introduced We show here the lower and upper limits of the first four insulating regions, corresponding to the CDW with n i = 1, n i+1 = 0 the Mott with one particle per site the CDW with n i = 2, n i+1 = 1 µ p (L + L/2) = E(L + L/2 + 1) − E(L + L/2) (37) and the MI with two particles per site The corresponding boundaries are plotted in Fig. 5. For small hopping amplitude, they match the values that are found later on with the numerical methods. But even for larger values, this approximation anticipates that the lobes are significantly larger for pair hopping t than for the correlated hopping j.

Phase model: analysis of the pair condensate
So far we have studied the many-body physics around the limit of strong interactions. However, the main goal of this work is to understand the effect of correlated hopping and the creation of a pair superfluid. In absence of a mean field theory, but still in the limit of dominant two-body hopping U, V ≪ t, we can use the number-phase representation, introduced in Ref. [27] for an ordinary BEC. Note, however, that the model in Ref. [27] cannot be directly applied here. Following that reference, one would assume a large number of particles per site, n i > 1, and introduce the basis of phase states | φ n| φ = (2π) −L/2 e i n· φ .
Using these states, one would then develop approximate representations for the operators a 2 i , a †2 i and n i , and diagonalize the resulting Hamiltonian in the limit of weak interactions. But after a few considerations one finds that the resulting phase model does not preserve an important symmetry of our system: if j = 0 particles can only move in pairs and the parity of each site, (−1) ni , is a conserved quantity.
To describe correlated hopping we must use a basis of states with fixed parity ν which is ν = 0 for the ground state we are interested in. As mentioned before, we now have to find expressions for the different operators, a 2 i , a †2 i and n i . We use the fact that our states will have a density close to the average valuen and approximate the action of the operators over an arbitrary state as a †2 i φ = (n + 1)(n + 2)e iφi φ .
Introducing the constant ρ 2 =n(n − 1)(n + 1)(n + 2) our Hamiltonian becomes similar to the quantum rotor model [27] For small U and V, the ground state of this model is concentrated around φ i −φ i+1 = 0. Expanding the Hamiltonian up to second order in the phase fluctuations around this equilibrium point, we obtain a model of coupled harmonic oscillators. This new problem can be diagonalized using normal modes that are characterized by a quasi- with normal frequencies and a global energy It is evident from Eq. (47) that our derivation is only self consistent for negative values of V. Otherwise, when V > 0 some of the frequencies become imaginary, signaling the existence of an unbounded spectrum of modes with |k| ≥ π/4 and that our ansatz becomes a bad approximation of the ground state. This strictly means that our choice φ i = φ i+1 only applies in the case of attractive nearest neighbor interactions, −U ≤ V ≤ 0, as we know that this interaction cannot destabilize a translational invariant solution such as the uniform Mott insulator. However, it does not mean by itself that the whole system becomes unstable for V > 0 -indeed, we will show numerically that it remains essentially in a similar phase for all values of V, but in the case of V > 0 the insulating phases are stable until values of the hopping slightly higher as in the V < 0 case.
If we focus on the regime of validity, we will find that the spectrum is very similar to that of a condensate. At small momenta the dispersion relation becomes linear, ω k ∝ v g k, with sound velocity v g = 4ρ √ 2U t/ , while at larger energies the spectrum becomes quadratic, corresponding to "free" excitations with some mass. This a consequence of the similarity between our approximate model for the pairs (45) and the phase model for a one-dimensional condensate. However, we can go a step further and conclude that the similarity extends also to the wavefunctions themselves, so that the state of a pair superfluid can be obtained from that of an ordinary BEC by the transformation n → 2n. This is indeed consistent with what we obtained for the diagonalization of a two-particle state in the limit j, U, V = 0 [See Eq. (15)].

Numerical methods
The previous sections draw a rather complete picture of the possible ground states in our model. In the limit of strong interactions we find both uniform insulators and a breakdown of translational invariance forming a CDW, while for dominant hopping we expect both single-particle superfluidity and a new phase, a pair superfluid. We now confirm these predictions using two different many-body variational methods.

Gutzwiller phase diagram
The first method that we use is a variational estimate of ground state properties based on a product state [28] |ψ GW = i ni Minimizing the expectation value of the free energy F = H − µN with respect to the variables f n , under the constrain of fixed norm n |f n | 2 = 1, we will obtain the phase diagram in the phase space of interactions and chemical potential (U, V, j, t, µ).
In our study we have made several simplifications. First of all, we assumed period-two translational invariance in the wavefunction, using only two different sets of variational parameters, f (2i+1) n = f 1 n and f (2i) n = f 0 n . In our experience, this is enough to reproduce effects such as the CDW. Next, since U ≥ 0 is required for the stability of the system, we have taken U = 1 as unit of energy. The limit U = 0 is approximated by the limits j, t ≫ 1 in our plots. Finally, in order to determine the roles of j and t, we have studied the cases j = 0 and t = 0 separately. The results are shown in Fig. 6 and Fig. 7 for V < 0 and V > 0, respectively.
The first interesting feature is that, as predicted by perturbation theory, we have large lobes both with integer 1, 2, . . . and with fractional 1/0, 2/1, . . . occupation numbers, forming uniform Mott insulators and CDW, respectively. The insulators are characterized by having a well defined number of particles per site, and thus no number fluctuations ∆n 2 = n 2 − n 2 = 0. While the size of the lobes does not depend dramatically on the sign of V, these are significantly larger for the pair hopping t than for the correlated hopping j, as already seen with perturbation theory.
The boundary of the insulating areas marks a second order phase transition to a superfluid regime, where we find number fluctuations ∆n = 0. In order to characterize these gapless phases we have computed the order parameter of a singleparticle condensate a , and two quantities that we use to detect pairing. The first one is a two-particle correlation that generalizes the order parameter of a BEC to the case of a pair-BEC a 2 . The second quantity ∆a 2 = | a 2 − a 2 | is used to correct the previous value eliminating the contribution that may come from a single-particle condensate coexisting with the pair-BEC.
When j = 0 we always find that a = 0, even outside the insulating lobes. This marks the absence of a single-particle BEC, which is expected since we do not have single-particle hopping. On the other hand, we now find long range coherence of the pairs and thus a 2 = 0 all over the non-insulating area, which we identify with the pair-superfluid regime.
The situation is slightly different for t = 0. The single-particle order parameter a no longer vanishes in the superfluid area, denoting the existence of single-particle coherence, but at the same time we find that the two-particle correlations exceed the contribution from the single-particle superfluid as ∆a 2 = 0, which we attribute to a coexistence of both a single-particle and a pair-superfluid, or a state with both features.
This picture does not change substantially when V is positive or negative. The only differences are in the insulating regions, where the CDW is either due to the incommensurability of the particle number (V < 0) or really gives rise to the separation of particles alternating holes and filled sites (V > 0). However, in the superfluid regime we find no significant changes and in particular we see no breaking of the translational invariance or modulation of the coherent phase.

Matrix Product States: long range pair correlations
The previous numerical simulations are very simple and cannot fully capture the single particle and two-particle correlators. To complete and verify the full picture we have searched the ground states of the full Hamiltonian using the so called iTEBD algorithm, which uses an infinite Matrix Product State ansatz together with imaginary time evolution [29].
Roughly, this ansatz is based on an infinite contraction of tensors that approximates the wavefunction of a translational invariant system in the limit of infinite size. Adapting the ansatz to our problem we write it as Here the Γ o and Γ e are matrices that depend on the state of the odd and even sites they represent, a dependence which is signaled by the n 2k+1 and n 2k+2 in the previous equation. These matrices are contracted with one-dimensional vectors of positive weights λ e,o α ≥ 0, which are related to the coefficients of the Schmidt decomposition. This variational ansatz is known to work well for states with fast decaying correlations, but it also gives a good qualitative description of the critical phases.
In order to optimize the iTEBD wavefunction we performed an approximate imaginary time evolution using a Trotter decomposition and local updates of the associated tensors, as described in Ref. [30]. Using the canonical forms for these tensors it is also straightforward to compute expectation values for different operators acting either on neighboring or separated sites.
In Fig. 8 we plot the most relevant results for three cuts across the phase diagram, µ = 0.5, 1.5 and 2.5, so that each line crosses both an insulating plateau and the superfluid region. We have used small tensor sizes from D = 16 up to 64, a value limited by the need of using large cutoffs for the site populations (n max = 8). As shown in the figures, when j = 0 the single-particle correlator is zero for distinct sites, and we are left only with two-particle correlations. In the MI case the pair correlations between neighboring sites decrease very quickly, while in the superfluid regime we see a critical behavior with an exponent that varies between α = 0.5 and α = 0.6, depending on the simulation parameters.

Detection
There are many ways of differentiating the phases we have found, each one having its own degree of difficulty. The simplest and best established detection methods are related to the insulating regimes. These phases, which involve both MI and CDW, are characterized by having a well defined number of particles at each lattice site, the lack of coherence and an energy gap that separates the insulator from other excitations. The energy gap in these insulators may be probed either by static or spectroscopic means as it has been done in experiments [7,31], determining that indeed the system is insulating. Second, the lack of coherence will translate into featureless time-of-flight images, having no interference fringes [7] at all. Even though there will be no fringes, the measured density will be affected by quantum noise. The analysis of the noise correlation will show peaks at certain momenta [32,33] that depend on the periodicity of the state, so that the number of peaks in the CDW phase will be twice those of the MI. In case of having access to the lattice sites, as in the experiments with electron microscopy [34], or in future experiments with large aperture microscopic objectives that can collect the fluorescence of individual lattice sites [35,36], the discrimination between the MI and the CDW should be even easier, since in one case we have a uniform density and in the other a periodic distribution of atoms.
When the system enters a superfluid phase, it becomes a perfect "conductor" with a gapless excitation spectrum. The lack of an energy gap, should be evident in the spectroscopic experiments suggested before. However, we are not only interested in the superfluid nature, but rather in the fact that this quantum phase is strongly paired. More precisely, we have found that for j = 0 the single-particle coherence is small or zero, and that the two-particle correlator decays slowly The first equation implies that the time of flight images will reveal no interference fringes and will exhibit noise correlations which will be similar to the MI. In order to probe C 2 ∆ and confirm the pairing of the particles, we suggest to use Raman photoassociation to build molecules out of pairs of atoms [37,38]. For an efficient conversion, it would be best to perform an adiabatic passage from the free atoms to the bound regime. As described in [39], we expect a mapping that goes from |2n → |n , where 2n is the number of bosons and n the number of molecules. More precisely, we expect the a †2 operator to be mapped into m † , so that the pair coherence of the original atoms translates into the equivalent of C 1 ∆ for the molecules. This order should reveal as an interference pattern in time-of-flight images of the molecules.
Finally, in cases with j = 0, we have found the coexistence of single-particle and two-particle coherences. This translates also into the coexistence of interference fringes with nonzero pair correlators.

Conclusions
Summing up, we have suggested a family of experiments with cold atoms that would produce correlated hopping of bosons. The mechanism for the correlated hopping is an asymmetry in the contact interaction between atoms. This asymmetry is exploited by trapping the atoms in dressed states, a configuration that gives rise to transport induced by collisions.
The main result of this paper is that there is a huge variety of interaction asymmetries that will give rise to long range pair correlations via interaction-induced transport. Formally, in the resulting effective models we recognize two dynamical behaviors. If we have a nonzero asymmetry in the interspecies interactions, the Hamiltonian will exhibit pair hopping, while an asymmetry in the intra-species scattering lengths gives rise to correlated hopping. However, we have given enough evidence that both Hamiltonian terms give rise to a novel quantum phase which we call pair superfluid. This phase is characterized by a gapless spectrum with a finite sound speed, zero single-particle correlations and long range pair coherence. All quantum phases are connected by second order quantum phase transitions. These phases can be produced and identified using variations of current experiments [40,41,33]. The nonperturbative nature of the effect should help in that respect.
Our ideas are not restricted to one dimension. It is possible to engineer also a twobody hopping using two-dimensional lattice potentials. Again, the basic ingredients would be atoms with an asymmetric interaction and an optical lattice that traps two states, |+ and |− with a relative displacement. Both in the one-and two-dimensional cases it is a valid question to ask whether the coupling between different trapped states, |± , can excite also transitions to higher bands, processes that have not been considered in the paper. Our answer here is no. There are only two sources of coupling to higher energy bands. One is the interaction, but we are already assuming that the interaction energies are much smaller than the band separation. Following the notation from Ref. [8], we have the constraint that the interaction energy should be smaller than the energy separation to the first excited state in a well of the periodic potential, n 2 U ≪ νn the same requisite as for ordinary Bose-Hubbard models [7]. The other source of coupling to higher bands would be single-particle hopping. However, unlike [8], here we are assuming that these terms are strongly suppressed compared with the interaction. In other words, realizing the models that we suggest in this paper, for realistic densities,n = 2, and simple potentials, imposes no further constraint in current experiments.
Finally, let us remark that transport-inducing collisons may be implemented using other kinds of spin-dependent interactions. For instance, correlated hopping appears naturally in state-dependent lattices loaded with spinor atoms, because their interactions can change the hyperfine state of the atoms while preserving total angular momentum [42].
We would like to thank Miguel Angel Martín-Delgado for useful discussions. M.E. acknowledges support from the CONQUEST project. J.J.G.R acknowledges financial support from the Ramon y Cajal Program of the Spanish M.E.C., from U.S. NSF

Appendix A.2. State-changing collisions
We will now express the interaction (4) in the basis of dressed states. We proceed using the change of variables in Eq. 6 to find the expression of the densities The first obvious conclusion is that the total density is independent of the basis on which it is written, Hence, the term of g 0 is insensitive to the state of the atoms. On the other hand, the asymmetric terms are not so simple. The g 1 interaction, which is proportional to the product of densities : ρ ↑ ρ ↓ : = 1 4 : (ρ + + ρ − ) 2 : − 1 4 : (ψ † + ψ − + ψ † − ψ + ) 2 : gives rise to a scattering that changes the state of interacting atoms from |− to |+ and viceversa, as in Fig. 1a. The term of g 2 has a lightly different effect, it gives rise to processes where one atom changes its state influenced by the surrounding environment. In the following subsections we will see what happens to the interaction terms (A.8), (A.9) and (A.10), when the atoms are confined in a lattice.

Appendix A.3. Final model
In this section we will put the previous results of this appendix together. We will take the tight-binding expansion of the field operators (A.3) and use it together with Eqs. (A.8), (A.9) and (A.10) to expand the interaction Hamiltonian (4). For convenience, we will rename the bosonic operators as c 2k = a k+ and c 2k+1 = a k− (A.11) according to the position at which their Wannier functions are centered (see Fig. 2c). Along the derivation, one obtains many integrals of ground state wavefunctions We will only keep those integrals with a separation smaller than a superlattice period. Taking Eq. (A.5), the expression for the superlattice localized states, one obtains where w(x) are the Wannier wavefunctions of the underlying sublattice. Using these tools, the symmetric interaction term becomes g 0 2 dx : (ρ ↑ (x) + ρ ↓ (x)) 2 := = g 0 2 N/2 k : n 2 2k C 2k,2k + n 2 2k+1 C 2k+1,2k+1 + 2n 2k n 2k+1 C 2k,2k+1 : and then finally the more complicated Eq. (A.10) Introducing constants that parameterize the on-site interactions and the strength of the underlying lattice (7) Completing terms and replacing the sum over k with a sum over nearest neighbors, we arrive at the desired model (1) with the parametrization given already in Eq. (8).