Quantum Simulations of Extended Hubbard Models with Dipolar Crystals

In this paper we study the realization of lattice models in mixtures of atomic and dipolar molecular quantum gases. We consider a situation where polar molecules form a self-assembled dipolar lattice, in which atoms or molecules of a second species can move and scatter. We describe the system dynamics in a master equation approach in the Brownian motion limit of slow particles and fast phonons, which we find appropriate for our system. In a wide regime of parameters, the reduced dynamics of the particles leads to physical realizations of extended Hubbard models with tuneable long-range interactions mediated by crystal phonons. This extends the notion of quantum simulation of strongly correlated systems with cold atoms and molecules to include phonon-dynamics, where all coupling parameters can be controlled by external fields.


Introduction
The recent preparation of cold ensembles of homonuclear and heteronuclear molecules in the electronic and vibrational ground state [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30] has opened the door to a new chapter in the theoretical and experimental study of trapped quantum degenerate gases [31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53]. Heteronuclear molecules, in particular, have large electric dipole moments associated with rotational excitations, and the new aspect in quantum gases of cold polar molecules is thus the large, anisotropic dipole-dipole interactions between the molecules which can be manipulated via external DC and AC fields in the microwave regime [31; 32; 33; 34; 52; 53; 54]. In combination with reduced trapping geometries, this promises the realization of novel quantum phases and quantum phase transitions, for example in the case of bosons the transition from a dipolar superfluid to a crystalline regime [52]. The theory of dipolar quantum gases, and various aspects of strongly correlated systems of polar molecules has been recently reviewed by Baranov [55] and by Pupillo et al. [56] Below we will extend these theoretical studies to mixtures of atomic and dipolar molecular quantum gases. More specifically, we will be interested in a situation where polar molecules form a self-assembled dipolar lattice, in which atoms can move and scatter. This scenario leads to a new physical realization of a Hubbard model where atoms see the periodic structure provided by the crystal formed by the polar molecules. In contrast to the familiar case of the realization of Hubbard models with cold atoms in optical lattices, where standing laser light waves produce a fixed periodic external lattice potential, self-assembled dipolar lattices have their own lattice dynamics represented by phonons. Thus atoms moving in dipolar crystals give rise to Hubbard models which include both (i) atom-phonon couplings, and (ii) atom-atom interactions. This extends the notion of quantum simulation of strongly correlated systems with cold atoms to include phonon-dynamics, where both the atomic and phonon coupling parameters can be controlled by external fields. We note that atommolecule mixtures arise naturally in photoassociation experiments, where a two-species atomic quantum degenerate gas is partially transferred to ground state molecules via formation of highlying Feshbach molecules, followed by a Raman transfer to the ground state. Similar Hubbard models result also in mixtures of two unbalanced species of polar molecules, where the first molecular species forms a crystal while the second species provides the extra particles hopping in the self-assembled dipolar lattice.
In a recent work [57], we have shown that the dynamics of atoms and molecules embedded in dipolar crystals is conveniently described by a polaronic picture where the particles are dressed by the crystal phonons. Standard treatments of polaron dynamics show a competition between coherent and incoherent hopping of a particle in the lattice. The former corresponds to tunneling of a particle from one site to to the next, while carrying the lattice distortion, without changing the phonon occupation. The latter corresponds to thermally activated particle hopping, related to incoherent hopping events where the number of phonons changes in the hop. That is, the polaron loses its phase coherence via the emission or absorption of phonons [58]. The physics of polarons has a long history, dating back to the seminal work of Landau [59]. Excellent reviews on the subject can be found e.g. in [60; 61; 62; 63; 64; 65; 66; 67]. Here we take the simple approach of placing these coherent and incoherent processes in the natural framework of a master equation treatment of the system dynamics, Figure 1. A dipolar crystal of polar molecules in 2D (a) and 1D (b,c) provides a periodic lattice Vcp for extra atoms or molecules giving rise to a lattice model with hoppingJ and long-range interactionsṼ i,j (see text). (a) In 2D a triangular lattice is formed by polar molecules with dipole moment dc perpendicular to the plane. A second molecular species with dipole moment dp ≪ dc moves in the honeycomb lattice Vcp (darker shading corresponds to deeper potentials). (b) A 1D setup with atoms scattering form a dipolar crystal with lattice spacing a. (c) A 1D dipolar crystal provides a periodic potential for a second molecular species moving in a parallel tube at distance b.
where the phonons are treated as a thermal heat bath. We extend our work [57] by re-deriving the results of [58] in the master equation context, and calculating additional corrections to the coherent and incoherent time evolution for our atomic and molecular mixtures. Since we find that in the latter the polaron dynamics is typically slow compared to the characteristic time evolution of the bath, we specialize to the Brownian motion limit of the master equation [68; 69]. We show that for the models of interest and low-enough temperatures, corrections to the coherent time evolution of the polaron system are small, and thus the dynamics of the dressed particles is well described by an effective extended Hubbard model in a wide range of realistic parameters.
The paper is organized as follows. In Section 2 below we derive the Brownian motion master equation for our system, providing explicit expressions for the coherent and incoherent polaron dynamics (details of the derivation can be found in Appendix C). In Section 3 we review the realization of dipolar crystals in two and one dimensions. In Section 4 we provide details of the implementation of polaronic extended Hubbard models for two one-dimensional configurations with atoms and molecules.

Dynamics of particles trapped in a crystal of polar molecules
We show below in Section 4 and Appendix A that in a wide range of system parameters the dynamics of atoms or molecules which move in a self-assembled lattice of dipoles is well described by the following Hamiltonian ω q,λ a † q,λ a q,λ + q,λ,j M q,λ e iqr 0 j c † j c j (a q,λ + a † −q,λ ), (1) where the brackets indicate sums over nearest neighbors. The latter is a singleband Hubbard model for the particles coupled to the acoustic phonons of the lattice. The first two terms on the r.h.s. of (1) describe the hopping of particles between neighboring sites of the lattice with a tunneling rate J, and the density-density interactions with strength V i,j for particles at site i and j, respectively. Here, c † i (c i ) denotes the creation (annihilation) operator for a particle at site i. These particles can be either fermions or bosons. The third term describes the excitations of the crystal given by acoustic phonons, where a † q,λ (a q,λ ) creates (destroys) a phonon with quasimomentum q in the mode λ, with dispersion relation ω q,λ . The last term in (1) is the particle-phonon coupling, which is of the density-displacement type, with coupling constant M q,λ . The microscopic derivation of (1) is detailed in Appendix A and Appendix B. For the models we consider (see Section 3 and Section 4 below), we find that: (i) The phonon coupling can largely exceed the hopping rate J, which precludes a naive treatment of the particle-phonon coupling as a (small) perturbation; (ii) We are generally interested in the so-called non-adiabatic regime, where the characteristic phonon frequency ω D (the Debye frequency) is typically (much) larger than the average kinetic energy of the particles ∼ J, that is ω D ≫ J, (see Section 4). This is due to the fact that, unlike e.g. the case of electrons in ionic crystals, in our system the mass of the particles and of the crystal dipoles are comparable, and the crystals can be made stiff.
Below, we derive a master equation for the dynamics of the particles only, while the crystal phonons are treated as a thermal heat bath. The master equation in the Markovian limit has the general forṁ where H S denotes a Hamiltonian term for the coherent time evolution of the reduced density matrix ρ S for the particles, while D describes dissipative processes responsible for the incoherent dynamics. The coherent processes correspond to hopping of a particle dressed by the crystal phonons (a polaron) and polaron-polaron interactions, while the latter are thermally activated incoherent hopping, where a particle loses its phase coherence by emitting or absorbing a phonon in the hopping process. As suggested by points i) and ii) above, we derive (2) in a perturbative, strongcoupling, approach where the hopping rate J in (1) acts as the small parameter. In addition, we work in the Brownian motion limit of the master equation, where the characteristic time evolution of the system τ S ∼ max(1/J, 1/V ij ) is (much) slower than the characteristic time evolution of the thermal heat bath τ B ∼ 1/ ω D . We obtain explicit expressions for the coherent and incoherent (dissipative) contributions to the system time evolution, and show that the latter can be made negligible in a wide range of realistic parameters for our models. We find that the effective Hamiltonian for the coherent dynamics of the particles only has the form which corresponds to an extended Hubbard model. HereJ andṼ ij are the tunneling rate for the particles and their mutual interactions, respectively, which are modified with respect to their bare values in (1) by the coupling to the crystal phonons.

Lang-Firsov Transformation (Polaronic picture)
In our system (see Section 4), we find that particles are slow and strongly coupled to the crystal phonons, with Hamiltonian (1). Following Lang and Firsov [61], it is convenient to change from this picture of bare particles strongly coupled to phonons to an equivalent one, where particles freely hop in the lattice while carrying the lattice distortion. This corresponds to dressing the particles with the lattice phonons, and the dressed particles are named polarons [58]. This is achieved by performing a unitary transformation of H asH = U HU † with and u q,λ = M q,λ / ω q,λ displacement amplitudes. In the new picture the Hamiltonian (1) reads [58 where now c j (c † j ) are the annihilation (creation) operators of a polaron located at site j. Here, is the displacement operator for the crystal molecules due to the presence of a particle located at site j. The energy E p in (5) is the polaron shift [58] and N p ≡ j c † j c j is the total number of particles. The quantitỹ is a modified particle-particle interaction, which comprises two terms: The first is the original (bare) interaction, while the latter is the phonon mediated particle-particle interactionṼ (1) ij , which in general (i) is long-ranged and (ii) can be comparable in strength to the bare interactions V ij , see Section 4.
For J = 0 the new Hamiltonian (5) is diagonal and describes interacting polarons and independent phonons. The latter are vibrations of the lattice molecules around new equilibrium positions with unchanged frequencies. For small finite J, we can treat the first term in (5) perturbatively, which is the starting point for the master equation approach detailed in the following sections.

Effective dynamics of polarons inside a thermal crystal
In this section we derive a master equation describing the effective dynamics of the polarons embedded in the crystal, which we assume to be in thermal equilibrium. That is, we consider the phonons to provide a heat bath at temperature T with a density matrix given by Heren q,λ (T ) denotes the (Bose-Einstein) phonon distribution at temperature T , where O ≡ tr B {Oρ 0 B } is the thermal expectation value of the bath operator O, with tr B denoting the trace over the bath degrees of freedom.
We split the Hamiltonian (5) into three parts asH = H S + H B + H I with Here, H S is the "reduced system" Hamiltonian describing the dynamics of polarons, with a hopping rate J X † i X j and interactionsṼ ij for two polarons at site i and j. The Hamiltonian H B is the Hamiltonian for the bath, and Hamiltonian H I gives the interactions between the system and the bath. In writing (10a)-(10c) we have conveniently added a term −J i,j c † i c j X † i X j to H S and subtracted it from H I . This ensures that the thermal average over the interaction Hamiltonian vanishes, H I = 0 [69].
The expectation value X † i X j in equation(10a) reads (see Appendix C.1) Equation (10a) shows that J ≡ J X † i X j = Je −ST plays the role of a phonon-modified tunneling rate, which is exponentially suppressed for S T > 0. In the following we will often distinguish between two regimes, i.e. a weak coupling regime where S T ≪ 1 (andJ ≃ J) and a strong coupling regime where S T ≫ 1 (andJ ≪ J).

The Quantum Brownian Motion Master Equation (QBMME)
In this section we derive a master equation for the coherent and incoherent dynamics of the reduced system of interacting polarons in the Brownian motion limit, where the system time evolution, τ S , is slow compared to the characteristic evolution time of the bath, τ B . That this approach may provide a reasonable description of the system is suggested by the fact that in a wide range of parameters for atoms and molecules trapped in dipolar crystals we find τ S ∼ max(1/J, 1/Ṽ ij ) ≫ τ B ∼ 1/ω D . We will thus derive the coherent time evolution (3), and provide explicit analytic expressions for the corrections both to the coherent and incoherent dynamics. In Section 4 we show that for the models of interest these corrections are in fact negligible in a wide regime of parameters, which provides an a posteriori self-consistency check for the approximations made here.
The time evolution for the density matrix of the entire system ρ(t) comprising the (polaronic) reduced system and the heat bath is dictated by the Liouville-von Neumann equatioṅ whereÃ(t) ≡ e iH0t/ Ae −iH0t/ denotes an operator A in the interaction picture with respect to H 0 = H S + H B . We assume that the reduced system and the heat bath are uncoupled at time t = 0, so that ρ(t = 0) can be written as the tensor product ρ(0) = ρ S (0) ⊗ ρ 0 B , with ρ S (0) and ρ 0 B the density matrices of the reduced system and the bath, respectively. The condition J ≪ ω D , which states that the interaction is a weak perturbation, forms the basis for a Born-Markov approximation with the phonons a finite temperature heat bath, providing the following master equation for the reduced density operator of the polaronsρ where ξ kl ij (τ, T ) are the (complex) bath correlation functions, with the quantities ∆ kl ij (T ), Γ kl ij (T ), δ kl ij (T ) and γ kl ij (T ) given by Here Re[f (x)] and Im[f (x)] denote real and imaginary part of f (x). The first term inside the commutator on the r.h.s. of (16) describes the coherent time evolution for the reduced system, with H S the effective Hubbard Hamiltonian (10a). The terms proportional to ∆ kl ij (T ) in (16) are self-energies, which in the single-polaron limit provide both a shift to the ground-state energy, and next-nearestneighbor hopping terms [71]. In addition, in the many-polaron problem they can provide offsite polaron-polaron interactions, whose strength will be evaluated in Section 2.2.3 below. The terms on the r.h.s. of (16) proportional to Γ kl ij (T ) are related to incoherent hopping events where the number of phonons changes in the hop. That is, the polaron loses its phase coherence via the emission or absorption of phonons. These processes are thermally activated and can dominate over the coherent hopping rateJ for large enough temperatures T [58].
The terms proportional to γ kl ij (T ) and δ kl ij (T ) are (small) corrections to the coherent and incoherent time evolution, respectively, which derive from the term proportional to τ in the expansion (15), and thus correspond to the Brownian motion corrections to the time evolution of the reduced system.
We notice that the terms proportional to (17a)-(17d) in (16) do not identify directly the corrections to the coherent time evolution given by (10a), because equation (16) is not diagonal. Instead these terms are elements of matrices whose eigenvalues are the corrections. The diagonalization of the master equation can be done analytically in the single polaron limit and is shown below in Section 2.2.4.
In Section 2.2.3 below we provide analytic expressions for the coherent and incoherent corrections to the coherent time evolution given by H S , cf. (10a). In Section 4 we show that these corrections are in fact negligible in a wide range of realistic parameters for atoms and molecules, and thus H S properly describes the dynamics of polarons inside the dipolar crystal.

Correlation Functions:
The bath correlation functions ξ kl ij (τ, T ) appearing in the master equation (16) are computed in Appendix C.1 and read Here we have introduced the spectral density where V BZ is the volume of the Brillouin zone, q d−1 denotes all components of the quasi-momentum vector except q x , andḡ kl ij (q) reads ij is a decaying function of the time τ with max(|Φ kl ij |) = 2S T , the actual decay rate depending strongly on the spectral density of the model [58]. In Sections 2.2.3 below we show that at small finite temperatures this decay rate is fast enough to ensure that the corrections to the coherent time evolution in the QBMME (16) for our 1D models are finite and small. This provides for an a posteriori check of the applicability of the QBMME to the polaron problem.

Self energies and dissipation rates in the strong and weak coupling limits
In this subsection we provide analytical approximate expressions for the matrix elements ∆ kl ij , Γ kl ij , γ kl ij and δ kl ij in the limits of strong and weak coupling S T ≫ 1 and S T ≪ 1, respectively. The details of the performed approximations are given in Appendix C, while explicit results for our 1D models are shown below in Section 4. Here we concentrate in particular on the leading self-energy corrections to the coherent time evolution determined by H S in (10a) in the two regimes.
In accordance with literature [61], we find that in the strong coupling limit, S T ≫ 1, the self-energies are suppressed by a factor of the order of (J/E p ) 2 , with E p the polaron shift (7). In addition we find that in the weak coupling limit, S T ≪ 1, these corrections are strongly suppressed by a factor (J/ ω D ) 2 and as a consequence, in Section 4 we show that they are negligible in a wide range of realistic parameters for our models.
It is shown in Appendix C that all matrix elements ∆ kl ij , Γ kl ij , γ kl ij and δ kl ij in the strong coupling limit are strongly suppressed by an exponential factor of the order e −2ST , unless i = l and j = k, [58; 61]. This makes sense, since processes for which i = l and j = k correspond to a double hop of a polaron, each one being suppressed by the factor e −ST . The conditions i = l and j = k describe a "swap"-process, which corresponds, e.g. in the case of ∆ kl ij , to the virtual hop of a polaron from its current position to a neighboring site and back.
The various corrections read where Erfi denotes the error function and where we have introduced the quantities Equations (22a)-(22d) result from an expansion of the function −Φ kl ij (τ, T ) appearing in the exponent of ξ kl ij (τ, T ), cf. (18), up to second order in the time τ around its maximum. For a vanishing A T , which corresponds to neglecting second order terms in the expansion of −Φ kl ij (τ, T ), we find A simple estimate of B, see Appendix C.2, shows that it is of the order of B ∼ 2E p / ω D for a sufficiently strong coupling and therefore ∆ 10 01 /E p ∝ (J/E p ) 2 and γ 10 01 ∝ (J/E p ) 2 . This dependence is known in literature as the "1/λ" strong coupling expansion, where λ = E p /J, see [61].
In the weak coupling limit, S T ≪ 1, the correlation functions ξ kl ij can be expanded to first order in Φ kl ij (τ, T ), which leads to the following expressions for the matrix elements as detailed in Appendix C.2. Here, P dx denotes the Cauchy principal value integral. In contrast to the strong coupling regime, as shown in Appendix C here the (small) parameter that gives the approximate size of the corrections is J/ ω D and not J/E p . Explicit results for the corrections in this limit are given in Section 4 below. The ratio (J/ ω D ) 2 also determines the size of the incoherent processes Γ kl ij . This can be seen by performing a low temperature approximation of expression (27b) above, for which we find Γ kl ij (T ) We notice that, in addition to being proportional to the small factor (J/ ω D ) 2 , these corrections depend linearly on temperature, a result also found in [58]. Because dipolar crystals can have a large Debye frequency, in our models we find J/ ω D ≪ 1 and thus corrections to the coherent time evolution determined by H S are small.

Diagonalization of the master equation
As noted above the corrections to the master equation Γ kl ij (T ), ∆ kl ij (T ), γ kl ij (T ) and δ kl ij (T ) are just matrix elements, while the actual energies and rates are obtained by diagonalizing the various terms in (16). In the following we sketch how to perform this diagonalization in the simple case of a single polaron for ∆ kl ij and Γ kl ij , while similar computations for all coherent and incoherent corrections are shown in Appendix C.3. In the case of a single polaron the terms proportional to ∆ kl ij are easily diagonalized, since the eigenergies for a particle in a lattice are readily determined by the energies associated with the various quasimomenta. That is where the sums over m, n range over basis vectors in the lattice, and i c † The eigenvalues can now be directly read-off from (29) To determine the rate of the dissipative term in (16) we notice that this term can be written as 1 The amplitude of the dominant rate can be now estimated from Calculations similar to the ones above lead to the eigenvalues γ q and δ q for the Brownian motion corrections (see Appendix C.3). In one dimension, which is relevant for the models of Section 4 below, we find In conclusion, we note that the corrections to the coherent time evolution given by ∆ kl ij , Γ kl ij , γ kl ij and δ kl ij can be made small both in the strong and weak coupling regimes, by ensuring that the ratios J/E p and J/ ω D are small, respectively. The smallness of these corrections provides for an a posteriori check of the applicability of the Brownian motion master equation approach to the polaron problem.

Crystals of polar molecules
In this section we briefly review how to realize self-assembled crystals of polar molecules. Following [52; 57], here we focus on crystals in two and one dimensions.  [72]. The crossover to the unstable regime for small replusion and finite confinement ω ⊥ is indicated (hatched region).
However, self-assembled crystals in three-dimensions can also be realized as detailed in [53].
Two-dimensional crystals: We consider a system of cold polar molecules in the presence of a DC electric field under strong transverse confinement, as illustrated in figure 2(a). A weak DC field along the z-direction induces a dipole moment d c in the ground state of each molecule. Thus, the molecules interact via the dipole-dipole This interaction is purely repulsive for molecules confined to the x, y-plane, while it is attractive for z > r/ √ 3, leading to an instability towards collapse in the many body system. In reference [52] it is shown that this instability can be suppressed for a sufficiently strong 2D confinement along z, as provided, for example, by a deep optical lattice with frequency ω ⊥ (blue arrows in the figure). In fact, a strong confinement with ω ⊥ D/a 3 , with a the mean interparticle distance, confines the molecules to distances larger than where V 3D dd (r) is purely repulsive, and thus the system is collisionally stable. Here, m c is the mass of the molecules. The 2D dynamics in this pancake configuration is described by the Hamiltonian which is obtained by integrating out the fast z-motion. Equation (37) is the sum of the 2D kinetic energy in the x,y-plane and an effective repulsive 2D dipolar interaction with ρ ij ≡ (x j − x i , y j − y i ) a vector in the x, y-plane. Tuning the induced dipole moment d c drives the system from a weakly interacting gas (a 2D superfluid in the case of bosons), to a crystalline phase in the limit of strong repulsive dipole-dipole interactions. This crystalline phase corresponds to the limit of strong repulsion where particles undergo small oscillations around their equilibrium positions. The strength of the interactions is characterized by the ratio r d of the interaction energy over the kinetic energy at the mean interparticle distance a This parameter is tunable as a function of d c from small to large r d . A crystal forms for where the interactions are dominant [52; 73]. For a dipolar crystal, this is the limit of large densities, as opposed to Wigner crystals with 1/r-Coulomb interactions. Figure 2(b) shows a schematic phase diagram for a dipolar gas of bosonic molecules in 2D as a function of r d and temperature T . In the limit of strong interactions r d > r c the polar molecules are in a crystalline phase for temperatures [72] the crystal recoil energy, typically a few to tens of kHz. The configuration with minimal energy is a triangular lattice with spacing a L = (4/3) 1/4 a see [52]. Excitations of the crystal are acoustic phonons with Hamiltonian given by equation (10b), and characteristic Debye frequency One-dimensional crystals: One dimensional analogues of the 2D crystals can be realized by adding an additional in-plane optical confinement to the configuration of figure 2(a) [57; 74; 75]. For large enough interactions r d ≫ 1, the phonon frequencies have the simple form while the classical melting temperature can be estimated to be of the order of T m ≃ 0.2r d E R,c /k B , see [75].
Finally, for a given induced dipole d c the ground-state of an ensemble of polar molecules is a crystal for mean interparticle distances a min a a max , where a max ≡ d 2 c m/ 2 r c corresponds to the distance at which the crystal melts into a superfluid. For SrO (RbCs) molecules with the permanent dipole moment d c = 8.9D (d c = 1.25D), a min ∼ 200nm(100nm), while a max can be several µm. Since for large enough interactions the melting temperature T m can be of order of several µK, the self-assembled crystalline phase should be accessible for reasonable experimental parameters using cold polar molecules.

Specific implementations with atoms and molecules in dipolar crystals
In this section we consider a mixture of two species of particles confined to one dimension. The first species of particles comprises (strongly interacting) molecules with dipole moment d c , forming a one-dimensional crystal. The second species of particles can be either atoms [see figure 1(b)] or molecules of a second species, with dipole moment d p ≪ d c [see figure 1(c)]. The former interact with the crystal molecules via a short range pseudopotential proportional to an elastic scattering length a cp , while the latter interact with the crystal molecules via long-range dipole-dipole interactions. For both configurations, we obtain explicit expressions for all parameters characterizing the coherent and the incoherent dynamics in the system. These one dimensional setups can be readily generalized to two dimensions. . We show the dispersion for a 1D and 2D dipolar crystal as a function of the quasimomentum in units of the crystal recoil energy E R,c . In 1D (a) the dispersion is peaked at the zone border and tends to zero ∝ q for small momenta. The 2D dispersion has two acoustic branches, a longitudinal and a transversal one. We plot the two dispersions as a function of the quasimomentum, choosing a path in the first Brillouin zone that is outlined in the inset of the figure.

Neutral atoms moving inside a crystal tube
As a first realization, we consider a setup, where an ensemble of neutral atoms is confined by an optical trap to the same 1D tube as the dipolar crystal, see Figure 1(b), say along x. For simplicity we assume the trap for the neutral atoms and the polar molecules to have the same harmonic oscillator frequency ω ⊥ .
An atom and a molecule inside the tube interact via a short range potential, which we model in the form of an effective 1D zero range potential, Assuming that the 3D scattering length a cp is (much) smaller than the harmonic oscillator length of the traps, a cp ≪ a p,⊥ = ( /m p ω ⊥ ) 1/2 , the effective 1D coupling strength is given by g cp ≈ 2 ω ⊥ a cp . In the following we focus on positive scattering lengths, a cp > 0, corresponding to a situation where the atoms and molecules effectively repel each other, cf. g cp > 0.

Tight binding limit and Hubbard models.
For a "frozen" crystal, i.e. without phonons, the molecules are at their equilibirum positions, ja + a/2, which provide for a static periodic potential for the atoms, The dynamics for a single neutral atom is then determined from the static Hamiltonian H p = p 2 /2m p + V p (x), corresponding to the Kronig-Penney model with a potential comb of strength g cp . Its energy spectrum is given in the form of a bandstructure, E n,q (with band-index n = 0, 1, . . .), which is obtained from where E R,p ≡ 2 π 2 /2m p a 2 denotes the recoil energy of an atom. We denote the bandwidth of the lowest band by 4J ≡ E 0,π/a − E 0,0 = E R,p − E 0,π/a , and the gap to the first band by ∆ ≡ E 1,π/a − E 0,π/a = E 1,π/a − E R,p . The latter are both shown in figure 4(a) as a function of the coupling strength g cp . We notice that for g cp E R,p a/2 (indicated by a vertical dashed line), the gap exceeds the band-width, ∆ > 4J. In the tight binding limit, cf. g cp ≫ E R,p a, the dispersion relation in the lowest band becomes E 0,q ≈ 4J sin 2 (qa/2), with J ≈ 2E 2 R,p /π 2 g cp + O(aE R,p /g cp ) and ∆ ≈ 3E R,p . The Wannier-functions for a particle become localized at site j and are approximated by w j (x) ≈ cos[π(x − ja)/a]/ a/2 for |x − ja| < a/2 and zero otherwise.
Obtaining the tight-binding limit requires that the ratio, g cp /aE R,p = 2a cp a/π 2 a 2 p,⊥ (largely) exceed the value ≈ 1/2. We notice that for current state-of the art optical traps, one can achieve harmonic oscillator lengths as small as a p,⊥ ∼ 20 nm, and taking a "typical" 3D scattering length of a cp ∼ 100a 0 ≈ 5 nm, we get that g cp /aE R,p 1/2 is attained for lattice spacings a π 2 a 2 ⊥ /a cp ∼ 200nm.
Analogous to the atom-molecule interactions, we model the interactions between two neutral atoms by a contact potential with a coupling strength given by g pp ≈ ω ⊥ a pp for the 3D atom-atom scattering length a pp ≪ a ⊥ . Then in the tight-binding limit the atom-atom interactions reduce to repulsive onsite energy shifts only The dynamics for an ensemble of bosonic atoms in the crystal is then described by a single band Bose-Hubbard model with hopping rate J and onsite repulsion V ii , provided that V ii ≪ ∆. On the other hand for an ensemble of (spin-polarized) fermionic atoms, we notice that the system reduces to a lattice model with hopping rate J and no interactions. The atoms couple to the crystal molecules via a density-displacement interaction [58]. To first order in the displacement (see Appendix B.2), the coupling constant reads whereṼ cp is the Fourier transform of the atom-molecule interaction potential, and β q is the Fourier transform of the square of the Wannier-functions. For g cp ≫ aE R,p the latter is The latter decreases with increasing q from β 0 ≡ 1 to β π/a = 8/3π ≈ 0.85. The (monotonical) dependence of the coupling constant M q on the quasimomentum is shown in figure 4(b). In particular, it has a maximum, M π/2 ≈ (8g cp /3a 2 )(2 /N m c ω D ) 1/2 , at the band-edges, while for small quasi-momenta q it vanishes as Finally, let us address the validity of the single band approximation in the Hubbard model (1), when coupled to phonons. For simplicity let us consider the limit of vanishing interactions V ii = 0 and a weak coupling M q : We notice that, for ω D < 4J + ∆ the second band is gapped from the branch of acoustic phonons, and therefore higher band excitations are (strongly) suppressed and can be neglected. Since E R,p ≤ 4J + ∆ ≤ 3E R,p , this requires a mass ratio m p /m c 3/ √ 2r d . While this for a soft crystal with r d ∼ 1 merely implies that m p < m c , for a stiff crystal with r d ∼ 200 this requires m p 0.15m c , which is quite restrictive. In the latter case, that is for a stiff crystal and comparable masses, cf.
√ 2r d m p /3m c > 1, we notice that the first excited band would cut the phonon branch at a frequency ω ⋆ ∼ 4J + ∆. However, in this regime a single band model may still hold, if one restricts the initial phonon population to sufficiently low temperatures, i.e. for k B T ≪ ∆. In the tight binding limit this requires temperatures T ≪ (3m c / √ 2r d m p ) × ω D /k B which even for a stiff crystal with r d ∼ 200 and comparable mass ratio m p ∼ m c yields temperatures on the order of 0.1 ω D /k B . These are smaller than the melting temperature of the crystal, T m ≃ 0.1 √ 2r d ω D /k B (in 1D), and are thus reasonable, even for a "soft" crystal with r d ∼ 1.

Extended Hubbard model for atomic polarons inside a dipolar crystal
As we have seen in Section 2.1, it is convenient to change from a picture of bare atoms and crystal phonons, to one of atoms dressed by their surrounding crystal displacements, i.e. polarons. The corresponding displacement amplitudes u q for the dressing then are which at the band-edge attain their minimum, u π/a ≈ 1.02U 2 0 /N 1/2 , while they diverge at small quasimomenta as ∼ 4.94|qa| −1/2 . In (44) we introduced the coupling ratio Thus we switch between the weak (S T ≪ 1) and strong (S T ≫ 1) coupling regimes by choice of the dimensionless prefactor U 0 . (b) The full particleparticle interactionṼ ij between two extra particles at interparticle distances |i − j| = 0, 1, 2, 3 for r d = 100 and mp/mc = 0.1 as a function of the coupling gcp in units of the particle recoil energy. Notice that the sign of the interaction alternates with every site and that the total interaction strength decreases with the interparticle distance as ∝ 1/|i − j| 2 .
which increases linearly with the atom-molecule coupling constant g cp /aE R,p and the mass ratio of m c /m p , but is inversely proportional to the "stiffness" of the crystal r d . The dependence of U 0 × m p /m c = g cp /aE cp r 3/4 d on the coupling g cp /aE cp and r d is also shown as dashed contour lines in figure 6(a). We see that for, e.g., a mass ratio of m p /m c ∼ 1 (e.g. for a gas of Cs atoms in a crystal of LiCs molecules ) U 0 can take values ranging from U 0 ≈ 0.01 (at g cp /aE R,p = 1/2, r d = 200) to U 0 ≈ 5 (at g cp /aE R,p = 5, r d = 1), while still having a crystal and being in the tight binding limit. For a mixture with a mass ratio of m p /m c ∼ 1/20 the coupling U 0 is significantly larger, i.e. U 0 ≈ 50 (at g cp /aE R,p = 5, r d = 1 for m p /m c ∼ 1/20).
Due to the dragging of the surrounding phonon cloud, the hopping rate for the polarons,J, compared to the (bare) hopping rate, J, is suppressed by the factor J/J = e −ST , where the exponent is, cf. (11), We notice that S T increases quadratically with U 0 (cf. the coupling g cp ), whereas the ratio S T /U 2 0 depends only on the temperature T (in units of the Debye frequency), which is shown in figure 5(a). In particular we remark that for T = 0 the exponent scales with the coupling ratio U 0 as S T =0 ≈ 0.9U  The phonon-coupling provides phonon mediated particle particle interactions of strengthṼ for two polarons at sites i and j, respectively. We notice thatṼ (1) ij are temperatureindependent and vary in sign and magnitude with the separation i − j, i.e. they are attractive for even i − j and repulsive for odd i − j, while their absolute value decreases with increasing inter-polaron separation |i − j|. In figure 5(b) we show the leading contributions for the total off-site shifts, which are purely induced by the phonons, for i = j, as a function of g cp /aE R,p , cf. the leading off-site terms decay as V i,j /2E p ≈ 0.16/|i − j| 2 . In addition we also plot the corresponding (modified) hopping rateJ for zero temperature T = 0. We notice that near g cp /aE R,p ∼ 1.6 the nearest neighbor interactions become comparable with the effective hopping rate, For bosonic particles, the phonon-mediated interactions also include an attractve on-site shift, the value of which is exactly twice the polaron shift, This can (in principle) lead to a collapse of the system, as it favors the piling up of polarons at a single site. However, since the full interactions comprise the bare and the phonon-mediated interaction, the total onsite shift,Ṽ ii = V ii − 2E p , is positive for V ii /2 > E p , which we require to ensure the stability of the bosonic system. We remark that by resorting to tune g pp via a Feshbach resonance, one can tune the onsite-shift up to V ii ∼ ∆, without breaking the single-band approximation in the Hubbard model (1). Thus we notice that in principle the stability of the system can be guaranteed, provided E p < (V ii /2) < ∆/2.

Corrections to the extended Hubbard model
In the following we are interested in higher-order corrections to the effective Hubbard model H S of (3), which we derived in Section 2 in terms of the spectral densities J kl ij (ω) for correlated nearest-neighbor hopping events i → j and k → l. For our atomic-crystalline mixture we find from (20) that the latter are given by where we took β q ≈ 1 and ω q ≈ ω D sin(qa/2) and q ω ≡ arcsin(ω/ω D )/a. The spectral density for the "swapping" of two particles on neighboring sites, J 10 01 (ω) is shown in Figure 4(c) and shows a Van Hove singularity at ω → ω D , due to the 1/ √ ω D − ω divergence of the density of states for the crystal phonons.
Strong coupling limit S T ≫ 1: The value S T is determined from equation (45). Values S T ≫ 1 are obtained for large coupling ratios U 0 and/or high temperatures T . However, already for a (reasonably small) ratio U 0 > 1.1 we have S T > 1 at T = 0 and thus we are in the strong-coupling regime for all temperatures T ≥ 0. In this limit the main corrections to the Hubbard model (3) are due to "swap" processes. The corresponding rates and coefficients for S T ≫ 1 are well approximated by the expressions (22b)-(22d), which are shown in figure 7 and figure 8 as a function of the coupling ratio U 0 and the temperature T . Notice that in (22a)-(22d) the parameter B = 4π 6 U 2 0 /3(186ζ(5)) 3/2 ≈ 0.48U 2 0 while the ratio A T /B exceeds (k B T / ω D ) tanh( ω D /2k B T ).
In figure 7(a-d) we show the leading contributions in the strong-coupling regime as a function of the ratio U 0 and the temperature T . In particular, panels (a), (b), (c), and (d) are contour plots of the incoherent rate Γ 10 01 /J, the coherent shift ∆ 10 01 /J, γ and δ, respectively, as obtained from the strong coupling approximation (22a)-(22d). The two dashed lines in each panel signal where S T = 10 and S T = 1, and the strongcoupling approximation is valid (S T > 1). We remark that in panels (a,b) Γ/J, ∆/J are divided by the (small) ratio J/ ω D , while in (c,d) γ and δ are divided by the even smaller ratio (J/ ω D ) 2 . We notice that in the spirit of the master equation approach, all quantities in the figures are plotted at finite temperature. The figure shows that, in a wide range of parameters, corrections to the coherent time evolution determined by H S can be made small in the strong coupling regime.  while U 0 1/2 for T = ω D ). We notice that in the area where S T 0.1 the ratios shown in panels (a-d) are smaller than ≈ 1, and thus all corrections are strongly suppressed compared toJ, provided J ≪ ω.

Polar molecules interacting with a one-dimensional crystal
As a second configuration, we consider a setup where polar molecules of a second species are trapped at a distance b from the crystal tube, under one-dimensional trapping conditions [see figure 1(c)]. An external electric field aligns all dipoles in the direction perpendicular to the plane containing the two tubes. Molecules trapped in the two different tubes interact via long-range dipole-dipole interactions.

Tight binding limit and Hubbard models
For crystal molecules fixed at the equilibrium positions with lattice spacing a, the particles (that is, the molecules of the second species) feel the following periodic potential where d p is the induced dipole moment of the second-species molecules. The potential above has a depth which determines the band-structure for the particles, with and E R,p = 2 π 2 /2m p a 2 the particle recoil energy. The lattice depth V 0 is shown in Fig. figure 9(a) to have a comb-like structure for b/a < 1/4, since the particles resolve the individual molecules forming the crystal, while it is sinusoidal for b/a 1/4. Figure 9(b) shows the width 4J of the lowest-energy band, with J/E R,p ∼ (V 0 /E R,p ) 3/4 e −2 √ V0/ER,p for b/a 1/4, together with the energy gap ∆ ≃ (4V 0 E R,p ) 1/2 , as a function of b/a and forv 0 = 1 and 50. For a single particle, the single-band model is valid for 4J < ∆. Figure 9(c) is a contour plot of the regimes of validity of the single-band model as a function of b/a andv 0 .
When more particles are considered, the strong dipole-dipole repulsion between the particles acts as an effective hard-core constraint [37]. We find that for 4J < ∆ and d p ≪ d c the bare off-site interactions satisfy V ij ∼ d 2 p /(a|i − j|) 3 < ∆, and thus the single-band model is still valid. give an idea about the extra-particle crystal interaction potential which determines the bandstructure. The hopping amplitude 4J and the gap to the first excited band ∆ are shown in panel (b) as a function of b/a in units of the particle recoil energy. The bandstructure strongly depends on the ratio of the particle-crystal interaction energy over the kinetic particle energyv 0 = (dp/dc)r d /(mp/mc). The hopping amplitude decreases while the gap increases rapidly with increasing coupling strength andv 0 . The single band model is only valid where the gap exceeds the bandwidth ∆ > 4J.
In this configuration, the particle-phonon coupling as obtained from equation (A.13) is given by where K 1 denotes the modified Bessel function of the second kind, and β q = dxe iqx |w 0 (x)| 2 , with w 0 (x) the lowest-band Wannier functions. Figure 10(a) shows that for b/a small-enough, such that the single-band approximation is fulfilled for all v 0 [see figure 9(c) above], M q becomes peaked at large quasimomenta q. We notice that consistency with the requirement of a stable crystal implies that the variance of the fluctuations of the crystal molecules around their equilibrium positions induced by the presence of a particle localized at a site j, δv ij , be small compared to the interparticle distance a, that is δv ij /a < 1. For a given ratio d p /d c , this limits how small the ratio b/a can realistically be, in order to avoid that the inter-species interactions destroy the crystalline structure [57]. For example, for a ratio d p /d c ≈ 0.1 the ratio b/a can be as small as b/a ≈ 0.2.

Extended Hubbard model for molecular polarons inside a dipolar crystal
We continue by determining the modified Hubbard parametersJ andṼ ij for this configuration. Here, the parameter S T , which determines the regime of interactions, is given by The latter depends strongly on the ratio b/a and is proportional to r For a given r d and d p /d c ratio, the regimes of weak and strong coupling, S T ≪ 1 and S T ≫ 1, respectively, can be directly determined from figure 11, which is a contour plot of S T as a function of the dimensionless temperature k B T / ω D and the ratio b/a. As in the previous model S T increases with increasing temperature and particlephonon coupling (that is, with decreasing ratio b/a).
The phonon mediated interaction as determined from equation (8) is given bỹ As in the previous model, the phonon mediated interactions show oscillations which for b/a 1/4 decay slowly as 1/|i − j| 2 and are thus long-ranged. Depending on their sign, they can enhance or reduce the bare dipole-dipole repulsions between two second-species molecules. The phonon-mediated interaction is strong and dominates the full particle-particle interactionṼ ij for small b/a. This is shown in figure 11(b) which is a plot ofṼ ij /V ij as a function of b/a, where for b/a 0.4 the value ofṼ ij /V ij can even change sign. With increasing intertube distance, the phonon-mediated term becomes small compared to the bare value V ij .

Corrections to the extended Hubbard model
In the following we are interested in coherent and incoherent corrections to the time evolution determined by the effective Hubbard Hamiltonian H S of (3). These can be calculated in terms of the spectral density, which has the form Figure 11. The quantity S T depends on the temperature k B T , the ratio b/a and is proportional to r 1/2 d (dp/dc) 2 . In figure (a) we show how S T behaves with the dimensionless temperature k B T / ω D and the ratio b/a. By fixing the dipole ratio and r d we can determine where the weak and strong coupling regimes are valid. Figure (b) shows the full particle-particle interactionṼ i,j , the sum of the bare dipole-dipole repulsion with the phonon-mediated interaction, in units of the bare nearest neighbour interaction V i,i+1 . The full interaction decays with increasing distance and shows an alternating sign for small ratios b/a where the phononmediated interaction V where we have used ω q ≈ ω D sin(qa/2), and q w a = 2 arcsin(w/ω D ). The spectral density tends to zero linearly for small w and shows an integrable (van Hove) singularity at w = ω D . The spectral density J 10 01 (w) is plotted in figure 10 for a few values of b/a.
General expressions for the corrections ∆ q (T ), Γ q (T ), γ q (T ) and δ q (T ) are given in (17a)-(17d). For the configuration that we consider here, they depend on b/a, the temperature k B T and the dimensionless parameter (d p /d c ) 2 r As explained in Section 2, these quantities correspond to the corrections for the "swap" process (i = l and k = j), which is not suppressed exponentially by a factor ∝ exp(−2S T ), and is therefore the dominant correction in the strong-coupling limit S T ≫ 1 [57; 61]. In particular, we can estimate ∆ q (T ) ≈ 2∆ 10 01 (T ) and Γ q (T ) ≈ 2Γ 10 01 (T ), while upper bounds for δ q (T ) and γ q (T ) can be estimated as δ max (T ) < 12δ 10 01 (T ) and γ max (T ) < 12γ 10 01 (T ).  2)] can exceed the effective tunneling rateJ when the latter is strongly suppressed for strong couplings S T ≫ 1. While large Γ 10 01 can in principle lead to significant decoherence, and thus a transition from coherent hopping to thermally-activated hopping for large enough temperatures, we find that for reasonable temperatures k B T / ω D 0.1 these processes are strongly suppressed. On the other hand, the self-energies ∆ 10 01 (T ) can lead to significant modifications to the coherent-time evolution determined by H S by providing next-nearest neighbor hopping and sizeable off-site interactions in the strong coupling regime S T ≫ 1.
Weak coupling limit S T ≪ 1: Figure 13  2)], in the regime of parameters where the single-band approximation is valid. We find that the weak coupling limit covers most of the accessible parameter regime for b/a. Figure 13 shows the maximal eigenvalues Γ max , ∆ max , γ max and δ max for all four corrections. The central result here is that we find that the latter are negligible compared to the coherent hoppingJ in the region where the weak coupling expansion is valid. In the figure, as a reference to identify how good the "weak-coupling" expansion is we also plot the values of S T for the various regimes of parameters.
Finally, we conclude this section by providing an example of the regime of validity of our model configuration. We consider a crystal of SrO where second-species molecules are KRb. The dipole and mass ratios are d p /d c ≈ 0.08 and m p /m c ≈ 1.2, respectively. We find that our treatment properly accounts for the system dynamics for separations 0.2 b/a 0.7, provided r d 80.

Conclusion
In this work we studied the realization of lattice models in mixtures of cold atoms and polar molecules, where a first molecular species is in a crystalline configuration and provides a periodic trapping potential for the second (atomic or molecular) species. We have treated the system dynamics in a master equation formalism in the Brownian motion limit for slow, massive, particles embedded in the molecular crystal with fast phonons. In a wide regime of parameters the reduced system dynamics corresponds to coherent evolution for particles dressed by lattice phonons, which is well described by extended Hubbard models. For two realistic one-dimensional setups with atoms and molecules these lattice models display phonon-mediated interactions which are strong and long-ranged (decaying as 1/|i − j| 2 ). The sign of interactions can vary with distance from repulsive to attractive, which can possibly lead to the realization of interesting phases of interacting polarons in one dimensional dipolar crystals. This study, and extensions to two dimensions will be the subject of future work.

Acknowledgments
We acknowledge funding from the European Union through the STREP FP7-ICT-2007-C project NAME-QUAM (Nanodesigning of Atomic and MolEcular QUAntum Matter) and the Austrian Science Foundation (FWF).

Appendix A. Hamiltonian for particles moving in a crystal
It is the aim of this section to derive the Hamiltonian (1) starting from a generic mixture of two interacting species of atoms or molecules.

Appendix A.1. Effective continuum Hamiltonian
A mixture of two interacting species is confined to one or two dimensions by a strong optical trapping potential. In [52] it is shown how such a trapped system can be reduced to an effective lower dimensional model in the low energy limit. The effective one-or two-dimensional Hamiltonian reads Here H c describes the effective motion of the crystal particles, H p the dynamics of the extra particles with an extra particle-crystal particle interaction given by H cp . From Hamiltonian (A.1) we identify where we denote p i (r i ) and P i (R i ) the momentum (position) of extra particles and crystal particles with masses m p and m c , respectively. The sums range over all the respective particles in the mixture. The interaction between two particles from the species α, β ∈ {c, p} in a distance r from each other is denoted by V αβ (r). We take the continuum Hamiltonian (A.1) as the starting point for the following discussion and derive our lattice model from it.

Appendix A.2. Crystal Hamiltonian
In the crystalline phase the particles are characterized by small fluctuations of their positions R j around their equilibrium positions R 0 j . Thus we expand the potential in a Taylor series to second order in the displacements u j = R j −R 0 j about the equilibrium positions [58] as which is readily diagonalized and provides the dispersion relation ω q,λ for the phonons with quasimomentum q and polarization e λ . We express the displacement u i in terms of the bosonic creation (annihilation) operators a † q,λ (a q,λ ) for the corresponding phonons as where N denotes the total number of crystal particles. We write the Hamiltonian as and can identify m c ω 2 q,λ from the eigenvalues of D ij from which we are able to read off the phononic dispersion relation, see Appendix B.

Appendix A.3. Extra particles inside the crystal
The dynamics of extra particles inside the crystal and their interaction with the crystal is described by H p + H cp .
In analogy to the proceedings of the last section we are interested in the dynamics of a stiff crystal and therefore we expand the particle-crystal potential in the small displacements u i of crystal constituents about their equilibrium positions, as Here we have introduced the potential V p (r) that an extra particle at position r feels from the entire crystal. The lowest order provides a static periodic potential We note that in writing equation (A.9) we only retain the first non-vanishing correction to the static trapping potential, i.e. the term linear in the displacement u i . Thereby we neglect terms of second (and higher) order in u j , which are expected to (merely) provide a renormalization of the phonon-spectrum, i.e. when including those terms in equation (A.6).
The extra particles that we consider are confined to a plane/tube parallel to a crystal plane/tube, as pictured in figure 1. Hamiltonian (A.3) together with the zeroth order contribution to the particle-crystal interaction, describes interacting particles for which the entire crystal provides a static periodic trapping potential, the ingredients for a simple Hubbard model. The particle crystal interaction then determines the band structure.
We consider a single band model where extra particles cannot be excited to the second band that is separated from the lowest band by an energy gap ∆. Then in the low energy limit the extra particles localize at single sites of position r 0 i and their wave functions become Wannier functions of the lowest band w 0 (r − r i ). In second quantization the field operator can then be expanded as where c † i and c i denote the creation and annihilation operators of extra particles at site i, which obey the canonical bosonic (fermionic) commutation (anticommutation) relations for bosons (fermions).
For such a model, see [76], the dynamics of the extra particles is described by the hopping amplitude J that calculates from the overlap of the wavefunctions at two neighbouring sites and the interaction between two particles located at sites i and j is given by Hamiltonian (A.10) can then be written as The single band approximation is valid if all particle energies are (much) smaller than the gap, cf. J, V ij ≪ ∆. We remark, that the Debye frequency ω D in our models is typically (much)larger than the gap while the particle phonon coupling is dominated at high frequencies, see discussion below. In order to avoid excitations beyond the gap by the coupling we put a constraint on the temperature in a way that all phonon modes with energies larger than ∆, are essentially unoccupied.
Let us now focus on the higher order terms of the interaction in equation (A.9) which correspond to the backaction of the crystal on the extra particles. The remaining part of the Hamiltonian is given by and describes a dynamic coupling of the particles to the vibrations of the crystal. In second quantization for the extra particles the remaining Hamiltonian is obtained (cf. Appendix B.2) as were we have introduced the particle-phonon coupling M q,λ given by Here β q denotes the Fourier transform of the modulus square of the Wannier function The function β q accounts for the localization of the particles with a finite width in the static trapping potential and thus, for small q, approaches 1. Similarly, for the class of potentials we consider (and discuss) in Section 4, the Fourier transform of the particlecrystal potentials approaches a finite value at q = 0. Therefore the particle-phonon coupling behaves like q 1/2 for small momenta, while it is large for q ∼ π/a. The Hamiltonian of the entire model is given by H eff = H c + H ′ p + H ′ cp and reads 14) The crystal motion given by Hamiltonian (A.14) corresponds to a set of uncoupled harmonic oscillators under the influence of an extra particle density dependent force.
Similarly to the problem of a charged harmonic oscillator in a constant electric field [58] such a force displaces the crystal molecules from their original position proportional to its strength while keeping the oscillation frequency fixed. A crystal molecule at position R i is therefore displaced by with the Fourier transform of extra particle density η q = j e −iqr 0 j c † j c j . If the relative displacement between two crystal molecules δv ij = v i − v j at neighbouring sites i, j is small on the scale of the lattice constant δv ij ≪ a this effect can be neglected.

Appendix B. Phonon spectrum and particle-phonon coupling
In this section we derive the specific form of the dispersion relation for one-and two-dimensional dipolar crystals as well as the crystal-particle interaction.

Appendix B.1. The phonon spectrum
The crystalline phase is characterized by small displacements of the molecules from their equilibrium positions R i = R 0 i + u i . A dipolar crystal that is trapped in one or two dimensions by an optical trap with a trapping frequency ω ⊥ is described by the following Hamiltonian where P i denote the momenta of the molecules. The dispersion relation is found from the second order correction of the expansion of the molecule-molecule interaction potential We make an ansatz for the displacement as and can identify a † q,λ and a q,λ as creation and annihilation operators of phonons with the polarization e λ and quasi momentum q iff the matrix D ij is diagonal. Then the polarization vectors for phonons with a dispersion ω q,λ are the eigenvectors of D ij .
With the discrete Fourier transform we can write Hamiltonian (B.1) as ω q,λ a † q,λ a q,λ + 1 2 and find the dispersion relation from calculating the eigenvalues of D ij . The +1/2 contribution in the last equation is the vacuum zero point energy, a constant energy shift that is ommitted henceforth.
Appendix B.1.1. Phonon spectrum in 1D: In a 1D crystal tube the second order of the expansion in the displacement of the interaction potential is given by In q-space we find

Phonon spectrum in 2D:
We consider a dipolar crystal that forms in the xy-plane. The expansion of the molecule molecule potential to second order in the displacement gives can write It is difficult to diagonalize the matrix in the last equation. We can however use the block diagonal form of the matrix in the last equation and, including the optical trapping potential as we have done in the 1D case above, write the dispersion relations as where f ± q denote the two eigenvalues of the xy-block of the matrix and f q is given by equation (B.8). The eigenvalues f ± q give a longitudinal and a transversal acoustic dispersion relation. As noted above we cannot write them in a closed form but an approximate result may be obtained by including only particles in the interaction that are within some finite range of each other. We show a numeric evaluation of the dispersion relation in figure 3.

Appendix B.2. The particle-phonon coupling
We denote the crystal-particle interaction potential by V cp (r − R j ) where r and R j denote the position of an extra particle and a crystal particle respectively. The full potential is expanded up to first order in the displacement R j = R 0 j + u j as where d denotes the dimension of the setup and we have used the discrete Fourier transform of the interaction potential The Hamiltonian is found by integration of the extra particle density ρ(r) over the interaction H I = drρ(r)V cp (r) which gives The extra particle density ρ(r) is defined through the single particle operator, which reads in site representation ψ(r) = m c m φ m (r) with φ m (r) the wavefunction of an extra particle at site m, and the annihilation operator c m as ρ(r) = ψ † (r)ψ(r). In the tight binding limit the particles are strongly localized at sites and the overlap of the wavefunctions of particles at different sites is neglected. This approximation gives We can work out β q = dr|φ(r)| 2 e iqr by making a separation ansatz for the wavefunction φ(r) = φ x (x)φ y (y)φ z (z) and thus β q = β qx β qy β qz . In the directions of confinement the wavefunction is taken to be Gaussian φ z (z) = e −z 2 /2a 2 is taken a Wannier function in all other directions thus β q = dr|w 0 (r)| 2 e iqr with w 0 (r) the Wannier function of the lowest Bloch band. The interaction Hamiltonian is where we can use (A.7) to identify Extra particles interacting with crystal molecules will force the latter to new equilibrium positions displaced from their original positions by an extra particle density dependent displacement The impact of the displacement is neglected if the relative shift between two neighbouring molecules |v i − v i+δ | = v i − v i+δ is small compared to the lattice constant.
we find where we have introduced Z q,λ kl = u q,λ e −iqr 0 l −e −iqr 0 k and the displacement operator D(a, α) ≡ e αa † −α * a . In Section 2.2 we have introduced the (thermal) expectation value of a bath operator O by O ≡ tr B {Oρ 0 B } where tr B denotes the trace over the bath degrees of freedom. The expectation value of the displacement operator is given by D(a q,λ , α) = exp[−(n q,λ (T ) + 1 2 )|α| 2 ] wheren q,λ (T ) = 1/(exp[ ω q,λ /k B T ] − 1) denotes the thermal expectation value of the phonon number operator. Therefore we find Notice that the link r 0 i − r 0 j is always a symmetry axis of the integration interval, the first Brillouin zone. Since u q,λ andn q,λ (T ) are invariant under a rotation with the symmetry of the Brillouin zone the sum over q in equation (C.4) becomes independent of the orientation of that link.
Expectation values of time dependent displacements can we worked out in a similar fashion.
With the explicit form of the bath Hamiltonian H B = q,λ ω q,λ a † q,λ a q,λ we can write a bath operator in the interaction picture asã q,λ (t) = a q,λ e iω q,λ t and we find D(a q,λ , Z q,λ kl e −iω q,λ t ). (C.5) In the following we introduceZ q,λ ij (t) = Z q,λ ij e iω q,λ t and use the identity for displacement operators q,p D(a q , α)D(a p , β) = q D(β, α/2)D(a q , α + β) to calculate Here Φ kl ij (τ, T ) is found by inserting forZ(t) as Φ kl ij (τ, T ) = q,λ u 2 q,λ n q,λ + 1 g kl ij e −iω q,λ τ +n q,λ g kl ij * e iω q,λ τ (C.7) with g kl ij = (e −iqr 0 j − e −iqr 0 i )(e iqr 0 l − e iqr 0 k ). The quantity Φ kl ij (τ, T ) depends only on the relative time τ = t − t ′ and the two links r 0 i − r 0 j and r 0 k − r 0 l . A symmetry argument that relies on the fact that (qr 0 i ) is a projection on a symmetry axis of the integration interval, the first Brillouin zone, while everything else is an even function in q allows us to write (C.7) as From this form one can immediately read off the important relations Φ kl ij (τ, ij * (τ, T ) and |Φ kl ij (τ, T )| ≤ 2S T . Equation (C.8) is in accordance with the literature [58].
For many sites we take the continuum limit and replace the summation over q by an integration over the first Brillouin zone of volume V BZ . A variable transformation where q x is replaced by f (q d−1 , w), a function of the remaining momenta q d−1 and w ≡ ω q,λ , allows us to introduce the spectral density where the q d−1 integration interval depends strongly on w and the form of the Brillouin zone while w ranges over all frequencies. In terms of the spectral density equation (C.8) reads Φ kl ij (τ, T ) = dwJ kl ij (w) coth w 2k B T cos(wτ ) − i sin(wτ ) . (C.10) The bath correlation functions ξ kl ij (τ, T ) that enter the master equation are thus given by ξ kl ij (τ, T ) = X † i (t)X j (t)X † k (t ′ )X l (t ′ ) − e −2ST (C.11) = e −2ST (e −Φ kl ij (τ,T ) − 1) (C. 12) and satisfy ξ ij kl (−τ, T ) = ξ kl ij * (τ, T ).

Appendix C.2. Corrections to the Master Equation
The We can find explicit approximate results to these corrections in the two limiting cases of a weak, S T ≪ 1, and a strong, S T ≫ 1, coupling. Figure C1. Real part of (a) the function Φ kl ij (τ, T ) and of (b) the correlation function ξ kl ij (τ, T ) for a 1D model with Mq ∝ q/ √ ωq and ωq = ω D sin(q/2) as a function of the (dimensionless) time τ ω D . Shown are the respective functions for the swapping of a particle, i = l and k = j (solid lines), and the hopping of particles into the same direction, j = k (dashed lines). Notice the overall minimum of the functions |Φ kl ij (τ, T )| is attained at τ = 0 for the swap process, i.e. i = l and k = j. In panel (b) we chose S T = 3/4, while for larger values of S T the tail of Re[ξ 10 01 (τ, T )], together with all other processes, is strongly suppressed. In the strong coupling limit (S T ≫ 1) the latter become of the order of O(e −2S T ). (C.25) The quantity B for this process reads B = dwJ 10 01 (w)w = 2 π dq M q ω q 2 sin(q/2) 3 (C.26) 2π ω D .
(C. 27) In the two models in Section 4 we find that the onsite phonon mediated interactioñ V (1) i,i = 2E p dominates over all the others. Therefore we can approximately write B ≈ E p /π ω D with which we find where the w integration ranges over all frequencies. This way the Debye frequency of our model functions as a natural cutoff (compare [68]). The functional identities where we have written out the Cauchy principal value integrals. We remark that the integration (C.35) becomes infinite in the limit of zero temperature.

Appendix C.3. Diagonalization of the master equation
As noted above the corrections to the master equation Γ kl ij (T ), ∆ kl ij (T ), γ kl ij (T ) and δ kl ij (T ) are just matrix entries and have little physical meaning by themselves.
In this section we show how to estimate the energy of the corrections to the master equation (16) for a single particle at many sites. We start by treating the self energies, The dissipative term proportional to γ kl ij (T ) in the single particle limit is given by where k ′ and l ′ denote the nearest neighbours of k and l. We can write the first two terms as ij kl In the second step we have used the fact that the origin of the indices i, j, k, l is the function g kl ij (q) and therefore the corrections show the same symmetries with respect to the indices as g kl ij (q) itself. In the same way we find the other eigenvalues as