3 gives the interaction picture Hamiltonian where we define the laser detuning Δ ≡ ω$_{0}$ − ω$_{L}$ (i.e. red detuning is positive) and the Rabi frequency We choose the relative phase of | 0 〉 , | 1 〉 so as to make Ω real. However, unlike the usual Bose-Hubbard model [10], we do not assume that the atoms are in the lowest motional band of the lattice, but allow them to be in higher mo- tional states. The excitations to these higher motional states are the phonons, with creation and annihilation operators B † $_{iα}$, B$_{iα}$ for a motional excitation of atom i in the direction α . For simplicity we consider the case of a single plane po- larized cooling laser, with peak electric field E$_{L}$ , wavevec- tor k$_{L}$ , and frequency ω$_{L}$ . To describe the spontaneous decay and dipole-dipole interaction, we also need to in- clude the free radiation field, with photon annihilation operators a$_{k}$$_{λ}$ labelled by a wavevector k and a polariza- tion index λ = 1 , 2. We treat the atom-field interactions in the electric dipole approximation, i.e. as the dot product of the atomic dipole moment and the electric field. The classical laser field is while the free photons have electric field operator where L 3 is the quantization volume and ϵ$_{k}$$_{λ}$ the polar- ization vector. Putting all these together, the Hamiltonian of the sys- tem in the Schr¨ odinger picture is As our system is not in an optical cavity, the free pho- tons a$_{k}$$_{λ}$ quickly leave the system, and are absorbed by the walls of the apparatus on a time scale Δ t ∼ L/c , mak- ing our system an open quantum system. We hence need to use a master equation rather than a Hamiltonian to describe its evolution. For realistic sizes, ω$_{0}$ ≫ 1 / Δ t ≫ ν, Δ , Ω , Γ (where Γ is the decay rate derived below), which allows us to use the method of Refs. [16, 17, 25] to derive a Markovian (memoryless) master equation. This method approximates the photon loss process by alternating free evolution under H$_{I}$ for time Δ t with the instantaneous loss of all photons, where ρ ( t ) is the (electronic and motional) density matrix of the atoms at time t , | 0 〉$_{pht}$ is the vacuum state of the photon field, and U$_{I}$ ( t$_{0}$ + Δ t, t$_{0}$ ) is the unitary evolution operator under H$_{I}$ from time t$_{0}$ to t$_{0}$ + Δ t . Since Δ t is short compared to the dynamical timescales of H$_{I}$ , we can approximate this by expand- ing the time ordered exponential U$_{I}$ ( t$_{0}$ + Δ t, t$_{0}$ ) = where x$_{i}$ is the position of atom i . We next eliminate the time dependence of the laser term, by moving to the interaction picture with respect to and applying the rotating wave approximation with re- spect to ω$_{0}$ (i.e. neglecting terms e i ωt with ω ≳ ω$_{0}$ ). This T exp$^{{}$ − i ℏ t$_{0}$ +Δ t ∫ t$_{0}$ dtH$_{I}$ ( t )$^{}}$ up to second order in Δ t : Δ ρ Δ t = Tr$_{photons}$$^{}$ $_{}$− i ℏ    1 Δ t t$_{0}$ +Δ t ∫ t$_{0}$ dtH$_{I}$ ( t ) (10) − i $_{ℏ}$Δ t t$_{0}$ +Δ t ∫ t$_{0}$ dt t ∫ t$_{0}$ dt $^{′}$H$_{I}$ ( t ) H$_{I}$ ( t $^{′}$)$^{}$   ρ ( t$_{0}$ ) ⊗ | 0 〉 〈 0 |$_{pht}$ + H . c . + 1 ℏ $^{2}$Δ t t$_{0}$ +Δ t ∫ t$_{0}$ dt t$_{0}$ +Δ t ∫ t$_{0}$ dt $^{′}$H$_{I}$ ( t ) ρ ( t$_{0}$ ) ⊗ | 0 〉 〈 0 |$_{pht}$ H$_{I}$ ( t $^{′}$)$^{}$  , H$_{0}$ = N ∑ i =1 $_{ℏ}$ω$_{L}$ S † $_{i}$S$_{i}$ + ∑ k λ $_{ℏ}$ω$_{k}$ a † k $_{λ}$a$_{k}$$_{λ}$ , (6) B. The interaction picture H = N ∑ i =1 $_{ℏ}$ω$_{0}$ S † $_{i}$S$_{i}$ + N ∑ i =1 e D$_{i}$ · $^{[}$E$_{L}$ ( x$_{i}$ , t ) + E ( x$_{i}$ ) ] + ∑ k λ $_{ℏ}$ω$_{k}$ a † k $_{λ}$a$_{k}$$_{λ}$ + N ∑ i =1 ∑ α = x,y,z $_{ℏ}$ν$_{α}$ B † $_{iα}$B$_{iα}$ , (5) ρ ( t$_{0}$ + Δ t ) = Tr$_{photons}$ [ U$_{I}$ ( t$_{0}$ + Δ t, t$_{0}$ ) × $^{{}$ρ ( t$_{0}$ ) ⊗ | 0 〉 〈 0 |$_{pht}$ $^{}}$U † I ( t$_{0}$ + Δ t, t$_{0}$ ) $^{]}$, (9) E ( x ) = ∑ k λ i$^{√}$ ω$_{k}$ 2 ϵ$_{0}$$_{ℏ}$L 3 ϵ$_{k}$$_{λ}$ e i k · $^{x}$a$_{k}$$_{λ}$ + H . c . , (4) E$_{L}$ ( x , t ) = E$_{L}$ cos( k$_{L}$ · x − ω$_{L}$t ) , (3) D$_{i}$ = D$_{01}$ S$_{i}$ + H . c . (2) C. The master equation Ω ≡ 2 e ℏ D$_{01}$ · E ∗ L . (8) H$_{I}$ = N ∑ i =1 $_{ℏ}$Δ S † $_{i}$S$_{i}$ + N ∑ i =1 $^{[}$1 2 $_{ℏ}$Ω e − i k$_{L}$ · x$_{i}$ S$_{i}$ − ∑ k λ i e$^{√}$ ω$_{k}$ 2 ϵ$_{0}$$_{ℏ}$L $_{3}$D$_{01}$ · ϵ ∗ k $_{λ}$e − i k · x$_{i}$ e i( ω$_{k}$ − ω$_{L}$ ) $^{t}$S$_{i}$ a † k λ ] + H . c . + N ∑ i =1 ∑ α = x,y,z $_{ℏ}$ν$_{α}$ B † $_{iα}$B$_{iα}$ , (7)