Simulating quantum-optical phenomena with cold atoms in optical lattices

We propose a scheme involving cold atoms trapped in optical lattices to observe different phenomena traditionally linked to quantum-optical systems. The basic idea consists of connecting the trapped atomic state to a non-trapped state through a Raman scheme. The coupling between these two types of atoms (trapped and free) turns out to be similar to that describing light-matter interaction within the rotating-wave approximation, the role of matter and photons being played by the trapped and free atoms, respectively. We explain in particular how to observe phenomena arising from the collective spontaneous emission of atomic and harmonic oscillator samples such as superradiance and directional emission. We also show how the same setup can simulate Bose-Hubbard Hamiltonians with extended hopping as well as Ising models with long-range interactions. We believe that this system can be realized with state of the art technology.


Introduction
Cold atoms trapped in optical lattices have been proved to be very versatile quantum systems in which a large class of many-body condensed-matter Hamiltonians can be simulated (see [1,2,3] for extensive reviews on this subject). One of the first proposals in this direction was the article by Jacksch et al. [4], where it was proved that the dynamics of cold atoms trapped in optical lattices is described by the Bose-Hubbard Hamiltonian provided that certain conditions are satisfied; shortly after, the superfluidto-Mott insulator phase transition characteristic of this Hamiltonian was observed in the laboratory [5]. Since then, the broad tunability of the lattice parameters, and the increasing ability to trap different kind of particles (like bosonic and fermionic atoms with arbitrary spin or polar molecules), has allowed theoreticians to propose optical lattices as promising simulators for different types of generalized Bose-Hubbard and spin models which are in close relation to important condensed-matter phenomena [1,2,3]. Recent experiments have shown that optical lattices can be used to address open problems in physics like, e.g., high-T c superconductivity [6], to study phenomena in low dimensions such as the Berezinskii-Kosterlitz-Thouless transition [7], or to implement quantum computation schemes [8].
In this article we keep digging into the capabilities of optical lattices as simulators, showing how they can also be a powerful tool for simulating quantum-optical phenomena. In particular, following our previous proposal [9], we introduce various schemes in which superradiance-like phenomena can be observed. In addition, we will show that these same schemes realize Bose-Hubbard and spin models with extended interactions (in particular, models where hopping between sites and ferromagnetic interactions between spins are long range instead of nearest-neighbors).
The concept of superradiance was introduced by Dicke in 1954 when studying the spontaneous emission of a collection of two-level atoms [10] (see [11] for a review). He showed that certain collective states where the excitations are distributed symmetrically over the whole sample have enhanced emission rates. Probably the most stunning example is the single-excitation symmetric state (now known as the symmetric Dicke state), which instead of decaying with the single-atom decay rate Γ 0 , was shown to decay with NΓ 0 , N being the number of atoms. He also suggested that the emission rate of the state having all the atoms excited should be enhanced at the initial steps of the decay process, which was a most interesting prediction from the experimental point of view, as this state is in general easier to prepare. However, Dicke used a very simplified model in which all the atoms interact with a common radiation field within the dipolar approximation; almost 15 years later, and motivated by the new atom-inversion techniques, several authors showed that dipolar interactions impose a threshold value for the atom density, that is, for the number of interacting atoms, in order for superradiance to appear [12,13,14].
It is worth noting that together with atomic ensembles, spontaneous emission of collections of harmonic oscillators were also studied at that time [13]. Two interesting features of this system were reported: (i) an initial state with all the harmonic oscillators excited does not show rate enhancement, on the contrary, most excitations remain within the sample in the steady state, while (ii) the state having all the oscillators in the same coherent state has a superradiant rate.
In the last years, there has been a renewed interest in this phenomenology because of its potential for quantum technologies [15,16,17,18,19,20]. In these new approaches, superradiant effects differ from that found by Dicke in that they are mediated by the initial entanglement of the emitters, not by a build up of their coherence. Of particular interest to our current work is the analysis performed in [18] (see also [17]) of the spontaneous emission by regular arrays of atoms separated a distance d 0 comparable to the wavelength λ of the emitted radiation. It was shown that directional emission can be obtained for particular single-excitation entangled states when λ > 2d 0 , the emission rate being enhanced by a factor χ ∝ (λ/d 0 ) 2 N 1/3 . Also of interest for our purposes is the work carried by Sajeev John and collaborators on the collective emission of atoms embedded in photonic band-gap materials [21,22]. In this kind of materials the density of states of the electromagnetic field is zero for frequencies laying within the gap, what gives rise to the phenomenon of light localization [23]. As for collective emission, it was shown by using a simplified Dicke-like model [21,22] that when the atomic transition lies close to the gap, radiation is emitted in the form of evanescent waves and the single-excitation symmetric state has a rate proportional to N 2 Γ 0 , which is even larger than the one predicted by Dicke in free space.
We will show that all this phenomenology can be observed in optical lattices with the setup already presented in [9], which is depicted in figure 1. Consider a collection of bosonic atoms with two relevant internal states labeled by a and b (which may correspond to hyperfine ground-state levels, see for example [1]). Atoms in state a are trapped by a deep optical lattice in which the localized wavefunction of traps at different lattice sites do not overlap (preventing hopping of atoms between sites), while atoms in state b are not affected by the lattice, and hence behave as free particles. A pair of lasers forming a Raman scheme drive the atoms from the trapped state to the free one [1], providing an effective interaction between the two types of particles. We consider the situation of having non-interacting bosons in the lattice [24], as well as hard-core bosons in the collisional blockade regime, where only one or zero atoms can be in a given lattice site [25]. In the first regime, the lattice consists of a collection of harmonic oscillators placed at the nodes of the lattice; in the second regime, two-level systems replace the harmonic oscillators, the two levels corresponding to the absence or presence of an atom in the lattice site. Therefore, it is apparent that this system is equivalent to a collection of independent emitters (harmonic oscillators or atoms) connected only through a common radiation field, the role of this radiation field being played by the free atoms. This system is therefore the cold-atom analog of the quantum-optical systems considered above, with the difference that the radiated particles are massive, and hence have a different dispersion relation than that of photons in vacuum. Moreover, we show that the field of free atoms can be characterized by a dispersion relation which is similar Figure 1. Scheme of our proposed setup. Atoms in state a are trapped in an optical lattice, while atoms in state b are free, and can thus have any momentum. An external pair of Raman lasers connect the two levels with some detuning ∆. We will show that above some critical value ∆ = ∆ c the atoms in state b are able to leave the trap (left), while below this, they are trapped forming a bound state with atoms in state a (right). This behavior is typical of band-gap systems where the density of states is zero between the connected states.
to the one obtained for photons within a photonic band gap material [9,21,22]. It is then to be expected, and so we will prove, that this system will show the same kind of phenomenology as its quantum-optical counterparts.
To conclude this introduction let us explain how the article is organized. We first introduce the model and find its associated Hamiltonian (Section 2). Then we study the radiative properties of the system when the lattice sites emit independently (Section 3); we will show in this section that our system shows localization of the free atoms, in analogy to the localization of light occurring in a photonic band gap material, and we identify as well the Markovian and non-Markovian regimes of the emission. In Section 4 we study collective effects. We first deduce a reduced master equation for the lattice atoms, explaining under which conditions it is valid. We then show how by changing the system parameters extended Bose-Hubbard and spin models (Section 4.1), Dicke superradiance of atomic (Section 4.2.1) and harmonic oscillator (Section 4.2.2) samples, and directional superradiance (Section 4.3) can be observed. We finally talk about the effect of restricting the motion of the free atoms to 2D or 1D traps in Section 5, and give some conclusions in Section 6.

The model
Our starting point is the Hamiltonian of the system in second quantization [26]. We will denote by |a and |b the trapped and free atomic states, respectively (having internal energies ω 0 a and ω 0 b ). Two-body interactions for the trapped atoms are included with the usual contact-like pseudopotential [3,24], but we neglect the collisions for the free atoms. The Hamiltonian is then written asĤ =Ĥ 0 +Ĥ a−b , witĥ H 0 contains the individual dynamics of the atoms, H j being the first-quantized motion Hamiltonian of the atom in the corresponding state, and g = 4π 2 a s /m, where a s is the s-wave scattering length of the trapped atoms (which have mass m).Ĥ a−b contains the Raman coupling between the atomic states, k L = k 1 − k 2 and ω L = ω 1 − ω 2 (laser wave vector and frequency in the following) being the relative wave vector and frequency of the two lasers involved in the Raman scheme (see figure 1), with Ω the corresponding two-photon Rabi frequency. For atoms in state |a , H a = − ( 2 /2m) ∇ 2 + V opt (r), where V opt (r) corresponds to a 3-dimensional optical lattice with cubic geometry and lattice period d 0 . We work with ultracold atoms under conditions such that their wavefunctions can be described by the set of first-band Wannier functions localized around the nodes of the lattice [1]. The traps of the optical lattice are approximated by isotropic harmonic potentials [1], what allows us to write the Wannier functions as where r j = d 0 j is the position of the j ∈ Z 3 lattice site (we consider M sites in each orthogonal direction defining the cubic lattice), and X 2 0 = /mω 0 , being ω 0 the frequency of the harmonic trap. The energy associated to these wave functions is E 0 = 3 ω 0 /2. We will assume that Wannier functions localized at different lattice sites do not overlap, hence preventing tunneling between sites.
On the other hand, atoms in state |b can move freely in every direction of space according to H b = − ( 2 /2m) ∇ 2 , and hence the plane waves ψ k (r) = e ik·r / √ V , with energy E k = 2 k 2 /2m, are their motion eigenfunctions (V is the total available volume for the free atoms, which we might take as infinite for calculations).
We consider two opposite regimes for the interaction between trapped atoms (remember that in (1a) any two-body interaction involving the free atoms has been neglected). The first limit consists in neglecting the interactions, which might be accomplished by, e.g., tuning the scattering length with an additional magnetic field through a Feshbach resonance [24]. Expanding the quantum fields onto their corresponding motion eigenstates, we get within the interaction picture: with being the detuning of the laser frequency respect to the |a ⇋ |b transition (ω a = ω 0 a + 3ω 0 /2 and ω b = ω 0 b ). The operators {â j ,â † j } and {b k ,b † k } satisfy canonical bosonic commutation relations, and create or annihilate an atom at lattice site j and a free atom with momentum k, respectively.
As for the second limit, we assume that the on-site repulsive atom-atom interaction is the dominant energy scale, and hence the trapped atoms behave as hard-core bosons in the collisional blockade regime, what prevents the presence of two atoms in the same lattice site [25]; this means that the spectrum ofâ † jâ j can be restricted to the first two states |0 j , |1 j , having 0 or 1 atoms at site j, and then the boson operators â † j ,â j can be changed by spin-like ladder operators σ † j ,σ j = {|1 j 0|, |0 j 1|}. In this second limit the Hamiltonian readŝ Hamiltonians (3) and (5) show explicitly how this system mimics the dynamics of collections of harmonic oscillators or atoms, respectively, interacting with a common radiation field. Note that, particularly, the dispersion relation appearing in (4) is similar to that of the radiation field within an anisotropic 3D photonic band-gap material, where photons acquire an effective mass close to the gap [21,22]. We will show how by varying the parameters of these Hamiltonians, the system has access to regimes showing many different phenomena, as explained in the Introduction.
Note finally that in order to satisfy that the trapped atoms are within the first Bloch band, it is required that ω 0 ≫ ∆, Ω.

Emission of an atom from a single site
Before studying the collective behavior of the atoms in the lattice, it is convenient to understand the different regimes of emission when the sites emit independently. To this aim we first study the emission properties of one atom in a single site following the analysis performed in [21,22] for an atom embedded in a photonic band-gap material. We are going to show that there exists a critical value of the detuning above which the atom is emitted, while below which the radiated atom gets bound to the trapped atom, and there is a nonzero probability for the atom to remain in the trap (see figure 1). This second regime is the analog of the photon-atom bound state predicted more than 25 years ago [23]. We will also identify the Markovian and non-Markovian regimes of the emission.
Following the Weisskopf-Wigner procedure [27,28], we write the state of an atom in a single site (which we take as j = 0) as where |1, {0} refers to the state with an atom in the trap and no free atoms, and |0, 1 k to the state with no trapped atoms and one free atom with momentum k. We can use the Schrödinger equation to write the evolution equations of the coefficients A and B k ; then, the equation of B k can be formally integrated arriving to a single equation for A given byȦ where we have assumed that the atom is initially in the trap, that is, A (0) = 1 and B k (0) = 0. This is a non-Markovian equation where the free atoms enter the dynamics of the trapped atoms as a reservoir with correlation function where we have taken the continuous limit for the momenta in the last equality and k L = 0 for simplicity. Note that this correlation function coincides exactly with that of atoms in anisotropic 3D photonic band-gap materials [21,22]. In Appendix A we show how to handle equation (7) by using Laplace transform techniques. In particular, in Appendix A.1 we show that in the strong confinement regime ω 0 → ∞, this equation can be solved as with Solution (9) has a very suggestive form. The first term has exponential evolution as expected for a Markovian radiation process, while the integral term has a nontrivial time dependence which always goes to zero for sufficiently large times. The time while the integral term is comparable to the exponential one defines the non-Markovian part of the evolution. This result also shows that the behavior of the emission is highly dependent on the parameters of the system: • In the∆ < 0 region Im b 2 + = 0, and hence the steady state population is nonzero and given by is radiated during the non-Markovian period which lasts longer as∆/α 2 goes to zero (see figure 2a). In this region the radiated atoms are emitted in the form of evanescent waves, and the equivalent of a photon-atom bound state [23] is formed.
• In the 0 <∆ < πα 2 region (figure 2b) the steady state population is zero, i.e., the atom leaves the trap eventually, and the evolution is dictated solely by the integral part of the solution (hence it is pure non-Markovian radiation).
• Finally, in the∆ > πα 2 region the atom is radiated in a Markovian fashion only for ∆ ≫ α 2 , with a decay rate given by Γ 0 = Im b 2 − ≈ 4Ω 2 2π∆/ω 3 0 ; if this is not the case (∆ ∼ α 2 ), the integral part is comparable to the exponential part during most of the evolution time, and hence the radiation process is non-Markovian (figure 2c).
Therefore, there exists a phase transition at ∆ = 4Ω 2 /ω 0 in close analogy to that of spontaneous emission in photonic band-gap materials [21,22]: Above this value the atom is forced to leave the trap, while below it there is a nonzero trapped population left in the steady state. Using the techniques we explain in Appendix A.2, it is possible to show that this phase transition is also present for finite ω 0 (see figure 2d).

Collective dynamics
Let us study now the collective emission of atoms from many sites. We analyze this through the master equation for the reduced density operator of the lattice within the Born-Markov approximation [29]. We introduce the Born approximation by fixing the state of the reservoir to a nonevolving vacuum. This leads to an effective interaction between different lattice sites mediated by the two-point correlation function being . The Born approximation yields the dynamics of the system to the leading order in the coupling parameters g k [29]. On the other hand, as commented at the end of the previous section, the Markov approximation is valid only if |∆| ≫ Γ, where Γ is a typical evolution rate for trapped atoms, which can differ from Γ 0 owed to a renormalization due to collective effects as we show below. We also include the single-site energy shift that we found in the previous section (the analogous of the Lamb shift in the optical case), and substitute ∆ by∆ in the two-point correlation function. The master equation for the reduced density operator of the trapped atomsρ is found by using standard techniques from the theory of open quantum systems [29], and reads dρ dt = j,l a similar equation is obtained for the hard-core bosons but replacing the boson operators by the corresponding spin operators. The effective interaction between different lattice sites is mediated by the Markov couplings (see Appendix B) . This expression has been evaluated in the strong confinement regime ω 0 → ∞ and considering k L ≪ X −1 0 (see Appendix B). Note that the term '1 − erf (d 0 |j − l| /2X 0 )' is basically zero for j = l, and therefore the ξ parameter dictates the spatial range of the interactions as can be appreciated in figure 3. Finally, we remind that Γ 0 = 4Ω 2 2π|∆|/ω 3 0 is the single-emitter decay rate.

Extended Bose-Hubbard and spin models
For k L = 0 and∆ < 0 the Markov couplings are purely imaginary and have negative imaginary parts. Under this conditions, the master equation takes a Hamiltonian form with effective Hamiltonianŝ for non-interacting and hard-core bosons, respectively. The first Hamiltonian describes extended hopping in the lattice, a feature that has been recently shown to be helpful for achieving true incompressible Mott phases which otherwise can be blurred because of the additional slowly varying harmonic trap used to confine the atoms within the lattice [30]. The second Hamiltonian is equivalent to an extended ferromagnetic Isinglike Hamiltonian, and connects our system to extended spin models. Note that for j = l the couplings have a Yukawa form and can even take a Coulomb form if ξ is large enough. These kind of interactions (specially the Coulomb-like) are difficult to obtain with other techniques. These results show that in the trapping region the system under study could be useful for quantum simulation of condensed-matter phenomena requiring extended interactions.

Dicke superradiance
In this section we connect our system to superradiant Dicke-like phenomena both in atomic and harmonic oscillator samples (we will still take k L = 0).
For∆ > 0, the Markov couplings are complex in general, and therefore the master equation of the system takes the form dρ dt with a dissipation term given by having collective decay rates and a reversible term corresponding to inhomogeneous dephasing with Hamiltonian being The same holds for hard-core bosons but replacing the boson operators by the corresponding spin operators. For ξ ≪ 1 the Markov couplings do not connect different lattice sites, that is Γ j−l ≃ Γ 0 δ j,l , and the sites emit independently. On the other hand, when ξ ≫ M the collective decay rates become homogeneous, γ j−l ≃ Γ 0 , and we enter the Dicke regime. Hence, we expect to observe the superradiant phase-transition in our system by varying the parameter ξ.
Note that in the Dicke regime the dephasing term cannot be neglected and connects the sites inhomogeneously with Λ j =l ≃ Γ 0 ξ/ |j − l|. This term appears in the optical case too, although it was inappropriately neglected in the original work by Dicke [10] when assuming the dipolar approximation in his initial Hamiltonian, and slightly changes his original predictions as pointed out in [31,32] (see also [11]).
In the following we analyze the superradiant behavior of our system by studying the evolution of the total number of particles in the lattice n T = j â † jâ j and the rate of emitted atoms in the last equality we have used that the total number operator kb † kb k + jâ † jâ j is a constant of motion.

4.2.1.
Hard-core bosons: Atomic superradiance. Let us start by analyzing the case of a lattice in an initial Mott phase having one atom per site in the collisional blockade regime, which is the analog of an ensemble of excited atoms [9]. As explained in the Introduction, Dicke predicted that superradiance should appear in this system as an enhancement of the emission rate at early times [10], although this was later proved to happen only if the effective number of interacting emitters exceeds some threshold value [12,13,14]: This is the superradiant phase transition. (dashed-dotted, green), and 10 (dotted, yellow) are considered. As expected from (a), the maximum of the rate is delayed above some critical ξ value. Note that both the rate and its derivative have been normalized to the values expected for independent emitters, which are 4M 3 Γ 2 0 and 2M 3 Γ 0 , respectively.
In our system, the number of interacting spins is governed by the parameter ξ (see figure 3), and the simplest way to show that the superradiant phase transition appears by varying it, is by evaluating the initial slope of the rate which can be written as ‡ This expression has a very suggestive form: The term corresponding to the rate associated to independent emitters is balanced by a collective contribution arising from the interactions between them. In figure 4a we show the dependence of this derivative with ξ for various values of the number of sites M 3 . It can be appreciated that there exists a critical value of ξ above which the sign of the derivative is reversed; hence, the rate increases at the initial time and we expect its maximum to be no longer at t = 0, which is a signature of superradiance.
The time evolution of the rate for a cubic lattice with M 3 = 27 sites is shown in figure 4b for different values of ξ. We can appreciate how above some critical ξ value the maximum rate of emission is delayed as expected. In order to find R(t) we have simulated the evolution equations for the coherences c jl = σ † jσ l and ‡ Note that the evolution equation of the expectation value of any operatorÔ can be written as and similarly the hard-core bosons in terms of the spin ladder operators. the populations s j = σ 3 j , which we close by using the semiclassical approximation σ † mσ 3 We have checked the validity of these semiclassical equations by comparing them with a direct simulation of the master equation for small number of sites in 1D and 2D geometries; except for small quantitative deviations, they offer the same results.

Non-interacting bosons: Harmonic oscillators superradiance.
Let us analyze now the case of having non-interacting bosons in the lattice, which is equivalent to a collection of harmonic oscillators as discussed before. In previous works on superradiance this system was studied in parallel to its atomic counterpart [13], and here we show how our system offers a physical realization of it. We will show that superradiant effects can be observed in the evolution of the total number of atoms in the lattice, both for initial Mott and superfluid phases §. § Let us note that a superfluid state with N excitations distributed over the entire lattice is more easily defined in the discrete Fourier-transform (DFT) basiŝ with q = (q x , q y , q z ) and q x,y,z = 0, 1, 2, ...M − 1, as the state having N excitations in the zeromomentum mode, that is, The evolution of the total number of atoms in the lattice is given by -see (22)- and hence depends only on the real part of the Markov couplings. Therefore, we restrict our analysis to the dissipative term D [ρ] of the master equation (16). By diagonalizing the real, symmetric collective decay rates with an orthogonal matrix S such that jl S pj γ j−l S ql =γ p δ pq , one can find a set of modes ĉ p = j S pjâj with definite decay properties. Then, it is completely straightforward to show that the total number of atoms can be written as a function of time as In general, γ j−l requires numerical diagonalization. However, in the limiting cases ξ ≪ 1 and ξ ≫ M, its spectrum becomes quite simple. Following the discussion after (20), in the ξ ≪ 1 limit γ j−l = Γ 0 δ jl is already diagonal and proportional to the identity. Hence, any orthogonal matrix S defines an equally suited set of modes all decaying with rate Γ 0 . Therefore, if the initial number of atoms in the lattice is N, this will evolve as irrespective of the particular initial state of the lattice (e.g., Mott or superfluid). The emission rate R = 2Γ 0 N exp (−2Γ 0 t), corresponds to the independent decay of the N atoms as expected in this regime having no interaction between the emitters. Let us consider now the opposite limit ξ ≫ M; in this case γ j−l = Γ 0 ∀ (j, l), and the dissipative term can be written in terms of the symmetrical DFT mode only as (25). Hence, the DFT basis diagonalizes the problem, and shows that all the modes have zero decay rate except the symmetrical one, which has an enhanced rate proportional to the number of emitters. Therefore, starting with a superfluid state, the N initial atoms will decay exponentially with initial rate NM 3 Γ 0 , that is On the other hand, if the initial state corresponds to a Mott phase, most of the atoms will remain in the lattice, as only the component which projects onto the symmetric mode will be emitted; concretely, from (27) and (25) the number of atoms in the lattice will evolve for this particular initial state as Hence, according to this picture, superradiant collective effects can be observed in our system by two different means. Calling t 0 the time needed to radiate the atoms in the absence of collective effects, one could start with a superfluid phase and measure this 'evaporation time' as a function of ξ; this should go from t 0 for ξ ≪ 1, to a much shorter time t 0 /M 3 for ξ ≫ M (see figure 5a). Alternatively, one could start with a Mott phase, and measure the number of atoms left in the lattice in the steady state as a function of ξ; in this case, it should go from n T,steady = 0 after a time t 0 for ξ ≪ 1, to n T,steady = N − N/M 3 after a time t 0 /M 3 for ξ ≫ M (see figure 5b).
In order to find the evolution of n T we have simulated the equations satisfied by the coherences c jl = â † jâ l , which read -see (22)- Note that in this case the equations are closed without the need of a semiclassical approximation, and hence they are exact. Note that they reproduce the analytic evolution of n T as given by (28), (30), and (29) in the corresponding limits (see figure  5).
Our results connect directly to those found by Agarwal some decades ago [13]. Working with a Dicke-like model, he showed that if all the oscillators start in the same coherent state |α , the initial number of excitations, which in that case is given by N = M 3 |α| 2 , decays following (29). This is not a coincidence, but rather a consequence that, if N is large enough, a multi-coherent state of that kind is a good approximation of a superfluid state with that number of excitations. He also predicted that if the oscillators start in a number state, most of the excitations would remain in the steady state as follows from (30).

Directional superradiance
The last phenomenon connected to quantum-optics that we want to discuss is how allowing k L = 0 the direction of emission of atoms can be controlled.
This follows from the fact that by tuning the laser wave vector such that k L = k 0 while working in the∆ > 0 regime, the master equation of our system (12) has exactly the same form as that of [18] (see also [17]), where the properties of light emitted by regular arrays of atoms were analyzed in terms of the rate and the direction of emission. Hence, we could expect the same phenomena to appear in our system, with the difference that in our case the emitted particles are atoms and not photons. Let us explain the results found in [18] and translate them to our system.
We are particularly interested in the analysis performed in [18] of the emission properties of an initial symmetric Dicke state (a symmetrical spin wave in its notation). It was shown that the average number of photons emitted in the direction of the vector u Ω = (sin θ cos φ, sin θ sin φ, cos θ), θ and φ being the polar and azimuthal spherical angles, respectively, can be written as where the total emission rate Γ is found from the normalization condition dΩI (Ω) = 1 (Ω is the solid angle). This expression was deduced in [18] in the M ≫ 1 limit by diagonalizing the master equation with the use of periodic boundary conditions. The limit Γ ≪ c/L was considered (c is the speed of light and L = Md 0 is the length of the lattice in a given direction), that is, it was assumed that the emission time is larger than the time taken by a photon to leave the sample. Following the same steps, but taking into account the different dispersion relation of the emitted particles in our system, the same expression (32) can be proved but now under the condition Γ ≪ v 0 /L, where v 0 = k 0 /m is the 'resonant' speed of the emitted atoms. We will explain latter what this condition means in terms of the system parameters.
In the optical case treated in [18], the angular distribution (32) applies strictly only to the single excitation sector. In our case, this result is specially relevant for the case of non-interacting atoms trapped in the optical lattice in an initial superfluid state: since a superfluid state can be described as a state of N atoms in the completely symmetric state, see (25), the results of [18] imply here that those N trapped atoms will emit N free atoms with an angular distribution given by (32).
The sums in (32) can be performed analytically, leading to the result where we have denoted byk L = k L /k L the unit vector in the 'laser' direction. This function has a diffraction maximum whenever the denominator vanishes, that is, for vectors m ∈ Z 3 such that u Ω =k L + 2πξm. In order to analyze directionality it is better to rewrite this condition as πξ|m| = |k L ·m|, wherem = m/|m|. Then we see that for πξ > 1 only m = 0 can satisfy this condition, and hence the atoms are mostly emitted in the laser direction. This is the directional regime. The emission rate Γ and the angular width of the atomic beam ∆θ must be evaluated numerically in general. However, approximating the diffraction peaks by Gaussian functions, a limit which is justified if the peaks are narrow enough (ξ/M ≪ 1), we can get the following expression for them [18] Γ = χΓ 0 with χ = π 3/2 Mξ 2 , and ∆θ = ξ/M.
Hence we see that in addition to directionality, the rate is enhanced by a factor χ. Let us finally explain what the initial approximation Γ ≪ v 0 /L means in terms of the system parameters. Using (34), this can be written as which can be always satisfied by a proper tuning of the two-photon Rabi frequency. This analysis shows that our scheme can be used to observe directional, superradiant emission by starting with a superfluid phase and tuning the laser wave vector so that k L = k 0 .

Considerations about dimensionality
So far we have allowed atoms in state |b to move freely in 3D. However, the free motion of these atoms can be restricted to 2D or 1D by coupling them to an additional 1D or 2D deep optical lattice, respectively. Under these conditions the momentum of the free atoms has components only along the non-trapped directions, say k 2D = k xx + k yŷ and k 1D = k xx ; from the point of view of the atoms in state |a , the dimensionality of the reservoir they are interacting through is therefore reduced to 2D or 1D. Hence, only sites laying in the same z = const plane (for the 2D reservoir) or (y, z) = const tube (for the 1D reservoir) would be connected through a two-point correlation function given by (11), but replacing G(τ ) by the correlation functions G 1D coincides with the correlation function of atoms embedded in isotropic band-gap materials [21,22], while G 2D has no quantum-optical analog to our knowledge. Hence, with an additional trapping of the free atoms, our scheme could be useful to study low-dimensional reservoirs described by these correlation functions.

Conclusions
To conclude, let us briefly summarize what we have shown in this article.
We have studied a model consisting of atoms trapped in an optical lattice (with trapping frequency ω 0 , lattice period d 0 , and cubic geometry with M 3 sites) which are coherently driven to a non-trapped state through a Raman scheme having detuning ∆ respect to the atomic transition, relative wave vector between the Raman lasers k L , and two-photon Rabi frequency Ω. The limits of non-interacting bosons and hardcore bosons -where either one or zero atoms can be in a given lattice site-have been considered for the trapped atoms, while interactions involving the free atoms have been neglected.
These are the main results that we have found (numbering follows the sections of the article): (2) We have first deduced the Hamiltonian of the system, showing that it is equivalent to that of a collection of harmonic oscillators or two-level systems for noninteracting atoms and hard-core bosons, respectively, interacting with a common radiation field consisting of massive particles.
(3) We have then studied the emission of free atoms when the lattice sites emit independently. We have found a phase transition at∆ = ∆ − 4Ω 2 /ω 0 = 0 from a regime (∆ < 0) where the trapped and free atoms form a bound state and there is a nonzero steady state population left in the lattice, to a radiative regime (∆ > 0) in which all the atoms leave the lattice eventually. We have identified the Markovian and non-Markovian regimes of the emission, showing that the condition for Markovianity is |∆| ≫ 8πΩ 4 /ω 3 0 . This behavior is analogous to the one predicted in [21,22] for the spontaneous emission of atoms embedded in photonic band-gap materials.
(4) As for collective effects, we have studied them by using a reduced master equation for the lattice within the Born-Markov approximation, explaining under which conditions this may hold in our system. This master equation has shown that the number of interacting sites is controlled by a single parameter, namely ξ = 1/d 0 k 0 with k 0 = X −1 0 2|∆|/ω 0 . Now let us summarize the main collective effects that we have predicted to appear.
(4.1) For∆ < 0 and k L = 0, the evolution of the lattice follows an effective Hamiltonian: Extended hopping for non-interacting atoms and extended ferromagnetic interactions for hard-core bosons. The coupling between lattice sites has a Yukawa form, but approaches a Coulomb form as ξ increases. This could be helpful for the simulation of condensed-matter phenomena requiring extended interactions.
(4.2) For∆ > 0 and k L = 0, the master equation has a collective dissipation term in addition to the Hamiltonian part. In the ξ ≪ 1 limit this dissipation term corresponds to the lattice sites emitting independently with the same rate, while in the ξ ≫ M limit it takes a Dicke-like form, and therefore superradiant effects are expected to appear. Hence, we have argued that the superradiant phase-transition could be observed in our system by varying the ξ parameter. Indeed, we have shown that both atomic and harmonic oscillator superradiance can be observed: (4.2.1) When starting with a Mott phase having one atom per site in the collisional blockade regime (which is the equivalent of a collection of excited atoms), there exists a critical value of ξ above which the emission rate is enhanced at the initial stage of the emission, and its maximum occurs for t = 0. This delay of the maximum rate was suggested by Dicke in his seminal work [10], and was shown to be a signature of the superradiant phase-transition for atomic samples [12,13,14].
(4.2.2) If the lattice atoms are non-interacting, there are two equivalent ways of observing superradiant behavior. One can start with a superfluid phase, and measure the time needed to radiate the lattice atoms, which should be reduced by a factor M 3 when going from the ξ ≪ 1 to the ξ ≫ M limit. Alternatively, one can start with a Mott phase having N atoms in the lattice, and measure the number of trapped atoms left in the steady state; this should go from zero in the ξ ≪ 1 limit, to N − N/M 3 in the ξ ≫ M limit. This behavior is consistent with the results found in [13] for the spontaneous emission of a collection of harmonic oscillators, and hence our scheme provides a physical realization of this system. (4. 3) The last phenomenon linked to quantum-optics that we have shown is that superradiant, directional emission can be obtained in our system. In particular, we have proved that by starting from a superfluid phase distributed over enough lattice sites (M ≫ 1), matching k L to k 0 , and working in the πξ > 1 regime, the atoms are emitted mostly within the k L direction. If the condition ξ/M ≪ 1 is also satisfied, the beam has been shown to have angular width ∆θ = ξ/M and a rate enhanced by a factor χ = π 3/2 Mξ 2 . These results are in correspondence with those of [17,18], where the emission of photons by regular arrays of atoms was studied.
(5) We have finally discussed what would change if the motion of the free atoms is restricted to 2D or 1D by using an additional optical lattice. We have shown that in this conditions the correlation function of the reservoir through which the lattice sites are interacting takes the form G(τ ) ∝ (1 + iω 0 τ /2) d/2 , where d is the dimensionality (d = 1 and 2 for 1D and 2D, respectively). For 1D this corresponds to the correlation function of the radiation field in an isotropic photonic band-gap material [21,22], while for 2D it has no quantum-optical counterpart.
All the proposed scenarios can be reached with currently available technology, and hence we believe that this work provides a solid base for considering optical lattices as quantum simulators of quantum-optical systems.
of independent lattice sites.
Making the Laplace transform of (7) the solution for the amplitude A in Laplace space is found as On the other hand, the Laplace transform of the correlation function (8) is or to first order in (∆ + is) /ω 0 Note that as we assume that s/ω 0 ≪ 1, this transform cannot describe properly times below ω −1 0 , and this is why sometimes this limit receives the name 'strong confinement', as it needs ω 0 → ∞ or otherwise the initial steps of the evolution are lost in the description. Unlike in the general case where the transform of the correlation function has a really complicated form (A.4), in this limit it is possible to make the inverse Laplace transform ofÃ (s). Hence, in the following we consider the two cases separately, as the first case requires special attention in order to extract results.

Appendix A.1. Evolution of the population in the strong confinement limit
The Laplace transform of A can be inverted as [33]   being R j the residues of the kernel K (s) =Ã (s + i∆) exp (st) at its poles [33]. The integrals over the paths C u , C d , and C c are zero when r → 0 and R → ∞.
Making the variable change s = e iπ x and s = e −iπ x in the integrals along R u and R d , respectively, we get which is not analytical, but can be performed numerically very efficiently.
On the other hand, the poles of K (s) are easily evaluated by making the change s = ix 2 , hence turning its denominator into x 2 + 2 √ παx +∆ whose roots are therefore the poles are simple and read s ± = ix 2 ± . In order to evaluate the residues, it is convenient to write the Kernel of (A.6) as from which the residues are easily found as The problem now is that the phase of s ± must be defined in the interval ]−π, π[ (remember the branch cut in figure A1), and hence we have to analyze their phases as a function of the phases of x ± , which we take in the interval [0, 2π] as a convention.
In terms of the parameters∆ and α, we can define 3 regions: • In the region∆ < 0 we have ϕ + = 0 and ϕ − = π. Hence, we have to write s + = e iπ/2 |x + | 2 and s − = e iπ/2 |x − | 2 e 2iπ × e −2iπ , (A.12) that is, we have to decrease the phase of s − by 2π to keep it within the interval ]−π, π[; with s ± written in this way, the residues read • In the region 0 <∆ < πα 2 , we have ϕ ± = π, and hence, we have to decrease the phase of both s ± by 2π to keep them within the interval ]−π, π[; therefore both residues are zero R ± = 0.

Appendix A.2. Steady state population in the general case
From the last discussion it should be clear that the long term population is dictated solely by the residues of K (s) =Ã (s + i∆) exp (st). Moreover, only the residues corresponding to pure imaginary poles will remain in the t → ∞ limit, and hence the steady state population will be given by (A. 16) Then, in order to know the steady state population in the finite trap case, we just need to find the pure imaginary poles of K (s) with the transform of the correlation function given by (A.4). By writing s = ix with x ∈ R, this amounts for solving the equation This equation cannot be solved analytically, but it can be efficiently done numerically. On the other hand, once we know the solutions x j of this equation, we can use them to evaluate their associated residues as [33] R j = e ix j t d ds s + i∆ +G (s + i∆) This is the method that we have used to compute the results shown in figure 2d.
Appendix B. Evaluation of the Markov couplings (13) In this appendix we show that in the strong confinement regime the Markov couplings have the form (13). After the Born-Markov approximation is carried, the effective coupling between different sites has the form ¶ In the limit k L ≪ X −1 0 (in which we can remove k L from g k -see (4)-), taking the continuous limit for the momentum, and performing the angular part of the integral over momenta by assuming that r j−l is within the z-axis, this expression can be reduced to Γ j−l = 8X 0 Ω 2 iπ 1/2 ω 0 e −ik L ·r j−l The second integral is not analytical in general. Nevertheless, in the strong confinement regime both |∆|/ω 0 and X 0 /d 0 are small, and hence the Gaussian function in the Kernel can be approximated by one, so that the remaining integral reads Introducing these integrals in (B.2) and talking the ε → 0 limit we arrive to expression (13) for the Markov couplings. ¶ Note that ε is more than just a mathematical artefact for regularizing the integral. It accounts for the decoherence channels not taken into account in the model, which make the correlation between lattice sites decay as G j−l (τ ) exp(−ετ ). Even if these channels appear at a time scale larger than the relevant times (ε ≪ 1), they remove the zeros of the kernel's denominator, showing that its singularities are not physical.