1.1 Optical Lattices

Optical lattices are periodic potentials created by light-matter interactions. When an atom interacts with an electromagnetic field, the energy of its internal states depends on the light intensity. Therefore, a spatially dependent intensity induces a spatially dependent potential energy. Optical lattices have been widely used in atomic physics as a way to trap and cool atoms.

Neutral atoms interact with light in both dissipative and conservative way.

The conservative interaction comes from the interaction of the light field with the induced dipole moment of the atom which causes a shift in the potential energy called AC Stark shift.

The dissipative interaction comes due to the absorption of photons followed by spontaneous emission. Laser cooling techniques make use of this dissipative process.

Single particle physics

Bloch functions

One of the most important characteristics of a periodic potential is the emergence of a band structure. Consider a one dimensional particle described by the Hamiltonian H = + V _lat (x), where V _lat (x) = V _lat (x + a). Bloch's theorem (^{Ashcroft & Mermin, 1976}; ^{Ziman, 1964}) states that the eigenstates Φ _q ⁽ⁿ⁾ (x) can be chosen to have the form of a plane wave times a function with the periodicity of the potential:

Using this ansatz into the Schrodinger equation, H Φ _q ⁽ⁿ⁾ (x)= E _q ⁽ⁿ⁾ (x) Φ _q ⁽ⁿ⁾ (x), yields an equation for u _q ⁽ⁿ⁾ (x) given by:

Bloch's theorem introduces a wave vector q. The quantity q should be viewed as a quantum number characteristic of the translational symmetry of the periodic potential, just as the momentum is a quantum number characteristic of the full translational symmetry of the free space. Even though is not the same, it turns out that hq plays the same fundamental role in the dynamics in a periodic potential than the momentum in the absence of the lattice. To emphasize this similarity ħq is called the quasimomentum or crystal momentum . In general the wave vector q is confined to the first Brilloiun zone, i.e. -Π/α < q ≤ Π/α.

The index n appears in Bloch's theorem because for a given q there are many solutions to the Schrodinger equation. Eq. (10) can be seen as a set of eigenvalue problems in a fixed volume, one eigenvalue problem for each q. Therefore, each of them, has an infinite family of solutions with a discretely spaced spectrum with modes labelled by the band index n. On the other hand, because the wave vector q appears only as a parameter in Eq. (10), the energy levels for a fixed n has to vary continuously as q varies. The description of energy levels in a periodic potentials in terms of a family of continuous functions Eq' each with the periodicity of a reciprocal lattice vector, 2Π/α, is referred as the band structure. The lattice potential is generally expressed in units of the atom recoil energy, E _R = ħ ² k ² /(2m), with k = Π/α.

Fig. 2 shows the band structure of a sinusoidal potential. For V₀ = 0, the particles are free so the spectrum is quadratic in q. As the potential is increased the band structure appears.The band gap increases and the band width decreases with increasing lattice depth.

In the absence of a lattice the eigenfunctions of the free system are plane waves (see Fig. 3). As the lattice depth is increased the barrier height between adjacent lattices sites increases and the eigenstates of the system tend to get localized at each lattice site (regions around the potential minima).

1.1.4 Wannier orbitais

Wannier orbitals are a set of orthogonal normalized wave functions that fully describe particles in a band and that are localized at the lattice sites. They are defined as:

where the sum is over the first Brillouin zone, M is the total number of lattice sites and x _i is the position of the i ^th lattice site. Wannier orbitals are thus a unitary transformation of the Bloch functions and are formally an equivalent representation to describe the periodic system.

The actual form of a Wannier function may be understood if we Eq.(8) to be approximately the same for all Bloch states in a band. Under this approximation the Wannier function at the origin can be shown to be

They look likes u ⁽ⁿ⁾ (x) at the site center, but it spreads out with gradually decreasing oscillations with are needed to ensure orhogonality. In Fig. 4 we show a Wannier orbital centered at the origin site.

1.1.5 Tight binding approximation

The tight-binding approximation deals with the case in which the overlap between Wannier orbitals at different sites is enough to require corrections to the picture of isolated particles but not too much as to render the picture of localized wave functions completely irrelevant.

In this regime to a very good approximation one can only take into account overlap between Wannier orbitals in nearest neighbor sites. Furthermore, if initially the atoms are prepared in the lowest band, the dynamics can be represented quite satisfactorily by using Wannier functions of only the lowest band. The tight-binding model is commonly used to solve the problem of a particle in a periodic potential when also an external potential V(x), not strong enough or sharp enough to induce interband transitions, is applied.

By expanding the wave function in lowest band Wannier orbitals,

we get the following equations of motion

with

J is the tunneling matrix element between nearest neighboring lattice sites and ε ₀ is the unperturbed on site energy shift. Eq. (14) is known as the discrete Schrodinger equation (DSE) or tight binding Schrodinger equation.

For the special case of a sinusoidal potential V _lat (x)= V ₀ sin² (kx) for which Mathieu functions are the exact solutions an analytic expression for J can be obtained:

with α =1.39666, β =1.051 and ү =2.12104 and Ṽ ₀ = V ₀ /E _R . Note that J decreases exponentially with the lattice depth.

1.1.6 Semiclassical dynamics

Solving Eq. (14) is presumably not possible for an arbitrary V(x_i). However, we get general insight to the nature of solutions by an application of the correspondence principle. It is well known that wave-packet solutions of the Schrodinger equation behave like classical particles obeying the equations of motion derived from the same classical Hamiltonian. The classical Hamilton equations are

Using the correspondence principle and the resemblance of the crystal momentum or quasi-momentum to real momentum, the semiclassical equations of a wave packet in the first band of a lattice can be written as (^{Ashcroft & Mermin, 1976}; ^{Ziman, 1964})

The semiclassical equations of motion describe how the position and wave vector of a particle evolve in the presence of an external potential entirely in terms of the band structure of the lattice. If we compare the acceleration predicted by the model with the conventional newtonian equation, m = -dV(x)/dx, we can associate an effective mass induced by the presence of the lattice, m*. This is given by

If no external potential is applied, V (x) = 0, Bloch waves are the solutions of tight binding equation: with M is the total number of lattice sites. If also periodic boundary conditions are assumed, the quasimomentum q is restricted to be an integer multiple of . The lowest energy band dispersion relation in this case is given by

From the above equation it is possible to see the connection between J and the band width

The effective mass becomes and close to the bottom of the band q → 0, m*→ . Since the tunneling decreases exponentially with the lattice depth, the effective mass grows exponentially. In the presence of interactions, the large effective mass manifest itself in a substantial enhancement of the interaction to kinetic energy ratio, in comparison to the free particle case. This is the reason why atoms in optical lattice can reach the strongly interaction regime.

1.2.1 Superfluid Mott insulator transition

At zero temperature the physics described by the Bose Hubbard Hamiltonian can be divided in two different regimes. The interaction dominated regime when J is much smaller than U, and the system is in the Mott insulator phase, and the kinetic energy dominated regime, when tunneling overwhelms the repulsion and the system exhibits superfluid properties. The transition between the two regimes is a consequence of the competition between the kinetic energy which tries to delocalize the particles and the interaction energy which tries to localize them and to make the number fluctuations small.

In the superfluid regime, the kinetic energy term dominates the Hamiltonian and the system behaves as a weakly interacting Bose gas. Quantum fluctuations can be neglected and the system can be described by a macroscopic wave function. For the translationally invariant case the ground state consists of all atoms in the 0-quasimomentum mode

where N is the total number of atoms and is the anhihilation operator of an atom with quasimomentum. Remember α is the lattice spacing.

Since the many body state is almost a product over identical single particle wave functions the system has a well defined phase and large number fluctuations ~ . The macroscopic occupation of a single mode implies that the commutator of the annihilation and creation operators that create and destroy a particle in the condensate mode can be neglected with respect to the total number of atoms in the mode. Therefore, to a good approximation the field operator can be replaced by a c-number,

This procedure can be interpreted as giving to the field operator a non zero average, 〈 ₀〉 = , and thus a well defined phase to system. Because the original Hamiltonian is invariant under global phase transformation, this definition of BEC corresponds to an spontaneous symmetry breaking (^{Bogoliubov, 1947}; Dalfovo, Giorgini, & Stringari, 1999; P. W. Anderson, 1966; Leggett, 2001; Leggett & Sols, 1991). The wave function associated with the condensate mode, in this case the zero quasi-momentum eigenstate is often called "condensate wave function". It characterizes the off-diagonal long range order present in the one particle density matrix (^{Ginzburg & Landau, 1950}; Penrose & On-sager, 1956).

Strictly speaking, in finite size systems with well defined number of particles neither the concept of broken symmetry, nor the one of off-diagonal long range order can be applied. However the condensate wave function can still be determined by diagonalization of the one particle density matrix. The eigenstate with larger eigenvalue corresponds to the condensate wave function.

As interaction increases the average kinetic energy required for an atom to hop from one site to the next becomes insufficient to overcome the energy cost. Atoms tend to get localize at individual lattice sites and number fluctuations are reduced. Deep in the Mott insulator phase the ground state of the system instead consist of localized atomic wave functions with a fixed number of atoms per site.

where n ₀ = N/M is the filling factor or mean number of particles per site. We have assumed that the number of particles is commensurate with lattice sites, i.e. that n ₀ is an integer.

The lowest lying excitations that conserve particle number are particle-hole excitations (adding plus removing a particle from the system). The Mott phase is characterized by the existence of an energy gap. The gap is determined by the energy necessary to create one particle-hole pair.

1.2.2 Mean Field Phase diagram

The zero temperature Bose-Hubbard phase diagram (M. P. A. Fisher et al., 1989; Krauth et al., 1992) of a translational invariant system Vj = 0 exhibits lobe-like Mott insulating phases in the J-μ plane, see Fig. 5. Each Mott lobe is characterized by having a fixed integer density. Inside these lobes the compressibility, ∂ρ/∂μ with ρ the average density of the system, vanishes.

The physics behind this diagram can be understood as follows: If we start at some point in the Mott insulating phase and increase μ while keeping J fixed, there is going to be a point at which the kinetic energy of adding an extra particle and letting it hop around will balance the interaction energy cost. With an extra particle free to move around the lattice phase coherence is recovered and the system enters the superfluid regime. Similarly, by decreasing μ from a point in the Mott phase, at some point eventually it will be energetically favorable to remove one particle from the system. The extra mobile hole will induce also phase coherence and the system will condense in a superfluid state. Since the kinetic energy of the system increases with J, the width of the lobes decreases with J. The distance in the μ direction at fixed J between the upper and lower part of the lobe is the energy gap. At J = 0 the gap is just equal to U. Also at J = 0 the intersection points between the lobes with density p = n ₀ and the ones with density ρ = n ₀ + 1 are degenerated. Because there is not energy barrier to the addition of an extra particle at this generated points the system remain superfluid. Mott insulator phases occur only at integer densities, non-integer densities contours lie entirely in the superfluid phase because there is always an extra particle that can hop without energy cost.

The phase diagram includes two different types of phase transition. One type takes place at any generic point of the phase boundary, and it is driven by the energy cost to add or subtract small numbers of particles to the incompressible Mott state as explained above. On the other hand, the other type only occurs at fixed integer density and takes place at the tip of the lobes. This transition is driven at fixed density by decreasing U/J and enabling the bosons to overcome the on site repulsion. The two kinds of phase transition belong to different universality classes. In the generic one, the parameter equivalent to the reduced temperature δ = T - T _c for finite temperature transitions that measures the distance from the transition is δ ~ μ- μ _c , with μ _c the chemical potential at the phase boundary. For the special fixed density, on the other hand, one must take δ ~ (J/U ) - (J/U)_c. As shown in M. P. A. Fisher et al. (1989) provided δ ~ μ - μ _c , the compressibility has a singular part scaling near the transition as δ ^-α , α > 0. This scaling does not apply at the special fixed density Mott-Superfluid transition for which δ ~ (J/U ) - (J/U ) _c and differentiation with respect to delta becomes inequivalent to differentiation with respect to the chemical potential.

Most of the features of the phase diagram discussed above can be verified by simple calculations done using a decoupling mean field approximation. Under the decoupling approximation one introduces a real superfluid order parameter ψ = = , and constructs a consistent mean-field theory by substituting the kinetic energy term by

into Eq. 25. This yields

where z = 2d is the number of nearest-neighbor sites. This Hamiltonian is diagonal with respect to the site index j, so we can use an effective onsite Hamiltonian. We will therefore drop the subscript j in the following.

By treating ψ as the perturbation parameter one can use perturbation theory to calculate the mean field energy. Assuming commensurate filling N/M = n ₀ = 1,2,... , to zero order in ψ the ground state of the system will have energy E ⁽⁰⁾ _no and exactly n ₀ particles per site (Eq. 32). In this unperturbed occupation number basis the odd powers of the expansion of the energy in ψ will always be zero. Up to four order the ground state energy becomes

If we now follow the usual Landau procedure for second order phase transitions and minimize Eq. 34 with respect to the superfluid order parameter ψ, we find that ψ ≠ 0 when E ⁽²⁾ < 0 and that ψ = 0 when E ⁽²⁾ < 0. This means that E ⁽²⁾ = 0 signifies the boundary between the superfluid and the insulator phases. The solution E ⁽²⁾ = 0 occurs for

By equating μ ⁺ = μ- one finds the point of smallest U for each lobe which corresponds to the critical value at the tip:

1.2.3 The Bose Hubbard model and Superfluidity

The concept of superfluidity is closely related to the existence of a condensate in the interacting many-body system. Formally, the one-body density matrix ρ ⁽¹⁾ (x, x') has to have exactly one macroscopic eigenvalue whichs defines the number of particles in the condensate; the corresponding eigenvector describes the condensate wave function Φ₀ (x) = . A spatially varying condensate phase, Θ (x), is associated with a velocity field for the condensate by

This irrotational velocity field is identified with the velocity of the superfluid flow, v _s (x) = v₀ (x) (Leggett, 1999; ^{Lifshitz & Pitaevskii, 1980}) and enables us to derive an expression for the superfluid fraction, f _s . Consider a system with a finite linear dimension, L along the e₁-direction, and a ground-state energy, E₀, calculated with periodic boundary conditions. Now we impose a linear phase variation, Θ (x) = θx ₁ /L with a total twist angle θ over the length of the system in the e₁-direction. The resulting ground-state energy, E _θ will depend on the phase twist. For very small twist angles, θ«Π, the energy difference, E ₉ - E₀, can be attributed to the kinetic energy, T _s , of the superflow generated by the phase gradient. Thus,

where m is the mass of a single particle and N is the total number of particles so that mN f _s is the total mass of the superfluid component. Replacing the superfluid velocity, v_s with the phase gradient according to Eq. (37) leads to a fundamental relation for the superfluid fraction.

where the second equality applies to a one dimensional lattice system on which a linear phase variation has been imposed. Here the distance between sites is α, the phase variation over this distance is ∆θ, and the number of sites is M. In this case, J = ħ ² /(2ma²).

Technically the phase variation can be imposed through so-called twisted boundary conditions (M. E. Fisher, Barber,& Jasnow, 1973). In the context of the discrete Bose-Hubbard model it is, however, more convenient to map the phase variation by means of a unitary transformation onto the Hamiltonian (Roth & Burnett, 2003; Rey et al., 2003). The resulting "twisted" Hamiltonian in this case

exhibits additional phase factors e^±i∆θ - the so-called Peierls phase factors - in the hopping term. We calculate the change in energy E _θ - E _θ under the assumption that the phase change ∆θ is small so that we can write:

Using this expansion the twisted Hamiltonian (40) takes the following form:

where we retain terms up to second order in ∆θ. The current operator (Note that the physical current is given by this expression multiplied by ) and the hopping operator are given by:

The change in the energy E _θ - E ₀ due to the imposed phase twist can now be evaluated in second order perturbation theory

The first order contribution to the energy change is proportional to the expectation value of the hopping operator

Here │ψ₀〉 is the ground state of the original Bose-Hubbard Hamiltonian. The second order term is related to the matrix elements of the current operator involving the excited states │Φ₀〉 (v = i, 2,...) of the original Hamiltonian

Thus we obtain for the energy change up to second order in ∆θ

The quantity D, defined above, is formally equivalent to the Drude weight used to specify the DC conductivity of charged fermionic systems (Shastry & Sutherland, i990). The superfluid fraction is then given by the contribution of both the first and second order term:

In general both, the first and the second order term contribute. The superfluid fraction is 1 for the noninteracting system and decreases slowly for small U/J. In the region of the Mott transition f _s drops rapidly and goes to zero in the Mott-insulator phase. The rapid decrease of the total superfluid fraction is largely due to the second order contribution.

1.2.4 Bogoliuov treatment of the superfluid to Mott transition

As mentioned before the Bose Hubbard Hamiltonian is the simples model that describes interacting bosons in a lattice. Even though the Hamiltonian looks simple in general it is very difficult to deal with exactly and perturbative or approximated treatments are in order. A systematic approach to deal with the weakly interacting regime uses the ideas proposed by ^{Bogoliubov in 1947} (^{Bogoliubov, 1947}). The Bogoliubov approximation starts by treating the field operator as a c-number plus a small fluctuating term. The c-number describes the condensate or the coherent part of the matter field and the fluctuating term accounts for quantum correlations. Since in a weakly interacting gas, quantum fluctuations are small, one can neglect third and higher terms in the fluctuating field and derive an effective quadratic Hamiltonian in the fluctuating field which can be exactly diagonalized.

In the following, we will apply the Bogoliubov approximation to the Bose Hubbard Hamiltonian and derive the correspondent Bogoliubov-de Genes(BdG) equations as developed in Refs. (van Oosten et al., 2001; Roth & Burnett, 2003; Rey et al., 2003). We will restrict our analysis to zero temperature T = 0.

In the very weakly interacting regime, to a good approximation, the field operator can be written in terms of a c-number plus a fluctuation operator:

Replacing this expression for â _n in the Bose Hubbard Hamiltonian leads to :

with

Where < n, n > restricts the sum over m to nearest neighbors and h.c. stands for the hermitian adjoint. The terms of the Hamiltonian have been grouped in equations according to the number of non condensate operators which they contain.

The first step is the minimization of the energy functional H ₀ . This requires the condensate amplitudes ψ _n to be a solution of the following non-linear equation,

This equation is known as the discrete nonlinear Schrödinger equation (DNLSE). It is a discrete version of the so called Gross-Pitaevskii equation (GPE) (E. P. Gross, 1961; Pitaevskii, 2003). It can be also obtained by spanning the GPE equation in a Wannier orbital basis and then using the tight binding approximation (see Eq. 14.) If ψ _n is a solution of Eq. 57 it means that the linear Hamiltonian ₁ vanishes.

The quadratic Hamiltonian ₂ allows the leading order effects of the non-condensate atoms to be taken into account. ₂ can be diagonalized by finding a basis of so-called quasiparticle states which do not interact which each other. Mathematically this involves a linear canonical transformation of the single-particle creation and annihilation operators and â _n into quasiparticle operators The transformation is what is known as Bogoliubov transformation and is given by

In general the spectrum of fluctuations includes a zero mode. This mode is the Goldstone boson associated with the breaking of global phase invariance by the condensate. The zero mode is essentially non perturbative and it introduces an artificial infrared divergence in low dimensional models. For this reason quadratic approximations are actually improved if the contribution from this mode is neglected all together (Al Khawaja, Andersen, Proukakis, & Stoof, 2002). A different way to deal with the zero mode has been proposed by Refs. (^{Castin & Dum, 1998}; Gardiner, 1997; ^{Morgan, 2000a}, 2000b). For simplicity we are going to ignore zero mode fluctuations and restrict the fluctuation operators to act only on the excited states.

The Bogoliubov transformation Eq. 58 is required to be canonical, which means that it preserves the commutation relations and leads to bosonic quasiparticles. To satisfy this, the amplitudes {u ^s _n , v ^s _n } are constrained by the conditions:

The necessary and sufficient conditions that the quasiparticle amplitudes have to fulfill to diagonalize the Hamiltonian are provided by the so called Bogoliubov-de Gennes (BdG) equations

with u ^s = (u ^s ₁ , u ^s ₂ ,...) , v ^s = (v ^s ₁ , v ^s ₂ ,...) and ψ =(ψ₁ ψ₂,.-). The matrices and are given by

with δ _n1 the Kronecker delta which is one if n = k and zero otherwise. If the BdG equations are satisfied, the quadratic Hamiltonian takes the form

The quasi-particle energies ħω ^B _s come in pairs: if ħω ^B _s is a solution for the amplitudes (u^s, v ^s ) then -hmf is also a solution for the amplitudes (v ^s *,u ^s *). The solution with zero energy is always a solution and in this case the amplitudes must be proportional to the condensate (v⁰,u⁰) ∝ (ψ, ψ). We explicitly exclude the zero mode solution to guarantee that the excitations are orthogonal to the condensate.

To have a complete description of the quadratic Hamiltonian, we define the one-body density fluctuation matrix , ρ _nm, and the anomalous average matrix, m _nm , as:

The averages denote an ordinary quantum expectation value. In the zero temperature situation, the average is over the ground quasiparticle state given by â _s │0〉 = 0 and the p _nm and m _nm ^can be written in terms of quasiparticle amplitudes as:

The quantity ρ _nn represents the non-condensate population at position n or depletion at position n. If the average number of atoms is N, the depletion and condensate atoms are related by the equation:

To understand many-body effects included in the quadratic approximations we study the case where no external confinement is present, V _n = 0. We assume a d dimensional separable square optical lattice with equal tunneling matrix element J in all directions and periodic boundary conditions. The total number of wells is M.

Due to the translational symmetry of the system the condensate amplitudes are constant over the lattice, ψ _n = . The quasiparticle modes have a plane wave character and therefore can be related to quasimomentum modes:

Here the vector q denotes the quasi-momentum whose components assume discrete values which are integer multiples of m with a is the lattice spacing. The amplitudes u _q and v _q

must satisfy the condition │u_q│² -│v_q│² = 1 and can all be chosen to be real and to depend only on the modulus of the wave vector (u_q = u_-q, v_q = v _-q ).

In the translationally invariant system the DNLSE reduces to

where z is the number of nearest neighbors z = 2d.

The BdG equations become the following 2 x 2 eigenvalue problem

with

Here we have introduced the definition ε _q =

The quasiparticle energies ħω _q and modes are found by diagonalizing Eq.(73):

And

with n the total density, n = N/M. The constrain that fixes the number of particles can be written as

As opposed to the free particle system where for high momentum always the single particle energy (which grows as q ² ) is dominant, in the presence of the lattice the single particle excitations are always bounded by 4Jd. Therefore, in the regime Un _c /J > 1 the interaction term dominates for all quasimomentum and so ħω _q ~ . On the other hand, in the weakly interacting regime, the most important contribution to the sum comes from the low lying modes. For these modes is also true that ħω _q ~ . Thus, to a good approximation the condensate fraction can be written as:

with

In Eq.(82) the term is a finite size effect term which is important to keep for low density systems. α is a dimensionless quantity which depends only on the dimensionality of the system. In the thermodynamics limit, M → ∞ the sum can be replaced by an integral and we find

The infrared divergency in the one dimensional thermodynamic limit is a consequence of the importance of long wave length correlations in low dimensional systems. For finite one dimensional systems however α (1, M) has a finite value.

Because Eqs. ((74)-(75)) are completely determined if n _c is known, by solving Eq. (80) we obtain all necessary information. The solution of the algebraic equation is:

Eq. (86) tells us that for very weak interactions almost all the atoms are in the condensate. As interaction increases the condensate fraction decreases but it only vanishes when U/J → ∞. The BdG equations therefore fail to predict any superfluid to Mott insulator phase transition. Since for non integer fillings the system is never a Mott insulator the BdG equations describe much better incommensurate systems.

1.2.5 External trapping potential

Up to this point we have discussed the basic aspects of the superfluid-Mott insulator transition in a translationally invariant system. The situation is different for a inhomogeneous system with a fixed total number of atoms and external confinement. This is the case realized in experiments, where besides the lattice there is a harmonic trap that collects the atoms at the center. In this case the density of atoms is not fixed since the atoms can redistribute over the lattice and change the local filling factor.

To deal with the inhomogeneous case, it is possible to define an effective local chemical potential, μj = μ , - Vj, at each lattice site j. If the change in the mean number of atoms between neighboring sites is small, the system can be treated locally as an homogeneous system. Because in the inhomogeneous case, the local chemical potential is fixed by the density, as the ratio U/J is changed the system can cross locally the boundary between the superfluid and Mott insulator phases. For example if the chemical potential at the trap center falls into the n ₀ = 2 Mott lobe, one obtains a series of Mott domains separated by a superfluid by moving to the boundary of the trap mimicking a wedding cake. In this manner, all different phases that exist for given J/U are present simultaneously.

The existence of such wedding-cake-like density profiles of a Mott insulator is supported by Monte Carlo calculations in one, two, and three dimensions. It has been pointed out recently nevertheless that at the quantitative level the local density approximation is not fully correct (Rigol, Batrouni, Rousseau, & Scalettar, 2009) and corrections of the local density approximation have been measured in 2D systems (Jimenez-Garcia et al., 2010).

The shell structure has been experimentally measured by using indirect methods such as mean field shifts (G. K. ^{Campbell et al., 2006}) and spatially selective microwave transitions together with spin changing collisions (Fölling, Widera, Müller, Gerbier, & Bloch, 2006). Just recently direct insituo imaging methods with single lattice site addressability have successfully resolved the shell structure in 2D systems (Gemelke, Zhang, Hung, & Chin, 2009; ^{Bakr et al., 2010}; Sherson, Weitenberg, Endres, Cheneau, & Bloch, 2010)

So far we have restricted our analysis to spinless bosons. In the following we will consider the more general situation in which the bosons carry an additional spin degree of freedom and discuss the application of spinfull bosons towards the investigation of quantum magnetism.

1.3.1 Exchange interactions

As initially proposed for electrons by ^{Dirac and Fowler (1926}, 1929) and ^{Heisenberg (1926,} 1928), effective spin-spin interactions can arise due to the interplay between the spin-independent Coulomb repulsion and exchange symmetry and do not require any direct coupling between the spins of the particles. The nature of such spin-exchange interactions is typically short-ranged, since it is governed by the wave function overlap of the underlying electronic orbitals.

The importance of exchange interactions can be already seen in two electron systems such as the excited states of the helium atom with one electron in the hydrogenic ground state Φ ₁ and the other in an excited state Φ ₀ . Since Coulomb repulsion is spin independent the ground state should be four fold degenerate, i.e. the states │↑, ↑〉, │ ↓, ↓〉, │↑, ↓〉, │ ↓, ↑〉 should have in principle the same energy. However, the Pauli exclusion principle enforces the wave function of fermionic particles to be antisymmetric under exchange and thus only the following configurations are allowed = (│Φ ₁ Φ ₂ 〉 - (│Φ ₁ Φ ₂ 〉) x │S = 1,M = ±1,0〉 also called orthohelium and (│Φ ₁ Φ ₂ 〉 + (│Φ ₁ Φ ₂ 〉) x │S = 0,M = 0〉 the so called parahelium. With │S = 1,M = 1〉 = │ ↑,↑〉, │S = 1,M = -1〉 = │↓,↓〉 and │S = 1,M = 0〉 = (│ ↑,↓〉 + │ ↓,↑〉), the spin triplet states and │S = 0,M = 0〉 = (│ ↑ ↓ 〉,- │↓ ,↑ 〉) the singlet state. Since in the singlet spin configuration the electrons tend to be closer due to the spatially symmetric wave functions, the singlet experiences higher interaction energy than the triplets which on the contrast minimize their interaction energy by avoiding each other. Consequently quantum statistics breaks the spin degeneracy and makes the orthohelium (triplets) the lower energy state.

The singlet-triplet energy splitting induced by exchange can be modeled by the following spin Hamiltonian

with V _ex a coupling constant proportional to the orbital overlap.

While exchange interactions have been extensively discussed in the context of fermions, similar exchange effects also apply to bosons, such as ⁸⁷Rb with nearly state-independent interactions in the lowest hyperfine manifold F = 1. Except here the particle wavefunctions are symmetrised rather than anti-symmetrised and the singlet state becomes the low energy state. Exchange interactions between pairs of ⁸⁷Rb atoms were measured and controlled in Anderlini et al. (2007) by initially preparing a Mott insulator state. Actually in the experiment the exchange interactions were used to create entanglement between pair of atoms.

Arrays of paired atoms were isolated and manipulated in the experiment via a double well lattice. This "superlattice" potential can be obtained by superimposing two standing light fields one with twice the periodicity of the other. Within a well the lowest and first vibrational states were used as the Φ ₀ and Φ ₁ orbitals, and the hyperfine states │F = 1, m _F = -1〉, │F = 1, m _F = 0〉 as the effective spin states ↑,↓

The steps followed in the experiment are summarized in Fig. 6

In a recent experiment exchange interactions were also probed and utilized for entanglement generation using optical tweezers instead of optical superlattices. See Kaufman et al. (2015) for details.

1.3.2 Superexchange interactions

In several topical insulators, such as ionic solids like e.g. CuO and MnO, however, an-tiferromagnetic order arises even though the wave function overlap between the magnetic ions is practically zero. In this case a "superexchange" interaction can be effective over large distance, as introduced in the seminal works of ^{Kramers and Anderson (Kramers, 1934}; P. W. Anderson, 1950). Here, the spin-spin interactions are mediated by higher-order virtual hopping processes, which in the case of bosons (fermions) leads to an (anti)-ferromagnetic coupling between particles on neighboring lattice sites (^{Auerbach, 1994}). Such superexchange interactions are believed to play an important role in the context of high-T _c superconductivity (Lee, Nagaosa, & Wen, 2006).

The first time-resolved experimental observation of superexchange interactions in double-well superlattices was reported in Ref. (^{Trotzky et al., 2008}), which we proceed to describe. It used isolated system of two coupled potential wells, which constitutes the simplest concept for the investigation of superexchange-mediated spin dynamics between neighboring atoms. The experiment used a single double well potential occupied by a pair of bosonic atoms with two different spin-states denoted by t and 4-. If the vibrational level splitting in each well is much larger than all other relevant energy scales and intersite interactions are neglected, the system can be described in a two-mode approximation by the Bose-Hubbard-type Hamiltonian which can be written as:

where the operators and â _σL ,_R create and annihilate an atom with spin σ in the left and right well respectively, count the number of atoms per spin-state and well. J _σ are the tunnel matrix elements (the subscript σ accounts for the fact that in principle the two spin states can see different lattice potentials) and are the onsite interaction energies between an atom in spin σ and other in σ'. The state of the system can be described as σ superposition of Fock states {│σ, σ'〉, │σσ', 0〉│0, σσ'〉}, where the left and right side in the notation represent the occupation of the left and right well, respectively, and the states │ σσ ',0〉 and │0, σσ '〉 are spin triplet states (the singlet states require the population of the second band to satisfy the symmetrization of the bosonic wave function and are consequently energetically forbidden).

In the limit of dominating interactions (U_σ , _σ ' » J), when starting in the subspace of singly occupied wells, the energetically high lying doubly occupied states can only be reached as "virtual" intermediate states in second order tunneling processes. Such processes lead to a nearest neighbor spin-exchange interaction, which couples the states │σ, σ'〉:

These type of processes can be recast in an effective spin Hamiltonian valid for an arbitrary spin configuration (^{Auerbach, 1994}; Duan, Demler, & Lukin, 2003; Kuklov & Svis-tunov, 2003):

where the effective coupling constants are given by

and acts as an effective magnetic field. Here constant terms have been dropped. Note that for fermionic particles the same Hamiltonian can be derived but one needs to replace → - and set the U _σ,σ terms to infinity, to account for the fermionic statistics.

For the generic case of spin independent interactions U _σ,σ = U (with a very good approximation the case for ⁸⁷Rb atoms in the hyperfine level F = 1, which was the system used in ^{Trotzky et al. (2008)}), and spin independent tunnelings J _σ,σ ' = J, one recovers the Heisen-berg Hamiltonian

with J_ex = . The negative sign is for bosons and favors ferromagnetic correlations. The positive sign is for fermions and favors antiferromagnetism.

Let's now allow for the presence of a bias or a tilt along the double well axis. This is capture by adding to the Hamiltonian a term

When a potential bias ∆ > 0 is applied, the degeneracy of the two intermediate states in the superexchange process is lifted. For J, ∆ «.U this leads to a modification of the effective superexchange coupling with now J_ex = J² / ( U + ∆ )+J ² /(U - ∆ ) = 2J² U /(U ² - ∆² ). By tuning the bias to A > U, it is possible to change the sign of J_ex and therefore to switch between ferromagnetic and antiferromagnetic superexchange interactions. For J << │U - ∆│ the picture of an effective coupling via two virtual intermediate states is again valid.

The measured dynamical observable was time-dependent Neel order parameter or spin imbalance N _z - ( n↑ _L + n↓_R - n↑ _R - n↓ _R )/2 starting with double wells initially prepared in │↑,↓). Here n↑,↓ _;L ,_R -〈 t,;_;í,_R〉 denote the corresponding quantum mechanical expectation values.

Since the total magnetization operator, Ŝ^z _L + Ŝ^z _R, commutes with the Hamiltonian, for the initial conditions in consideration we can restrict the Hilbert space to those states with zero net total magnetization i.e. │↑,↓〉, │↑, ↓〉,│0, ↑↓〉 and │ ↑↓,0〉 . In this case for symmetric double wells (∆ = 0), the Hamiltonian Eq. 90 can be diagonalized analytically to give a valid picture for all values of J and U within the single band Bose-Hubbard model (BHM). A convenient basis is given by the spin triplet and singlet state │t/s〉 = (│↑,↓〉 ± │↑,↓〉) and the states │±〉 ≡ (│↑↓,0〉 ±|0, ↑↓〉)/ . Two of the eigenstates are linear combinations of │t〉 and │+〉, where the one having the larger overlap with │t〉 is the ground state. The spin singlet │s〉 and the state │ - 〉 are already eigenstates themselves with energy 0 and U respectively. As a direct consequence, │ - 〉 cannot be reached from the initial state │↑,↓〉 = (│s〉 + │s〉)/ . Therefore, the dynamical evolution of the spin imbalance contains only two frequencies

The extraction of these frequencies from time-resolved measurements allows for the determination of 2J = ħ and U = ħ(ω ₁ - ω ₂ ) within the BHM.

In order to investigate the spin dynamics atoms in ^{Trotzky et al. (2008)} atoms were prepared in a 3D array of double wells with a Neel-type antiferromagnetic order, │↑,↓〉, in each double-well. After rapidly ramping down the short lattice and thereby the double well barrier the system was let to evolve for a hold time t, after which the spin-configuration was frozen out by ramping up the barrier and then the ensemble average N _z (t) measured.

Two typical time traces characterized the dynamics. For low barrier depths (J/U > 1), there was a pronounced time evolution of the spin imbalance N _z (t) consisting of two frequency components with comparable amplitudes and frequencies. With increasing interaction energy U relative to J, the frequency ratio increased, leaving a slow component with almost full amplitude and an additional high-frequency modulation with small amplitude (Fig. 1.3.2). For J/U«1, the fast modulation is completely suppressed and the only process visible was the superexchange oscillation.

1.3.3 Extended Bose-Hubbard Model

The inclusion of next-neighbor interactions to the simple BHM Hamiltonian defined in Eq. 90 leads to the extended Bose-Hubbard-type Hamiltonian. For simplicity let's focus on the experimentally relevant case of two coupled wells occupied by exactly two atoms in the two different spin-states │↑〉 and │↓〉

where U_LR corresponds to the inter-well interaction energy U_LR = and AJ = and ∆J= . It is apparent that those terms in Eq. 96 proportional to AJ act in the same way on an arbitrary Fock state in the system as the tunneling operator does and therefore modify J to J' = J + ∆J, while the terms proportional to U_LR lead to an energy shift of the states with exactly one atom on each site with respect to the states with double occupancy in a single well. The states │s〉 and │ - 〉 stay eigenstates of the system and the states │t〉 and │+〉 are now coupled via the matrix

with the eigenvalues

In the regime of strong interactions, the frequency ω ₂ is assigned to the superexchange process and by means of perturbation theory up to second order, one finds the effective coupling parameter to be J'_ex = 2J' /U - U_LR. The direct spin-exchange term U_LR favors an antiferromagnetic ground state for repulsively interacting bosons. However the case U_LR > 2J' ² /U is never reached and the ground state remains always ferromagnetic for bosons (anti-ferromagnetic for fermions). This is in agreement with the Lieb-Mattis theorem, which states that the ground state of two bosons in a spin independent potential is a spin triplet state. Experimentally deviations from the simple BHM were observed in ^{Trotzky et al. (2008)}. In fact, the EBHM description yields ħ(ω ₁ -ω ₂ ) = U + 3U_LR which better described the experimentally measured frequency difference for small short-lattice depths.