Dissipative dynamics of atomic Hubbard models coupled to a phonon bath: dark state cooling of atoms within a Bloch band of an optical lattice

We analyse a laser assisted sympathetic cooling scheme for atoms within the lowest Bloch band of an optical lattice. This scheme borrows ideas from sub-recoil laser cooling, implementing them in a new context in which the atoms in the lattice are coupled to a Bose–Einstein condensate (BEC) reservoir. In this scheme, excitation of atoms between Bloch bands replaces the internal structure of atoms in normal laser cooling, and spontaneous emission of photons is replaced by creation of excitations in the BEC reservoir. We analyse the cooling process for many bosons and fermions, and obtain possible temperatures corresponding to a small fraction of the Bloch band width within our model. This system can be seen as a novel realisation of a many-body open quantum system.


Introduction
New frontiers in atomic physics have often been enabled by the development of new cooling techniques. Examples are provided by laser cooling and evaporative cooling [1]- [4], which underly the exciting experimental advances to realize Bose-Einstein condensates (BEC) and quantum degenerate Fermi gases of atoms and molecules [5]- [15]. 4 Recently, quantum degenerate gases of bosons and fermions have been loaded into optical lattices in one, two-and three-dimensional (2D, 3D) configurations [16]- [25]. This makes it possible to realize atomic Hubbard dynamics with controllable parameters, opening the door to the study (and simulation) of strongly correlated systems with cold atoms. Control via external fields allows the engineering of atomic lattice Hamiltonians for boson, fermion and spin models [26]- [28], which for a long time have been the focus of research in theoretical condensed matter physics. However, one of the most challenging obstacles in the realization of some of the most interesting condensed matter systems in current experiments is the need for lower temperatures of the atoms within a Bloch band of the optical lattice [29]. In this paper, we will analyse a configuration where atoms moving in the lowest Bloch band of an optical lattice are cooled via laser assisted sympathetic cooling with a heatbath represented by a BEC of atoms. A unique feature of the present scheme is that the achievable temperatures of the atoms within a single Bloch band are significantly lower than those of the cooling reservoir, reaching a temperature of a small fraction of the Bloch band width. From a physics point of view, a guiding idea of the present study is formal analogies with laser cooling [1]- [4]. In laser cooling the motional degrees of an atom are coupled via laser excitation of the electrons to the effective zero temperature (photon) reservoir, i.e., the vacuum modes of the optical light field. The spontaneous emission of photons carries away the entropy allowing a purification (cooling) of the motional state of the atom. In the present study, we will follow a similar scenario by coupling atoms moving in the lowest Bloch band via laser assisted processes to an excited band, which can decay back to the lowest band by spontaneous emission of phonon (or particle-like) excitations in the BEC. We note that this cooling process is atom number conserving, in contrast to filtering and evaporative cooling technquies [7,30,31]. The present paper focuses on a specific protocol for cooling within a Bloch band which follows analogies with (subrecoil) Raman laser cooling, as developed in a seminal study by Kasevich and Chu [32].
We see an important aspect of the present study in pointing out the formal analogies between open quantum systems familiar from quantum optics with light fields and the present set-up of atoms in optical lattices coupled to a phonon bath. These analogies not only stimulate the transfer of well established ideas of, e.g., laser cooling, to a new context, but also provide a new realization of an open quantum optical system with dissipative dynamics of cold atoms in optical lattices. The present paper extends and gives details of the study published in [33], in which we assumed a set-up with a homogeneous optical lattice, without an additional trapping potential for lattice atoms. As an additional remark, we also present an extension of these ideas to a form of spatial sideband cooling in a harmonic trapping potential.
We consider a set-up as shown in figure 1, where the internal electronic states of atoms familiar from laser cooling, are replaced with the two lowest Bloch bands of an optical lattice, trapping atoms of a species a. In analogy to laser cooling the energy of particles in the lowest Bloch band is upconverted by transferring lattice atoms to the first excited band via a Raman 4 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT transition. We immerse the atoms in the optical lattice in a BEC [34], which serves as a reservoir for the lattice atoms, and should be seen as the counterpart of the vacuum modes of the radiation field in laser cooling, carrying away motional energy from the system. In our scheme, cooling is achieved in two steps: (i) we design a sequence of excitation pulses which efficiently excites atoms with high quasi-momenta in the lattice, but do not couple to atoms with q = 0. (ii) In a second step, the coupling to the BEC reservoir leads to decay of excited lattice atoms back to the ground state via emission of a Bogoliubov excitation into the BEC, thereby randomizing the quasi-momentum. Repeating these two processes leads to an accumulation of the atoms in a narrow region around q = 0, i.e., to cooling, in analogy to the Kasevich-Chu scheme of Raman laser cooling for free atoms [32,35]. We will show below that our method works away from the limit of unit filling of the lattice, and is capable of cooling single atoms and many non-interacting bosons and fermions. The method can be utilized in a scheme to create strongly interacting phases of atoms in the optical lattice, by first cooling non-interacting atoms (thereby exploiting the tunability of interactions, e.g., via Feshbach resonances) and subsequently ramping up the interaction.
The paper is organized as follows. In the following section 2, we will present a short overview of laser cooling as is relevant for understanding the physical analogies with our cooling method outlined in later sections. Section 3 presents our model. We will then illustrate the details of the cooling protocol for the case of a single particle in section 4. In section 5, we will discuss how the system can be adapted to achieve cooling for many non-interacting bosons and fermions and we investigate the effects of interaction on the cooling scheme for many bosons. In section 6, we will show how the interaction strength between already cooled lattice atoms can be ramped up to achieve cold strongly correlated gases. In section 7, we discuss how similar ideas could be used to compact fermions in a harmonic trap, and we then conclude and summarize the ideas and main results in section 8.

Laser cooling
In this section, we will outline the main ingredients needed for the description of laser cooling and give a short overview of the development of the different laser cooling schemes. The structure presented here is directly related to the cooling scheme we present in section 3, as are the concepts of subrecoil laser cooling, which we discuss in more detail at the end of this section.
As described in the introduction, in a laser cooling scheme, an atom interacts with a laser and is coupled to the electromagnetic radiation field, i.e., one considers a system of the form Here, the atomic HamiltonianĤ a , includes the atom's motion, possibly an external trapping potential and the internal structure, usually modelled as a two or a three level system. The atom interacts with the electromagnetic radiation field, described byĤ b , and this interaction is described byĤ int . In addition, a laser set-up couples the atom's internal states, which is described byĤ LS . The two steps of laser cooling, i.e., the upconversion of the energy due to the laser and the removal of system energy due to spontaneous emission, can be described in terms of a 5 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT master equation in the Born-Markov approximation for the reduced system density operatorρ (see, e.g., [36], and taking units whereh = 1 throughout this article) where we have specialized to atomic motion in 1D and to a two level system for sake of simplicity. The atom's position operator is denotedx,σ + (σ − ) is the raising (lowering) operator corresponding to the electronic transition, denotes the linewidth of the excited state, and k l the wavenumber of the laser light. Note that photons can be spontaneously emitted in all 3D, and the angular dependence of the spontaneous emission is accounted for by the normalized dipole distribution N(u) in this 1D model. These basic ingredients have been utilized in the development of various different laser cooling schemes over the last years. This began with Doppler cooling [1], where two counterpropagating laser beams incident on an atom, so that the radiation pressure of the two beams compensates for an atom at rest, but is Doppler shifted towards resonance and thus enhanced for the counter-propagating beam in the case of a moving atom. In Doppler cooling, the final temperature is limited by the linewidth of the excited state. Lower temperatures were obtained in various schemes, including polarization gradient cooling and Sisyphus cooling [37], where the limiting temperature is given by the recoil energy an atom receives during the emission of a single photon. Cooling of atoms even below this single photon recoil limit was proposed and observed in the form of Raman laser cooling and velocity selective coherent population trapping (VSCPT) [32,35].

Subrecoil laser cooling
The basic idea in subrecoil laser cooling is to make the photon absorption rate velocity dependent and, in particular, vanishing for a dark state. The atoms will consequently accumulate in this state during the cooling process. In the following, we will describe the examples of Raman cooling and VSCPT.
2.1.1. Raman cooling. In Raman cooling [32] one considers a -system as shown in figure 2, consisting of two (hyperfine) ground states g 1 and g 2 and one electronically excited state e. The energy of the moving atom is given by the sum of its internal energy ω i , where i = g 1 , g 2 , e, and its kinetic energy q 2 /2m. The state |g 1 , q is coupled to |g 2 , q + (k 1 + k 2 ) via a Raman procesŝ H LS , which is far detuned from the excited state. Here, q denotes the momentum of the atom, k 1 and k 2 are the wavenumbers of the two Raman laser beams.
Cooling is achieved in a two step process. In the first step a set of Raman laser pulses is designed that efficiently transfer atoms with high momenta in the ground state g 1 to g 2 , but do not couple atoms with zero momentum. In the second step, which is applied after each single excitation pulse, atoms are optically pumped back to the state g 1 via the excited state e (see figure 2). The spontaneous emission process randomizes the momentum of the atom, which leads to a finite probability of falling into a region near the dark state, which in this case is the state |g 1 , q = 0 . Repeating the two steps leads to accumulation of atoms in state |g 1 with momentum near zero, i.e., to cooling of the atom. Schematic picture of one excitation and decay step in subrecoil Raman laser cooling. (a) Atoms from a region with high momentum |q| > 0 are transferred from one ground state g 1 to a second ground state g 2 . (b) The atoms in state g 2 are optically pumped back to the initial state g 1 via the electronically excited state e, thereby randomizing the momentum due to the recoil kick. Cooling is achieved by designing a sequence of excitation pulses to efficiently transfer all atoms, except for those in the dark state with q = 0, to the second internal state, each pulse is followed by the optical pumping process, and repeating the sequence leads to accumulation of atoms in the dark state, i.e., to cooling. (a) Schematic picture of the setup in VSCPT. Two circularly polarized laser beams incident on the atom, couple the two ground states |g −1 and |g 1 via the excited state |e 0 . A dark state (see text) appears due to quantum interference, where atoms accumulate during the cooling process as they decay there via spontaneous emissions from the excited state. (b) Theoretical model for the excitation rate in the Lévy statistics analysis. For momenta |q| < q 0 the excitation rate follows a power law, outside this region it is assumed constant.

VSCPT.
In VSCPT (see figure 3(a)) [35], two counter-propagating circularly polarized laser beams couple the two ground states |g −1 and |g 1 via the excited state |e 0 . In the setup states with angular momentum J = 1 are used in the ground and excited state manifold, so that only the three states |g −1 , |g 1 and |e 0 are coupled via the lasers, and spontaneous emission only leads to decay back to the states g −1 and g 1 since the transition e 0 → g 0 is forbidden. A dark state forms due to quantum interference and is given by |D = (|g −1 , k l + |g 1 , −k l )/ √ 2, where k l is the wavenumber of the two lasers. This dark state does not couple to the Raman process and will be increasingly populated during the cooling process due to the decay of atoms in the excited state via spontaneous emission. This leads to cooling of the atom, characterized by two 7 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT narrow peaks at k l and −k l in the final momentum distribution of the atom, as the dark state |D is a superposition of these two momentum eigenstates.

Lévy statistics
Analytical calculations based on Lévy statistics [38] have been shown to be a very powerful and accurate tool in the context of subrecoil laser cooling. They will be applied below to describe our optical lattice Raman cooling process. In this section, we briefly review the underlying model and some results in the context of Raman laser cooling for the most important physical parameters, such as the temperature and the fraction of the momentum distribution in the dark state.
In Lévy statistics [38] a trapping region near q = 0, with |q| q trap and a recycling region with |q| > q trap are defined. The parameter q trap is an auxiliary variable and the results for real physical quantities do not depend on it. During the cooling process, the atom will undergo a random walk in momentum space, with N trapping periods of duration τ i , alternating with N recycling periods of durationτ i , during which the atom is inside and outside the trapping region respectively. The total time of the cooling process can be written as the sum of the total trapping time T N and the total recycling timeT N where the τ i and theτ i are independent random variables. The efficiency of the cooling process can be quantified by the statistics of the total trapping and recycling times, which are determined by the excitation rate R(q) of atoms with momentum q from the state g 1 to the state g 2 . The rate R(q) is modelled as as schematically plotted in figure 3(b). The values of τ 0 , q 0 and λ depend on the details of the excitation pulses and are mainly determined by the duration, momentum transfer and shape of the last two Raman pulses. In subrecoil cooling, typically λ > 1 (especially λ = 2 for the case of time square excitation pulses in Raman laser cooling) and thus the probability distribution for the total trapping times T N is a broad distribution, i.e., the expectation values T N and T 2 N diverge. In this case, a generalized central limit theorem (for details see [38]) predicts a Lévy distribution for the probability distribution for the total trapping time. This distribution is the starting point for a calculation of the relevant physical parameters, like the width of the momentum distribution or the height of the central peak at q = 0 (for details again see [38]), with the following results. We define the temperature in terms of the half width q at e −1/2 of the maximum of the velocity distribution and find The collisional interaction with the BEC atoms is switched on, the resulting decay of the excited lattice atoms leads to a randomization of the quasi-momentum. Sequences of pulses, each one followed by a decay time τ c , efficiently excite all atoms outside a narrow region around q = 0. Repeating the sequence leads to accumulation of atoms in the dark state region around q = 0, i.e., to cooling. especially, we find T ∝ −1 for the case of time square excitation pulses. Similarly we derive an expression for the population density at q = 0 and time , n 0 ( ), which is given by where in this equation, γ(x) denotes the Euler γ function.

Raman cooling in an optical lattice
In this section, we will describe our Raman cooling scheme for atoms in an optical lattice, therefore we will again briefly describe the idea of the cooling process and then analyse the different parts of the set-up in detail. In our scheme, cooling of the atoms is achieved in a two steps (cf figures 4(a) and (b)) and can be seen in analogy to free space Raman cooling (figure 2). (i) Raman laser pulses are designed to excite atoms with large quasi-momenta |q| > 0 from the lowest band to the first excited band, whilst not exciting atoms in the dark state (with q = 0) (see figure 4(a)). (ii) The decay of excited lattice atoms goes along with the emission of a phonon (Bogoliubov excitation) into the BEC reservoir (see figure 4(b)). Assuming a BEC temperature k B T b ω the Bogoliubov modes with an energy corresponding to the separation of the Bloch bands are essentially in the vacuum state and the BEC effectively acts as a T = 0 reservoir. This is in analogy to the coupling of the atoms to the vacuum modes of the radiation field, giving rise to spontaneous emission into a T = 0 reservoir in the case of laser cooling. Repeating these two processes leads to cooling of the atoms, as they accumulate in a region near the dark state with q = 0.
Note that in contrast to the free space version of Raman cooling, in our set-up the upconversion of energy is already performed in the excitation step. This is only possible, because in our set-up the spontaneous decay back to the lowest band is due to collisional interaction with the BEC reservoir which can be switched off during the excitation step, e.g., via Feshbach 9 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT resonances.A direct analogue to the cooling protocol in free space Raman cooling could, however, also be achieved with a slightly different set-up: for example, one could use a spin-dependent lattice, and perform the excitation step from the lowest Bloch band for lattice atoms a to the lowest Bloch band of a second species a , followed by an optical pumping step via an excited Bloch band, as in free space Raman cooling.
In our set-up, we assume a homogeneous optical lattice, i.e., no additional external (harmonic) trapping potential for the lattice atoms. This assumption is an experimental requirement inherent in the realization of many strongly correlated phases, and a topic of significant current interest experimentally (see, e.g., [39]). Our analysis will be performed for the ideal case of a homogeneous system, and inhomogeneities appearing in a real experiment will limit the achievable temperatures. In practice, we expect that advances in homogeneous traps and cooling methods will occur somewhat iteratively, in that clean flat-bottomed traps allowing better cooling, which will provide the opportunity to cool atoms to even lower temperatures where the remaining imperfections become more noticeable.

Lattice atoms and laser set-up
The HamiltonianĤ a for the lattice atoms can be written asĤ a =Ĥ 0 +Ĥ I , wherê describes the kinetic energy of the atoms and the optical lattice potential V a (x) = V a,x sin 2 (πx/d) + V a,y sin 2 (πy/d) + V a,z sin 2 (πz/d). Here, d = λ/2 is the lattice spacing with λ the wavelength of the lasers generating the lattice potential, and V a,j is the strength of the optical lattice potential in j = x-, y-, z-direction, m a is the mass of the lattice atoms andψ † a (x) andψ a (x) are the field operators for the lattice atoms, which will satisfy (anti-)commutation relations in the case of (fermions) bosons. Onsite interactions between lattice atoms are represented bŷ where g aa = 4πa aa /m a with a aa the s-wave scattering length. The interaction of the two laser beams generating the Raman transition with the lattice atoms in the rotating wave approximation can be written aŝ where the Raman detuning δ = ω 1 − ω 2 is given by the frequency difference of the two lasers, is the detuning from the excited state, and 1 and 2 denote the single photon Rabi frequencies of the two lasers. Both lasers are far detuned from resonance with the transition between ground and excited state, which can thus be adiabatically eliminated (see appendix A).
For simplicity, we will assume in the following an anisotropic lattice potential, with V a,y = V a,z ≡ V a,⊥ V a,x so that the atoms effectively move in 1D along the x-direction, whereas the transverse hopping is suppressed due to the large trapping potential. Such an analysis is readily extended to higher dimensions. In the case where only the two lowest bands play a role in the dynamics of the system, we can express the field operators in terms of the Bloch functions for the lowest two bands in the x-direction, φ α q (x), where the band index α ∈ {0, 1}, and localized wavefunctions w 0 y (y) and w 0 z (z) representing the confinement in the transverse directions. This imposes requirements on the laser coupling strength and detuning, and also on the interaction strength between atoms in the lattice (see below). We obtain where (A α q ) † creates a particle with quasi-momentum q in Bloch band α. The operators (A α q ) † and (A α q ) will again satisfy (anti-)commutation relations for (fermions) bosons. Inserting this into equations (8) and (9), we obtain where ω is the energy separation of the Bloch bands, and the kinetic energy of the lattice atoms In a set-up where the wavenumber difference of the two running wave laser beams inducing the Raman process is parallel to the x-direction and of magnitude δq, we can define the effective Rabi frequency as Here, the Wannier functions w α (x) are defined as Similarly, we find that the terms describing onsite interactions between lattice atoms can be expressed aŝ with Here, g ⊥ = 2ω ⊥ a aa and ω ⊥ = 4V a,⊥ ω R is the oscillation frequency of the transverse confinement, ω R is the recoil frequency. Note that this expression is valid for U αα ω.
In order to obtain the expressions for equations (16) and (17) we have inserted equation (11) into equation (9) and performed the integrals over the transverse directions, where we have approximated the lattice potential by harmonic oscillators with frequency ω ⊥ . The interaction strengths (17) can be explicitly calculated if we also approximate the Wannier functions in x-direction with harmonic oscillator wavefunctions, and we find and U 10 = U 00 /2, U 11 = 3U 00 /4. In summary, the above model requires J α , U α,α , ω, ω ω ⊥ .

BEC reservoir and interaction with the lattice atoms
The BEC reservoir is described as a 3D homogeneous quantum gas consisting of N b particles of mass m b in a volume V bŷ where g bb = 4πa bb /m b , with a bb the scattering length for the interaction of reservoir atoms, and the operatorb † k creates a Bogoliubov excitation with momentum k and energy in the reservoir. The sound velocity in the BEC is given by c = The interaction of the lattice atoms with the superfluid reservoir is modelled by the densitydensity interaction Hamiltonian with the interaction strength g ab = 4πa ab /2m r , where a ab is the inter-species s-wave scattering length and m r = m a m b /(m a + m b ) the reduced mass. The field operators for the BEC can be expressed asψ and 12 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT in terms of creation and annihilation operatorsb † k andb k for a Bogoliubov excitation with momentum k = (k, k y , k z ). The coefficients u k and v k can be written as where We neglect the terms proportional to δψ † b (x)δψ b (x) [34] and using equation (11) and leaving out the constant mean field shift, we obtain Here, the coupling is where S(k) is the static structure factor In equation (27) we have again neglected the overlap of Wannier functions in different lattice sites and performed the integration over transverse lattice directions, where the wavefunctions are again approximated with ground state harmonic oscillator wavefunctions and result in factors of one, provided m b ω/2m a ω ⊥ 1, which is usually fulfilled for our 1D lattice potential. The interaction Hamiltonian equation (26) describes scattering processes, where a Bogoliubov excitation with momentum k in the (3D) BEC reservoir and a lattice atom with momentum q − k (with k the component of the momentum of the Bogoliubov excitation along the x-direction) in Bloch band α are annihilated (created) and a lattice atom with momentum q in Bloch band α is created (annihilated).
The structure factor S(k) → 1 for energies much larger than the chemical potential µ = m b c 2 = g bb ρ b , where the corresponding Bogoliubov excitations are particle-like, whereas for energies much less than µ the excitations are phonons and S(k) ∝ |k| → 0. In our set-up, we will choose 4J 0 µ ω (see figure 5), which corresponds to typical experimental parameters. In this case, interband transitions will involve absorption and emission of Bogoliubov excitations in the particle branch (E k ∼ ω). These will have much larger coupling strengths (S(k) ∼ 1) than intraband transitions, which involve absorption and emission of Bogoliubov excitations in the phonon branch (E k ∼ J α , S(k) ∝ |k| → 0) (cf figure 5).
In the derivation of the master equation (subsection 3.3) we will restrict ourselves to interband processes, i.e., to the decay from the first excited Bloch band back to the lowest band. In the following we will comment on intraband processes, which can lead to sympathetic heating or cooling of atoms within a Bloch band due to collisions with Bogoliubov excitations in the BEC reservoir. Excitations corresponding to the band separation ω are particle-like, small excitations with energies on the order of a Bloch band width are phonon-like. The chemical potential is larger than the width 4J 0 of the lowest band, but less than the band separation ω.
The heating/cooling process is described by the interaction Hamiltonian (26) and the corresponding rate can be estimated, e.g., with Fermi's golden rule. The heating/cooling processes involve the scattering of a lattice atom with a single Bogoliubov excitation, and Fermi's golden rule implies both energy and momentum conservation. We will illustrate this for the case of sympathetic heating within a Bloch band, as this represents a possible imperfection for the cooling process (whereas sympathetic cooling would only speed up the cooling) and the same arguments are valid for intraband cooling. The typical case of such a heating process in our scheme would be a scattering of an atom with momentum q ≈ 0 and energy ε 0 q≈0 to a momentum q with energy ε 0 q + c|k| (excitations within the Bloch band are sound waves), as typically most of the atoms are within a narrow region around q = 0 after a few iterations of the cooling process. This process is however only allowed if energy conservation c|k| = ε 0 q − ε 0 q and momentum conservation along the lattice axis k = q − q are fulfilled. As c|k| = c| k 2 + k 2 y + k 2 z | c|k| we find that this process is forbidden if the condition is fulfilled (see figure 6). For the typical parameters in our set-up (see figure 5), where µ, ω R J 0 , this condition will always be fulfilled, and thus sympathetic heating and cooling within a Bloch band are forbidden by energy and momentum conservation. Higher order processes arising from scattering with two or more BEC excitations will be small [34].

Master equation for decay
The collisional interaction (26) of atoms in the excited Bloch band of the optical lattice with the BEC reservoir leads to a decay of excited lattice atoms back to the lowest band, in analogy with spontaneous emission (see section 2). For typical BEC temperatures k B T b ω, the Bogoliubov modes corresponding to the band separation ω will initially be in the vacuum state, and we can derive an effective T = 0 master equation in the Born-Markov approximation for the reduced system density operatorρ, which describes atoms moving in the lattice, while the BEC is treated as a reservoir of Bogoliubov excitations. This is done in close analogy to [40] and [34] in the context of lattice loading of fermions from an external reservoir, and cooling of single atoms in a harmonic trap immersed in a BEC, respectively. We find (see appendix B)ρ where the 1D momentum k along the lattice axis is bounded by |k| k max = √ 2m b ω due to energy conservation, and the jump operatorsĉ k are defined aŝ with the position space operatorsÂ Note that an operator written A α q−k should always be understood asÂ α q , where q is a quasi-momentum in the first Brillouin zone, found by subtracting an integer multiple of the reciprocal lattice vector from q − k, i.e., q = q − k − zG, z ∈ Z, G = 2π/d. Similarly, the operatorĉ k ≡ĉ k+zG .
If a single atom is present in the excited Bloch band, the total decay rate via creation of excitations with all possible values of k is given by = k k . We can then define the distribution of emitted excitations, d /dk, which for deep lattices (where we can approximate Wannier functions by harmonic oscillator ground states, and ω |J 1 |, J 0 ) can be written explicitly as d dk= Here, a 0 = √ 1/m a ω is the size of the ground state of the harmonic oscillator in the x-direction and L = Md is the length of the lattice along the x-direction.
The total decay rate can be explicitly calculated in the harmonic oscillator approximation as where erf(x) denotes the error function. The value of can be tuned by changing the scattering length, the density of the BEC reservoir or the depth of the optical lattice.
If we compute the distribution for the change in quasi-momentum, q − q of the lattice atoms, we must translate values for k outside the first Brillouin zone, |k| > π/d back into the first Brillouin zone, as q = q − k − zG, z ∈ Z, G = 2π/d. This is indicated by the regions between the dashed lines in figure 7. We can distinguish two interesting limits for our parameters based on the ratio of the upper bound on |k|, k max = √ 2m b ω and the extent of the first Brillouin zone, π/d = √ 2m a ω R .
1. k max π/d. Here, the distribution d /dk extends over k values much larger than the first Brillouin zone, as depicted in figure 7. When we compute the distribution of changes in quasi-momentum for the lattice atoms, this will be approximately uniform over the first Brillouin zone. Note that this limit also corresponds to a situation in which the wavelength of emitted excitations, λ k = 2π/k is typically much shorter than the separation between lattice sites d. Thus collective effects of decay on different lattice sites (analogous to superradiance and sub-radiance in atomic decay) are suppressed. 2. k max < π/d. Here, the distribution d /dk is cut off before it reaches the edge of the first Brillouin zone. As a result, the distribution of q is localized at low values, peaking at the cutoff value k = ±k max . This can be used to target the decay to one area within the first Brillouin zone. For example, we can choose to excite atoms to a quasi-momentum value in the first band from which decay into the dark state will be strongly favoured.
Note that this limit also corresponds to a situation in which the wavelength of emitted excitations, λ k > d. Thus collective effects involving decay on different lattice sites occur, and this decreased spatial resolution of the decay process corresponds to the increased resolution that we observe in momentum space. These effects are properly accounted for in our calculations.

Single particle cooling
In this section, we will analyse the cooling process consisting of the two steps (i) the Raman laser excitation and (ii) the decay of excited lattice atoms. We will describe how efficient excitation laser pulses can be designed and present the results obtained from both numerical simulations and from analytical calculations based on Lévy statistics. In the excitation step, which represents the first part of our cooling protocol, we will assume that the interaction with the BEC reservoir can be switched off (e.g., via a Feshbach resonance [41]- [43]), whereas the Raman coupling is switched off during the decay step.

Designing the required laser pulses
We define the probability P j (q) that the j-th pulse (with j = 0..N p − 1 and N p the number of pulses) excites an atom with initial quasi-momentum q from the lowest band to the excited band and require P j (q) = 0 for q ≈ 0 and P j (q) → 1 for states with higher quasi-momentum (cf figure  4(a)). The probability P j (q) can be obtained by solving the Heisenberg equations of motion for the system operatorṡÂ 0 where the effective detuning δ q+δq j ≡ ω + ε 1 q+δq j − ε 0 q − δ. These equations are valid in the subspace of a single lattice atom, where the interaction HamiltonianĤ I = 0, and can be analytically solved for two simple cases. In the case of weak excitation, j (t) dt 1, we find the probability in terms of the Fourier transform of the Raman Rabi frequency j (t). In the case of a time square pulse ( j (t) = j for 0 t τ j and j (t) = 0 otherwise), we have The goal of the excitation step is to design efficient laser pulses, which excite atoms with large quasimomentum |q| > 0 but do not couple the atoms in the dark state q = 0. We would also like to do this on a fast timescale and therefore require a π-pulse, i.e., τ j = π/ j in the case of a time square pulse, for the resonant transition and adjust the parameters to always keep the first node of the sinc-function at q = 0. Such a pulse sequence consisting of N p = 3 pulses, is shown in figure 8. We start with an intense laser pulse to excite atoms with momentum qd ∼ π around the edges of the Brillouin zone, and then move the resonance closer to q = 0 by adjusting the Raman detuning δ j and momentum kick δq j and at the same time reducing the intensity of the laser beams. To be able to resolve the band structure with our Raman pulses we always have to fulfill 8|J 1 | and consequently τ π/8|J 1 |. Note, that the relevant energy scale is the hopping |J 1 | in the upper band, which is typically an order of magnitude larger than the hopping J 0 in the lower band for lattice strengths V a,x ∼ 10ω R . The parameters here have also been carefully chosen to avoid unwanted excitation to higher bands.

Results
In this section, we will quantitatively analyse the cooling process, computing the final temperature and cooling timescales. We make use of both numerical and analytical methods, and compare the results we obtain in each case.
In the numerical analysis we simulate the time evolution of the system density operator using a Monte Carlo wavefunction method [44,45]. In the simulations we start with an initial mixed state according to a thermal distribution of atoms in the lowest Bloch band with a typical temperature 4J 0 k B T ω. To obtain good approximations of the system density operator at all times we evolve the state according to the master equation (30) and typically take the statistical average over ∼ 10 5 trajectories in the simulations, for a 1D optical lattice, with M = 101 lattice sites and two bands.
The analytical calculations make use of Lévy statistics, similar to the corresponding calculations for freespace subrecoil laser cooling (see subsection 2.2). We define the temperature of the atom in terms of the width of the quasi-momentum distribution, and find with q denoting the half width of the momentum distribution at e −1/2 of the maximum. This is in close analogy to free space Raman cooling [38], with the mass in free space replaced with the effective mass in the optical lattice, and the quasi-momentum now playing the role of momentum in free space. We again find T ∝ −2/λ (see subsection 2.2, and consequently time square excitation pulses again lead to efficient cooling, as λ = 2 for time square pulses also in the presence of the optical lattice. Similarly we again find n 0 ( ) ∝ 1/λ . In figure 9, we compare the analytical and the numerical results for the temperature and the fraction of atoms in the dark state. In the simulations, repeated application of the N p excitation and decay steps leads to the development of a sharp peak in the momentum distribution n 0 (q) in the lowest Bloch band already after a few iterations, as can be seen in figure 9 (a).
In figures 9(b) and (c), we plot the temperature of the lattice atom, as defined in equation (37) and the height of the peak at q = 0 against time in units of J 0 , where circles denote numerical results from the quantum trajectory simulations, and solid lines are the analytical results from Lévy statistics. In both cases we find excellent agreement of the numerical results with the analytical predictions. Typical temperatures of k B T/4J 0 ∼ 2 × 10 −3 and a typical fraction on the order of a few tens of percent in the central peak can be reached in t f J 0 ∼ 50 for the pulse sequence shown in figure 8. In the simulations we have furthermore used = 1ω R , which can be obtained from equation (33) for ρ b = 5 × 10 14 cm −3 and a s = 350a b (e.g., via a Feshbach resonance) for 87 Rb in the lattice and 23 Na in the reservoir, a b is the Bohr radius. Note, that for the parameters we use (see caption in figure 8), the cooling timescale is mainly determined by the duration of the last two excitation pulses. The cooling timescale is typically much faster than the excitation pulses, and thus smaller values of does not significantly slow down the cooling process, unless the cooling timescale becomes comparable to the duration of the excitation pulses. In our cooling scheme we assume that the interaction of lattice atoms with the reservoir can be switched off during the excitation process. In a real experiment this can be done, e.g., via optical or magnetic Feshbach resonances [41]- [43], eventual residual finite interactions (e.g., due to magnetic field fluctuations for magnetic Feshbach resonances) can lead to a change of the excitation profile and thus represent a possible source of imperfection. From numerical simulations, however, we find that for J 1 , the final temperatures (see figure 10(a)) and the fraction of atoms in the dark state (see figure 10 (b)) do not change significantly. From equation (33), we find that ∼ 10 −3 J 1 can be achieved for ρ b = 5 × 10 14 cm −3 and a s = 20a b for 87 Rb in the lattice and 23 Na in the reservoir.

Many particle cooling
In this section, we will adapt the cooling scheme to many quantum degenerate bosons or fermions. We will assume that no interactions between the lattice atoms are present during the cooling process. In the case of bosons this means that the interaction Hamiltonian H I = 0 (see equation (16)), which can be achieved experimentally, e.g., via an optical or magnetic Feshbach resonance ([41]- [43]). We consider a single species of fermions, and thus s-wave interactions are forbidden due to Pauli blocking.
The cooling process again follows the protocol as in the case of a single lattice atom: (i) in the excitation step Raman laser pulses are designed to excite atoms outside the dark state region to the first excited Bloch band, the coupling of the lattice atoms to the reservoir atoms are switched off during this step. (ii) in the decay step the dissipative coupling to the BEC reservoir randomizes the quasi-momentum of lattice atoms decaying to the lowest band and consequently to a finite probability of falling into the dark state region.
One additional question in the context of cooling many atoms is that of re-absorption of emitted excitations, processes in which the reservoir excitation created during the decay of one atom is 'absorbed' by another atom, transferring it to the excited Bloch band. This question has also been addressed in the context of laser cooling (see, e.g., [46]). Importantly, re-absorption effects will be very substantially reduced for cooling 1D or 2D systems, as reservoir excitations will generally be emitted away from axis or plane of the system being cooled. In the present context, where we are often interested in systems of reduced dimension, this should restrict any adverse effects from re-absorption. In order to minimise re-absorption in the present context, the BEC reservoir could also be evaporatively cooled in the sense that particle excitations are allowed to quickly leave the trap, and do not interact multiple times with the lattice atoms.
For non-interacting lattice atoms, the dynamics of the cooling process is again described by the master equation (30). However, exact numerical simulations of the master equation are impractical, as the discretisation of the momentum space grid must be very fine to make possible calculation of the low final temperatures. Therefore, we perform the analysis of the cooling process by numerical simulations of a quantum Boltzman master equation (QBME) [45,47]. The QBME is one of the simplest versions of the more general quantum kinetic master equation (QKME) [45], which represents a fully quantum mechanical kinetic theory for the time evolution of the system density operator, and was originally developed to analyse formation of BEC in atomic gases. The QBME is an equation for the diagonal elements w m ≡ m|ρ|m of the reduced system density operator, which describes the time evolution of the system in terms of classical configurations w m of atoms occupying momentum states m = [{m 0 q } q , {m 1 q } q ] in the two Bloch bands. Here, m α q denotes the occupation of momentum state q in Bloch band α. In addition to the Born approximation and the Markov approximation made in the derivation of the master equation, the QBME neglects the off-diagonal coherences between different classical configurations contained in the QKME.
In our numerical simulations we require a full QBME only for the decay step (ii), as for non-interacting atoms the excitation step (i) can be exactly calculated from the Heisenberg equations (35), i.e., from the excitation probability P j (q) as The QBME for the second step (ii), the decay of excited lattice atoms back to the lowest Bloch band can be written as (see appendix C) where m = m − e q−k,q and the upper (lower) signs are for bosons (fermions). Finally, we want to remark that the coherences which are neglected in the description of the decay step in terms of a QBME could, in principle, even be destroyed in a real experiment. For example, this could be done, by modulating the lattice depth and thereby randomizing the off diagonal elements, similar to the twirl in state purification protocols [48].
In the following, we will present the results obtained from our numerical simulations first for the case of non-interacting bosons, then we will describe how finite interactions affect the excitation steps and finally we will present the results for spin polarized fermions.

Non-interacting bosons.
For N non-interacting bosons, the T = 0 ground state is the fully occupied q = 0 momentum state, as in the case of a single atom. As a consequence we can use the same excitation pulse sequences as for a single atom (cf figure 8). In figure 11(a), we plot temperature against time which we obtain from numerical simulations [45] of the QBME equation (39). Temperature here is calculated by a Gaussian fit to the momentum distribution, excluding the central peak. Cooling to similar temperatures as in the single particle case can be obtained on shorter timescales for many atoms. This is due to the bosonic enhancement factor, which appears as the factor (1 + m 0 q−k ) in the QBME (38). In figure 11 (b), we show the increase of the central peak against time.

Excitation profile for bosons with interactions.
Until now, we have assumed that the interactions between the Bosons a is negligible. In this section, we investigate how small finite interactions will alter the excitation profile, and give approximate values for interaction strengths that can be safely neglected. Here, we must compute the time evolution of the many-body system during the excitation step, which is described by the full two-band Bose-Hubbard-Hamiltonian for interacting atoms given in subsection 3.1 asĤ 0 +Ĥ LS +Ĥ I . This is possible in 1D at temperature T = 0 using time-dependent DMRG methods [49]- [51]. These numerical methods allow us to compute the time evolution of a 1D manybody system where the Hilbert space can be expressed as the product of local Hilbert spaces of dimension d forming a 1D chain. In these methods, the state of the system is written as a truncated matrix product state representation [49], in which χ states are retained in each Schmidt decomposition of the system. The dynamical evolution for 1D systems with sizes similar those seen in experiments can then be computed, starting both from the ground state and from weakly excited states.
We simulate the time evolution in this way for the duration of one coherent Raman pulse. We consider the situation with all the interactions equal, i.e., U 00 = U 01 = U 11 . To minimize the influence of box boundary conditions, we used M = 41 sites and an initial Fock state of the form |0 · · · 0111110 · · · 0 with N = 5 atoms located at the centre of the system in the lower band. Since Fock states have a flat equally occupied momentum distribution they allow the excitation profile to be determined by examining the final momentum distribution. A χ = 50 and an occupancy cutoff of n 0 max 4 and n 1 max 2 atoms per site for the lower and upper band respectively (equivalent to s = 15) was found to be sufficient. In figure 12(a), the momentum distribution for the lowest band n 0 (q) after the pulse has been applied is displayed for a sequence of interaction strengths U αα /J 0 . For the parameters chosen when U αα = 0 the initially flat momentum distribution has a hole carved out at q = 4 q, where q = 2π/Md. This demonstrates the selective excitation of this momentum and crucially the dark state at q = 0 is left unchanged. As interactions are switched on the pulse begins to distort the final momentum distribution from this characteristic shape with a peak in the population of momenta either side of q = 4 q as well as an increase in the population remaining at q = 4 q itself. This latter quantity provides a useful measure of the shape of the excitation profile and is plotted in figure 12(b) and displays a linear increase only after U αα /J 0 > 0.2. The most relevant quantity for the cooling scheme is n 0 (0), the population at q = 0, which is seen to increase linearly for small interaction strengths in figure 12 (c). In total this reveals that for U αα 1/τ |J 1 | the excitation profile remains close to that assumed for U αα = 0. If we include the additional constraint that U αα |J 0 |, which ensures that interactions do not substantially alter the ground state, the conclusions for the excitation pulse used earlier should not change.

Fermions.
In contrast to bosons, in the case of N non-interacting (spin-polarized) fermions the T = 0 ground state is characterized by the T = 0 Fermi distribution, i.e., a step function n 0 (q) = (q − q F ) for the momentum distribution in the lowest band, where the Fermi momentum q F = π(N − 1)/Md. The excitation pulse sequence thus has to be changed in order to create a dark state region for atoms with momenta |q| q F . As a consequence the use of time square pulses is no longer advantageous, due to the large sidelobes they create in the  excitation probability. We thus use a sequence of Blackman pulses [32] where these sidelobes are suppressed, as shown in figure 13, where we compare typical excitation probabilities for a Blackman pulse and a time square pulse. The sequence of excitation pulses is now designed as in the previous cases, however the excitation probability can only be calculated numerically and the π-pulse condition changes to τ = π/0.42 for Blackman pulses [38].
In figure 14(a) we plot the momentum distribution (again obtained from Monte Carlo simulations of the QBME) after j = 0, 1, 2, 20 cooling cycles and find that the expected shape close to the expected T = 0 Fermi distribution appears after a few cooling cycles. Temperature is now obtained by fitting a Fermi distribution to occupation of the momentum states in the lowest Bloch band. In figure 14(b) we plot temperature against time in units of J 0 and find that temperatures k B T/4J 0 ∼ 10 −2 can be obtained in tJ 0 ∼ 500 for the parameters used (see caption of figure 14). The timescales for many fermions are slower than those of bosons due to Pauli blocking, which increasingly slows down the decay into an increasingly filled Fermi sea.

Realizing cold strongly-interacting gases
In this section, we investigate how strongly correlated regimes can be realized where the cooling scheme cannot be applied directly. This can be done by decoupling the optical lattice from the reservoir b and adiabatically ramping up the interaction strength. In this case our attention is restricted to the evolution governing byĤ 0 +Ĥ I with no Raman coupling and α = 0, and so describes the Bose-Hubbard model of the lowest band only. We analyse this by again making use of time-dependent DMRG methods [49]- [51].

Ramping the ground state
In principle adiabaticity requires an infinitely slow ramping. Here, we determine a finite timescale in which near-adiabatic ramping can be achieved. We consider a 1D lattice of M = 21 sites containing N = 10 atoms initially prepared in the ground state with U 00 = 0 and then raise the interaction strength over a time τ r according to The constant C is fixed so as to make U 00 (t = 0) = 0 and we choose the parameter w = 18 τ r so that U 00 (t = τ r ) U max . As this simulation begins from a non-interacting limit we took the occupancy cut-off to be n 0 max six atoms per site (equivalent to s = 7). We found that for the computation of groundstates retaining χ = 50 states in the matrix product decomposition was sufficient for this system, however, to accurately describe the dynamical ramping a χ = 100 was required. For the ramping itself we took U max = 20 J 0 and so approached the hard-core Bose lattice gas in 1D (Tonks gas).
To quantify the achievement of near-adiabatic ramping we computed the energy deposited into the system E(τ r ) = [E(τ r ) − E 0 ]/NJ 0 and the many-body overlap F(τ r ) = | ψ(τ r )|ψ 0 | as a function of the ramping time τ r of the final ramped state |ψ(τ r ) with energy E(τ r ) and the groundstate |ψ 0 with energy E 0 for U 00 (τ r ). Note that we express all energy differences E as a fraction of NJ 0 . This is a useful energy scale since it is (up to a constant prefactor of order one) the maximum energy difference which N non-interacting atoms inside the lowest band could have. Both E(τ r ) and F(τ r ) are plotted together in figure 15 (a). We observe that for τ r > 10/J 0 there is an exponential decrease in E(τ r ) with τ r , and thus for τ r of order U 00 (τ r )/J 0 there is negligible heating within the Bloch band. Additionally F(τ r ) can be seen to rapidly approach unity on the same timescale rigourously verifying that this final state is converging to the groundstate in the strongly correlated regime. We note, however, that since these calculations are based on a single system size they do not address the issue of scaling with M. Despite this, we expect the results presented to be applicable to systems which are currently realized in experiments since their size is typically of the same order as considered here.

Ramping weakly excited states
In the previous section the degree of excitation induced by the ramping process on the ground state was quantified. Here, we confirm that the near-adiabatic timescale determined for the ground state also applies to good approximation to low-lying excited states. This is done by performing the ramping with a weakly excited initial state and demonstrating that the resulting final state remains weakly excited in the strongly interacting regime. Specifically, we generate an excited state by evolving the U 00 = 0 ground state for a time t ex in the presence of a spatially homogeneous on-site interaction which is varying randomly in time in the range U 00 (t) ∈ [0, 1 2 J 0 ]. With this method we constructed three excited initial states using J 0 t ex equal to (i) 1, (ii) 5 and (iii) 9. We characterize these states by their normalized energy difference E, many-body overlap F with the ground state, quantum depletion D (equal to N minus the largest eigenvalue of the singleparticle density matrix), and energy spread defined by ( E/NJ 0 ) 2 = Ĥ 2 − Ĥ 2 (see caption of figure 15 for specific numbers). Despite this being classical noise which coherently excites the ground state it does produce excited states with features similar to those found in experiments. In particular the characteristics for state (i) may be representative of a weakly excited superfluid state resulting from the cooling scheme outlined.
The three initial states were then ramped in an identical way to subsection 6.1 using τ r = 100/J 0 . In figure 15(b) the evolution of energy difference E(t) as a function of the time t during a ramp is shown for each of them. From the two most excited of these states, (ii) and (iii), E(t = τ r ) is seen to level-off at around four times their initial values giving on the order of 10% heating within the Bloch band, and have E(t = τ r ) 0.17, while the many-body overlaps of their final ramped states with the strongly-correlated ground state reduce to F = 0.55 and F = 0.41 respectively. This indicates that the ramping is not entirely adiabatic for these states. For the least excited initial state (i) we find that both E(t = τ r ) and E(t = τ r ) have approximately doubled, but crucially the heating E(t = τ r ) is still less than the 4%. Finally, for (i) the many-body overlap F = 0.89 suffers a smaller reduction and importantly remains sizable showing that the final ramped state is still weakly excited in the strongly-correlated regime.

Spatial sideband cooling in a harmonic trap
Up until now, we have discussed the implementation of a Raman cooling scheme in an optical lattice with a flat external potential, based on the coupling of atoms in a lattice to a reservoir gas. As an additional remark, we describe in this section how these ideas can be applied in a different way to cool atoms in a harmonic trap to sites of lower energy, in analogy with sideband cooling to lower motional states in an ion trap. [36]. This example serves to strengthen the formal Figure 16. (a) Techniques involving spontentous creation of excitions in the reservoir can also be used to cool fermions in an external trapping potential to sites of lower energy, compacting them so that in the centre of the trap we observe unit occupancy. (b) Atoms are excited to the first excited motional level in each well by a Raman transition (with two-photon Rabi frequency ) detuned to energies below resonance with the excited state. Atoms can then tunnel to neighbouring sites with tunnelling rate J 1 , and can decay back to the ground motional state by interaction with the external reservoir gas. The dashed arrows show the effective result of the two-step excitation and tunnelling process, and we see that transitions to sites of lower energy are favoured because these are closer to resonance. (c) Results for a single particle being cooled in this manner: plot of the mean position of the atom in a single quantum trajectory as a function of time, beginning with an atom displaced 12 sites from the centre of a harmonic potential, analogies between open quantum systems encountered in quantum optics, and systems of atoms in optical lattices coupled to a phonon bath. It can also be applied in the context of preparation of a quantum register with fermions [40,52], in this case compacting the cloud of atoms as they collect in the centre of the trap, at the sites of lower energy [see figure 16 (a)]. Here, we outline how single atoms can be cooled to sites of lower energy, and then briefly discuss the application to many fermions.
We consider the same set-up discussed in section 1, with an identical Hamiltonian, equation (12), except that here we impose an additional external trapping potential (normally a harmonic trap),Ĥ E = i ini , wheren i is the number operator on site i. In this scheme, Raman coupling of atoms from the ground motional state in each lattice site to the excited motional state is not pulsed, but continuously switched on. The interaction with the external reservoir gas b, giving rise to the master equation (30) is also identical, and is also continuously switched on, giving rise to onsite decay. We assume here that the lattice depth is large, so that we are in a limit where collective processes (analogous to superradiance and subradiance in phonon emission) can be neglected (see subsection 3.3).

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
The basic concept of this process is shown in figure 16(b). We choose the detuning of the Raman excitation so that it is below resonance with the excited motional level. If the atom subsequently tunnels (with amplitude J 1 ) to a neighbouring lattice site, and then undergoes a decay, it can be transferred to neighbouring lattice sites. (Note that tunnelling in the lowest band J 0 J 1 , and is further suppressed by the potential offset between the lattice sites, due to which tunnelling in the lowest band is no longer a resonant process.) Transfer to sites of lower energy is favoured, however, because of our choice of detuning for the Raman excitation. In figure 16(b), dashed arrows indicate an effective two-step process including the Raman excitation and tunnelling to neighbouring sites. We note that this process can be made resonant for a site of lower energy, whilst far detuned from a site of higher energy. This is similar to the concept of sideband cooling in ion traps. There, the goal is to cool the motional state of an atom in a single trap, by coupling to states in which the electronic state is excited, but the motional state is reduced by one quantum (this is called the red sideband). We can draw a figure analogous to figure 16(b), so that in comparison with sideband cooling, we have replaced excited and ground electronic states with excited and ground motional states in each well, and the coupling to lower motional states in sideband cooling is replaced by coupling to neighbouring sites of lower energy. In ion traps the goal is to cool an ion to the lowest motional state, whereas in this proposal, sideband cooling leads to a spatial redistribution of atoms towards the centre of the harmonic trap.
In analogy with sideband cooling, the cooling steps in this scheme (coupling to sites of lower energy) will more strongly dominate heating steps (coupling to sites of higher energy) when the energy offset between the sites is made larger, and when the detuning δ ∼ − , where is the energy offset between sites. Note that in a harmonic trapping potential, where the potential difference between neighbouring sites varies, the implementation should involve a sweep in the detuning, making the cooling more efficient in different parts of the trap as a function of time.
In figure 16(c), we plot example results for a single atom being cooled in a harmonic trapping potential, obtained from Monte Carlo wavefunction simulations. The solid line in the figure shows the mean position of the atom as a function of time for a single quantum trajectory, which illustrates the cooling and heating transitions as the atom moves closer to and further from the centre of the trap respectively. We also show the root-mean-square displacement from the centre of the trap [i.e., ( j j 2n j )/M ] when we begin with a uniform distribution for the initial position. The finite final temperature is determined by competition between heating and cooling steps, as is clearly illustrated by the transitions in the example trajectory.
This method can be simply extended to fermions, where double occupation of the ground motional state in any lattice site is prevented by Pauli blocking. For this case it is possible to derive a QBME in analogy to that in section 5, and to simulate these equations. Compaction of the fermions into the centre of the trap is observed, producing a quantum register with one fermion on each lattice site [53]. An additional step will, in general, be required to remove any remaining atoms from the upper motional level at the end of the cooling process. In the present form, this scheme is not well suited for bosons, as it would be difficult to prevent many atoms from collecting on a single site, even in the presence of interactions. However, the existence of this analogy is again a strong demonstration that ideas from quantum optics can be used in the context of this new type of open quantum system, an idea that has many possibilities for future exploration.

Conclusions
We have analysed a new cooling scheme based on ideas borrowed from sub-recoil Raman cooling schemes, but where these ideas are now placed in the context of a new form of open quantum system, where atoms in an optical lattice are coupled to a BEC reservoir. In our case, this setup provides laser assisted sympathetic cooling, in which the final temperature of atoms in the lattice is not limited by the temperature of the reservoir. This is motivated by the practical requirement of achieving low temperatures in an optical lattice, which is important for simulation of many strongly correlated quantum systems. Our scheme indeed provides a possibility to achieve low temperatures, and from our model we predict that temperatures a small fraction of the Bloch band width can be achieved in this way.
Equally importantly, however, this study opens questions and possibilities on the implementation of open quantum systems in this new context. Here, we have demonstrated how ideas from quantum optics can be applied in a many-body system where we have strong control over many parameters, especially interactions between atoms, and thus our system-reservoir coupling. The possibility exists to extend these ideas to different types of dissipative Hubbard models, and quantum reservoir engineering more generally. These ideas could be used both in the context of state preparation and assist the engineering of new models with systems of many atoms in optical lattices.
In the opposite direction, such open quantum systems also offer possibilities for investigating effects that are of fundamental interest in quantum optics, here in parameter regimes that can be very different from regular quantum optics experiments. One example that is clear in the current context would be the possibility to study the analogue of superradiant or subradiant decay, using the lattice and BEC parameters to engineer the wavelength of the excitations created in the BEC, and therefore the strength of these effects.
after integrating the resulting equation, and assuming 1 (t), 2 (t) and A α q (t) to be slowly varying on the timescale 1/ .
As the detuning of the two Raman lasers from the excited state is the largest frequency, i.e., ω α , δ, ε α q , we can neglect all terms rotating with , and thuŝ Here, we have assumed two running wave Raman lasers with relative momentum δq and neglected the overlap of Wannier functions in different lattice sites. Transforming back to the Schrödinger picture with respect to the lattice and only taking into account the lowest two Bloch bands we obtain the Hamiltonian as given in equation (12).

Appendix B. Derivation of the master equation
In this section, we will derive the master equation for the reduced system density operatorρ for the decay of lattice atoms from the first excited Bloch band back to the lowest band as given in equation (30) . As described above, the BEC can be treated as a (3D) T = 0 reservoir, and the Born-Markov master equation in the interaction picture is given by [44]ρ whereρ R is the density operator for the BEC reservoir, and Tr R { } denotes the trace over the reservoir states. In a standard way we change to the variable τ ≡ t − t and extend the integration 31 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT equation (30). We project on the diagonal elements of the density operator and finḋ The expectation values in the first part (C.1) of this equation can be calculated as where the upper (lower) signs are for bosons (fermions), the product of the two δ functions appearing in these expressions can be evaluated as δ m,n+e q −k,q δ m,n+e q−k,q = δ m,n+e q−k,q δ q,q . For the second and third line (C.2) and (C.3) we find the identical expressions