Transport of strong-coupling polarons in optical lattices

We study the transport of ultracold impurity atoms immersed in a Bose-Einstein condensate (BEC) and trapped in a tight optical lattice. Within the strong-coupling regime, we derive an extended Hubbard model describing the dynamics of the impurities in terms of polarons, i.e. impurities dressed by a coherent state of Bogoliubov phonons. Using a generalized master equation based on this microscopic model we show that inelastic and dissipative phonon scattering results in (i) a crossover from coherent to incoherent transport of impurities with increasing BEC temperature and (ii) the emergence of a net atomic current across a tilted optical lattice. The dependence of the atomic current on the lattice tilt changes from ohmic conductance to negative differential conductance within an experimentally accessible parameter regime. This transition is accurately described by an Esaki-Tsu-type relation with the effective relaxation time of the impurities as a temperature-dependent parameter.


Introduction
The study of impurities immersed in liquid helium in the 1960s opened a new chapter in the understanding of the structure and dynamics of a Bose-condensed fluid [1]. Two of the most prominent examples include dilute solutions of 3 He in superfluid 4 He [1,2] and the measurement of ionic mobilities in superfluid 4 He [1,3]. More recently, the experimental realization of impurities in a Bose-Einstein condensate (BEC) [4,5] and the possibility to produce quantum degenerate atomic mixtures [6][7][8][9][10][11] have generated renewed interest in the physics of impurities. Of particular importance, collisionally induced transport of impurities in an ultra-cold bosonic bath has been recently observed [12] in an experimental setup of a similar type as to the one considered in this paper.
In the context of liquid helium and Bose-Einstein condensates, a plethora of theoretical results have been obtained. Notably, the effective interaction between impurities [13][14][15][16] and the effective mass of impurities [17][18][19][20][21], both of which are induced by the condensate background, have been studied in detail. Moreover, the problem of self-trapping of static impurities has been addressed [22][23][24] in the framework of Gross-Pitaevskii (GP) theory [25,26], and quantized excitations around the ground state of the self-trapped impurity and the distorted BEC have been investigated in [27]. In the present paper we study the dynamics of impurities immersed in a nearly uniform BEC, which in addition are trapped in an optical lattice potential [28,29]. This setup has been considered in [30] by the present authors, however, we here motivate and develop the small polaron formalism in more detail and extend the results to the case of a tilted optical lattice potential, as illustrated in Fig. 1. In contrast to similar previously-studied systems, the presence of an optical lattice allows us to completely control the kinetic energy of the impurities. Crucially, for impurities in the lowest Bloch band of the optical lattice the kinetic energy is considerably reduced by increasing the depth of the lattice potential [28,29]. As a consequence, it is possible to access the so-called strong-coupling regime [31], where the interaction energy due to the coupling between the impurity and the BEC is much larger than the kinetic energy of the impurity. This regime is particularly interesting since it involves incoherent and dissipative multiphonon processes which have a profound effect on the transport properties of the impurities.
As a starting point we investigate the problem of static impurities based on the GP approximation and by means of a Bogoliubov description of the quantized excitations [26,32]. The main part of the paper is based on small polaron theory [33][34][35], where the polaron is composed of an impurity dressed by a coherent state of Bogoliubov phonons. This formalism accounts in a natural way for the effects of the BEC background on the impurities and allows us to describe the BEC deformation around the impurity in terms of displacement operators. As a principal result, we obtain an extended Hubbard model [36,37] for the impurities, which includes the hopping of polarons and the effective impurity-impurity interaction mediated by the BEC. Moreover, the microscopic model is particularly suitable for the investigation of the transport of impurities from first principles within the framework of a generalized master equation (GME) approach [38][39][40][41].
The transport properties of the impurities are considered in two distinct setups. We first show that in a non-tilted optical lattice the Bogoliubov phonons induce a crossover from coherent to incoherent hopping as the BEC temperature increases [30]. In addition, we extend the results in [30] to the case of a tilted optical lattice (see Fig. 1). We demonstrate that the inelastic phonon scattering responsible for the incoherent hopping of the impurities also provides the necessary relaxation process required for the emergence of a net atomic current across the lattice [42]. In particular, we find that the dependence of the current on the lattice tilt changes from ohmic conductance to negative differential conductance (NDC) for sufficiently low BEC temperatures. So far, to our knowledge, NDC has only been observed in a non-degenerate mixture of ultra-cold 40 K and 87 Rb atoms [12] and in semiconductor superlattices [43,44], where the voltagecurrent dependence is known to obey the Esaki-Tsu relation [45]. By exploiting the analogy with a solid state system we show that the transition from ohmic conductance to NDC is accurately described by a similar relation, which allows us to determine the effective relaxation time of the impurities as a function of the BEC temperature.
The implementation of our setup relies on recent experimental progress in the production of quantum degenerate atomic mixtures [6][7][8][9] and their confinement in optical lattice potentials [10,11]. As proposed in [46], the impurities in the optical lattice can be cooled to extremely low temperatures by exploiting the collisional interaction with the surrounding BEC. The tilt imposed on the lattice potential may be achieved either by accelerating the lattice [47] or superimposing an additional harmonic potential [12,48].
An essential requirement for our model is that neither interactions with impurities nor the trapping potential confining the impurities impairs the ability of the BEC to sustain phonon-like excitations. The first condition limits the number of impurity atoms [10,11,49], and thus we assume that the filling factor of the lattice is much lower than one, whereas the second requirement can be met by using a species-specific optical lattice potential [50]. Unlike in the case of self-trapped impurities [22][23][24], we assume that the one-particle states of the impurities are not modified by the BEC, which can be achieved by sufficiently tight impurity trapping as shown in Appendix C.
The paper is organized as follows. In Section 2 we present our model. In Section 3, we motivate the small polaron formalism by considering the problem of static impurities based on the GP equation and within the framework of Bogoliubov theory. In Section 4, the full dynamics of the impurities is included. We develop the small polaron formalism in detail and derive an effective Hamiltonian for hopping impurities. In Section 5, the GME is used to investigate the transport properties of the impurity atoms. First, we discuss the crossover from coherent to incoherent hopping and subsequently study the dependence of a net atomic current across a tilted optical lattice on the system parameters. We discuss possible extensions of our model and conclude in Section 6. Throughout the paper we consider impurity hopping only along one single direction for notational convenience. However, the generalization to hopping along more than one direction is straightforward.

Model
The Hamiltonian of the system is composed of three partsĤ =Ĥ B +Ĥ χ +Ĥ I , wherê H B is the Hamiltonian of the Bose-gas,Ĥ χ governs the dynamics of the impurity atoms andĤ I describes the interactions between the impurities and the condensate atoms.
The Bose-gas is assumed to be at a temperature T well below the critical temperature T c , so that most of the atoms are in the Bose-condensed state. The boson-boson interaction is represented by a pseudo-potential gδ(r − r ′ ), where the coupling constant g > 0 depends on the s-wave scattering length of the condensate atoms. Accordingly, the grand canonical HamiltonianĤ B of the Bose-gas readŝ where µ b is the chemical potential and is the single-particle Hamiltonian of the non-interacting gas, with m b the mass of a boson and V ext (r) is a shallow external trapping potential. The boson field operatorsψ † (r) and ψ(r) create and annihilate an atom at position r, respectively, and satisfy the bosonic commutation relations [ψ(r),ψ † (r ′ )] = δ(r − r ′ ) and [ψ(r),ψ(r ′ )] = [ψ † (r),ψ † (r ′ )] = 0.
To be specific we consider bosonic impurity atoms with mass m a loaded into an optical lattice, which are described by the impurity field operatorχ(r). The optical lattice potential for the impurities is of the form V L (r) = V x sin 2 (kx) + V y sin 2 (ky) + V z sin 2 (kz), with the wave-vector k = 2π/λ and λ the wavelength of the laser beams. The depth of the lattice V ℓ (ℓ = x, y, z) is determined by the intensity of the corresponding laser beam and conveniently measured in units of the recoil energy E R =h 2 k 2 /2m a . Since we are interested in the strong-coupling regime we assume that the depth of the lattice exceeds several recoil energies. As a consequence, the impurity dynamics is accurately described in the tight-binding approximation by the Bose-Hubbard model [28,29], where the mode-functions of the impurities χ j (r) are Wannier functions of the lowest Bloch band localized at site j, satisfying the normalization condition dr|χ j (r)| 2 = 1. Accordingly, the impurity field operator is expanded aŝ χ(r) = j χ j (r)â j and the Bose-Hubbard Hamiltonian is given bŷ where µ a is the energy offset, U is the on-site interaction strength, J is the hopping matrix element between adjacent sites and i, j denotes the sum over nearest neighbours. The operatorsâ † j (â j ) create (annihilate) an impurity at lattice site j and satisfy the bosonic commutation relations jâj is the number operator. The last term in Eq. (3) was added to allow for a tilted lattice with the energy levels of the Wannier states χ j (r) separated by Bloch frequency ω B , as shown in Fig. (1).
For the lattice depths considered in this paper it is possible to expand the potential wells about their minima and approximate the Wannier functions by harmonic oscillator ground states. The oscillation frequencies of the wells are ω ℓ = 2 √ E R V ℓ /h and the spread of the mode-functions χ j (r) is given by the harmonic oscillator length σ ℓ = h/m a ω ℓ . Explicitly, the Gaussian mode-function in one dimension reads where x j = aj with a = λ/2 being the lattice spacing. We note that the mode-functions χ j,σ (r) are highly localized on a scale short compared to the lattice spacing a since σ ℓ /a = (V ℓ /E R ) −1/4 /π ≪ 1 for sufficiently deep lattices. Consistency of the tight-binging model requires that excitations into the first excited Bloch band due to Landau-Zener tunneling or on-site interactions are negligible, i.e. ω B ≪ ω ℓ and 1/2 Un j (n j − 1) ≪hω ℓ , with n j the occupation numbers. These inequalities are readily satisfied in practice. The density-density interaction between the impurities and the Bose-gas atoms is again described by a pseudo-potential κδ(r − r ′ ), with κ the coupling constant, and the interaction Hamiltonian is of the form We note that the effect ofĤ B andĤ I dominate the dynamics of the system in the strongcoupling regime and hence the kinetic part inĤ χ can be treated as a perturbation.

Static impurities
In a first step we investigate the interaction of static impurities with the BEC, i.e. we neglect the hopping term inĤ χ . We first derive a mean-field Hamiltonian based on the GP approximation and subsequently quantize the small-amplitude oscillations about the GP ground state using standard Bogoliubov theory [26,32].
In the GP approximation, the bosonic field operatorψ(r) is replaced by the classical order parameter of the condensate ψ(r). Accordingly, the impurities are described by the impurity density ρ χ (r) = χ † (r)χ(r) = j n j |χ j (r)| 2 , where the average is taken over a fixed product state |Υ representing a set of static impurities, i.e. |Υ = {j}â † j |0 , with {j} a set of occupied sites and |0 the impurity vacuum. We linearize the equation for ψ(r) by considering small deviations ϑ(r) = ψ(r) − ψ 0 (r) from the ground state of the condensate ψ 0 (r) in absence of impurities. This approach is valid provided that |ϑ(r)|/ψ 0 (r) ≪ 1 and corresponds to an expansion ofĤ B +Ĥ I to second order in κ. Explicitly, by replacingψ(r) with ψ 0 (r)+ϑ(r) inĤ B +Ĥ I we obtain the GP Hamiltonian In absence of impurities, i.e. for κ = 0 and ϑ(r) identically zero, the stationary condition δH GP /δψ * 0 (r) = 0 implies that ψ 0 (r) satisfies the time-independent GP equation The deformed ground state of H GP due to the presence of impurities is found by imposing the stationary condition δH GP /δϑ * (r) = 0 or equivalently which is the Bogoliubov-de Gennes equation [26,32] for the zero-frequency mode.
For the case of a homogeneous BEC, the effect of the impurities on the ground state of the condensate can be expressed in terms of Green's functions [22,23,27]. If we assume for simplicity that both ψ 0 (r) and ϑ(r) are real then Eq. (8) reduces to the modified Helmholtz equation Figure 2. The order parameter ψ(r) of a two-dimensional homogeneous condensate deformed by two highly localized impurities with κ > 0 according to Eq. (10). The deformation of the BEC results in an effective interaction potential V i,j between the impurities at lattice sites i and j. The range of the potential is characterized by the healing length ξ, which is comparable to the lattice spacing a for realistic experimental parameters.
where ξ =h/ √ m b gn 0 is the healing length, D is the number of spatial dimensions, and where the Green's functions are with K 0 (x) the modified Bessel function of the second kind. Thus, independent of the dimension D of the system, the BEC deformation induced by the impurities falls off exponentially on a length scale set by ξ, as illustrated in Fig. (2). It follows from Eq. (10) that for occupation numbers n j ∼ 1 the condition with d = n −1/D 0 the average separation of the condensate atoms. The generalization of this condition to the non-homogeneous case is α(r) = (|κ|/g) [d(r)/ξ(r)] D ≪ 1 under the assumption that the BEC is nearly uniform. We note that α(r) ∝ n 0 (r) (D−2)/2 , and hence α(r) diverges in 1D as n 0 (r) → 0, e.g. near the boundary of the condensate. However, this is consistent with the GP approach, which is no longer applicable in the dilute limit ξ(r)/d(r) ≪ 1 of a 1D Bose-gas [26].
The change in energy of the BEC due to the impurities is found by inserting the formal solution for ϑ(r) in Eq. (10) into H GP . Provided that ϑ(r) satisfies Eq. (8) the identity H lin + 2H ϑ = 0, or equivalently H GP = H ψ 0 + 1 2 H lin , holds. Using the latter expression for H GP we find whereĒ = κn 0 is the first order contribution to the mean-field shift. In Eq. (13) we have neglected constant terms which do not depend on the impurity configuration, i.e. contributions containing ψ 0 (r) only. The off-site interaction potential between the impurities is given by and E j = 1 2 V j,j is the potential energy of an impurity, both resulting from the deformation of the condensate, as illustrated in Fig. (2).
Consequently, the total HamiltonianĤ stat of the BEC and the static impurities is composed of three partŝ whereĤ ζ governs the dynamics of the Bogoliubov quasi-particles and H GP is the GP ground state energy, which for the case of a homogeneous BEC takes the simple form in Eq. (13). The third term Ĥ χ represents the average value ofĤ χ with respect to the fixed product state |Υ introduced earlier with n j impurities at each site j, giving explicitly The ground state of the system corresponds to the Bogoliubov vacuum defined bŷ β ν |vac = 0.

Hopping impurities in the polaron picture
Given the HamiltonianĤ stat derived in the previous section it is straightforward to determine the ground state energy for a set of static impurities. However, an analysis of the full dynamics of the impurities requires an alternative description in terms of polarons, i.e. impurities dressed by a coherent state of Bogoliubov quasi-particles. In this picture small polaron theory allows us, for example, to calculate the effective hopping matrix element for the impurities, which takes the BEC background into account. Our approach is based on the observation that the HamiltonianĤ ζ of the quantized excitations and consequently the Bogoliubov-de Gennes equations (17) are independent of ρ χ (r). In other words, the effect of the impurities is only to shift the equilibrium position of the modesβ ν without changing the spectrumhω ν . Therefore it is possible to first expand the bosonic field operator asψ(r) = ψ 0 (r) +θ(r), withθ(r) = ϑ(r) +ζ(r), subsequently expressθ(r) in terms of Bogoliubov modes about the state ψ 0 (r) and finally shift their equilibrium positions in order to (approximately) minimize the total energy of the system.
Similarly to the case of static impurities, we first replaceψ(r) with ψ 0 (r) +θ(r) in H B +Ĥ I to obtain the HamiltonianĤ GP =Ĥ ψ 0 +Ĥ ϑ +Ĥ lin witĥ We note that nowĤ GP contains the density operatorχ † (r)χ(r) instead of ρ χ (r) in order to take into account the full dynamics of the impurities. The expansion ofθ(r) in terms of Bogoliubov modes readŝ where the spectrumhω ν and coefficients u ν (r) and v ν (r) are determined by Eqs. (17). The bosonic operatorsb † ν (b ν ) create (annihilate) a Bogoliubov excitation around the ground state of the condensate ψ 0 (r) in absence of impurities, and thus, importantly, do not annihilate the vacuum |vac defined in the previous section, i.e.b ν |vac = 0. By substituting the expansion in Eq. (21) forθ(r) in the HamiltonianĤ GP and using the identityχ † (r)χ(r) = i,j χ * i (r)χ j (r)â † iâj we find that up to constant terms § with the matrix elements where we assumed for simplicity that ψ 0 (r) is real. The non-local couplings M i,j,ν andĒ i,j with i = j resulting from the off-diagonal elements inχ † (r)χ(r) are highly suppressed because the product of two mode-functions χ i (r) and χ j (r) with i = j is exponentially small. As a consequence, the evolution of the BEC and the impurities, including the full impurity HamiltonianĤ χ , is accurately described by the Hubbard-Holstein Hamiltonian [33][34][35] H hol =Ĥ χ + j,νh where we discarded the non-local couplings and introducedĒ j =Ē j,j and M j,ν = M j,j,ν .
Since we intend to treat the dynamics of the impurities as a perturbation we shift the operatorsb ν andb † ν in such a way thatĤ hol in Eq. (24) is diagonal, i.e. the total energy of the system is exactly minimized, for the case J = 0. This is achieved by the unitary Lang-Firsov transformation [34,35] from which it follows that by applying the Baker-Campbell-Hausdorff formula [52] Ub † § As forĤ ζ , the details of establishing the above form ofĤ ϑ are contained in Ref. [32].
The operatorX † j is the sought after displacement operator that creates a coherent state of Bogoliubov quasi-particles, i.e. a condensate deformation around the impurity. As a result, the transformed HamiltonianĤ LF =ÛĤ holÛ † is given bŷ with the effective energy offsetμ j = µ a +Ē j − E j , the effective on-site interaction , and E j = 1 2 V j,j the so-called polaronic level shift [34,35]. In the case of static impurities and, more generally, in the limit ζ = J/E j ≪ 1, the polarons created byâ † jX † j are the appropriate quasi-particles. Consequently, the Hamiltonian in Eq. (27) describes the dynamics of hopping polarons according to an extended Hubbard model [36,37] with an non-retarded interaction potential V i,j . The contributions toĤ LF are qualitatively the same as for static impurities, except for the additional hopping term. In the thermodynamic limit, the interaction potential V i,j and the polaronic level shift E j are identical, respectively, to their GP counterparts V i,j and E j . The reason for the simple connection between GP theory and small polaron results is that both involve a linearization of the equations describing the condensate. It should be noted that the interaction potential V i,j and the polaronic level shift E j can also be obtained exactly fromĤ hol by applying standard Rayleigh-Schrödinger perturbation theory up to second order in κ since all higher order terms vanish. However, the merit of using the Lang-Firsov transformation lies in the fact that it yields a true many-body description of the state of the system, which would require a summation of perturbation terms to all orders in κ [53].
We gain qualitative and quantitative insight into the dependence of the quantities in Eq. (27) on the system parameters by considering the specific case of a homogeneous BEC with the total Hamiltonian with the energy offsetμ = µ a + κn 0 − E p , the on-site interaction strengthŨ = U − 2E p , the polaronic level shift E j ≡ E p and the phonon momentum q. For a homogeneous BEC the Bogoliubov coefficients are of the form u q (r) = u q exp(iq·r) and v q (r) = v q exp(iq·r) with [26,32] u q = 1 √ 2Ω ε q + gn 0 hω q + 1 Here, ε q = (hq) 2 /2m b is the free particle energy,hω q = ε q (ε q + 2gn 0 ) is the Bogoliubov dispersion relation and Ω is the quantization volume. The corresponding matrix elements are found to be where For the Gaussian mode-function χ j,σ (r) in Eq. (4) the factor f j (q) takes the simple form The q-dependence of the impurity-phonon coupling is M j,q ∝ f j (q)/ |q| in the long wavelength limit |q| ≪ 1/ξ, which corresponds to a coupling to acoustic phonons in a solid state system [35], whereas in the free particle regime |q| ≫ 1/ξ we have M j,q ∝ f j (q)/q 2 . The potential V i,j and the polaronic level shift E p reduce to a sum over all momenta q and can be evaluated in the thermodynamic limit Ω −1 q → (2π) −D dq where V i,j → V i,j . In particular, for the Gaussian mode-functions χ j,σ (r) one finds where the functions G(r, σ) are defined in Appendix A and plotted in Fig. 3a for a three-dimensional system. It follows from the definition of G(r, 0) that G(r) ≡ G(r, 0), and hence the shape of the potential is determined by the Green's function G(r) in the limit σ ≪ ξ. The σ-dependence of G(0, σ), which determines the depth of the interaction potential V i,j and the level shift E p = 1 2 V j,j , is shown in Fig. 3b. The range of the potential, characterized by the healing length ξ, is comparable to the lattice spacing a for realistic experimental parameters, and hence the off-site terms V j,j+1 are non-negligible. In particular, as shown in [30,54], the off-site interactions can lead to the aggregation of polarons on adjacent lattice sites into stable clusters, which are not prone to loss from three-body inelastic collisions.

Transport
The coupling to the Bogoliubov phonons via the operatorsX † j andX j changes the transport properties of the impurities notably. Since we assume the filling factor of the lattice to be much lower than one we can investigate the transport properties by considering a single polaron with the Hamiltonianĥ 0 +ĥ I given bŷ To start we investigate the crossover from coherent to diffusive hopping [33,35] in a non-tilted lattice (ω B = 0) and then extend the result to a tilted lattice (ω B = 0) to demonstrate the emergence of a net atomic current across the lattice [42] due to energy dissipation into the BEC.

Coherent versus incoherent transport
We first consider coherent hopping of polarons at small BEC temperatures k B T ≪ E p , where incoherent phonon scattering is highly suppressed. In the strong-coupling regime ζ = J/E j ≪ 1, and hence the hopping term in Eq. (35) can be treated as a perturbation. The degeneracy of the Wannier states requires a change into the Bloch basis with k the quasi-momentum and N the number of lattice sites. Applying standard perturbation theory in the Bloch basis and describing the state of the system by |k, {N q } , where {N q } is the phonon configuration with phonon occupation numbers N q , we find the polaron energy up to first order in ζ where we defined the effective hoppingJ = J q N q |X † jXj+1 |N q and a is the position vector connecting two nearest neighbor sites. In particular, for the case of a thermal phonon distribution with occupation numbers N q = (eh ωq/k B T − 1) −1 the effective hopping is [33][34][35] Thus, the hopping bandwidth of the polaron band is highly suppressed with increasing coupling κ and temperature T . At high temperatures E p ≪ k B T ≪ k B T c inelastic scattering, in which phonons are emitted and absorbed, becomes dominant, and thus the transport of impurities through the lattice changes from being purely coherent to incoherent. While matrix elements {N q }|X † jXj+1 |{N q } ′ involving two different phonon configurations {N q } and {N q } ′ vanish at zero temperature they can take non-zero values for T > 0. The condition for energy and momentum conservation during a hopping event implies that incoherent hopping is dominated by a three-phonon process that involves phonons with a linear dispersionhω q =hc|q|, where c ∼ gn 0 /m b is the speed of sound. This process is reminiscent of the well-known Beliaev decay of phonons [26].
We investigate the incoherent transport properties by using a generalized master equation (GME) [38] for the site occupation probabilities P j (t) of the impurity. This formalism has been applied to the transfer of excitons in the presence of electron-phonon coupling [39,40] and is based on the Nakajima-Zwanzig projection method [41]. The generalized master equation is of the form where the effect of the condensate is encoded in the memory functions W i,j (s). As shown in Appendix B the memory function to second order in ζ is given by In the case ω B = 0, the nontrivial part of W i,j (s) takes the values 2(J/h) 2 at s = 0 and 2(J/h) 2 in the limit s → ∞ due to the cancellation of highly oscillating terms, as shown in Fig. 4a.
In the regime k B T ≪ E p , the effective hoppingJ is comparable to J and the memory function W i,j is well approximated by 2(J/h) 2 Θ(s), with Θ(s) the Heaviside step function, and thus describes purely coherent hopping. In the regime E p ≪ k B T ≪ k B T c , processes involving thermal phonons become dominant and coherent hopping is highly suppressed, i.e.J ≪ J. In this case the memory function W i,j drops off sufficiently fast for the Markov approximation to be valid, as illustrated in Fig. 4a. More precisely, one can replace P j (t − s) by P j (t) in Eq. (39) and after intergration over s the GME reduces to the standard Pauli master equation where the hopping rates w i,j are given by The Pauli master equation describes purely incoherent hopping with a thermally activated hopping rate w i,j [33,35]. The evolution of an initially localized impurity at different temperatures T for a one-dimensional 41 K -87 Rb system [8] is shown in Figs. 4b -4d, which were obtained by numerically solving the GME with the memory function W i,j (s) in Eq. (40). It can be seen that for small temperatures k B T ≪ E p the hopping is coherent with two wavepackets moving away from the initial position of the impurity atom. In contrast, for high temperatures k B T ≫ E p inelastic scattering of phonons results in a diffusive motion of the impurity, where the probability P j (t) remains peaked at the initial position of the impurity. Figure 5. Analysis of the mean-square displacement l 2 (t) = l l 2 P l (t), obtained from the evolution of an initially localized impurity in a one-dimensional system for the time t evol = 10h/J according to the GME. The mean-squared displacement was assumed to be of the form l 2 (t) = A t α . (a) The exponent α versus temperature T for the full evolution time (solid line) and the period t evol /2 to t evol (dashed line). The drop from α ≈ 2 at zero temperature to α ≈ 1 at high temperatures k B T ≫ E p clearly indicates the crossover from coherent to diffusive transport. (b) The prefactor A in units of (J/h) α versus temperature T for the full evolution time (solid line) and the period t evol /2 to t evol (dashed line). The system parameters are the same as for Fig. 4.
This crossover from coherent to diffusive hopping can be quantitatively analyzed by considering the mean-squared displacement of the impurity, l 2 (t) = l l 2 P l (t), which we assume to be of the form l 2 (t) = A t α . The exponent α takes the value α = 2 for a purely coherent process, whereas α = 1 for diffusive hopping. Figure (5) shows α and A in units of (J/h) α as functions of the BEC temperature T , which were obtained from the evolution (according to the GME) of an impurity initially localized at j = 0. We see that the exponent α drops from α ≈ 2 at zero temperature to α ≈ 1 at high temperatures k B T ≫ E p , thereby clearly indicating the crossover from coherent to diffusive transport. We note that the transition from coherent to diffusive hopping takes place in a temperature regime accessible to experimental study and therefore, importantly, may be observable.

Atomic current across a tilted lattice
The inelastic phonon scattering responsible for the incoherent hopping of the impurities also provides the necessary relaxation process required for the emergence of a net atomic current across a tilted optical lattice. This is in contrast to coherent Bloch oscillations, which occur in an optical lattice system in absence of incoherent relaxation effects or dephasing [47]. As pointed out in [42], the dependence of the atomic current on the lattice tilthω B changes from ohmic conductance to NDC in agreement with the theoretical model for electron transport introduced by Esaki and Tsu [45].
To demonstrate the emergence of a net atomic current, and, in particular, to show that impurities exhibit NDC, we consider the evolution of a localized impurity atom in a one-dimensional system. With the impurity initially at site j = 0 we determine its average position x d = j aj P j (t d ) after a fixed drift time t d of the order ofh/J. This allows us to determine the drift velocity v d = x d /t d as a function of the lattice tilt hω B and the temperature T of the BEC. In analogy with a solid state system, the drift velocity v d and the lattice tilthω B correspond to the current and voltage, respectively. Figure 6a shows the voltage-current relation at different temperatures T , which was obtained by numerically solving the GME with the memory function in Eq. (40). We see that for a small lattice tilt the system exhibits ohmic behavior v d ∼hω B , whereas for a large lattice tilt the current decreases with increasing voltage as v d ∼ 1/(hω B ), i.e. the impurities feature NDC.
Following [42] we describe the voltage-current relation for the impurities by an Esaki-Tsu-type relation whereṽ 0 =Ja/h the characteristic drift velocity, τ the effective relaxation time of the impurities and γ a dimensionless prefactor. Fitting the Esaki-Tsu-type relation in Eq. (43) to the numerical data allows us to extract the parameters τ and γ, both depending on the BEC temperature. As can be seen in Fig. 6b, the effective relaxation time τ decreases with increasing BEC temperature T and is of the order of τ 0 =h/(gn 0 ). The significance of 1/τ is that of an average collision rate between the impurities and Bogoliubov excitations, which would allow us, for example, to formulate the problem of transport in terms of a classical Boltzmann equation for the distribution function of the impurities [35]. The prefactor γ varies only slightly since the exponential dependence of the current on the temperature is accounted for byṽ 0 . We note that independently of ω B τ the maximum drift velocity is given by γṽ 0 . The Esaki-Tsu-type relation in Eq. (43) reflects the competition between coherent and incoherent, dissipative processes. In the collision dominated regime ω B τ ≪ 1, inelastic scattering with phonons destroys Bloch oscillations, whereas in the collisionless regime ω B τ ≫ 1 the evolution of the impurities is mainly coherent, i.e. Bloch oscillations of the impurities lead to a suppression of the net current. The crossover between the two regimes is most pronounced at zero temperature, whereJ is comparable to J. However, the change from ohmic to negative differential conductance is identifiable even at finite temperatures and thus should be observable in an experimental setup similar to the one used in [12].

Conclusion
We have studied the transport of impurity atoms in the strong-coupling regime, where the interaction energy due to the coupling between the BEC and the impurities dominates their dynamics. Within this regime, we have formulated an extended Hubbard model describing the impurities in terms of polarons, i.e. impurities dressed by a coherent state of Bogoliubov phonons. The model accommodates hopping of polarons and the effective off-site impurity-impurity interaction mediated by the BEC.
Based on the extended Hubbard model we have shown from first principles that inelastic phonon scattering results in a crossover from coherent to incoherent hopping and leads to the emergence of a net atomic current across a tilted optical lattice. In particular, we have found that the dependence of the current on the lattice tilt changes from ohmic conductance to negative differential conductance for sufficiently low BEC temperatures. Notably, this transition is accurately described by an Esaki-Tsu-type relation with the effective relaxation time of the impurities as a temperature-dependent parameter.
Using the techniques introduced in this paper, qualitatively similar phenomena can also be shown to occur for fermionic impurities and, moreover, for impurities of different species [55]. For instance, in the case of two impurity species A and B with the couplings κ A > 0 and κ B < 0, respectively, the effective off-site impurity-impurity interaction is attractive for the same species, but repulsive for different species. In either case, observation of the phenomena reported in this paper lies within the reach of current experiments, which may give new insight into the interplay between coherent, incoherent and dissipative processes in many-body systems.
Starting point of the derivation of the GME is the Liouville-von Neumann equation [38] ∂ρ(t) ∂t = Lρ(t) , (B.1) whereρ(t) is the density matrix of the system expressed in the eigenbasis ofĥ 0 . The Liouville operator L is defined by L = L 0 + L I , with L 0 = −i/h [ĥ 0 , · ] and where P is the projection operator on the relevant part ofρ(t) and Q = 1 − P the complementary projection operator to P. The time evolution of the relevant part Pρ(t) is governed by the Nakajima-Zwanzig equation [41] ∂ ∂t Specifically, the projection operator P employed in the derivation of the GME is defined by Pρ(t) =ρ B ⊗ D Tr Bρ (t) [38,39]. Here,ρ B is the density matrix of the condensate in thermal equilibrium, Tr B is the trace over the condensate degrees of freedom and D is the projection operator on the diagonal part. Thus̺ A (t) = D Tr Bρ (t) is the diagonal part of the reduced density matrix. We note that the trace over the Bogoliubov phonon states in the definition of P introduces irreversibility into the system. The Nakajima-Zwanzig equation can be simplified given that PL 0 = L 0 P = 0 and PL I P = 0 for the definitions of L 0 , L I and P above. In addition, we assume that the initial density matrix of the total system has the formρ(0) =ρ B ⊗̺ A (0) so that Qρ(0) = 0, and hence the inhomogeneous term in Eq. (B.2) vanishes. Taking these simplification into account and using Tr B Pρ =̺ A we find from Eq. (B.2) that̺ A (t) evolves according to [38] ∂ ∂t̺ with the memory kernel The memory kernel to second order in ζ, i.e. dropping L I in the exponent in Eq. (B.4), can be expressed in tetradic form as [39,41] with the partition function of the condensate and wherehΩ i,j is the energy difference between the impurity configurations i and j, andhω {Nq} is the energy of the phonon configuration {N q }.
The explicit expression for the memory function W i,j (s) = K 0 ii,jj (s) for the Hamiltonianĥ 0 +ĥ I has been evaluated in [39,40] based on Eq. (B.4), however, we here give an alternative derivation of W i,j (s) starting from Eq. (B.5). The evaluation of K 0 ii,jj (s) in Eq. (B.5) can be separated into a phonon part and an impurity part, where the latter is given by Thus, for the phonon part we only have to consider operators of the formX † jXj±1 , which can be written in terms of displacement operators D(β) = exp(βb † − β * b ) aŝ with β j,q = M * j,q − M * j±1,q and Φ j,q the corresponding phase. This allows us to treat each phonon mode in Eq. (B.5) separately, and the problem reduces to the summation To evaluate the sum over n we use the fact that the following relation for generalized Laguerre polynomials holds [57] ∞ n=0 n! L γ n (x)L γ n (y)z n Γ(n + γ + 1) provided that |z| < 1. Here, Γ(x) is the Gamma function and I γ (x) are modified Bessel functions of the first kind. Using the relation in Eq. (B.12) we find that the sum over n yields with (x) n = x(x + 1)(x + 2) · · · (x + n − 1) the Pochhammer symbol, we find that which is identical to expression (B.11) except for the upper limit in the sum over l.
Using relation (B.12) again and discarding the double counting of l = 0 we find that the total sum in Eq. (B.9) is given by To evaluate the sum over l we use the identity [56] ∞ l=−∞ I l (x)t l = exp and find that Eq. (B.17) equals which can be written as where we used Z B = q Z q .
Appendix C. Self-trapping Impurities immersed in a BEC get self-trapped for sufficiently strong impurity-BEC interactions, even in the absence of an additional trapping potential [22][23][24]. Based on the results for static impurities in Section 3 we now show that for the parameter regime considered in this paper self-trapping effects can be neglected.
As pointed out by Gross [58] the coupled equations describing the impurity and the condensate in the Hartree approximation are given by µ b ψ(r) = −h 2 2m b ∇ 2 ψ(r) + g|ψ(r)| 2 ψ(r) + κ|χ(r)| 2 ψ(r) (C.1) εχ(r) = −h 2 2m a ∇ 2 χ(r) + κ|ψ(r)| 2 χ(r) , where χ(r) is the wavefunction of the impurity, m a is the impurity mass and ε the impurity energy. To determine whether the impurity localizes for given experimental parameters one has, in principle, to solve the coupled equations for ψ(r) and χ(r) [24]. Alternatively, as suggested in [23], we use the Gaussian mode-function χ j,σ (r) in Eq. (4) as a variational wavefunction for the impurity, with the spread σ as a free parameter, and minimize the total energy of the system in the regime α ≪ 1, where the linearization of the GP equation is valid. For a homogenous condensate, the potential energy of the impurity and the BEC is (−E p ), and thus adding the kinetic energy of the impurity yields the total energy The impurity localizes if E(σ) has a minimum for a finite value of σ, which depends on the dimensionless quantity In one dimension, we find by asymptotically expanding E(σ) in the limit σ/ξ ≫ 1 that there exists a self-trapping solution for arbitrarily small α ′ and that For the two-dimensional case, asymptotically expanding E(σ) in the limit σ/ξ ≫ 1 yields a critical value α ′ c = 2π, above which self-trapping occurs. The corresponding spread of the self-trapping solution σ 2D diverges close to α ′ c , which validates the asymptotic expansion of E(σ). In three dimensions, numerical minimization of E(σ) shows that the critical value is α ′ c ≈ 31.7, and the corresponding self-trapped state is highly localized with σ 3D ≈ 0.87 ξ.
According to Eq. (C.5) the spread of the self-trapping solution exceeds several lattice spacings for the parameter regime considered in this paper. In other words, the spread σ 1D is much larger than the harmonic oscillator length σ ℓ in practice, and hence self-trapping effects are indeed small.