Physical replicas and the Bose-glass in cold atomic gases

We study cold atomic gases in a disorder potential and analyze the correlations between different systems subjected to the same disorder landscape. Such independent copies with the same disorder landscape are known as replicas. While in general these are not accessible experimentally in condensed matter systems, they can be realized using standard tools for controlling cold atomic gases in an optical lattice. Of special interest is the overlap function which represents a natural order parameter for disordered systems and is a correlation function between the atoms of two independent replicas with the same disorder. We demonstrate an efficient measurement scheme for the determination of this disorder-induced correlation function. As an application, we focus on the disordered Bose-Hubbard model and determine the overlap function within perturbation theory and a numerical analysis. We find that the measurement of the overlap function allows for the identification of the Bose-glass phase in certain parameter regimes.


I. 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,2,3], the superconducting-toinsulator 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,8,9,10,11,12,13,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,17,18,19,20,21,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 phase, the Mottinsulator, and a Bose-glass phase. While the superfluid phase exhibits a finite condensate fraction and is characterized by off-diagonal long-range order, the Bose glass is only distinguished from the incompressible Mottinsulator 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 al-α β FIG. 1: (Colour online) Schematic representation of the implementation of disorder replicas with the cold atomic gases toolbox: Two planes α and β with equal disorder realisations, illustrated here for the case of a disorder landscape introduced using a second particle species (black dots), in which additional probe atoms (blue/lighter dots) evolve. The planes can be combined to measure correlation functions, such as the Edwards-Anderson overlap function.
lowing 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 ex-pressed as the correlation between independently evolving systems with the same disorder landscape (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 these.
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,25,26,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 superfluid-Mott-insulator quantum phase transition [28].
Note that this measurement scheme for the overlap between system with the same disorder can be applied to any disordered one-or two-dimensional 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 one-dimensional numerical simulations. We find that the overlap function q exhibits a sharp crossover from the Mott-insulating phase with q ≈ 0 [29] to the Bose-glass phase with q ≈ p(1 − p), with p ∈ (0, 1), and thus makes it possible to distinguish the two phases. Because the superfluid 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 experi-mentally 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 [30,31,32,33,34] were performed with laser speckles that had a disorder correlation length larger than the condensate healing length. Localisation effects were thus washed out by interactions. As an alternative, superlattice techniqueswhich 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 [35]. Only very recently the Palaiseau group has managed to create speckle potentials on the sub-micron scale and directly observe Anderson localisations effects in a BEC released in a onedimensional waveguide [36] based on the theoretical predictions of [37,38]. The Florence group reported observations of localisations phenomena in quasi-disordered potentials in BEC of K 39 , which allows for complete control of the strength atomic interactions using Feshbach resonances [39]. 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 [40,41]. 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 [42,43]. While production of 2D planes with equal disorder realisations 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 II, 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 III, we describe in detail the preparation of the system and the subsequent measurement scheme for the overlap function. In section IV, we determine the overlap for the disordered Bose-Hubbard model via perturbation theory for physicallyrelevant limits and compare the result with numerical simulations of the one-dimensional system. Details of the calculations are presented in the appendices. 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 [30,31,32,35,42,43] 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 short-range disorder with the disorder in different wells independent of each other. The probability distribution P (δ) and its correlation length depend 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 [30,31], 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 [40,41]. Alternatively, a disorder potential can also be created in mixtures of cold atomic gases in optical lattices [45,46] 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, whilst 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 give 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 superfluid, the Mott-insulator, and the Bose-glass phases. In the twodimensional regime of interest here, the superfluid phase appears for large hopping energies J U, ∆ and is characterized by off-diagonal (quasi) long-range order (finite superfluid 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 Mott-insulator phase with an excitation gap and a commensurate filling factor, i.e., the averaged particle number in a single well n i ∈ N is an integer. While the quantum phase transition from the superfluid to the Mott-insulator has been experimentally identified [47], the disorder potential gives rise to an additional Bose-glass 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 [43,48,49,50] and in other regimes via numerical methods such as quantum Monte Carlo [51], density matrix renormalization group (DMRG) [52], dynamical mean-field theory [53], and analytical methods [54,55]. Furthermore, recent interests focused also 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 [56] while random onsite interactions can give rise to a Lifshits glass [57]. The superfluid 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 Mott-insulating phase is characterized by the disappearance of these interference patterns and as well as a change in the behaviour of excitations [24,47]. However, the identification of the Bose-glass phase in the presence of disorder requires an additional observable allowing for the distinction between the Mottinsulator and the Bose-glass, where the condensate fraction vanishes. Such an additional property is known from 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 [58]. 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 and 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 [58].

III. 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: (1) initial preparation of disorder replicas, (2) introduction of probe atoms and (3) 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.

A. Initial preparation of disorder replicas
We consider a situation in which atoms are confined to a three-dimensional 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 tunneling 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 realisation.
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 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 [59]. 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 realisations.

B. 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 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][4][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 III B 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 that 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 III B 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.

Quantities that need to be measured
We need to measure two quantities to determine the overlap q, namely [ q 1 ] av , wherê with α = β and L is the total number of lattice sites, corresponding to the density-density correlations between the two layers, and 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 realisations of the disorder [60]. In general, one can expect that a measurement on 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 [61]. 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. First we focus on measurement of [ q α 2 ] av [ q β 2 ] av , i.e., the second term of the overlap in (4). If the layers are prepared as discussed above, the symmetry (1/L) i n α i = (1/L) i n β i is preserved and thus we can measure the average over two replicas, [ q 2 ] 2 where [62] 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 [ q 1 ] i.e., the first part of the overlap in (4). In this case the densitydensity 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 III B 2) and subsequently applying a similar scheme to the one presented in [28] (described in detail in section III B 3). Using the following equations that relate both p n and p α,β n to p nm which in conjunction with p 3 = p 12 + p 21 and p 4 = p 22 gives all necessary terms. Then, (12) simplifies to Finally, note thatq 2 as given in (10), can also be expressed as q 2 = 1 2 n np n . 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 III B 2 we describe the process of combining the two layers into a single one whilst in section III B 3 we describe a similar scheme to the one presented in [28] that can be used to distinguish different occupation numbers.

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 U kl = 2ω ⊥ a dz|w k (z)| 2 |w l (z)| 2 , 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 0 (x) and w 0 (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 nondegenerate 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.

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 III B 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 observed the Mott-insulating shell structure that arises due to a harmonic trapping potential commonly present in cold atomic gases experiments [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 "final" state (one particle of species b and all others 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 the different particle numbers. Since particles in the combined layer can be spread across multiple motional states (see section III B 2 above), different energy shifts occur for different configurations of the same number of particles. Whilst 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 Hamil-tonian (within a rotating-wave approximation) (15) 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 Ω(t)dt π, 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 a k = a † k a k and n b k = b † k b k . We further make the assumption 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 01 = U ij 00 and U ij 11 = (3/4)U ij 00 . 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 0 , N a 1 ) the number of particles of species a, 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 00 = − 1 and ∆E/U aa 00 = 3 4 ( − 1), where we have defined ≡ U ab 00 /U aa 00 . 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 |1, 1 ± ≡ (1/ √ 2) (|0, 1; 1, 0 ± |1, 0; 0, 1 ), with eigenenergies E + /U aa 00 = and E − /U aa 00 = 0. The associated energy shifts for these two possible final states are ∆E/U aa 00 = − 1 and ∆E/U aa 00 = −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 el- where Ω f i ≡ f |H RC |i , while for the state |1, 1 − we find Ω f i = 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 Ω f i . We note that the energy shifts, ∆E, for the different final states to which the initial states couple (i.e., Ω f i = 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, 0 or |0, N a 1 ; 0, 0 , the corresponding final states are given by |N a 0 −1, 0; 1, 0 and |0, N a 1 − 1; 0, 1 . The corresponding energy shifts are then given by ∆E/U aa 00 = (N a 0 − 1)( − 1) and ∆E/U aa 00 = 3 4 (N a 1 − 1)( − 1), and the matrix coupling element in each case is Ω f i = N a 0 and Ω f i = N a 1 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 Ω f i can be directly determined. For a total initial number of three and four particles the expressions for |m, n ± , E (m,n) ± , ∆E and Ω f i are given in appendix A 1 and A 2, respectively, for all the different initial configurations. In table I 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 Ω f i and the third column lists the associated energy shifts ∆E.
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) [63]. 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 Ω f i and the energy shifts ∆E in columns four and five, respectively, in table I. 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 t n = ∆t/ √ n at all resonant frequencies in table I with total number of particles n (ignoring the resonances with small and vanishing couplings Ω f i ). The total number of transferred particles is then proportional to p n and efficiently avoids problems arising from nonadiabatic combination of the layers.

IV. 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 distribu- I: 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 column we list the matrix coupling elements Ω f i 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 Ω f i and ∆E for the specific value = 0.98. tion 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 [64], 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 [58], or quantum systems with extended interactions [65]). For the latter type of system, measurement of the overlap function using physical replicas (see section III) is thought to yield different values for the overlap depending on α and β, the defining signature of replica symmetry-breaking.
In section IV A we analytically calculate the overlap function, for the various phases of the disorder Bose-Hubbard model (1), whilst in section IV B we present results for the overlap function from numerical simulation of the one-dimensional system.

A. Analytical determination of the overlap function
Starting with the limit of vanishing hopping J = 0, (1) reduces to an onsite Hamiltonian and the ground state is given by |Ω = Π i |Ω i with |Ω i the lowest energy state within each site i. This lowest energy states 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. First, for a chemical potential within the range U (n 0 − 1) + ∆ < µ < U n 0 − ∆ the above constraint in (18) is fulfilled independently of the site parameter i for the integer value n 0 , i.e., at each site the ground state is characterized by the integer particle density [ n i ] av = n 0 and we obtain a Mott-insulator (MI) phase (see figure 2). Furthermore, the overlap parameter vanishes identically in this limit, i.e., q = 0 for a Mott-insulator at J = 0.
Second, for chemical potentials in the range U n 0 −∆ < µ < U n 0 + ∆, 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 0 , while at sites with disorder potential ∆ i = −∆ the lowest energy state is characterized by n 0 + 1 particles. Therefore, the averaged particle density and the overlap parameter take the form This situation corresponds to the disorder-induced Boseglass (BG) phase (see figure 2) and the ground state is characterized by nonhomogeneous 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 determine the overlap function using 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 perturbation theory by computing Ω|ñ i |Ω whereñ i ≡ U † n i U . Details are given in appendix B.
We start with the Mott-insulator 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 0 (using |Ω i = |n 0 in the MI). To lowest order in perturbation theory the overlap is then given by with and |φ jk = |n 0 + 1 j |n 0 − 1 k Π l =ij |n 0 l . Using these expressions it is straightforward to compute the expectation value of the density correction Ã 2 = n 0 (n 0 + 1) where E 0 −E ij = −U +∆ j −∆ i . It is now straightforward to calculate the disorder averages and the overlap is given by q = 64(n 0 (n 0 + 1)) 2 p(1 − p)z(z + 1) where z is the number of nearest neighbours.
Next we turn to the Bose-glass 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 0 = p(1 − p) (the overlap at zero hopping),Ã 2 and S 1 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 0 ≥ 1 we obtain where and Once again a straightforward calculation of the disorder averages gives the following expression (n 0 ≥ 1) for the overlap: An analogous calculation for n 0 = 0 gives which corresponds to the U → 0 limit of (30). When the effects of hopping dominate, the bosonic atoms are in the superfluid phase (see figure 2). Then, the influence of weak disorder on the superfluid 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 ψ i = e iφ √ n 0 + δn i where n 0 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 0 +δn. We substitute this ansatz into H(ψ i ), expand ψ i for small fluctuations, δn i n 0 (see appendix C) and minimize the resulting expression in (C1) at fixed particle number, which is equivalent to requiring [δn i ] av = 0. Note that for large dimension 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 n 0 = (µ + Jz)/U + 1 2 .
In summary, we find that in the limit of small hopping J/U 1, the overlap signals a sharp crossover between the Mott-insulator with q ≈ 0 and the Bose-glass with q ≈ p(1 − q), see figure 2. In turn, the superfluid 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.

B. 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 [66,67]. 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 replicasymmetry 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 n is the density of the system. Computation of the overlap function can be simplified in this way as L −1 i n i self-averages to [ n i ] av for sufficiently large systems, which in turn self-averages to n.
The ground states used in (34) are computed by imaginary-time evolution, c.f. [67], 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 (Mott-insulator or Bose-glass regime). Regarding the number of disorder realizations for each value of ∆, we find that whilst 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 is figure 6 would suggest.
As shown in figure 6, we generally found very good agreement between numerical simulations and perturbation theory in its region of validity. For the Mottinsulator 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 Bose-glass regime, breakdown of perturbation theory sets in for low ∆, as degeneracy occurs there for ∆ = 0.
FIG. 6: Overlap function q as a function of impurity strength ∆ in the Mott-insulating phase at filling factor n0 = 1 (black lines) and in the Bose-glass phase for n0 = 0 (grey lines). In both cases, the impurity probability is p = 0.5 and U/J = 10, for L = 64. Simulation results (solid lines) are shown versus results from perturbation theory (dashed lines). Much of the variation in the simulated overlap function is due to the fact that we allow the number of impurities to fluctuate around the value given by Lp in each particular disorder realization. Note that the perturbation result in the Mott-insulator diverges as ∆ approaches U due to energy degeneracy [c.f. (24)]. In the Bose-glass phase, second-order perturbation theory yields a divergence as ∆ approaches 0, again due to energy degeneracy [c.f. (31)].

V. CONCLUSION
We have shown how the Bose-glass 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 explained in detail how to determine the occupation number distribution in this situation. Using perturbation theory we have calculated the overlap for weak hopping and have shown the different behaviours exhibited in the Mott-insulator and Bose glass phases; these results are in good agreement with the results from a one-dimensional 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 [58,65] would provide interesting insights to 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 [58], i.e., the results for the overlap depend on the particular replica. Our method could also give novel insights into physics of systems that possess multiple metastable states, such as ultracold dipolar gases [68] or quantum emulsions [69]; in particular statistics of overlaps of such metastable states could be measured. In addition to equilibrium properties considered here, the dynamics of the overlap might also posses an interesting behaviour, and could be measured in a similar setup. In one-dimension the time evolution of the overlap could be determined numerically [67]. 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 III B 3 and will not be discussed here further.
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 Ω f i = √ 3, whilst for the state |2, 1 − we find Ω f i = 0. The energy shifts for the two different possible processes are given by ∆E/U aa 00 = 2( − 1) and ∆E/U aa 00 = 2 − 2. 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 00 = E (1,2) ± − 3.

APPENDIX B: PERTURBATION THEORY
In this section we derive the corrections to the overlap function due to small hopping in perturbation theory following closely the method presented in [70]. We start with the zero hopping Hamiltonian H 0 with ground state |Ω and eigenenergy E 0 . We treat the hopping Hamiltonian, which has eigenstates denoted by |φ ij and eigenenergies denoted by E ij in perturbation theory. We proceed by determining a unitary transformation that transforms the total Hamiltonian H = H 0 + JH 1 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 1 + J 2 S 2 + · · · . 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 Ω|ñ i |Ω whereñ i ≡ U † n i U . 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 1 , are given by

Mott-insulator
In this section we derive the leading order correction to the overlap function for small hopping in perturbation theory for the MI phase. We begin by expanding S to fourth order and calculate the average corrected density to be Note that because n i = n 0 at all lattice sites (independent of the disorder) the second-and third-order contributions cancel.
Using these expressions it is straightforward to compute the expectation value of the density correction A 2 = n 0 (n 0 + 1)

Bose-glass
In this section we derive the leading order correction to the overlap function for small hopping in perturbation theory for the BG phase. Again we consider the corrected density ñ i as given in (B2) 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 0 = [ n i 2 ] av − [ n i ] 2 av is the overlap at zero hopping derived in section IV A. Since n i is now site dependent the second-order contribution does not vanish.
To determine the required commutators [S 1 , [S 1 , n i ]] we again write S 1 and n i as projectors 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 0 ≥ 1 we obtain and

APPENDIX C: MEAN FIELD THEORY
We replace the annihilation operators b i in the Hamiltonian in (1) by the mean field ψ i = e iφ √ n 0 + δn i and expand the fluctuations δn i n 0 up to second order to obtain Next, we minimize (C1) at fixed particle number which is equivalent to requiring [δn i ] av = 0. Using (32) and consistently solving gives Finally, then the overlap can be expressed in terms of these fluctuations as follows Substituting (C2) into (C3) and performing the disorder average for a bimodal distribution gives the result in (33) for the overlap function q.