Physical replicas and the Bose glass in cold atomic gases

S Morrison^1,2,7, A Kantian^1,2, A J Daley^1,2, H G Katzgraber³, M Lewenstein^4,5, H P Büchler⁶ and P Zoller^1,2

¹ Institute for Theoretical Physics, University of Innsbruck, Technikerstr. 25, A-6020 Innsbruck, Austria
² Institute for Quantum Optics and Quantum Information of the Austrian Academy of Sciences, A-6020 Innsbruck, Austria
³ Theoretische Physik, ETH Zurich, CH-8093 Zürich, Switzerland
⁴ ICAO-Institut de Ciéncies Fotóniques, Parc Mediterrani de la Tecnologia, E-08860 Castelldefels, Barcelona, Spain
⁵ ICREA—Instituciò Catala de Ricerca i Estudis Avançats, 08010 Barcelona, Spain
⁶ Institute for Theoretical Physics III, University of Stuttgart, Pfaffenwaldring 57, 70550 Stuttgart, Germany

⁷ Author to whom any correspondence should be addressed.

E-mail: sarah.morrison@uibk.ac.at

Received 5 May 2008
Published 16 July 2008

Contents

• 1. Introduction
• 2. Bose-Hubbard model with disorder
• 3. Measurement of the overlap function
  • 3.1. Initial preparation of disorder replicas
  • 3.2. Measurement process
    • 3.2.1. Quantities that need to be measured
    • 3.2.2. Combining layers
    • 3.2.3. Number distribution characterization
• 4. Overlap function in the disordered Bose-Hubbard model
  • 4.1. Analytical determination of the overlap function
  • 4.2. Numerical results
• 5. Conclusion
• Acknowledgments
• Appendix A. Energy shifts and couplings for three and four particles
  • A.1. Three particles
  • A.2. Four particles
• Appendix B. Perturbation theory
  • B. 1. MI
  • B.2. BG
• Appendix C. Mean field theory
• References

1. Introduction

The interplay between disorder and interaction gives rise to a plethora of fundamental phenomena in condensed matter systems. The most predominant examples include spin glasses [1]–[3], the superconducting-to-insulator transition in thin superconducting films [4, 5], and localization phenomena in fermionic systems, such as weak localization and the metal–insulator transition [6]. While the nature of order in spin glasses and its theoretical description is still highly debated [1], [7]–[14], a substantial contribution to the understanding of bosons in a disordered medium has been provided by work on the disordered Bose–Hubbard model introduced by Fisher et al  [15]. The disordered Bose–Hubbard model has recently attracted considerable attention due to its potential realization with cold atomic gases in an optical lattice [16]–[22].

A general challenge in the description of glass phases in disordered media is the absence of a simple order parameter distinguishing the different ground states. This problem becomes evident in the disordered Bose–Hubbard model, where the phase diagram is determined by the competition between a superfluid (SF) phase, the Mott-insulator (MI) and a Bose-glass (BG) phase. While the SF phase exhibits a finite condensate fraction and is characterized by off-diagonal long-range order, the Bose glass is only distinguished from the incompressible MI by a vanishing excitation gap and a finite compressibility. However, in experiments on cold atomic gases the excitation gap and the compressibility are difficult to determine accurately and are also obscured by the finite harmonic trapping potential. The challenge is therefore to develop measurement schemes allowing for the experimental determination of these observables or to develop new observables to characterize the glass phase. While inaccessible in `real materials', the Edwards–Anderson order parameter [1, 23] is commonly studied analytically and measured numerically to quantify the `order' in a disordered system. It can be expressed as the correlation between independently evolving systems with the same disorder landscape (the so-called replicas), or as a correlation of the same system at different temporal measurements. Because the latter depends on an extra variable (temporal measurement window) it is advantageous to study two replicas of the system with the same disorder. Thus, the measurement of this order parameter requires the preparation of several samples with exactly the same disorder landscape, and the subsequent measurement of correlations between them.

We demonstrate that such a procedure is feasible in cold atomic gases allowing one to gain access to characteristic properties of disordered systems naturally hidden in condensed matter systems. The basic idea is to focus on cold atomic gases in an optical lattice: along one direction, the optical lattice is very strong and divides the system into independent two-dimensional (2D) planes [24]–[27]. In addition, the system is subjected to a disorder potential and we are interested in the situation, where in each plane the same disorder landscape is realized (see figure 1). Although the atoms in different planes are decoupled, the presence of the same disorder landscape within each plane induces a correlation between the different realizations: this correlation is of special interest in the replica theory of spin glasses and allows one to measure the overlap function, a characteristic property of spin glasses. Furthermore, we show that this correlation function is accessible in experiments on cold gases: the main idea is to quench the motion of the atoms by a strong optical lattice, and combine the different planes into a single one. Subsequently, the particle number occupation within each well carries the information about the correlation function. The possibility for the accurate determination of the particle number within each well of an optical lattice has recently been demonstrated experimentally for the SF–MI quantum phase transition [28].

Note that this measurement scheme for the overlap between systems with the same disorder can be applied to any disordered 1D or 2D system realized with cold atomic gases in an optical lattice. Of special interest is the realization of spin glasses and their study in the quantum regime. As an application, we focus here on the disordered Bose–Hubbard model and calculate the overlap function analytically in different regimes and compare it with 1D numerical simulations. We find that the overlap function q exhibits a sharp crossover from the Mott-insulating phase with q≈0 ^Note8  to the BG phase with q≈p(1–p), with p(0, 1), and thus makes it possible to distinguish the two phases. Because the SF phase can be detected via the interference peaks in a time-of-flight measurement, we propose that this novel measurement scheme for the overlap function can be used for a qualitative experimental verification of the phase diagram for the disordered Bose–Hubbard model.

Different implementations of disorder or quasi-disorder in cold atomic gases have been discussed and experimentally realized. Several groups attempted to search for traces of Anderson localization and the interplay of disorder–nonlinear interactions in Bose–Einstein condensates (BECs). The first experiments  [29]–[33] were performed with laser speckles that had a disorder correlation length larger than the condensate healing length. Localization effects were thus washed out by interactions. As an alternative, superlattice techniques—which are combinations of several optical potentials with incommensurable spatial periods—were used to produce quasi-disordered potentials with short correlations. This approach allowed the observation of some signatures of the Bose glass by the Florence group [34]. Only very recently, the Palaiseau group has managed to create speckle potentials on the sub-micron scale and directly observe Anderson localizations effects in a BEC released in a 1D waveguide [35] based on the theoretical predictions of [36, 37]. The Florence group reported observations of localizations phenomena in quasi-disordered potentials in BEC of K³⁹, which allows for complete control of the strength atomic interactions using Feshbach resonances [38]. It is worth noting that experimental attempts for local addressability in an optical lattice also offer the possibility to create disorder with correlation lengths comparable to the lattice spacing [39, 40]. Mixtures of cold gases, on the other hand, provide an alternative approach for the realization of disorder by quenching the motion of one species by a strong optical lattice, which then provides local impurities for the second species [41, 42]. While production of 2D planes with equal disorder realizations follows naturally if the disorder is implemented with laser speckles, we also show how such a system can be realized in the case of disorder induced by a mixture of atomic gases.

In section 2, we start with a detailed description of the disordered Bose–Hubbard model and the different disorder realizations. After a short review of the phase diagram of the disordered Bose–Hubbard model, we provide the definition of the overlap function q characterizing the disorder-induced correlations between independent realizations of the system. In section 3, we describe in detail the preparation of the system and the subsequent measurement scheme for the overlap function. In section 4, we determine the overlap for the disordered Bose–Hubbard model via the perturbation theory for physically relevant limits and compare the result with numerical simulations of the 1D system. Details of the calculations are presented in the appendices.

2. Bose–Hubbard model with disorder

Cold atomic gases subjected to an optical lattice are well described by the Bose–Hubbard model [33, 43]

with b^†_i (b_i) the bosonic creation (annihilation) operator and n_i = b_i^†b_i the particle number operator at the lattice site i. The first term describes the kinetic energy of the atoms with hopping energy J between nearest-neighbour sites, the second term accounts for the onsite interaction between the atoms, and the third term describes the disorder potential with random-site off-sets Δ_i of the chemical potential μ.

The disorder potential {Δ_i} can be generated via different methods [29]–[31], [34, 41, 42] and is determined by a probability distribution P(δ) describing the probability of having an onsite shift with strength δ. The mean square of the disorder distribution

gives rise to a characteristic energy scale Δ of the disorder potential (note that the mean energy shift of the disorder potential is absorbed in the definition of the chemical potential). In addition, the disorder potential is characterized by a correlation length. Here, we are interested in the short-range disorder with the disorder in different wells independent of each other. The probability distribution P(δ) and its correlation length depends on the source of the disorder potential, which varies depending on the microscopic implementation in cold gases. The most promising possibilities in order to produce disorder with short correlation length involve the use of laser speckle patterns and mixtures of different cold atomic gases. While first experiments with laser speckles had a disorder correlation length larger than the lattice spacing [29, 30], experimental efforts towards local addressability in an optical lattice also offer the possibility for the creation of disorder with a correlation length comparable to the lattice spacing, as well as Gaussian probability distributions [39, 40]. Alternatively, a disorder potential can also be created in mixtures of cold atomic gases in optical lattices [44, 45] by quenching the motion of one species. Then, the disorder correlation length is of the order of the Wannier function, which can be smaller than the lattice spacing due to the strong lattice that quenches the motion of the disorder species, while the probability distribution becomes bimodal with Δ_i = ±Δ. Here, we are primarily interested in such a short-ranged disorder with Δ_i being independent in the different wells. Note that both Gaussian disorder, as well as bimodal disorder generally gives rise to different physical phenomena in glass physics and thus both are interesting in their own right.

The zero-temperature phase diagram of the Hamiltonian (1) was first studied by Fisher et al  [15], where three different phases were discussed :  the SF, the MI, and the BG phases. In the 2D regime of interest here, the SF phase appears for large hopping energies JU, Δ and is characterized by off-diagonal (quasi) long-range order (finite SF condensate) and a linear excitation spectrum giving rise to a finite compressibility. On the other hand, for dominating interaction energies U  Δ, J the ground state corresponds to an incompressible MI phase with an excitation gap and a commensurate filling factor, i.e. the averaged particle number in a single well n_iN is an integer. While the quantum phase transition from the SF to the MI has been experimentally identified [46], the disorder potential gives rise to an additional BG phase characterized by a vanishing excitation gap and finite compressibility; see figure 2 for a sketch of the phase diagram. The details of the phase diagram have been studied via analytical methods in the regime of Anderson localization [42], [47]–[49] and in other regimes via numerical methods, such as the quantum Monte Carlo [50], density matrix renormalization group (DMRG) [51], dynamical mean-field theory [52], and analytical methods [53, 54]. Furthermore, recent interests also focused on the appearance of alternative phases for different disorder types in the Bose–Hubbard model, e.g. off-diagonal disorder giving rise to a Mott-glass phase [55] whereas random onsite interactions can give rise to a Lifshits glass [56].

Figure 2. Sketch of the phase diagram for the disordered Bose–Hubbard model: the MI appears for integer filling and the overlap function q vanishes for small hopping J/U  1, whereas in the BG phase the overlap function approaches a finite value q = p(1–p) with p(0, 1) characterizing the bimodal disorder distribution (see section 4). Consequently, a measurement of the overlap allows for a clear distinction between the BG and the MI phase in this regime. In turn, the SF phase is characterized by off-diagonal long-range order resulting in coherence peaks in a time-of-flight experiment.

The SF phase in the absence of disorder is experimentally identified via a measurement of the coherence peaks characterizing the condensate fraction, while the transition to a MI phase is characterized by the disappearance of these interference patterns as well as a change in the behaviour of excitations  [24, 46]. However, the identification of the BG phase in the presence of disorder requires an additional observable allowing for the distinction between the MI and the BG, where the condensate fraction vanishes. Such an additional property is known from the spin-glass theory as the Edwards–Anderson order parameter [1, 23] which appears as an order parameter in the mean field theory on spin glasses [57]. In the present situation with a disordered Bose–Hubbard model, its generalization leads to

where n = [n_i]_av denotes the mean particle density in the sample. It is important to note the different averages involved in the definition of the order parameter q_EA: for a fixed disorder realization, the average  ... denotes the thermodynamic average over the ground state of the system, while [ ...]_av describes the disorder average over different disorder realizations. While the experimental determination of this quantity in the Bose–Hubbard model requires an exact state tomography of the system for each lattice site, a simpler and alternative route is obtained by preparing two systems with the same disorder landscape and measuring the correlations between these two decoupled systems: the remaining correlation is disorder induced and disappears for weak disorder. This so-called overlap function q between the two systems realized with the same disorder landscape takes the form

where α and β describe the two different systems with the same disorder. This definition follows in close analogy of the mean-field order parameter in the replica theory of spin glasses [57].

3. Measurement of the overlap function

In this section, we outline the key process of this work, i.e. the basic experimental procedure by which physical replicas can be prepared and used to measure the overlap function q defined in (4). This consists of three essential steps: (i) initial preparation of disorder replicas, (ii) introduction of probe atoms, and (iii) the measurement process involving recombination of the replicas and spectroscopy to measure the overlap itself. At no stage do we assume individual addressability of lattice sites, but rather consider global operations and the measurement of global quantities.

3.1. Initial preparation of disorder replicas

We consider a situation in which atoms are confined to a 3D optical lattice, which is particularly deep along the z-direction, so that an array of planes with 2D lattice potentials are formed in the x–y-plane, with no tunnelling between the planes. Each of these 2D planes of the potential can be divided into two subplanes, again along the z-direction, by adding a superlattice of half the original period which can be controllably switched on and off [25, 26]. The resulting two subplanes then constitute our replicas, and our goal is that these pieces should have the same disorder realization.

In the case of disorder generated optically (using, e.g. speckle patterns), the different layers exhibit the same disorder landscape for a suitable orientation of the lasers. However, for a disorder induced by a second species, the replicas must be produced in several steps. We start with a two-species mixture in the undivided plane. A fast ramping up of the optical lattice depth within the plane for the disorder species results in a Poissonian distribution of particles in the lattice. The next step is then a filtering procedure removing singly occupied sites and keeping only doubly occupied and unoccupied sites (which are distributed randomly). Such a filtering procedure has been achieved in recent experiments where atoms are combined to the Feshbach molecules, and atoms remaining unpaired are removed from the system using a combination of microwave and optical fields that is energy-selective based on the molecular binding energy [58]. Finally, the superlattice is applied adiabatically providing a double-layer structure, where exactly one atom of the disorder species is transferred to each of the two subplanes. The result is two planes with a random distribution of atoms of the `disorder species', with an exact copy in each subplane, i.e. a replica of the random disorder pattern. This is illustrated in figure 1.

Interactions between atoms of a probe species and the disorder species can now be adiabatically increased, or the probe species otherwise adiabatically introduced to the system (e.g. using spin-dependent lattices or superlattices). This results in two replicas of the same system, which we label α and β, with atoms in each replica evolving independently according to their system Hamiltonian in the same disorder potential. We again reiterate that the superlattice separating the replicas should be deep, so that there is no tunnelling between the replicas. The only correlations between the two should thus be determined by the identical disorder realizations.

3.2. Measurement process

To measure the overlap q as given by (4), we need to compute correlation functions that compare the state of probe atoms in the two replicas α and β. In each case, the correlation functions can be computed if we measure the joint probability distribution for occupation numbers in the two replicas (this is discussed in detail below). In particular, we need to determine the joint probability of having n particles in layer α and m particles in layer β at a given lattice site, denoted by p_nm, and the probability of having n particles in layer α, β at a given lattice site, denoted by p_n^{α, β} (both depicted in figure 3). We thus first freeze the state in each replica by ramping the lattice suddenly to deeper values in all the three directions, so that tunnelling no longer occurs. We can then measure the different possible configurations of particles occupying individual lattice sites in the layers α and β (see figures 3–5).

This is done in two steps: first, the layers are combined so that at each site each initial eigenstate involving different occupation numbers in the two replicas becomes a superposition of degenerate states involving the total number of atoms spread over motional states in a single well (see figure 4, section 3.2.2 below). The final state is, of course, strongly dependent on the combination process, as well as the initial state before the combination of the layers. However, we need only determine the fraction of lattice sites with a total of n particles, denoted p_n, which is unchanged provided the lattice is deep within each plane, to prevent tunnelling of atoms. The relative frequency with which each final configuration p_n occurs can then be measured based on the different spectroscopic energy shifts arising from different numbers of atoms in the final well, and different distributions of atoms over the motional states (see figure 5, section 3.2.3 below). This involves weak coupling of atoms to an auxiliary internal state, similar to that performed in [28]. We show that the set of possible final states given an initial configuration can be distinguished spectroscopically from the set of states corresponding to initial configurations giving different contributions to p_n, and thus we can reconstruct the original joint probability distribution p_nm prior to combining the layers. The reduced probabilities p_n^{α, β} can be determined from similar spectroscopy without first combining the layers, allowing the full overlap function q to be constructed. Below, we give further details of each step in this procedure.

3.2.1. Quantities that need to be measured We need to measure two quantities to determine the overlap q, namely , where

with α≠β and L is the total number of lattice sites, corresponding to the density–density correlations between the two layers, and , , where

corresponding to the average density of each layer. Then, we can express

The disorder average should be performed by repeated measurements of the quantum average for different realizations of the disorder^Note9 . In general, one can expect that a measurement of n_i and n_j for large separation between sites i and j is independent of each other, and consequently the summation over lattice sites automatically performs a statistical quantum average^Note10 . In the following, we are interested in small filling factors such that lattice sites with three or more particles are rare and can be neglected in all three phases of the system.

Firstly, we focus on measurement of , i.e. the second term of the overlap in (4). If the layers are prepared as discussed above, the symmetry is preserved and thus we can measure the average over two replicas, where ^Note11 

For the average density, we obtain:

where p_n^{α, β} = N_n^{α, β}/L with N_n^{α, β} the number of sites with n particles in layers α, β, and thus

Next, we consider the measurement of , i.e. the first part of the overlap in (4). In this case, the density–density correlations can be expressed as

where p_nm = N_nm/L with N_nm the total number of sites with n particles in layer α and m particles in layer β, and thus, keeping terms up to two particles in total over the replicas at a given site

Now consider p_n = N_n/L, with N_n the total number of particles from both layers, which can be measured by combining both layers (described in detail in section 3.2.2) and subsequently applying a similar scheme to the one presented in [28] (described in detail in section 3.2.3). Using the following equations that relate both p_n and p_n^{α, β} to p_nm

we can show that which in conjunction with p₃ = p₁₂ + p₂₁ and p₄ = p₂₂ gives all necessary terms. Then, equation (12) simplifies to

Finally, note that as given in (10), can also be expressed as .

In the following two sections, we show how to measure these conditional probabilities by combining the two replicas together, and performing spectroscopy where a single atom is coupled to a different internal state. In section 3.2.2, we describe the process of combining the two layers into a single one; whereas in section 3.2.3, we describe a similar scheme to the one presented in [28] that can be used to distinguish different occupation numbers.

3.2.2. Combining layers We now describe in more detail the process of combining the two layers into a single one and determine the resulting joint state for the different possible initial configurations. The combination of the layers is achieved by lowering the potential barrier that separates both layers (switching off the superlattice). Specifically, we now determine the final (after completely lowering the barrier) state for each of the different initial configurations of particles (prior to lowering the barrier).

After the superlattice is removed and layers are combined, the system is described by the Hamiltonian

where , with a the scattering length, ω_⊥ the transverse trapping frequency, w_k(z) the Wannier functions for band k, and _k the band energy (k,l{0, 1}). The expression for U_kl is a special case of the general 3D expression

where w₀(x) and w₀(y) are the Wannier functions for the lowest motional state in the 2D replica plane (x, y), with the additional assumption of an isotropic and deep lattice potential in the replica plane characterized by the frequency ω_⊥.

States which are initially non-degenerate will evolve adiabatically to the corresponding eigenstate of the Hamiltonian in (14), provided the barrier is lowered slowly with respect to the energy separation between the (nondegenerate) eigenstates. However, we note that some of the initial states are doubly degenerate due to the equal single particle energies in the two layers. Thus, when the barrier is lowered these degenerate states become coupled to each other (because the single particle states change time dependently) and the final state will evolve into a process-dependent superposition of the different corresponding eigenstates of the Hamiltonian in (14). However, this poses no problem for our measurement scheme since the states that are initially degenerate contribute equally to the probabilities we are trying to measure. Moreover, for a given initial configuration it is not necessary to know the full final state, in fact it is sufficient to know only the possible corresponding final eigenstates of the Hamiltonian in (14) (with the same number of particles as the initial configuration) that contribute to the final state. Figure 4 depicts all the different possible states before and after combining the layers for up to a maximum of four particles in total.

In the next section, we show how the different configurations for a given particle number can be measured and thus the total occupation numbers determined.

3.2.3. Number distribution characterization We now investigate a scheme similar to the one presented in  [28] that uses density-dependent energy shifts to spectroscopically distinguish different particle numbers at each site of the combined layer (see section 3.2.2 above) resulting in a measurement of p_n. The measurement of p_n^{α, β} follows from the same spectroscopic analysis, but performed on each individual layer α, β before combining the layers. Such a scheme has already been implemented experimentally and was utilized to observe the Mott-insulating shell structure that arises due to a harmonic trapping potential usually present in cold atomic gases [28].

We begin with all atoms in an internal state a, and investigate the weak coupling of particles to a second internal state b (with at most one particle transferred). Due to the difference in interaction energy between particles of species a and b, and particles solely of species a there is an energy shift between the `initial' state (all particles of species a) and the `final' state (one particle of species b and all other particles of species a). Thus, different initial particle numbers can be distinguished by the energy shift (corresponding to the Raman detuning at which the transfer is resonant), provided the shift is different for different particle numbers. Since particles in the combined layer can be spread across multiple motional states (see section 3.2.2 above), different energy shifts occur for different configurations of the same number of particles. While we find that the energy shifts of some of the configurations (of the same number of particles) are identical or very similar, there are other configurations having substantially different energy shifts. Thus, we have to check that the energy shifts for all configurations of a given number of particles can be distinguished from energy shifts of configurations corresponding to a different total number of particles. To this end, we now present a detailed calculation of the energy shifts arising from all the different initial configurations of particles and show in which cases the total number of particles can be distinguished.

Let a_k^† and b_k^† denote the creation operators for particles of species a and b, respectively, in band k, where k{0, 1}. Then the Raman coupling between the two internal states a and b is described by an effective Hamiltonian (within a rotating-wave approximation)

where Δ(t) is the Raman detuning, Ω(t) the effective two-photon Rabi frequency, and we have assumed that Ω(t)  Δ  ' where Δ  ' are the detunings of the lasers from the atomic excited state. In what follows, we assume that , i.e. weak coupling to internal state b and that the two lasers creating the Raman transition are running waves with equal wave vectors. Note that processes in which the internal state and the band would change do not occur provided the Raman detuning and effective two-photon Rabi frequency are smaller than the band separation.

The onsite interaction energy for atoms of species a and b in the lowest two bands is described by the following Hamiltonian

where

with a_ij the scattering length between two particles in internal states i and j (i,j{a,b}), ω_⊥ the transverse trapping frequency, and w_k(z) the Wannier functions corresponding to band k (k,l{0, 1}). Note we have also introduced the particle number operators n_k^a = a_k^†a_k and n_k^b = b_k^†b_k. We further assume that the lattice is very deep so that we can approximate all the Wannier functions by simple harmonic oscillator eigenfunctions. This assumption of a deep lattice gives the following simplified expressions 2U₀₁^ij = U₀₀^ij and U₁₁^ij = (3/4)U₀₀^ij.

Using (16), we now explicitly compute the energy shift ΔE = E_f–E_i associated with the processes |i→|f, where |i is the initial state corresponding to all particles of species a and having energy E_i, and |f is the final state corresponding to one particle transferred to internal state b and having energy E_f. Clearly, ΔE depends on how the particles are initially distributed in the lower and higher bands after combining the two layers (see the right-hand side of figure 4). There are two different cases to be considered; either all particles are initially in the same band (this is the situation in  [28]) or particles are initially in both bands. In the latter case, it is then possible that particles in either band are transferred to internal state b and since these different possible resulting states are coupled by the interaction Hamiltonian in (16) the final resulting state is a superposition of these, corresponding to an eigenstate of the Hamiltonian in (16).

For example, let us consider a total initial number of two particles as depicted in figure 5. The different initial and final states are most conveniently expressed in the occupation number basis |N₀^a, N₁^a; N₀^b, N₁^b with N₀^{a, b} (N₁^{a, b}) the number of particles of species a and b in the lower (higher) band at a given lattice site. First consider the case where both particles are initially in the lower or higher bands corresponding to the initial states |2, 0; 0, 0 or |0, 2; 0, 0 (see top and bottom of figure 5). The corresponding final states are given by |1, 0; 1, 0 and |0, 1; 0, 1, and the corresponding energy shifts are ΔE/U₀₀^aa = –1 and , where we have defined ≡U₀₀^ab/U₀₀^aa. Next, we consider the case where there is a particle in each band initially corresponding to the initial state |1, 1; 0, 0 (see middle of figure 5). Then the corresponding final states are the eigenstates of the Hamiltonian in (16) in the subspace corresponding to one particle of species a and one particle of species b given by , with eigenenergies E ₊ /U₀₀^aa =  and E_–/U₀₀^aa = 0. The associated energy shifts for these two possible final states are ΔE/U₀₀^aa = –1 and ΔE/U₀₀^aa = –1, respectively. Note that the Raman transition only couples the state |1, 1 ₊  to the initial state |1, 1; 0, 0 with a matrix coupling element , where Ω_fi≡f|H_RC|i, while for the state |1, 1_– we find Ω_fi = 0. Figure 5 summarizes the different possible configurations for a total of two particles together with the corresponding energy shift ΔE and matrix coupling element Ω_fi. We note that the energy shifts, ΔE, for the different final states to which the initial states couple (i.e. Ω_fi≠0) are very similar, this is of advantage since our scheme needs only to distinguish energy shifts for different total particle numbers.

Similarly, we can determine the energy shifts and coupling matrix elements for other numbers of particles. If all particles are initially in the lower or higher state corresponding to the initial states |N₀^a, 0; 0, 0 or |0, N₁^a; 0, 0, the corresponding final states are given by |N₀^a–1, 0; 1, 0 and |0, N₁^a–1; 0, 1. The corresponding energy shifts are then given by ΔE/U₀₀^aa = (N₀^a–1)(–1) and , and the matrix coupling element in each case is and , respectively. For other initial configurations with m particles in the lower band and n particles in the higher band the corresponding final states are found by diagonalizing the Hamiltonian in (16) in the appropriate particle subspace. This always results in a pair of eigenstates denoted |m,n_± with eigenenergy E_±^{(m, n)} and using these expressions the energy shift ΔE and coupling matrix elements Ω_fi can be directly determined. For a total initial number of three and four particles the expressions for |m,n_±, E_±^{(m, n)}, ΔE and Ω_fi are given in appendices A.1 and A.2, respectively, for all the different initial configurations. In table 1, we list the results for up to a maximum number of four particles; the first column lists the different initial configurations and corresponding final states, the second column lists the corresponding matrix coupling elements Ω_fi and the third column lists the associated energy shifts ΔE.

Table 1. Resonant energies and coupling strengths for the transfer of single atoms in the combined layer to a different internal state. Column one lists the different initial configurations and the corresponding final states up to a maximum number of four particles. In the second and third columns, we list the matrix coupling elements Ωfi and energy shifts ΔE, respectively, corresponding to the different initial configurations in the first column. The states |m, n± are listed in the text and appendices together with the corresponding eigenenergies E±(m, n) and coupling matrix elements η±(m, n). The last two columns give numerical values of Ωfi and ΔE for the specific value  = 0.98.

|i→|f	Ωfi	ΔE/U00aa	Ωfi	ΔE/U00aa
|1, 0; 0, 0→|0, 0; 1, 0	1	0	1	0
|0, 1; 0, 0→|0, 0; 0, 1	1	0	1	0
|2, 0; 0, 0→|1, 0; 1, 0			1.4	–0.018
|2, 0; 0, 0→|1, 0; 1, 0		(–1)	1.4	–0.024
|1, 1; 0, 0→|1, 1 +		(–1)	1.4	–0.024
|1, 1; 0, 0→|1, 1–	0	–1	0	–1
|3, 0; 0, 0→|2, 0; 1, 0		2(–1)	1.7	–0.049
|2, 1; 0, 0→|2, 1 +		2(–1)	1.7	–0.049
|2, 1; 0, 0→|2, 1–	0		0	–1.5
|1, 2; 0, 0→|1, 2 +	η + (1, 2)	E + (1, 2)–3	1.7	–0.29
|1, 2; 0, 0→|1, 2–	η–(1, 2)	E–(1, 2)–3	–0.0034	–1.8
|0, 3; 0, 0→|0, 2; 0, 1			1.7	–0.037
|4, 0; 0, 0→|3, 0; 1, 0		3(–1)	2	–0.073
|3, 1; 0, 0→|3, 1 +		3(–1)	2	–0.073
|3, 1; 0, 0→|3, 1–	0	–3	0	–2.0
|2, 2; 0, 0→|2, 2 +	η + (2, 2)		2.0	–0.070
|2, 2; 0, 0→|2, 2–	η–(2, 2)		–0.0031	–2.0
|1, 3; 0, 0→|1, 3 +	η + (1, 3)		2.0	1.4
|1, 3; 0, 0→|1, 3–	η–(1, 3)		–0.0054	–0.52
|0, 4; 0, 0→|0, 3; 0, 1			2	–0.055

It is possible to distinguish different occupation numbers for specific ranges of if the energy shifts for different total numbers of particles are distinct. The desired range of is then between points where any energy shifts for different total number of particles would coincide (e.g.  = 1)^Note12 . For a specific example, we consider the value  = 0.98 (corresponding to scattering length values in [28]) and give numerical values for the matrix coupling elements Ω_fi and the energy shifts ΔE in columns four and five, respectively, in table 1. Then, to a high accuracy, the determination of p_n is obtained by applying the coupling between the different internal states a and b for a short time at all resonant frequencies in table 1 with total number of particles n (ignoring the resonances with small and vanishing couplings Ω_fi). The total number of transferred particles is then proportional to p_n and efficiently avoids problems arising from non-adiabatic combination of the layers.

4. Overlap function in the disordered Bose–Hubbard model

In this section, we determine the overlap function q (see (4)) in the disorder Bose–Hubbard model, both analytically and numerically. We focus on a bimodal disorder potential Δ_i = ±Δ with the probability distribution P(Δ) = 1–P(–Δ) = p, which is the appropriate description for the chosen disorder implementation (see discussion in earlier sections). In the following, we present the modifications due to the presence of a weak (Δ < U/2) bimodal disorder potential, and the prediction for the overlap function within the different phases (see figure 2). It follows from this constraint that we do not enter the compressible disordered phase, and hence the system will always have a unique ground state. Thus, the ground states in two exact copies α, β of the system (i.e. same disorder distribution) are the same, and consequently n_i^α = n_i^β. We therefore expect a Sherrington–Kirkpatrick-type behaviour for the overlap [60], i.e. the overlap between any two copies will always be the same, and thus the overlap function q can in fact be calculated using (3).

While the above is a self-evident statement for a replica-symmetric system such as the one we examine here, it may no longer be true for systems thought to exhibit replica-symmetry breaking (e.g. a classical spin glass type model [57], or quantum systems with extended interactions [61]). For the latter type of system, measurement of the overlap function using physical replicas (see section 3) is thought to yield different values for the overlap depending on α and β, the defining signature of replica-symmetry breaking.

In section 4.1, we analytically calculate the overlap function for the various phases of the disorder Bose–Hubbard model (1); whereas in section 4.2, we present results for the overlap function from numerical simulation of the 1D system.

4.1. Analytical determination of the overlap function

Starting with the limit of vanishing hopping J = 0, equation (1) reduces to an onsite Hamiltonian and the ground state is given by with |Ω_i the lowest energy state within each site i. This lowest energy state takes the form |Ω_i = |n_i with |n_i being a state with fixed particle number n_i at site i, with n_i subject to the constraint

For a weak bimodal disorder potential (Δ < U/2), we have to distinguish two different cases.

Firstly, for a chemical potential within the range U(n₀–1) + Δ < μ < Un₀–Δ the above constraint in (18) is fulfilled independently of the site parameter i for the integer value n₀, i.e. at each site the ground state is characterized by the integer particle density [n_i]_av = n₀ and we obtain a MI phase (see figure 2). Furthermore, the overlap parameter vanishes identically in this limit, i.e. q = 0 for a MI at J = 0.

Secondly, for chemical potentials in the range Un₀–Δ < μ < Un₀ + Δ, the particle number at sites with different disorder potential Δ_i differs by a single particle, i.e. at each site with disorder potential Δ_i = Δ the particle number is determined by n₀, whereas at sites with disorder potential Δ_i = –Δ the lowest energy state is characterized by n₀ + 1 particles. Therefore, the averaged particle density and the overlap parameter take the form

This situation corresponds to the disorder-induced BG phase (see figure 2) and the ground state is characterized by non-homogeneous filling. Similarly to the case of the MI phase, the particle number is fixed over a finite range of the chemical potential and thus the BG is incompressible. In contrast to the MI phase, the BG phase is characterized by a non-vanishing overlap parameter dependent on p.

Next, we focus on small hopping, J  U, Δ and μ, and determine the overlap function using the perturbation theory. We seek a unitary transformation of the form U = e^–iS which transforms the total Hamiltonian in (1) to an effective Hamiltonian having the same eigenspectrum but acting only in the space of the unperturbed eigenstates. Using this unitary transformation, we can directly determine the corrections due to hopping in the perturbation theory by computing , where . Details are given in appendix B.

We start with the MI and consider the leading order (non-vanishing) correction to the overlap function which occurs at fourth order due to a site (disorder) independent density, Ω|n_i|Ω = n₀ (using |Ω_i = |n₀ in the MI). To the lowest order in the perturbation theory the overlap is then given by

where

with

and |_jk = |n₀ + 1_j|n₀–1_kΠ_l≠ij|n_0l. Using these expressions, it is straightforward to compute the expectation value of the density correction

where E₀–E_ij = –U + Δ_j–Δ_i. It is now straightforward to calculate the disorder averages and the overlap is given by

where z is the number of nearest neighbours.

Next, we turn to the BG phase and consider the leading order (non-vanishing) correction to the overlap function which occurs here at second order. The expression for the overlap in lowest-order perturbation theory is then given by

where q₀ = p(1–p) (the overlap at zero hopping), and S₁ are as given above for the MI phase but with

and

where d_{j, k, l}(1–d_{j, k, l}) is zero (one) for Δ_{j, k, l} = Δ (Δ_{j, k, l} = –Δ). Again, using these expressions it is straightforward to compute the expectation value of the density correction and for n₀≥1, we obtain

where

and

Once again a straightforward calculation of the disorder averages gives the following expression (n₀≥1) for the overlap:

An analogous calculation for n₀ = 0 gives

which corresponds to the U→0 limit of (31).

When the effects of hopping dominate, the bosonic atoms are in the SF phase (see figure 2). Then the influence of weak disorder on the SF phase can be studied within a generalized mean-field approach. In analogy to the Gross–Pitaevskii formalism, we replace the bosonic operators b_i in the Hamiltonian (1) by a local complex field ψ_i, and minimize the corresponding free energy functional H(ψ_i). The effect of disorder is to produce fluctuations in the local density. Hence, we make the ansatz , where n₀ is the homogeneous density in the absence of disorder and δn_iΔ_i is the disorder-induced density fluctuation. Introducing the notation δn = [δn_i]_av as the disorder-averaged density fluctuation, the particle density reduces to [n_i]_av = n₀ + δn. We substitute this ansatz into H(ψ_i), expand ψ_i for small fluctuations, δn_i  n₀ (see appendix C) and minimize the resulting expression in (C.1) at fixed particle number, which is equivalent to requiring [δn_i]_av = 0. Note that for large dimensions and thus a large number of nearest neighbours, we may replace

Using this relation and self-consistently solving for the density fluctuations δn_i (see appendix C) gives the overlap function

where .

In summary, we find that in the limit of small hopping J/U  1, the overlap signals a sharp crossover between the MI with q≈0 and the BG with q≈p(1–q), see figure 2. In turn, the SF phase is characterized by off-diagonal long-range order giving rise to coherence peaks in a time-of-flight picture. Consequently, the overlap and the coherence peaks in time-of-flight allow for a clear identification of the phases in the disordered Bose–Hubbard model. However, we would like to point out that the overlap function behaves smoothly across the phase transition and consequently is not suitable for the precise determination of the exact location of the phase transition.

4.2. Numerical results

We have simulated the disordered Bose–Hubbard model (1) in 1D using a number-conserving algorithm based on time-dependent matrix product states [62, 63]. In the matrix product state representations, we reduce the Hilbert space by retaining only the χ most significant components in a Schmidt decomposition of the state at each possible bipartite splitting of the system. States in 1D can be efficiently represented in this way, with relatively small χ giving essentially unity overlap between the represented state and the actual state. The simulations both serve as a check for the perturbation theory, and allow computation of the overlap function beyond the perturbation limit from the microscopic model.

Because we consider the regime 2Δ < U, the system remains incompressible and does not enter the disordered glass phase. Therefore, we do not expect replica-symmetry breaking to occur. For our simulations, the presence of replica-symmetry is equivalent to convergence to a unique ground state for any given disorder independent of initial conditions. The overlap function can then be computed as

where is the density of the system. Computation of the overlap function can be simplified in this way as self-averages to [n_i]_av for sufficiently large systems, which in turn self-averages to .

The ground states used in (35) are computed by imaginary-time evolution, cf [63], of the 1D Bose–Hubbard model given by (1) with a randomly chosen bimodal distribution, and different impurity strengths 2Δ. The lattice size is L = 64, representative of current experimental capabilities, and we compute the overlap as a function of Δ for several different values of U/J and the fraction of impurity occupied sites, p. In figure 6, we display results for U/J = 10 and p = 0.5. We have performed convergence tests in the time step, the duration of imaginary-time evolution, the number of disorder realizations and the number of retained states χ. We find that the overlap function converged even for surprisingly small values of χ~20 states in the large U/J limit (MI or BG regime). Regarding the number of disorder realizations for each value of Δ, we find that while 10 random realizations of bimodal disorder are not sufficient to narrow the error down, 20 realizations suffice. Much of the observed fluctuations in the value of the overlap function for a given Δ stem from the fact that we allow the number of impurities in each realization to fluctuate around the mean Lp. Conversely, in our simulation data, we observe different disorder realizations with the same number of impurities having overlaps closer to each other than the error bars in figure 6 would suggest.

As shown in figure 6, we generally found very good agreement between numerical simulations and the perturbation theory in its region of validity. For the MI regime (commensurate filling factor), this is the regime where 2Δ  U. We observe the breakdown of our second-order perturbation theory as Δ approaches U/2 = 5J, which is caused by energy-degeneracy of adjacent impurity- and non-impurity occupied sites when Δ = U/2. In the BG regime, breakdown of the perturbation theory sets in for low Δ, as degeneracy occurs there for Δ = 0.

5. Conclusion

We have shown how the BG phase can be identified by measuring the overlap function characterizing the correlations between disorder replicas, i.e. systems having identical disorder landscape. We have described a procedure to create disorder replicas using cold atomic gases in an optical lattice, focusing on the case of disorder induced by a second atomic species. A scheme to measure the overlap has been presented, which involves the characterization of the occupation numbers in both the individual replicas and the joint system, where both replicas have been combined. Specifically, we have shown that after combining the replicas together, particles can be distributed across multiple motional states and we have explained in detail how to determine the occupation number distribution in this situation. Using the perturbation theory, we have calculated the overlap for weak hopping and have shown the different behaviours exhibited in the MI and BG phases; these results are in good agreement with the results from a 1D numerical simulation. In the opposite limit of very large hopping, we also obtain analytical results for the overlap within a mean-field theory treatment.

Applying our proposed measurement scheme for the overlap to other types of disordered models [57, 61] would provide interesting insights into characteristic properties of disorder-induced quantum phases, such as the possible appearance of broken ergodicity and different degenerate ground states. Of particular interest would be the quantum spin glass, where the overlap corresponds to the Edwards–Anderson order parameter which has been conjectured to exhibit replica-symmetry breaking [57], i.e. the results for the overlap depend on the particular replica. Our method could also give novel insights into the physics of systems that possess multiple metastable states, such as ultracold dipolar gases [64] or quantum emulsions [65]; in particular, statistics of overlaps of such metastable states could be measured. In addition to the equilibrium properties considered here, the dynamics of the overlap might also possess interesting behaviour, and could be measured in a similar setup. In 1D, the time evolution of the overlap could be determined numerically [63].

Acknowledgments

Work in Innsbruck was supported by the Austrian Science foundation through SFB F15, the EU network SCALA and project I118_N16 (EuroQUAM_DQS). S M thanks the Department of Physics of the University of Auckland for their hospitality. H G K acknowledges support from the Swiss National Science Foundation under grant no. PP002-114713. M L acknowledges support of ESF PESC Programme `QUDEDIS' and Spanish MEC grants (FIS 2005-04627, Conslider Ingenio 2010 `QOIT').

Appendix A. Energy shifts and couplings for three and four particles

In the following two sections, we determine expressions for the eigenstates and eigenenergies of (16), i.e. |m,n_± and E_±^(mn) and associated energy shifts and couplings for all possible initial configurations of three and four particles. Note that the case where all particles are in the lowest or highest energy band has already been explained in section 3.2.3 and will not be discussed here further.

A.1. Three particles

We consider the two different configurations described by the initial states (i) |2, 1; 0, 0 and (ii) |1, 2; 0, 0.

For case (i), the two eigenstates are given by

with eigenenergies E ₊ ^{(2, 1)} = 2 + 1 and E_–^{(2, 1)} = /2 + 1, respectively. We see that only the state |2, 1 ₊  couples to the initial state |2, 1; 0, 0, with a matrix coupling element of , whereas for the state |2, 1_– we find Ω_fi = 0. The energy shifts for the two different possible processes are given by ΔE/U₀₀^aa = 2(–1) and .

For case (ii), the two eigenstates are given by

where

and with eigenenergies

In this case, both states |1, 2_± have nonzero coupling to the state |1, 2; 0, 0 given by

The energy shifts for the two different possible processes are given by ΔE/U₀₀^aa = E_±^{(1, 2)}–3.

A.2. Four particles

We consider the three different configurations described by the initial states (i) |3, 1; 0, 0, (ii) |2, 2; 0, 0 and (iii) |1, 3; 0, 0.

For case (i) the eigenstates are given by

with eigenenergies E ₊ ^{(3, 1)} = 3( + 1) and E_–^{(3, 1)} =  + 3, respectively. We see that only the state |3, 1 ₊  couples to the initial state |3, 1; 0, 0, with a matrix coupling element of , whereas for the state |3, 1_– we find Ω_fi = 0. The energy shifts for the two different possible processes are given by ΔE/U₀₀^aa = 3(–1) and ΔE/U₀₀^aa = –3.

For case (ii), the eigenstates are given by

where

and with eigenenergies

In this case, generally both states |2, 2_± have nonzero coupling to the state |2, 2; 0, 0 given by

The energy shifts for the two different possible processes are given by .

For case (iii), the eigenstates are given by

where

and with eigenenergies

In this case, generally both states |1, 3_± have nonzero coupling to the state |1, 3; 0, 0 given by

The energy shifts for the two different possible processes are given by .

Appendix B. Perturbation theory

In this section, we derive the corrections to the overlap function due to small hopping in the perturbation theory following closely the method presented in [66]. We start with the zero hopping Hamiltonian H₀ with ground state |Ω and eigenenergy E₀. We treat the hopping Hamiltonian, JH₁ with , which has eigenstates denoted by |_ij and eigenenergies denoted by E_ij in the perturbation theory. We proceed by determining a unitary transformation that transforms the total Hamiltonian H = H₀ + JH₁ into an effective Hamiltonian H ' with the identical eigenspectrum but acting only in the unperturbed Hilbert space of the state |Ω. The required unitary transformation may be written as U = e^–iS and can be determined consistently to any order by writing S = JS₁ + J²S₂ + ···. We assume that S has only non-diagonal matrix elements which connect the unperturbed ground state |Ω and the excited states |_ij. Using the unitary transformation, we can directly calculate the correction to the particle density due to the hopping by computing , where . We show in the following sections that it is sufficient to determine S to first order. Moreover, it is straightforward to show that to first order the matrix elements of S, i.e. S₁, are given by

B. 1. MI

In this section, we derive the leading order correction to the overlap function for small hopping in the perturbation theory for the MI phase. We begin by expanding S to fourth order and calculate the average corrected density to be

where

Firstly, note that since we have chosen S to have no diagonal matrix elements [S_n, n_i] = 0 to all orders.

We now calculate the first part of the overlap function as given in (3) and find

Similarly, for the second term of the overlap in (3), we obtain the following expression

Using n_i = n₀, we obtain the following expression for the overlap

Note that because n_i = n₀ at all lattice sites (independent of the disorder) the second- and third-order contributions cancel.

To determine the required commutators [S₁, [S₁, n_i]], we write S₁ and n_i as projectors

where

and

Using these expressions it is straightforward to compute the expectation value of the density correction

where E₀–E_ij = –U + Δ_j–Δ_i.

B.2. BG

In this section, we derive the leading order correction to the overlap function for small hopping in the perturbation theory for the BG phase. Again, we consider the corrected density as given in (B.2) but only up to second order in J. Again, we have [S_n, n_i] = 0 because we have chosen S to have no diagonal matrix elements. Thus, the overlap to second order is given by

where q₀ = [n_i²]_av–[n_i]_av² is the overlap at zero hopping derived in section 4.1. Since n_i is now site dependent, the second-order contribution does not vanish.

To determine the required commutators [S₁, [S₁, n_i]], we again write S₁ and n_i as projectors

where

with

and

where d_{j, k, l} (1–d_{j, k, l}) is zero (one) for Δ_{j, k, l} = Δ (Δ_{j, k, l} = –Δ). Again, using these expressions it is straightforward to compute the expectation value of the density correction and for n₀≥1, we obtain

where

and E₀–E_ij = –U[d_i + (1–d_j)] + Δ_j–Δ_i.

Appendix C. Mean field theory

We replace the annihilation operators b_i in the Hamiltonian in (1) by the mean field and expand the fluctuations δn_i  n₀ up to second order to obtain

Next, we minimize (C.1) at fixed particle number which is equivalent to requiring [δn_i]_av = 0. Using (33) and consistently solving gives

Finally, the overlap can be expressed in terms of these fluctuations as follows:

Substituting (C.2) into (C.3) and performing the disorder average for a bimodal distribution gives the result in (34) for the overlap function q.

References

[1]

Binder K and Young A P 1986 Rev. Mod. Phys. 58 801 
CrossRef

[2]

Mézard M, Parisi G and Virasoro M A 1987 Spin Glass Theory and Beyond  (Singapore: World Scientific) 

[3]

Young A P 1998 Spin Glasses and Random Fields  (Singapore: World Scientific) 

[4]

Finkel'Stein A M 1994 Physica B 197 636 
CrossRef

[5]

Fisher M P A 1990 Phys. Rev. Lett. 65 923 
CrossRefPubMed

[6]

Lee P A and Ramakrishnan T V 1985 Rev. Mod. Phys. 57 287 
CrossRef

[7]

Parisi G 1979 Phys. Rev. Lett. 43 1754 
CrossRef

[8]

Mézard M, Parisi G, Sourlas N, Toulouse G and Virasoro M 1984 Phys. Rev. Lett. 52 1156 
CrossRef

[9]

Bray A J and Moore M A 1986 Heidelberg Colloquium on Glassy Dynamics and Optimization ed L Van Hemmen and I Morgenstern (New York: Springer) p 121 

[10]

Fisher D S and Huse D A 1986 Phys. Rev. Lett. 56 1601 
CrossRefPubMed

[11]

Huse D A and Fisher D S 1987 J. Phys. A: Math. Gen. 20 L997 
IOPscience

[12]

Newman C and Stein D L 1996 Phys. Rev. Lett. 76 515 
CrossRefPubMed

[13]

Krzakala F and Martin O C 2000 Phys. Rev. Lett. 85 3013 
CrossRefPubMed

[14]

Palassini M and Young A P 2000 Phys. Rev. Lett. 85 3017 
CrossRefPubMed

[15]

Fisher M P A, Weichman P B, Grinstein G and Fisher D S 1989 Phys. Rev. B 40 546 
CrossRef

[16]

Damski B, Zakrzewski J, Santos L, Zoller P and Lewenstein M 2003 Phys. Rev. Lett. 91 080403 
CrossRefPubMed

[17]

Sanpera A, Kantian A, Sanchez-Palencia L, Zakrzewski J and Lewenstein M 2004 Phys. Rev. Lett. 93 040401 
CrossRefPubMed

[18]

Roth R and Burnett K 2003 J. Opt. B: Quantum Semiclass. Opt. 5 S50 
IOPscience

[19]

Roth R and Burnett K 2003 Phys. Rev. A 67 031602 
CrossRef

[20]

Bloch I, Dalibard J and Zwerger W 2007 Many-body physics with ultracold gasesPreprint 0704.3011
Preprint

[21]

Lewenstein M, Sanpera A, Ahufinger V, Damski B, Sen A and Sen U 2007 Adv. Phys. 56 243 
CrossRef

[22]

Jaksch D and Zoller P 2004 Ann. Phys. 315 52 

[23]

Edwards S F and Anderson P W 1975 J. Phys. F: Met. Phys. 5 965 
IOPscience

[24]

Stöferle T, Moritz H, Schori C, Köhl M and Esslinger T 2004 Phys. Rev. Lett. 92 130403 
CrossRefPubMed

[25]

Sebby-Strabley J, Anderlini M, Jessen P S and Porto J V 2006 Phys. Rev. A 73 033605 
CrossRef

[26]

Fölling S, Trotzky S, Cheinet P, Feld M, Saers R, Widera A, Müller T and Bloch I 2007 Nature 448 1029 
CrossRefPubMed

[27]

Spielman I B, Phillips W D and Porto J V 2007 Phys. Rev. Lett. 98 080404 
CrossRefPubMed

[28]

Campbell G K, Mun J, Boyd M, Medley P, Leanhardt A E, Marcassa L G, Pritchard D E and Ketterle W 2006 Science 313 649 
CrossRefPubMed

[29]

Lye J E, Fallani L, Modugno M, Wiersma D S, Fort C and Inguscio M 2005 Phys. Rev. Lett. 95 070401 
CrossRefPubMed

[30]

Fort C, Fallani L, Guarrera V, Lye J E, Modugno M, Wiersma D S and Inguscio M 2005 Phys. Rev. Lett. 95 170410 
CrossRefPubMed

[31]

Clément D, Varón A F, Hugbart M, Retter J A, Bouyer P, Sanchez-Palencia L, Gangardt D M, Shlyapnikov G V and Aspect A 2005 Phys. Rev. Lett. 95 170409 
CrossRefPubMed

[32]

Clément D, Varón A F, Retter J A, Sanchez-Palencia L, Aspect A and Bouyer P 2006 New J. Phys. 8 165 
IOPscience

[33]

Schulte T, Drenkelforth S, Kruse J, Ertmer W, Arlt J, Sacha K, Zakrzewski J and Lewenstein M 2005 Phys. Rev. Lett. 95 170411 
CrossRefPubMed

[34]

Fallani L, Lye J E, Guarrera V, Fort C and Inguscio M 2007 Phys. Rev. Lett. 98 130404 
CrossRefPubMed

[35]

Billy J, Josse V, Zuo Z, Bernard A, Hembrecht B, Lugan P, Clément D, Sanchez-Palencia L, Bouyer P and Aspect A 2008 Nature 453 891 
CrossRefPubMed

[36]

Sanchez-Palencia L, Clément D, Lugan P, Bouyer P, Shlyapnikov G V and Aspect A 2007 Phys. Rev. Lett. 98 210401 
CrossRefPubMed

[37]

Sanchez-Palencia L, Clément D, Lugan P, Bouyer P and Aspect A 2008 New J. Phys. 10 045019 
IOPscience

[38]

Roati G, D'Errico C, Fallani L, Fattori M, Fort C, Zaccanti M, Modugno G, Modugno M and Inguscio M 2008 Nature 453 895 
CrossRefPubMed

[39]

Lee P J, Anderlini M, Brown B L, Sebby-Strabley J, Phillips W D and Porto J V 2007 Phys. Rev. Lett. 99 020402 
CrossRefPubMed

[40]

Gorshkov A V, Jiang L, Greiner M, Zoller P and Lukin M D 2008 Phys. Rev. Lett. 100 093005 
CrossRefPubMed

[41]

Vignolo P, Akdeniz Z and Tosi M P 2003 J. Phys. B: At. Mol. Opt. Phys. 36 4535 
IOPscience

[42]

Gavish U and Castin Y 2005 Phys. Rev. Lett. 95 020401 
CrossRefPubMed

[43]

Jaksch D, Bruder C, Cirac J I, Gardiner C W and Zoller P 1998 Phys. Rev. Lett. 81 3108 
CrossRef

[44]

Günter K, Stöferle T, Moritz H, Köhl M and Esslinger T 2006 Phys. Rev. Lett. 96 180402 
CrossRefPubMed

[45]

Ospelkaus S, Ospelkaus C, Wille O, Succo M, Ernst P, Sengstock K and Bongs K 2006 Phys. Rev. Lett. 96 180403 
CrossRefPubMed

[46]

Greiner M, Mandel O, Esslinger T, Hänsch T W and Bloch I 2002 Nature 415 39 
CrossRefPubMed

[47]

Kuhn R C, Miniatura C, Delande D, Sigwarth O and Müller C A 2005 Phys. Rev. Lett. 95 250403 
CrossRefPubMed

[48]

de Martino A, Thorwart M, Egger R and Graham R 2005 Phys. Rev. Lett. 94 060402 
CrossRefPubMed

[49]

Lugan P, Clément D, Bouyer P, Aspect A and Sanchez-Palencia L 2007 Phys. Rev. Lett. 99 180402 
CrossRefPubMed

[50]

Scalettar R T, Batrouni G G and Zimanyi G T 1991 Phys. Rev. Lett. 66 3144 
CrossRefPubMed

[51]

Rapsch S, Schollwöck U and Zwerger W 1999 Europhys. Lett. 46 559 
IOPscience

[52]

Byczuk K, Hofstetter W and Vollhardt D 2005 Phys. Rev. Lett. 94 056404 
CrossRefPubMed

[53]

Yukalov V I, Yukalova E P, Krutitsky K V and Graham R 2007 Phys. Rev. A 76 053623 
CrossRef

[54]

Krutitsky K V, Thorwart M, Egger R and Graham R 2008 Phys. Rev. A 77 053609 
CrossRef

[55]

Sengupta P and Haas S 2007 Phys. Rev. Lett. 99 050403 
CrossRefPubMed

[56]

Lugan P, Clément D, Bouyer P, Aspect A, Lewenstein M and Sanchez-Palencia L 2007 Phys. Rev. Lett. 98 170403 
CrossRef

[57]

Parisi G 1983 Phys. Rev. Lett. 50 1946 
CrossRef

[58]

Winkler K, Thalhammer G, Lang F, Grimm R, Hecker Denschlag J, Daley A J, Kantian A, Büchler H P and Zoller P 2006 Nature 441 853 
CrossRefPubMed

[59]

Paredes B, Verstraete F and Cirac J I 2005 Phys. Rev. Lett. 95 140501 
CrossRefPubMed

[60]

Sherrington D and Kirkpatrick S 1975 Phys. Rev. Lett. 35 1792 
CrossRef

[61]

Giamarchi T, Le Doussal P and Orignac E 2001 Phys. Rev. B 64 245119 
CrossRef

[62]

Vidal G 2003 Phys. Rev. Lett. 91 147902 
CrossRefPubMed

[63]

Menotti C, Trefzger C and Lewenstein M 2007 Phys. Rev. Lett. 98 235301 
CrossRefPubMed

[64]

Roscilde T and Cirac J I 2007 Phys. Rev. Lett. 98 190401 
CrossRefPubMed

[65]

Vidal G 2004 Phys. Rev. Lett. 93 040502 
CrossRefPubMed

[66]

Cohen-Tannoudji C, Dupont-Roc J and Grynberg G 1992 Atom Photon Interactions: Basic Processes and Applications  (New York: Wiley) 

Notes

Note8  Note that our overlap function is defined as the averaged squared deviation from the mean. That is why, while the overlap of the two Mott replicas is indeed perfect, our overlap function, defined via deviations, is zero.

Note9  Note that in some situations, quantum averaging can be exploited to assist in producing disorder averages, see [59].

Note10  For small systems it may be necessary to improve the statistics by a repetition of the measurement.

Note11  In case the occupation number distribution in each replica is not equal, each should be individually addressed and measured, e.g. using a magnetic-field gradient during the measurement.

Note12  These points can easily be determined analytically by equating energy shifts for different total particle numbers (as given in table 1) and solving for .