text:
# Indication of a Co-Existing Phase of Quarks and Hadrons in Nucleus - Nucleus Collisions

## Abstract

The variation of average transverse mass of identified hadrons with charge multiplicity have been studied for AGS, SPS and RHIC energies. The observation of a plateau in the average transverse mass for multiplicities corresponding to SPS energies is attributed to the formation of a co-existence phase of quark gluon plasma and hadrons. A subsequent rise for RHIC energies may indicate a deconfined phase in the initial state. Several possibilities which can affect the average transverse mass are discussed. Constraints on the initial temperature and thermalization time have been put from the various experimental data available at SPS energies.

Results based on QCD (Quantum Chromodynamics) renormalization group approach predict that strongly interacting systems at high temperatures and/or densities are composed of weakly interacting quarks and gluons collins ; kisslin  due to asymptotic freedom and the Debye screening of color charge. Nucleus-Nucleus (A-A) collisions at very high energies may create situations conducive for the formation of a thermodynamic state where the properties of the system are governed by quarks and gluonic degrees of freedom. Such a state is called quark gluon plasma (QGP). QCD lattice gauge theory suggest that the critical temperature ($`T_c`$) for such a transition is $`170`$ MeV lattice . This has led to intense theoretical and experimental activities in this field of research qm .

Experimental detection of the QGP in A-A collisions is a non-trivial task because of the small space-time volume of the system. Among many signals of the QGP formation probe , one of the earliest is based on the relation of the thermodynamical variables, temperature and entropy to the average transverse momentum and multiplicity, respectively. This was originally proposed by Van Hove in the context of proton-proton collisions vanhove . It was argued that a plateau in the transverse momentum beyond a certain value of multiplicity will indicate the onset of the formation of mixed phase of QGP and hadrons; analogous to the plateau observed in the variation of temperature with entropy in a first order phase transition scenario. One hence looks for the variation of average transverse momentum ($`p_T`$) or transverse mass ($`m_T`$, $`m_T=\sqrt{p_T^2+m^2}`$) with respect to the total number of particles produced per unit rapidity ($`dN/dY`$) in A-A collisions at high energies. An increase in energy density should increase $`m_T`$ till it reaches a critical density for phase transition to occur, where it should then show a plateau due to formation of mixed phase. Further increase of energy density should again increase the $`m_T`$ of produced particles.

Several attempts have been made so far to look for such signals. While the data from cosmic ray experiments cosmic  are inconclusive, the data from A-A collisions in the laboratory for fixed center of mass energies ($`\sqrt{s}`$) with different centralities (impact parameters) wa98 ; rhic1 ; rhic2  have so far not shown this kind of behaviour. This implies that a mere change in the centrality of the collisions does not change the energy density of the system formed after the collisions required to create any change of phase. So it is imperative to study the variation of $`m_T`$ with $`dN_{ch}/dY`$ for a broad range of beam energies for fixed centrality.

In this letter the variation of $`m_T`$ with $`dN_{ch}/dY`$ is examined for AGS, SPS and RHIC energies spanning $`\sqrt{s}`$ from 2A GeV to 200A GeV. We then carry out explicit theoretical calculations to understand the observed behaviour in terms of the properties of the matter produced in the initial stages.

In fig. 1 the variation of $`m_T`$ with charge multiplicity is depicted for pions, kaons and protons at AGS, SPS and RHIC energies around mid-rapidity. The data shown here correspond to central events for different colliding systems, (produced) particle types and center of mass energies (see Table 1). The experimental data on transverse mass spectra can be parametrized as,

$$\frac{dN}{m_Tdm_T}C\mathrm{exp}\left(\frac{m_T}{T_{eff}}\right),$$
(1)
where the inverse slope parameter $`T_{eff}`$ is the effective temperature (effective because it includes the contribution from both the thermal and collective motion in the transverse direction). The variation of $`T_{eff}`$ with $`\sqrt{s}`$ for kaons has recently been shown in Ref marek . The average transverse mass of the particles obtained from Eq. 1 is:

$$m_T=T_{eff}+m+\frac{(T_{eff})^2}{m+T_{eff}}.$$
(2)
From the results shown in fig 1, one observes an increase in $`m_T`$ with $`dN_{ch}/dy`$ for AGS energies followed by a plateau for charge multiplicities corresponding to SPS energies for all the particle types, pions, kaons and protons. This may hint at the possible co-existence of the quark and hadron phases. For charge multiplicities corresponding to RHIC energies, the $`m_T`$ shows an increasing trend indicating the possibility of a pure QGP formation.

It is essential now to investigate whether the above experimental variation of $`m_T`$ can arise from various physical effects other than due to QGP formation. We will address this issue by trying to answer the following questions : (a) How does (3+1) dimensional hydrodynamical evolution of the system, with and without QCD phase transition, formed in heavy ion collisions affect the $`m_T`$ at freeze-out ? (b) Does the analysis of experimental data on hadrons, photons and dileptons indicate similar values of the initial temperature ($`T_i`$)? How does the value of $`T_i`$ compare with the $`T_c`$ predicted by lattice QCD lattice ? (c) How does the $`m_T`$ get affected by the gain in transverse momentum by hadrons through their successive collisions (Cronin effect cronin ) with the particles in the system ? (d) How does the $`m_T`$ of produced particles vary with $`\sqrt{s}`$ in a standard event generator based on the principle that nucleus-nucleus collision is a superposition of nucleon-nucleon collisions ?

(a) To understand the variation of $`m_T`$ with $`dN_{ch}/dY`$ for a system undergoing a (3+1) dimensional expansion we solve the hydrodynamical equations for various initial energy densities ($`ϵ_i`$hvg . The initial radial velocity is taken as zero. The momentum distribution of particles originating from a system undergoing transverse expansion hvg  with boost invariance along the longitudinal direction jdb , contains the effect of thermal motion as well as transverse flow. Hence it is imperative to examine how the value of the inverse slope or average transverse mass is affected by transverse motion. To achieve this we proceed as follows. We evaluate the $`m_T`$ spectra for various hadrons at the freeze-out point as a function of $`ϵ_i`$ for two different kind of space time evolution scenarios.

Firstly, we consider a hadronic gas at various values of $`ϵ_i`$ from $`0.05`$ to 4 GeV/fm3 without introducing a QCD phase transition. For all these initial energy densities the (3+1) dimensional hydrodynamic equations are solved to obtain the freeze-out surface at a temperature of 120 MeV. This is used as an input to evaluate the transverse mass spectra of hadrons. The hadronic equation of state (EOS) bm  contains hadronic degrees of freedom up to strange sector pdg . From the $`m_T`$ spectra of hadrons (pions, kaons and protons) corresponding to each $`ϵ_i`$ we evaluate the $`m_T`$. This is plotted in Fig. 2 (dashed line).

In the second scenario, we introduce the QGP (composed of up, down, strange quarks and gluons) at an energy density above $`ϵ_Q1.7`$ GeV/fm3 corresponding to critical temperature, $`T_c=170`$ MeV which converts to the hadronic phase at an energy density, $`ϵ_H0.56`$ GeV/fm3 at the same temperature. The effective degeneracy ($`g_{eff}`$) at $`T_c`$ for the hadronic phase is $`16`$ obtained from the study of hadronic EOS in bm . So in this scenario we have a mixed phase within the energy range $`0.56<ϵ`$(GeV/fm3)$`<\mathrm{\hspace{0.17em}1.7}`$. For $`ϵ>1.7`$GeV/fm3 a QGP bag model EOS has been used. The variation of $`m_T`$ as a function of $`ϵ_i`$ for this scenario is shown in Fig. 2 (solid line). It is clear that $`m_T`$ varies at a much slower rate for the second scenario. Although a complete plateau structure is not observed, it can be argued that the formation of a mixed phase slows down the rate of increase of $`m_T`$ with multiplicity. In fact, the slope of the curve in the case of pions for the first scenario is larger than the second by a factor of 2.5 at $`ϵ_i`$ = 1 GeV/$`fm^3`$. The $`m_T`$ for pions is about 400 MeV for values of $`ϵ_i`$ corresponding to mixed phase, which is similar to those obtained experimentally.

(b) The lower bound on the $`T_i`$ realized at SPS energies can be obtained from the measured $`m_T`$ provided we can estimate the effects of the flow. The inverse slope of the $`m_T`$ spectra contains the effect of flow and can be parametrized as uh ,

$$T_{eff}=T\sqrt{\frac{1+v_r}{1v_r}}$$
(3)
where $`v_r`$ is the average velocity and $`T=T_F`$ is the “true” freeze-out temperature. We start by noting that $`m_T`$ is actually the average energy, $`E`$ of the particle at zero rapidity ($`E=m_Tcoshy`$) and in a thermal system the average energy is proportional to temperature. Therefore, for a system which has zero transverse velocity but non-zero longitudinal velocity one must have

$$m_T_i>m_T_f,$$
(4)
because of the loss of thermal energy due to work done in longitudinal expansion. Here the subscript $`i`$ ($`f`$) refers to initial (freeze-out) state. Assuming $`m_T_f=m_T_{expt}`$, (i.e. all the hadrons are emitted from the freeze-out surface) and using Eq’s. 2 and 4 we get

$$T_i>m_T_{expt}/2.$$
(5)
A value of $`v_r0.4`$ reproduces the slope of the $`m_T`$ distribution of pions at SPS energies reasonably well. For this value of $`v_r`$ we get $`T_F120`$ MeV from Eq. 3, indicating, $`m_T=m+T_F+T_F^2/(m+T_F)315`$ MeV. It should be noted here that the effects of the transverse flow is “subtracted” out from $`m_T`$ now. Using Eq. 5 we obtain $`T_i>158`$ MeV.

Photons and dileptons emitted from the strongly interacting matter constitute one of the most efficient probes to measure the initial temperature of the system probe . It was shown earlier interf  that direct photons (measured by WA98 collaboration wa98\_photon ) having transverse momentum in the range $`1<p_T`$ (GeV)$`<\mathrm{\hspace{0.17em}2.5}`$ is expected to originate from a thermal source. The inverse slope of this spectra is $``$ 275 MeV. For $`v_r0.4`$ we obtain an average temperature $`T_{av}178`$ MeV by using Eq. 3. A similar value of $`T_{av}`$ is obtained from the $`p_T`$ distribution of $`e^+e^{}`$ ceres . Since photons and dileptons are emitted from all the stages of space time evolution, $`T_{av}`$ satisfies the condition, $`T_i>T_{av}>T_F`$. All these indicate that a value $`T_i170`$ MeV may have been realized at SPS energies.

The total multiplicity of hadrons is related to the initial temperature ($`T_i`$) and thermalization time ($`\tau _i`$) as:

$$\frac{dN}{dY}4\frac{\pi ^2}{90}g_{eff}\pi R^2T_i^3\tau _i/3.6,$$
(6)
$`R`$ is the radius of the colliding nuclei and $`g_{eff}`$ is the effective statistical degeneracy. Taking $`dN/dY700`$ corresponding to the highest SPS energies and using Eq. 5 we get $`\tau _i1.2`$ fm/c $`1/\mathrm{\Lambda }_{QCD}`$ where $`\mathrm{\Lambda }_{QCD}`$ is the QCD scale parameter.

(c) Recently, there have been some attempts to explain the transverse mass spectra of produced particles without introducing a transverse expansion. The hadrons gain transverse momentum through successive collisions with the particles in the system. This can be treated as a random walk of the hadrons hs  in the medium. The $`m_T`$ spectra of pions at SPS energies can be reproduced by this model, but it fails to reproduce the slope of kaon and proton distribution at the same colliding energies.

(d) As a last consideration, we investigate how the $`m_T`$ of produced particles vary with $`\sqrt{s}`$ of A-A systems in a standard event generator based on the principle that the A-A collision is a superposition of nucleon-nucleon collisions. This is done by using the event generator, HIJING hijing  which contains the effect of mini-jets. The results of this calculation do not show the variation of $`m_T`$ with $`\sqrt{s}`$ as has been observed experimentally.

In summary, a detailed analysis of the experimental data of the transverse mass spectra of identified hadrons, at AGS, SPS and RHIC energies, show a characteristic behaviour of $`m_T`$ as a function of multiplicity similar to that proposed by Van Hove. The plateau in $`m_T`$ corresponding to the multiplicity at SPS energies may hint at the formation of a mixed phase of QGP and hadrons in the initial stages. The average temperature obtained from the analysis of photons and dileptons lends support to this. We have shown by solving the hydrodynamical equations that the presence of a mixed phase really slows down the growth of $`m_T`$ with $`dN/dY`$ substantially. In the present work all the hadronic degrees of freedom up to mass 2.5 GeV has been included in the EOS; a fewer hadronic degrees of freedom will slow down this growth even more. Since a lower bound of the initial temperature is obtained here, the formation of QGP at a higher temperature cannot be ruled out. However, even if QGP phase is created at SPS, its space time volume could be too small to affect the dynamics of the system and hence our conclusions. The upward trend in the slope of $`m_T`$ at RHIC energies is indicative of a deconfined phase in the initial state.

The authors are thankful to Bikash Sinha for fruitful discussions. One of us (B.M.) is grateful to the Board of Research on Nuclear Science and Department of Atomic Energy, Government of India for financial support.


format: text

text:
# Phase Diagram of Bosonic Atoms in Two-Color Superlattices

## I Introduction

With the first experimental observation of the superfluid to Mott-insulator transition in ultracold atomic Bose gases trapped in optical lattice potentials GrMa02a ; GrMa02b  the unique prospects of this new class of systems became evident. Compared to experiments with Bose-Einstein condensates in shallow magnetic traps, the use of optical standing-wave lattices gives access to a new regime which is dominated by strong correlations. The most striking ramification of these correlations are quantum phase transitions Sach99  which are driven by the interplay of the different contributions to the Hamiltonian of the system.

The superfluid to Mott-insulator transition is governed by competition between kinetic energy and repulsive atom-atom interaction: in the presence of strong repulsive interactions the ground state is almost a pure Fock state with a definite occupation number at the individual lattice sites. This is associated with a vanishing of the superfluid density which is the natural order parameter for superfluid to insulator phase transitions JaBr98 ; RoBu03a .

New quantum phase transitions can be expected if other terms in the Hamiltonian come into play. One interesting case is the presence of some additional spatial structure of the lattice potential. This might range from a regular modulation of the well depths to random disorder. From several studies in the context of solid-state systems it is known that random disorder gives rise to new zero-temperature phases such as the Anderson localized phase and the Bose-glass phase ScBa91 ; KrTr91 .

Cold atoms in optical lattices offer unparalleled possibilities to study these disorder-related phenomena RoBu03b . By a superposition of multiple laser beams it is possible to generate a huge variety of different lattice topologies in a perfectly controlled manner. At the same time Feshbach resonances can be utilized to modify the strengths of the interatomic interaction. These powerful methods of tuning the relevant parameters are supplemented by versatile techniques to probe the state of the many-boson system, e.g. by observing the matter-wave interference pattern after the atoms were released from the lattice GrMa02a ; OrTu01  or by measuring the energy gap between the ground state and the first excited state GrMa02b . Atoms in an optical lattice are therefore an ideal model system for the experimental investigation of the fundamental questions associated with quantum phase transitions in lattice systems.

The aim of this paper is twofold: Firstly, we want to review the important physical quantities which describe the properties of Bose gases in optical lattices and investigate the relations between them, both, analytically and numerically. The emphasis in this part of the paper is on the superfluid fraction and its connection to directly observable quantities. Secondly, we want to systematically explore the zero-temperature phase diagram for bosons in two-color superlattices.

The theoretical framework for the description of strongly correlated Bose gases in lattices is the Bose-Hubbard model JaBr98 , which is summarized in Sec. II. We define the fundamental quantities, condensate fraction, quasimomentum distribution, and superfluid fraction, and discuss their relation on the formal level in Sec. III. The connection of these fundamental properties to the experimental observables, e.g., the matter-wave interference pattern, is investigated in Sec. IV. The formal definition shows that the superfluid fraction depends crucially on the excitation spectrum RoBu03a . The basic experimental observable, i.e. the interference pattern, only probes the ground state of the system and cannot therefore provide full information on the superfluid properties of the system. Using the example of the superfluid to Mott-insulator transition we illustrate these differences in the exact numerical solution of the Bose-Hubbard model for a one-dimensional uniform lattice. Finally, in a Sec. V we extend our studies to non-uniform lattice potentials. The emphasis is on so-called two-color superlattices which exhibit a sinusoidal variation of the well depths RoBu03b . We map out the zero-temperature phase diagram as function of the interaction strength and the amplitude of the modulation and compare the behavior of number fluctuations, condensate fraction, superfluid fraction and energy gap. This enables us to identify additional zero-temperature phases, i.e. a localized phase and a quasi Bose-glass phase, and to discuss possibilities for their experimental detection.

## II Bose-Hubbard model

Quantum phase transitions of Bose gases in optical lattices are associated with complex correlations in the many-body state which go far beyond a simple Gross-Pitaevskii-like description. A theoretical model which is capable of describing these dominating correlations is the Bose-Hubbard model. First introduced for model studies related to 4He liquids in porous media and granular superconductors FiWe89 ; BaSc90  it was recently applied to ultracold atomic Bose gases in optical lattices JaBr98 .

In the following we will briefly review the assumptions which lead to the Bose-Hubbard Hamiltonian and summarize solution methods and the basic observables. Because our emphasis is on conceptual and general aspects we restrict ourselves to one-dimensional lattices, the formal generalization to two or three dimensions is straightforward.

### II.1 Bose-Hubbard Hamiltonian

The basic assumption of the Bose-Hubbard model is that the lattice wells are sufficiently deep for the state of the system to be described using a basis of single-particle wavefunctions, localized at the individual lattice sites. Only these localized ground states are taken into account and excited vibrational states are neglected. In the language of band structure theory the model space of the Bose-Hubbard model comprises the lowest energy band only, all excited bands are excluded. The localized single-particles states are given by the Wannier functions for the lowest band. The energy gap between the ground state band and the first excited band has therefore to be sufficiently large for admixtures from the excited bands to be negligible.

A natural way to characterize the many-boson states in the model space of the Bose-Hubbard model is the occupation number representation. Let’s assume we are dealing with a system of $`N`$ bosons in a lattice composed of $`I`$ lattice sites. We introduce a set of occupation numbers $`\{n_1,\mathrm{},n_I\}`$ which specify the number of bosons in the localized single-particle state at the individual lattice sites. The set of Fock states $`|n_1,\mathrm{},n_I`$ for all possible compositions of the occupation numbers under the constraint $`_{i=1}^In_i=N`$ forms a complete basis the Bose-Hubbard model space.

We can define the associated creation and annihilation operators, $`\widehat{\mathrm{a}}_i^{}`$ and $`\widehat{\mathrm{a}}_i`$, for a boson localized at the $`i`$th lattice site:

$$\begin{array}{cc}\hfill \widehat{\mathrm{a}}_i^{}|n_1,\mathrm{},n_i,\mathrm{},n_I& =\sqrt{n_i+1}|n_1,\mathrm{},n_i+1,\mathrm{},n_I\hfill \\ \hfill \widehat{\mathrm{a}}_i|n_1,\mathrm{},n_i,\mathrm{},n_I& =\sqrt{n_i}|n_1,\mathrm{},n_i1,\mathrm{},n_I.\hfill \end{array}$$
(1)
The occupation number operator for the $`i`$th site is given by

$$\widehat{\mathrm{n}}_i=\widehat{\mathrm{a}}_i^{}\widehat{\mathrm{a}}_i.$$
(2)
Throughout this paper the indices $`i,j=1,\mathrm{},I`$ label the lattice sites. All operator-valued quantities are marked by a hat.

Using these creation and annihilation operators one can easily translate the many-body Hamiltonian of the system—consisting of the kinetic energy, the external lattice potential, and the two-body interaction—into its second quantized form. This procedure leads directly to the Bose-Hubbard Hamiltonian for the one-dimensional lattice system

$$\begin{array}{cc}\hfill \widehat{\mathrm{H}}_0=& J\underset{i=1}{\overset{I}{}}(\widehat{\mathrm{a}}_{i+1}^{}\widehat{\mathrm{a}}_i+\text{h.a.})\hfill \\ & +\underset{i=1}{\overset{I}{}}ϵ_i\widehat{\mathrm{n}}_i+\frac{V}{2}\underset{i=1}{\overset{I}{}}\widehat{\mathrm{n}}_i(\widehat{\mathrm{n}}_i1).\hfill \end{array}$$
(3)
The first term describes the tunneling between adjacent lattice sites characterized by a strength parameter $`J`$. It is associated with the kinetic energy part of the first quantized Hamiltonian. In general we will use periodic boundary conditions, i.e., tunneling between the first and the last lattice site is included (the site index $`I+1`$ which appears in the summation is implicitly replaced by $`1`$).

The second term in (3) gives an on-site single-particle energy which originates from the external potentials and the on-site part of the kinetic energy. The on-site energies $`ϵ_i`$ are constant for a translationally invariant lattice and might be set to zero. However, in the presence of an additional parabolic trapping potential or for irregular lattices the $`ϵ_i`$ vary with the site index $`i`$. Finally, the third term comprises the on-site two-body interaction characterized by the interaction strength $`V`$.

The parameters $`J`$, $`ϵ_i`$, and $`V`$ are given by the matrix elements of the terms of the first-quantized Hamiltonian the in the localized single particle states JaBr98 ; OoSt01 . In the following we will use the tunneling parameter $`J`$ as the unit of energy; the dimensionless ratios $`V/J`$ and $`ϵ_i/J`$ are varied to explore different regions of the zero-temperature phase diagram of the lattice system. As mentioned in the introduction one can envisage being able to tune the parameters of the Hamiltonian (3) experimentally, either by changing the lattice geometry or by exploiting a Feshbach resonance.

Besides the fundamental restriction of the Bose-Hubbard model to the ground state band some additional assumptions are used to construct the Hamiltonian (3): (*i*) Only tunneling between adjacent sites is included, long range hopping over several lattice sites is neglected. (*ii*) Only short-range on-site interactions are included; long-range interactions as they can occur for example in the presence of dipole-dipole forces are not considered here. (*iii*) The tunneling strength $`J`$ is assumed to be independent of the site index. This is obviously the case in regular lattice potentials, but it is an approximation for irregular lattices DaZa03 . If required all these simplifications can be easily abandoned.

### II.2 Exact solution

The most direct way to obtain the ground state and the excitation spectrum of a zero-temperature Bose gas in one-dimensional lattices is through the exact numerical diagonalization of the Bose-Hubbard Hamiltonian. Within the occupation number basis the construction of the matrix representation of the Bose-Hubbard Hamiltonian is trivial. The only substantial hurdle is the dimension of the number basis. For fixed number of sites $`I`$ and number of particles $`N`$ the dimension of the number basis is given by

$$D=\frac{(N+I1)!}{N!(I1)!}.$$
(4)
It grows dramatically with increasing system size: For fixed filling factor $`N/I=1`$ the basis dimension for $`I=8`$ is $`D=6435`$, for $`I=10`$ it grows to $`D=92378`$, and for $`I=12`$ it reaches $`D=1352078`$. The Hamilton matrix is, however, extremely sparse because only the tunneling term in (3) generates off-diagonal matrix elements. Efficient iterative Lanczos-type algorithms can therefore be utilized to determine the few lowest eigenvalues and the corresponding eigenvectors. This allows to treat systems with up to $`I=N=12`$ on a standard PC.

As the result of the numerical solution of the eigenvalue problem we obtain the energy eigenvalues $`E^{(\nu )}`$ and the corresponding eigenvectors $`C_\alpha ^{(\nu )}`$ for a few states ($`\nu =0,1,2\mathrm{}`$). The eigenvectors provide the coefficients of the expansion of the eigenstates in the number basis

$$|\mathrm{\Psi }_\nu =\underset{\alpha =1}{\overset{D}{}}C_\alpha ^{(\nu )}|\{n_1,\mathrm{},n_I\}_\alpha .$$
(5)
Here the index $`\alpha =1,\mathrm{},D`$ labels the different Fock states, i.e., the different sets $`\{n_1,\mathrm{},n_I\}_\alpha `$ of occupation numbers.

Besides the exact numerical diagonalization, which explicitly yields the interesting state vectors, it is possible to compute selected observables using Monte-Carlo simulations BaSc92  which allow to treat systems with several hundred lattices sites KaPr02 ; BaRo02 . Because of its simplicity and transparency we will, nevertheless, restrict ourselves to the exact diagonalization for relatively small lattices, which already exhibit the important universal properties of this class of systems.

Several approximation methods have been developed to solve the Bose-Hubbard model. One of them is the so-called mean-field approximation where the state is assumed to be a direct product of independent single-site states. Thereby complex correlations between the individual sites that can be of particular importance for the description of quantum phase transition cannot be described. On top of this independent-site approximation the Gutzwiller ansatz KrCa92  can be used to perform variational calculations of the ground state.

Another approximation scheme uses the Bogoliubov approach. Here the annihilation operators $`\widehat{\mathrm{a}}_i`$ are replaced by a complex amplitude $`z_i`$ plus a fluctuation operator $`\widehat{\delta }_i`$. The fluctuation part is assumed to be small such that terms involving squares of the fluctuation operators can be neglected. The application and the limitations of this approach are discussed in ReBu03 .

### II.3 Simple observables

From the ground state $`|\mathrm{\Psi }_0`$ and a few excited states obtained from the exact diagonalization of the Bose-Hubbard Hamiltonian one can extract several simple observables by computing expectation values.

The simplest observable is the mean occupation number of the different sites

$$\overline{n}_i=\mathrm{\Psi }_0|\widehat{\mathrm{n}}_i|\mathrm{\Psi }_0.$$
(6)
For a translationally invariant lattice $`\overline{n}_i`$ will be the same for all lattice sites and equal to the filling factor $`N/I`$ independent of the parameters $`V`$ and $`J`$. This, however, does not mean that the structure of the state does not change. A more telling quantity are the fluctuations around this mean occupation number given by

$$\sigma _i=\sqrt{\mathrm{\Psi }_0|\widehat{\mathrm{n}}_i^2|\mathrm{\Psi }_0\mathrm{\Psi }_0|\widehat{\mathrm{n}}_i|\mathrm{\Psi }_0^2}.$$
(7)
These number fluctuations provide direct information on the structure of the ground state (5). If $`|\mathrm{\Psi }_0`$ is a superposition of many Fock states $`|\{n_1,\mathrm{},n_I\}_\alpha `$ then many different occupation numbers occur at the individual sites and the number fluctuations $`\sigma _i`$ will be large. If, on the other hand, the ground state is described by a single Fock state, then the fluctuations will vanish.

As a complementary measure the magnitude of the largest coefficient in the Fock state expansion (5) of the ground state can be used

$$C_{\mathrm{max}}^2=\mathrm{max}(C_\alpha ^2).$$
(8)
Although this is not directly observable it yields sensitive information on the structure of the state. A small value of $`C_{\mathrm{max}}^2`$ indicates that $`|\mathrm{\Psi }_0`$ is a superposition of many number states; $`C_{\mathrm{max}}^21`$ means that the ground state is a pure Fock state.

To illustrate these quantities we solve the eigenvalue problem of the Bose-Hubbard Hamiltonian (3) for a regular lattice potential ($`ϵ_i=0`$) with filling factor $`N/I=1`$ for different ratios $`V/J`$ of the interaction strength and the tunneling coefficient. The results for the number fluctuations $`\sigma _i`$ and the largest coefficient $`C_{\mathrm{max}}^2`$ are depicted in Fig. 1. The vertical gray line indicates the critical interaction strength $`(V/J)_{\text{cr}}=4.65`$ for the superfluid to Mott-insulator transition in an infinite one-dimensional lattice extracted from a Monte Carlo calculation BaSc92  and from a strong coupling expansion FrMo96 . It marks the region where structural changes in the ground state are also expected in the finite size systems considered here.

Our first important observation is that the number fluctuations decrease rather slowly with increasing interaction strength. Moreover, there appears no clear signature of a phase transition in their variation. This is not an artifact of the limited lattice size, because the number fluctuation exhibit no noticeable size dependence as the three curves for $`I=12,10,`$ and $`8`$ in Fig. 1(a) illustrate. Thus we conclude that the number fluctuations and related quantities do not reveal the quantum phase transition from a superfluid to a Mott-insulator state and that the insulating state still has significant number fluctuations, i.e., is not a pure Fock state.

The last point is confirmed by the behavior of the largest coefficient $`C_{\mathrm{max}}^2`$ plotted in Fig. 1(b). This quantity is very small in the insulating phase and starts to grow in the region of the phase transition. However, even for large values of the ratio $`V/J`$ the largest coefficient in (5) remains smaller than $`1`$ indicating that there is a dominant Fock state (the one with $`n_i=N/I`$) but that other Fock states also contribute to the ground state.

Another important observable is the energy gap, i.e., the difference between the energies of the ground state and the first excited state,

$$E_{\text{gap}}=E^{(1)}E^{(0)}.$$
(9)
This quantity is of particular importance as the excited states are crucial for the superfluid properties of the system. Figure 2 depicts the behavior of the energy gap as function of the interaction strength. For small $`V/J`$ the energy gap is small and increases only slowly with increasing interaction strengths. The slope then starts to increase until, for $`V/J10`$, one observes a linear increase of the energy gap with a slope $`1`$, i.e., changing the interaction strength $`V`$ by a certain amount will change the energy gap by the same amount. Qualitatively, this behavior can be explained by assuming that the excited state is a superposition of different one-particle one-hole excitations of the ground state. The predominant Fock states in the expansion (5) of the excited state thus have one doubly occupied site which leads to the linear dependence of the energy gap $`E_{\text{gap}}`$ on $`V`$. These signatures of the Mott-insulator phase have been investigated experimentally GrMa02b  as well as theoretically SaSe02 ; BrDu02  by tilting the lattice potential.

The behavior of the energy gap indicates that the character of the energy spectrum changes in the region $`V/J5`$ to $`8`$. As we will see in Sec. III.2 this change is accompanied by a vanishing of superfluidity, i.e. the superfluid to Mott-insulator transition.

## III Condensate and superfluidity

Besides the simple quantities discussed so far there are two more complex physical quantities which reveal substantial information on fundamental physical properties of the Bose gas in the lattice. These are the condensate fraction and the superfluid fraction. This section is devoted to the definition and discussion of these non-trivial quantities.

### III.1 Condensate fraction

Unlike the Gross-Pitaevskii description of weakly interacting Bose gases DaGi99  the exact solution of the Bose-Hubbard model does not assume a perfect Bose-Einstein condensate from the outset. Moreover, due to the general representation of the ground state in terms of Fock states the presence or absence of a Bose-Einstein condensate is not obvious. For example, the exact ground state (5) of a noninteracting system at zero temperature, where one expects a perfect Bose-Einstein condensate, is a superposition of many Fock states with coefficients given by the multinomial distribution.

#### III.1.1 Definition

In order to define a Bose-Einstein condensate for the ground state $`|\mathrm{\Psi }_0`$ of a general interacting many-boson system in a lattice we consider the one-body density matrix associated with $`|\mathrm{\Psi }_0`$

$$\rho _{ij}^{(1)}=\mathrm{\Psi }_0|\widehat{\mathrm{a}}_j^{}\widehat{\mathrm{a}}_i|\mathrm{\Psi }_0.$$
(10)
The eigenvectors of the one-body density matrix describe the so-called natural orbitals and the eigenvalues the corresponding occupation numbers.

Following the formulation of Penrose and Onsager PeOn56  a Bose-Einstein condensate is present if one of the natural orbitals is *macroscopically* occupied. Its occupation number is just the number of condensate particles $`N_\text{c}`$ and the eigenvector constitutes the condensate wave function $`\chi _{\text{c},i}`$.

Strictly speaking, macroscopic occupation of an orbital means that the ratio of its occupation number and the total particle number, in the following called condensate fraction $`f_\text{c}=N_\text{c}/N`$, remains finite in the the thermodynamic limit ($`N,I\mathrm{}`$, $`N/I=\text{const.}`$). This implies that one cannot rigorously determine the *absence* of a Bose-Einstein condensate in finite-size systems. From the normalization of the one-body density matrix, $`\mathrm{Tr}\rho _{ij}^{(1)}=N`$, it follows immediately that for any system there is an eigenvalue larger or equal $`N/I`$. Thus the condensate fraction $`f_\text{c}`$ is always larger than $`1/I`$ and cannot vanish for finite-size systems. Despite of these formal complications the ratio of the largest eigenvalue of (10) and the total particle number provides important physical information even in finite lattices. For simplicity we will use the term condensate fraction for this quantity.

The presence of a Bose-Einstein condensate implies off-diagonal long range order Yang62 . That is, the matrix elements of the one-body density matrix do not go to zero far off the diagonal but remain finite:

$$\rho _{ij}^{(1)}0\text{for}|ij|\mathrm{}.$$
(11)
The prove of this connection is straightforward if one considers the spectral decomposition of the one-body density matrix which allows us to separate the condensate part from a residual positive semi-definite density matrix $`\stackrel{~}{\rho }_{ij}^{(1)}`$

$$\rho _{ij}^{(1)}=N_\text{c}\chi _{\text{c},i}^{}\chi _{\text{c},j}+\stackrel{~}{\rho }_{ij}^{(1)}.$$
(12)
It is important to realize that off-diagonal long range order is associated directly with the phenomenon of Bose-Einstein condensation whereas the connection with superfluidity is, at most, indirect.

The dependence of the condensate fraction, i.e., the relative occupation of the natural orbital with largest eigenvalue, on the interaction strength $`V/J`$ in a regular lattice is depicted in Fig. 3. In the noninteracting case $`V/J=0`$ all particles are in the condensate. As soon as a finite repulsive interaction is present ($`V/J>0`$) the condensate fraction is reduced, i.e., the condensate is depleted and non-condensate states are populated. For large interaction strengths the condensate fraction converges to $`1/I`$. This is the reason for the pronounced size dependence of the condensate fraction in the regime of strong depletion.

#### III.1.2 Quasimomentum Distribution

The eigenvalues of the one-body density matrix (10) provide more information than just the condensate fraction. For a translationally invariant lattice one can show that the natural orbitals satisfy the Bloch theorem and therefore are quasimomentum eigenstates JoMa85a . Thus the occupation numbers of the natural orbitals correspond to occupation numbers $`\stackrel{~}{n}_q`$ of the Bloch states with different quasimomenta $`q`$. Notice that for a finite lattice of length $`L`$ the quasimomenta can assume only discrete values which are integer multiples of $`2\pi /L`$. The Bloch function for quasimomentum $`q=0`$ describes the condensate state as defined above.

An alternative way to obtain the quasimomentum distribution for a regular lattice is to use the relation between the localized Wannier functions and delocalized Bloch functions, which are the quasimomentum eigenstates. The Bloch functions $`\psi _q(x)`$ of the lowest band can be decomposed in terms of Wannier functions $`w(x\xi _i)`$:

$$\psi _q(x)=\frac{1}{\sqrt{I}}\underset{i=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}q\xi _i}w(x\xi _i),$$
(13)
where $`\xi _i`$ is the coordinate of the $`i`$th lattice site. Using the fact that the creation operators $`\widehat{\mathrm{a}}_i^{}`$ of the Bose-Hubbard model create a boson in the Wannier state $`w(x\xi _i)`$ we can readily define creation operators $`\widehat{\mathrm{c}}_q^{}`$ for a boson in the Bloch state with quasimomentum $`q`$ OoSt01

$$\widehat{\mathrm{c}}_q^{}=\frac{1}{\sqrt{I}}\underset{i=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}q\xi _i}\widehat{\mathrm{a}}_i^{}.$$
(14)
The occupation numbers of the Bloch states with quasimomentum $`q`$ are thus given by

$$\begin{array}{cc}\hfill \stackrel{~}{n}_q=\mathrm{\Psi }_0|\widehat{\mathrm{c}}_q^{}\widehat{\mathrm{c}}_q|\mathrm{\Psi }_0& =\frac{1}{I}\underset{i,j=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}q(\xi _i\xi _j)}\mathrm{\Psi }_0|\widehat{\mathrm{a}}_j^{}\widehat{\mathrm{a}}_i|\mathrm{\Psi }_0\hfill \\ & =\frac{1}{I}\underset{i,j=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}qa(ij)}\rho _{ij}^{(1)},\hfill \end{array}$$
(15)
where $`a=\xi _{i+1}\xi _i`$ is the lattice spacing. As expected, the quasimomentum distribution is related to the Fourier transform of the one-body density matrix.

Figure 4 shows the occupation number for the three lowest quasimomentum states as function of the interaction strength $`V/J`$. The population for the $`q=0`$ state corresponds to the condensate fraction and is depleted as the interaction strength is increased. Simultaneously with this, the occupation numbers of Bloch states with higher quasimomenta, which vanish in the noninteracting system, increase successively. The repulsive interactions thus lead to a redistribution of the particles from the condensate to states with higher quasimomenta. In the limit of strong interactions all quasimomentum states exhibit the same occupation, i.e., the band is filled uniformly.

### III.2 Superfluid fraction

One of the most interesting physical quantities in an interacting many-boson system is the superfluid density. It is the natural order parameter for the various superfluid to insulator phase transitions, e.g., the superfluid to Mott-insulator transition observed experimentally GrMa02a .

#### III.2.1 Definition

The macroscopic, as well as the microscopic definition of superfluidity, are connected to flow properties. Macroscopically, one usually employs the two-fluid picture to describe systems which exhibit superfluidity. The superfluid and the normal fluid are distinguished through their behavior in the presence of moving boundaries. For example, in a narrow channel with moving side walls only the normal fluid is dragged along with the walls, whereas the superfluid stays at rest. In the rest frame of the walls only the superfluid moves. From the difference between the energies (in the rest frame of the walls) of a system with moving and one with stationary walls, one can directly extract the kinetic energy of the superfluid component and hence, via its velocity, the superfluid density PoCe87 ; LiSe02 .

A microscopic definition of superfluidity can be constructed in analogy to this simple macroscopic picture. First of all we consider the generalized condensate wave function $`\chi (\stackrel{}{x})=\mathrm{e}^{\mathrm{i}\vartheta (\stackrel{}{x})}|\chi (\stackrel{}{x})|`$ as defined in the previous section. A spatial variation of the phase $`\vartheta (\stackrel{}{x})`$ of the condensate wave function is connected to a velocity field Legg99

$$\stackrel{}{v}_\text{s}=\frac{\mathrm{}}{m}\stackrel{}{}\vartheta (\stackrel{}{x}).$$
(16)
The flow described by this velocity field is non-dissipative and irrotational, properties which conform with the macroscopic picture of superfluidity. This velocity is therefore identified with the flow field of a superfluid Legg99 . In principle Eq. (16) provides the basic microscopic distinction of the superfluid through its flow behavior.

As in the macroscopic case we can develop an explicit expression for the superfluid fraction by considering the total energy of the system. For simplicity we restrict ourselves to a one-dimensional system. Assume that we impose a linear phase variation with a total phase twist $`\mathrm{\Theta }`$ over the length $`L`$ of the system, i.e., $`\vartheta (x)=\mathrm{\Theta }x/L`$. One way to accomplished this is by imposing twisted boundary conditions on the many-body wave function

$$\mathrm{\Psi }(x_1,\mathrm{},x_k+L,\mathrm{},x_N)=\mathrm{e}^{\mathrm{i}\mathrm{\Theta }}\mathrm{\Psi }(x_1,\mathrm{},x_k,\mathrm{},x_N)$$
(17)
for all $`k`$ FiBa73 ; ShSu90 ; Poil91 . This imposed phase gradient generates a flow with a velocity $`v_\text{s}=\mathrm{}\mathrm{\Theta }/(mL)`$. The portion of the system which responds to the phase gradient by flowing with this velocity is called superfluid. The presence of the flow will increase the total energy of the system. As long as other excitations are absent, that is, as long as the imposed phase gradient is small, the increase in the total energy can be attributed solely to the kinetic energy of the superflow

$$E_\mathrm{\Theta }E_0=\frac{1}{2}M_\text{s}v_\text{s}^2.$$
(18)
Here $`E_0`$ is the ground state energy of the system without phase twist and $`E_\mathrm{\Theta }`$ is the ground state energy obtained with twisted boundary conditions. The mass $`M_\text{s}`$ corresponds to the total mass of the superfluid portion, by introducing the superfluid fraction $`f_\text{s}`$ is can be rewritten as $`M_\text{s}=mNf_\text{s}`$. From Eq. (18) we can extract a simple expression for the superfluid fraction FiBa73 ; Krau91 ; SiRo94

$$f_\text{s}=\frac{2mL^2}{\mathrm{}^2N}\frac{E_\mathrm{\Theta }E_0}{\mathrm{\Theta }^2}\text{for}\mathrm{\Theta }\pi .$$
(19)
Hence the superfluid fraction is determined by the stiffness of the system under phase variations. We have to stress that the twist angle $`\mathrm{\Theta }`$ has to be sufficiently small to avoid effects other than the collective flow of the superfluid component. We will return to this point at the end of this section.

The definition of the superfluid fraction through the energy change under imposed phase variations corresponds to the so-called helicity modulus introduced in Ref. FiBa73 . Moreover, it is equivalent to the concept of winding numbers which is used in path-integral Monte Carlo methods PoCe87 . All these measures predict the noninteracting Bose gas to be a perfect superfluid. This is in contradiction to the Landau picture of superfluidity which is based on the dispersion relation of elementary excitations which defines a critical velocity. In the Landau picture the ideal Bose gas is not considered a superfluid since its critical velocity is zero. A possible connection between these two pictures of superfluidity is still a matter of debate LiSe02 . We adopt a pragmatic view and use the definition (19) keeping in mind that it does not tell anything about the dynamical stability of the superfluid flow at finite velocities.

A second comment is appropriate concerning the relation between superfluidity and Bose-Einstein condensation. First of all one should realize that although the phase of the condensate wave function determines the superfluid velocity (16), the superfluid and the condensate fraction are not the same. A famous example is liquid 4He at zero temperature which is 100% superfluid but only 10% of the atoms are in the condensate Legg99 . From the formal point of view we used the condensate wave function only to motivate the connection between twisted boundary conditions (17) and the superfluid flow. On the level of the resulting Eqs. (17) and (19) the condensate does not appear. In general it is not strictly evident that superfluidity presupposes the presence of a condensate LiSe02 .

The relation (19) for the continuous system can be easily transfered to the Bose-Hubbard model. Basically one replaces the length $`L`$ of the continuous system by the number of lattice sites $`I`$ and the prefactor of the kinetic energy $`\mathrm{}^2/(2m)`$ by the tunneling strength $`J`$ of the Bose-Hubbard Hamiltonian. The superfluid fraction for the discrete lattice thus reads RoBu03a

$$f_\text{s}=\frac{I^2}{JN}\frac{E_\mathrm{\Theta }E_0}{\mathrm{\Theta }^2}\text{for}\mathrm{\Theta }\pi .$$
(20)
In order to compute the ground state energy of a Bose-Hubbard system with imposed phase gradient, the use of twisted boundary conditions is impracticable. By means of the local unitary transformation $`\widehat{\mathrm{U}}_\mathrm{\Theta }=_{k=1}^N\mathrm{exp}(\mathrm{i}\mathrm{\Theta }\widehat{\mathrm{x}}_k/L)`$ we can map the phase twist from the state onto the Hamiltonian. The eigenvalues of the resulting twisted Hamiltonian for periodic boundary conditions are identical to those of the original Hamiltonian with twisted boundary conditions. The twisted Bose-Hubbard Hamiltonian reads

$$\begin{array}{cc}\hfill \widehat{\mathrm{H}}_\mathrm{\Theta }=& J\underset{i=1}{\overset{I}{}}(\mathrm{e}^{\mathrm{i}\mathrm{\Theta }/I}\widehat{\mathrm{a}}_{i+1}^{}\widehat{\mathrm{a}}_i+\text{h.a.})\hfill \\ & +\underset{i=1}{\overset{I}{}}ϵ_i\widehat{\mathrm{n}}_i+\frac{V}{2}\underset{i=1}{\overset{I}{}}\widehat{\mathrm{n}}_i(\widehat{\mathrm{n}}_i1).\hfill \end{array}$$
(21)
Compared to the original Hamiltonian (3) only the off-diagonal hopping term is modified by an additional factor $`\mathrm{e}^{\mathrm{i}\mathrm{\Theta }/I}`$, the so-called Peierls phase factor ShSu90 ; Poil91 . Physically this term ensures that any particle which tunnels to an adjacent site acquires the correct phase in order to establish the linear phase variation across the system.

There are several ways to impose this phase factor experimentally. One possibility is to add a linear external potential or a homogeneous electric field in the case of charged particles. Another possibility is a constant acceleration of the lattice, which in the case of optical standing wave lattices can be realized by ramping the detuning between the two counter-propagating beams. One might be able to utilize this possibility to probe superfluidity directly.

A straightforward procedure to compute the superfluid fraction is now evident: We numerically solve the eigenvalue problems of the twisted and the non-twisted Hamiltonian, $`\widehat{\mathrm{H}}_\mathrm{\Theta }`$ and $`\widehat{\mathrm{H}}_0`$, with periodic boundary conditions and insert the resulting ground state energies $`E_\mathrm{\Theta }`$ and $`E_0`$ into (20).

Before doing so, it is useful to study how the energy difference $`E_\mathrm{\Theta }E_0`$ entering into Eq. (20) depends on the twist angle $`\mathrm{\Theta }`$. We have already assumed that $`\mathrm{\Theta }`$ has to be sufficiently small to avoid excitations other than the collective flow of the superfluid component. Figure 5 depicts the result of a numerical calculation for $`E_\mathrm{\Theta }E_0`$ as function of $`\mathrm{\Theta }`$ for different values of $`V/J`$. For the noninteracting system (black solid curve) the energy is proportional to $`\mathrm{\Theta }^2`$ up to large values of $`\mathrm{\Theta }`$, i.e., the superfluid fraction is independent of the actual value of $`\mathrm{\Theta }`$ used in the calculation. However, one should note that at $`\mathrm{\Theta }=\pi `$, which corresponds to antiperiodic boundary conditions, a level crossing between the ground state and the first excited state of the twisted system appears. The presence of interactions leads to a “level repulsion” between these two states (see dashed and dotted curves in Fig. 5). Hence the energy deviates significantly from the $`E_\mathrm{\Theta }E_0\mathrm{\Theta }^2`$ behavior for $`\mathrm{\Theta }\pi `$. If one were to use Eq. (20) for $`\mathrm{\Theta }\pi `$ one would underestimate the superfluid fraction significantly. A physically meaningful application of Eq. (20) to determine the superfluid fraction therefore requires $`\mathrm{\Theta }\pi `$. Formally one might even impose the limit $`\mathrm{\Theta }0`$. For all the following numerical calculations we will use $`\mathrm{\Theta }=0.1`$.

#### III.2.2 Direct numerical calculation

We now employ Eq. (20) to compute the superfluid fraction $`f_\text{s}`$ from the energies $`E_\mathrm{\Theta }`$ and $`E_0`$ obtained numerically. The results for systems with commensurate filling $`N/I=1`$ and lattice sizes up to $`I=12`$ are shown in Fig. 6(a).

In the range $`3<V/J<6`$ the superfluid fraction, i.e. the order parameter for the superfluid to Mott-insulator transition, declines rapidly with increasing interaction strength $`V/J`$. In these relatively small systems we already observe a rapid transition from the superfluid phase at small $`V/J`$ to an insulating phase with vanishing superfluid fraction at large $`V/J`$. The transition region is in good agreement with the critical point $`(V/J)_{\text{cr}}=4.65`$ for the infinite lattice determined by Monte Carlo simulations BaSc92  or a strong coupling expansion FrMo96 .

The dependence of the superfluid fraction on the lattice size shown in Fig. 6(a) is rather moderate. In the superfluid phase for up to $`V/J4`$ there is no noticeable size dependence. Only the shoulder of the curve exhibits a size dependence: With increasing $`I`$ the decrease in $`f_\text{s}`$ steepens and the tail of the curve is compressed and shifted towards lower $`V/J`$. From this observation we conclude that systems with $`I10`$ already allow very good qualitative, and even semi-quantitative conclusions, about the superfluid behavior of larger systems.

The rapid decrease of the superfluid fraction is in stark contrast to the behavior of the number fluctuations shown in Fig. 1(a). Evidently, the number fluctuations still exhibit values of up to $`\sigma _i0.3`$ in cases where the superfluid fraction has already vanished. This proves that the Mott-insulator state can still be a superposition of many Fock states and thus show considerable number fluctuations. Only in the limit $`V/J\mathrm{}`$ when $`\sigma _i`$ vanishes it becomes a pure Fock state.

#### III.2.3 Perturbative treatment

Apart from the direct solution of the eigenvalue problems of the twisted and the non-twisted Hamiltonian one can use a perturbative approach to determine the energy difference $`E_\mathrm{\Theta }E_0`$. As the superfluid fraction (20) is defined only for very small twist angles $`\mathrm{\Theta }`$ this perturbative treatment does not introduce any additional approximation. All terms of the expansion that contribute in the limit $`\mathrm{\Theta }0`$ will be included.

The derivation consists of two steps: We first expand the twisted Hamiltonian in a power series up to second order in $`\mathrm{\Theta }`$

$$\widehat{\mathrm{H}}_\mathrm{\Theta }\widehat{\mathrm{H}}_0+\frac{\mathrm{\Theta }}{I}\widehat{\mathrm{J}}\frac{\mathrm{\Theta }^2}{2I^2}\widehat{\mathrm{T}}=\widehat{\mathrm{H}}_0+\widehat{\mathrm{H}}_{\text{pert}},$$
(22)
where the current operator $`\widehat{\mathrm{J}}`$ and the kinetic energy operator $`\widehat{\mathrm{T}}`$ are given by:

$$\begin{array}{cc}\hfill \widehat{\mathrm{J}}& =\mathrm{i}J\underset{i=1}{\overset{I}{}}(\widehat{\mathrm{a}}_{i+1}^{}\widehat{\mathrm{a}}_i\text{h.a.})\hfill \\ \hfill \widehat{\mathrm{T}}& =J\underset{i=1}{\overset{I}{}}(\widehat{\mathrm{a}}_{i+1}^{}\widehat{\mathrm{a}}_i+\text{h.a.}).\hfill \end{array}$$
(23)
We then compute the energy shift $`E_\mathrm{\Theta }E_0`$ caused by the perturbation $`\widehat{\mathrm{H}}_{\text{pert}}`$ in (22) via second order perturbation theory. The lowest order contributions to the energy shift are quadratic in $`\mathrm{\Theta }`$. Inserting them into (20) leads to an expression for the superfluid fraction consisting of two terms

$$f_\text{s}=f_\text{s}^{(1)}f_\text{s}^{(2)}.$$
(24)
The first order term is proportional to the kinetic energy expectation value

$$f_\text{s}^{(1)}=\frac{1}{2NJ}\mathrm{\Psi }_0|\widehat{\mathrm{T}}|\mathrm{\Psi }_0=\frac{1}{N}\underset{i=1}{\overset{I}{}}\rho _{i,i+1}^{(1)},$$
(25)
where $`\mathrm{\Psi }_0`$ is the ground state of the original Hamiltonian $`\widehat{\mathrm{H}}_0`$ and $`\rho _{ij}^{(1)}`$ the associated one-body density matrix introduced in Eq. (10). The second order term involves matrix elements of the current operator between the ground state $`|\mathrm{\Psi }_0`$ and all the excited states $`|\mathrm{\Psi }_\nu `$ ($`\nu =1,2,\mathrm{}`$) of the original Hamiltonian $`\widehat{\mathrm{H}}_0`$,

$$f_\text{s}^{(2)}=\frac{1}{NJ}\underset{\nu 0}{}\frac{|\mathrm{\Psi }_\nu |\widehat{\mathrm{J}}|\mathrm{\Psi }_0|^2}{E^{(\nu )}E^{(0)}}.$$
(26)
All higher-order terms of the perturbative expansion of the energy shift do not contribute to the superfluid fraction in the limit $`\mathrm{\Theta }0`$. Thus Eq. (24) is an exact expression for $`f_\text{s}`$.

We note that this derivation and the resulting expression is closely related to the so-called Drude weight which characterizes the dc conductivity of charged fermionic systems FyMa91 .

This formulation of the superfluid fraction provides a detailed insight into the mechanisms that govern the superfluid properties of the system. First of all it shows that superfluidity is not a property of the ground state, but rather the response of the system to an external perturbation. The second order term (26) introduces an explicit dependence on the whole excitation spectrum. Since this second order contribution always lowers the total superfluid fraction we can conclude that the coupling to excited states leads to a suppression of superfluidity. The first order term (25), which is just the ground state expectation value of the kinetic energy and can be expressed in terms of the one-body density matrix of the ground state alone, gives only an upper bound for $`f_\text{s}`$ PaTr98 .

The behavior of the first and second order contributions to the superfluid fraction for a regular lattice are shown in Fig. 6(b). The first order term (25), i.e. the rescaled kinetic energy expectation value, is approximately one for small $`V/J`$ and decreases slowly with increasing interaction strength. Within the insulating phase, where the total superfluid fraction shown in Fig. 6(a) has vanished completely, the first order contribution $`f_\text{s}^{(1)}`$ still has values of up to $`0.3`$.

The vanishing of the total superfluid fraction is caused by a characteristic behavior of the second order contribution (26): In the transition region $`f_\text{s}^{(2)}`$ shows a threshold-like increase from zero to its maximum value. For values of $`V/J`$ beyond this maximum the second order term quickly converges to the value of $`f_\text{s}^{(1)}`$. Thus the vanishing of the superfluid fraction in insulating phases, e.g., the Mott-insulator, results from the cancellation between first and second order term which both have considerable size.

This clearly demonstrates that the coupling to the excited states and the structure of the excitation spectrum is essential for the superfluid properties of the system. The ground state alone, represented by the first order contribution $`f_\text{s}^{(1)}`$, provides only limited information—namely an upper bound—on the superfluid fraction. Therefore, any observable which is sensitive to the ground state only, such as number fluctuations, condensate fraction, and interference pattern, cannot provide clear information the superfluid properties of the system.

The only other observable which is sensitive to the properties of the excitation spectrum is the energy gap $`E_{\text{gap}}`$ between the ground state and the first excited state. As we have shown in Sec. II.3 the energy gap starts to grow linearly with $`V/J`$ in the Mott-insulator phase. This is another indicator for the change in the excitation spectrum that also leads to the vanishing of the superfluid fraction. However, as Eq. (26) illustrates, the dependence of the superfluid fraction on the excited states is more complex than just the energy gap.

## IV Experimental observables

The two quantities, condensate fraction and superfluid fraction, discussed in the preceeding section describe two important and fundamental properties of the system, however, they are not directly accessible to simple experiments. In this section we discuss two observables which are directly accessible to present day experiments and show their relation to the fundamental quantities discussed above.

### IV.1 Matter-wave interference pattern

The simplest experimental procedure to obtain information on the state of the Bose gas in the lattice is to release the atoms from the lattice and detect the matter-wave interference pattern after some time of flight OrTu01 ; GrMa02a ; GrMa02b . The crucial question is, how much the interference pattern can tell us about the superfluid or the condensed portion of the system.

Neglecting interactions during the expansion of the gas we can describe the time evolution of the single-particle states, which are initially localized at the different lattice sites, by Gaussian wave-packets. The matter-wave intensity observed at a point $`\stackrel{}{y}`$ after ballistic expansion for a time $`\tau `$ can be written as RoBu03a

$$(\stackrel{}{y})=\mathrm{\Psi }_0|\widehat{\mathrm{A}}^{}(\stackrel{}{y})\widehat{\mathrm{A}}(\stackrel{}{y})|\mathrm{\Psi }_0.$$
(27)
The amplitude operators are given by

$$\widehat{\mathrm{A}}(\stackrel{}{y})=\underset{i=1}{\overset{I}{}}G_i(\stackrel{}{y},\tau )\widehat{\mathrm{a}}_i,$$
(28)
where $`G_i(\stackrel{}{y},\tau )`$ is the amplitude of a Gaussian wave-packet that originates from site $`i`$ after free expansion for a time $`\tau `$. We are interested only in the generic structure of the intensity pattern, i.e., the presence or the absence of interference peaks. Therefore, we can simplify the treatment significantly by dropping all the terms related to the spatial envelope of the interference pattern. The Gaussian can thus be replaced by a phase factor and the amplitude operator reduces to

$$\widehat{\mathrm{A}}(\stackrel{}{y})=\underset{i=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}\varphi _i(\stackrel{}{y},\tau )}\widehat{\mathrm{a}}_i=\underset{i=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}\delta \varphi (\stackrel{}{y},\tau )i}\widehat{\mathrm{a}}_i,$$
(29)
Here $`\varphi _i(\stackrel{}{y},\tau )`$ is the phase accumulated on the path from the $`i`$th site to the observation point $`\stackrel{}{y}`$. In the far-field approximation this can be replaced by $`\delta \varphi (\stackrel{}{y},\tau )i`$, where $`\delta \varphi (\stackrel{}{y},\tau )`$ is the phase difference between paths originating from adjacent lattice sites. The matter-wave intensity (27) as function of $`\delta \varphi `$ now reads

$$\begin{array}{cc}\hfill (\delta \varphi )& =\frac{1}{I}\underset{i,j=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}\delta \varphi (ij)}\mathrm{\Psi }_0|\widehat{\mathrm{a}}_j^{}\widehat{\mathrm{a}}_i|\mathrm{\Psi }_0\hfill \\ & =\frac{1}{I}\underset{i,j=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}\delta \varphi (ij)}\rho _{ij}^{(1)}.\hfill \end{array}$$
(30)
This result proves the direct relation between the interference pattern and the quasimomentum distribution of the gas in the lattice given by Eq. (15) RoBu03a . The occupation number $`\stackrel{~}{n}_q`$ of the quasimomentum $`q`$ state is given by the intensity $`(\delta \varphi )`$ for a phase difference $`\delta \varphi =qa`$:

$$\stackrel{~}{n}_q=(\delta \varphi =qa).$$
(31)
As a special case of this result, the condensate fraction, i.e. the occupation of the $`q=0`$ state, is proportional to the intensity for $`\delta \varphi =0`$. Of course, this result agrees with the naive view that the interference pattern simply gives the momentum distribution of the trapped gas KaPr02 .

Furthermore Eq. (30) enables us to draw an important conclusion about the relation between the matter-wave interference pattern and the superfluid properties of the system. The interference pattern depends exclusively on the one-body density matrix $`\rho _{ij}^{(1)}`$ for the ground state. The superfluid fraction, in contrast, depends crucially on the excited states of the system as we have discussed in detail in Sec. III.2.3. Thus the interference pattern cannot give complete information on the superfluid properties of the system. It does not probe the physics that is crucial for superfluidity.

The matter-wave interference patterns which result from Eq. (30) for three different interaction strength $`V/J`$ are depicted in Fig. 7. For the noninteracting system (solid curve) it exhibits sharp interference peaks at $`\delta \varphi =0`$ and $`2\pi `$ which correspond to the principal peaks observed experimentally GrMa02a . Because we discarded the terms related to the envelope, the intensity $`(\delta \varphi )`$ is simply a $`2\pi `$-periodic function in $`\delta \varphi `$. With increasing interaction strength $`V/J`$ (dashed and dotted curves in Fig. 7) the amplitude of the principal peaks is reduced and an “incoherent background” emerges. Effectively this leads to a gradual broadening of the interference peaks which was observed experimentally GrMa02a .

The emergence of the background has a straightforward interpretation in terms of the quasimomentum distribution. The gray arrows in Fig. 7 mark the values of $`\delta \varphi `$ which correspond to the discrete quasimomenta $`q`$ allowed in the lattice. According to (31) the intensity at these points just gives the occupation numbers for the different quasimomentum states. Thus the depletion of the principal peak with increasing $`V/J`$ just indicates the depletion of the condensate. The emergence of a background shows that states with nonvanishing quasimomenta are successively populated—as we have already discussed in Sec. III.1.2. In the limit of strong interactions all quasimomentum states of the band are occupied uniformly and the corresponding interference pattern is perfectly flat.

The fact that the vanishing of the superfluid fraction is not associated with the vanishing of the interference fringes is highlighted by these numerical results. The values of $`V/J`$ used to compute the three intensity distributions in Fig. 7 correspond to superfluid fractions of $`f_\text{s}=1`$ (solid curve), $`f_\text{s}0.5`$ (dashed), and $`f_\text{s}0`$ (dotted). In the latter case, interference fringes are still clearly visible although the system is a perfect insulator, i.e., the superfluid fraction is zero. Thus, we see that with increasing $`V/J`$ the superfluid component vanishes much earlier than the interference fringes. In other words, the Mott insulator phase still exhibits a degree of phase coherence.

As a simple quantity to characterize the interference pattern one might introduce the visibility of the fringes. The standard definition of the fringe visibility

$$𝒱=\frac{_{\mathrm{max}}_{\mathrm{min}}}{_{\mathrm{max}}+_{\mathrm{min}}}$$
(32)
only relies on the maximum and minimum values of the intensity. For a regular lattice with even number of sites we immediately find $`_{\mathrm{max}}=(0)`$ and $`_{\mathrm{min}}=(\pi )`$. On the basis of (31) we can identify the maximum (minimum) intensity with the largest (smallest) quasimomentum occupation number. This definition of the visibility therefore provides a direct measure for the inhomogeneity of the quasimomentum distribution: $`𝒱`$ is $`1`$ as long as there are unoccupied quasimomentum states and it goes to zero if the band is filled uniformly.

The dependence of the visibility measure on the interaction strength is depicted in Fig. 8. The visibility decreases only very slowly with increasing $`V/J`$ and can have values up to $`𝒱0.7`$ in the insulating phase where $`f_\text{s}0`$. This confirms the observation that superfluidity vanishes much faster than the interference fringes and thus the coherence of the system RoBu03a . The fact that the curves in Fig. 8 are practically independent of the lattice size proves that this statement also holds for large lattices and is not an artifact of the relatively limited size. With regard to the quasimomentum distribution the slow decrease of $`𝒱`$ indicates that very large ratios $`V/J`$ are required to achieve a perfectly uniform occupation of the band.

One can define other measures for the fringe visibility which are more sensitive to particular changes in the interference patterns. However, the most important and decisive measure for the coherence properties of the system is the interference pattern itself.

### IV.2 Structure factor

In the presence of a periodic modulation of the lattice potential the atom density distribution will develop a spatial structure RoBu03b . Systems of this type will be discussed in detail in Sec. V. These periodic density modulations should be observable using Bragg diffraction of light from the trapped atomic gas. Experiments of this kind have already been performed using laser-cooled atomic gases in optical lattices BiGa95 ; WeHe95 .

The quantity directly accessible through Bragg diffraction experiments is the structure factor $`S(k)`$ which is given by OtWa95

$$S(k)=\frac{1}{I^2}\underset{i,j=1}{\overset{I}{}}\mathrm{e}^{\mathrm{i}k(\xi _i\xi _j)}\mathrm{\Psi }_0|\widehat{\mathrm{n}}_i\widehat{\mathrm{n}}_j|\mathrm{\Psi }_0.$$
(33)
This quantity is sensitive to density-density correlations and is therefore a suitable indicator for the presence of diagonal long-range order. For a homogeneous occupation of all lattice sites, present in the regular lattices discussed so far, the structure factor exhibits distinct peaks at integer values of $`ka/(2\pi )`$. The presence of any additional periodic structure in the atomic density distribution will lead to the appearance of new peaks at noninteger values of $`ka/(2\pi )`$. These signal the presence of diagonal long-range order. From the position and the relative amplitude of the additional peaks in the structure factor one can derive detailed information on the spatial density distribution. We will discuss this in detail in Sec. V.3.

## V Two-color superlattices

Having discussed the relevant observables on a formal level and compared their behavior for the superfluid to Mott-insulator transition in a regular lattice we now want to explore another aspect of the phase diagram of the zero temperature Bose gas. Optical lattices offer the unique possibility of modifying the topology of the lattice in a flexible and perfectly controlled manner. They thus open new possibilities for studying many-body physics in non-uniform lattice potentials. In contrast to most solid-state systems, where the precise topology of the irregular lattice is difficult to specify, the high degree of experimental control facilitates investigations into specific aspects of “disorder”-induced phenomena.

### V.1 Optical two-color lattices

The simplest way to experimentally generate a non-uniform lattice potential is through a superposition of two collinear optical standing waves with different wavelengths—a so called two-color lattice. The interference between these standing waves produces a superlattice with a sinusoidal modulation of the depth of the lattice wells. More complex topologies have already been realized experimentally by using a superposition of several laser beams with different wave vectors GuTr97 . Truly random lattices can be generated by superimposing a speckle pattern DaZa03 . Thus optical lattices cover the whole range from spatially modulated lattice potentials to disordered systems.

In the following we consider a superlattice generated by the interference of two standing waves with wavelengths $`\lambda _1`$ and $`\lambda _2=\frac{5}{7}\lambda _1`$. The intensity of the second laser field is assumed to be much smaller than the intensity of the first laser, in this example $`I_2=0.04I_1`$. This combination results in a superlattice consisting of identical unit cells each composed of $`5`$ elementary lattice sites. Figure 9(a) depicts the lattice potential $`U_{\text{lat}}(x)`$ for a certain value of the relative phase.

In the language of the Bose-Hubbard model this spatial modulation of the well depths corresponds to a modulation of the on-site energies $`ϵ_i`$. We assume a sinusoidal variation of the $`ϵ_i`$ in the range $`[\mathrm{\Delta },0]`$ as depicted in Fig. 9(b). The parameter $`\mathrm{\Delta }`$ characterizes the amplitude of the modulation, for $`\mathrm{\Delta }=0`$ we recover the regular lattice potential. The modulation of the lattice potential also influences the tunneling coefficient $`J`$ which, in general, varies for the different lattice sites. For the lattices considered here the influence of this variation can be neglected in good approximation DaZa03 .

### V.2 Localized phase

As a first step we consider a noninteracting Bose gas in the superlattice depicted in Fig. 9. We solve the eigenvalue problem for the Bose-Hubbard Hamiltonian including two unit cells with filling factor $`1`$, i.e., $`I=N=10`$. For $`V/J=0`$ we successively increase the amplitude $`\mathrm{\Delta }`$ of the spatial modulation of the on-site energies. The resulting distributions of mean occupation numbers $`\overline{n}_i`$ and number fluctuations $`\sigma _i`$ for the individual sites are depicted in Fig. 10. For a regular system, without any modulation of the $`ϵ_i`$, the lattice is occupied homogeneously. With increasing modulation amplitude $`\mathrm{\Delta }/J`$ the particles preferentially occupy those lattice sites with lower $`ϵ_i`$ as one would classically expect. The degree of localization is governed by the competition between the on-site energies and the tunneling term of the Bose-Hubbard Hamiltonian (3). For sufficiently large $`\mathrm{\Delta }/J`$ only the site with the lowest $`ϵ_i`$ in each unit cell is populated. The number fluctuations $`\sigma _i`$ behave accordingly: with increasing $`\mathrm{\Delta }/J`$ the number fluctuations at the deepest lattice well grow whereas the fluctuations at all other wells slowly decrease. In the limit of large $`\mathrm{\Delta }/J`$ the ground state is a superposition of all Fock states which have nonzero occupation numbers only at the deepest well of each unit cell.

There are several possibilities for detecting this localized phase in an experiment. The simplest signature appears in the interference pattern after the expansion of the gas. Figure 11 shows the intensity distribution that results from Eq. (30) for three different values of $`\mathrm{\Delta }/J`$. In the absence of the spatial modulation (upper panel) we observe the interference pattern of a regular lattice discussed in Sec. IV.1. With increasing $`\mathrm{\Delta }/J`$, that is increasing degree of localization, additional peaks emerge. First the peaks next to initial principal peaks appear, then the next-order peaks emerge. In the limit of perfect localization we obtain a regular interference pattern with a peak separation corresponding to the inverse of the number of elementary sites in a unit cell. In this manner the relative strength of the emerging peaks provides a direct experimental measure for the degree of localization in these superlattices.

Another experimental observable that is sensitive to the degree of localization is the structure factor $`S(ka)`$ introduced in Sec. IV.2. We will come back to this quantity in the next section.

### V.3 Quasi Bose-glass phase

It is intuitively clear that repulsive interactions will have a strong influence on the localized phase: A concentration of many particles at a single lattice site will generate a large interaction energy contribution. Thus the interaction enforces delocalization of the particles.

To illustrate the competition between lattice modulation and repulsive two-body interaction Fig. 12 depicts the change of the mean occupation numbers $`\overline{n}_i`$ and the number fluctuations $`\sigma _i`$ with increasing interaction strength $`V/J`$ for a fixed modulation amplitude $`\mathrm{\Delta }/J=50`$.

The perfectly localized configuration persists only at very small interaction strengths. Already for $`V/J2`$ it is energetically favorable to reduce the maximum occupation number and thus the interaction energy by redistributing particles to sites with higher on-site energies $`ϵ_i`$. This region of interaction-induced delocalization is sometimes termed Anderson-glass phase ScBa91 ; DaZa03 . Initially the changes in the mean occupation numbers happen continuously. However, for $`V/J10`$ this behavior changes. As shown in Fig. 12(a) there are extended intervals of interaction strengths where the mean occupation numbers are constant and approximately integer. Between these regions of “stable” configurations there are narrow windows in which some of the $`\overline{n}_i`$ change by roughly $`\pm 1`$. In the present example one of these almost steplike rearrangements happens at $`V/J30`$. This behavior of the mean occupation numbers is also reflected in the number fluctuations depicted in Fig. 12(b). Within the stable regions all $`\sigma _i`$ are rather small, i.e., the dominant contribution to the expansion (5) of the ground state comes from a single Fock state. The steplike rearrangements are accompanied by an increase of the number fluctuations of those sites that change their mean occupation number.

It is the subtle interplay between interaction and lattice irregularity that drives these successive rearrangements between configurations with integer occupation numbers. The sequence of rearrangements depends crucially on the topology of the lattice, i.e., the set of on-site energies. If the distribution of the $`ϵ_i`$ becomes more homogeneous (meaning that the number of different values of the $`ϵ_i`$ increases) then the number of possible rearrangements increases and the stable regions shrink. In the limit of a lattice with random disorder this leads to the so-called Bose-glass phase ScBa91 ; KrTr91 . For the superlattices with a rather small number of different on-site energies we will use the term quasi Bose-glass phase. The study of the superfluid fraction in Sec. V.4 shows that the quasi Bose-glass phase is insulating and that the energy gap is small.

If the interaction strength reaches a value $`V=\mathrm{\Delta }`$ a final and very sharp rearrangement happens and the system assumes a homogeneous occupation of all sites. This is the transition to the Mott-insulator phase where despite of the presence of the irregular lattice the ground state is dominated by the Fock state with uniform occupation numbers $`n_i=N/I`$. This transition between quasi Bose-glass and Mott-insulator phase happens independently of the lattice topology if the interaction strength $`V`$ exceeds the amplitude $`\mathrm{\Delta }`$ of the lattice modulation.

The sequence of rearrangements in the quasi Bose-glass phase can be directly observed in experiment. The observable which is most sensitive to the changes in the spatial distribution of the atoms in the superlattice is the structure factor $`S(ka)`$ introduced in Sec. IV.2. For a state with a homogeneous population of the lattice sites, e.g., in the Mott-insulator phase, $`S(ka)`$ shows peaks at integer values of $`ka/(2\pi )`$ only. If mean-occupation numbers, i.e. the diagonal elements of the one-body density matrix, exhibit some periodic long-range structure then new peaks will appear at noninteger values of $`ka/(2\pi )`$ from which information on the long-range order can be extracted.

The localized phase in a superlattice at small interaction strengths is the most pronounced realization of diagonal long-range order: For the superlattice discussed here each $`5`$th lattice site is occupied and all others are practically empty. The corresponding structure factor $`S(ka)`$ is shown in the upper panel of Fig. 13. Peaks of almost equal height emerge at $`ka/(2\pi )=\nu /5`$ ($`\nu `$ integer) indicating the perfect localization at a single site in each unit cell. With increasing interaction strength $`V/J`$ the particles are gradually redistributed to other lattice sites. This leads to a suppression of the peaks at non-integer values of $`ka/(2\pi )`$ as the two lower panels Fig. 13 illustrate. Thus the amplitude of the peaks in $`S(ka)`$ provides detailed information on the complicated spatial structure within the quasi Bose-glass phase.

To highlight this point further, Fig. 14 shows the dependence of the amplitude of the first peak of $`S(ka)`$ located around $`ka/(2\pi )1/5`$ as function of $`V/J`$ for fixed amplitude $`\mathrm{\Delta }/J`$ of the lattice modulation. The parameters correspond to those used in Fig. 12 and cover the full range from the localized phase to the Mott-insulator. Obviously each of the rearrangements of the mean occupation numbers within the quasi Bose-glass phase leaves a distinct signature in the peak amplitude. The regions of stable configurations as well as the steplike rearrangements that are visible in Fig. 12(a) directly show up in $`S(ka)`$. Thus the structure factor, which is experimentally accessible through Bragg scattering of light from the bosons in the lattice, provides a comprehensive insight into the spatial density structure and its complex dependency on the interaction strength and the modulation amplitude.

### V.4 Phase diagrams in the $`V`$-$`\mathrm{\Delta }`$ plane

We shall now illustrate the behavior of the superfluid fraction, condensate fraction, and energy gap using two-dimensional phase diagrams in the $`V`$-$`\mathrm{\Delta }`$ plane. This representation exposes the complicated interplay between lattice modulation and two-body interaction which leads to the rich phase structure.

The gross phase structure is already revealed in the behavior of the largest coefficient $`C_{\mathrm{max}}^2`$ in the expansion (5) of the ground state which is depicted in Fig. 15(a). The dark shadings correspond to values of $`C_{\mathrm{max}}^2`$ close to $`1`$, i.e. ground states which are almost pure Fock states. The prominent region of large $`C_{\mathrm{max}}^2`$ for values $`V>\mathrm{\Delta }`$ corresponds to the Mott-insulator phase. The lobe-like structures at large modulation amplitudes $`\mathrm{\Delta }`$ belong to the quasi Bose-glass phase. They correspond to the regions of stable configurations discussed in connection with Fig. 12 in which the mean occupation numbers are constant and approximately integral. The valleys of small $`C_{\mathrm{max}}^2`$ which separate the lobes are connected to the rearrangements between different stable configurations. As mentioned earlier the number and position of these lobes depends on the topology of the lattice potential.

In close connection with the behavior of $`C_{\mathrm{max}}^2`$ are the maximum number fluctuations $`\sigma _{\mathrm{max}}=\mathrm{max}(\sigma _i)`$ depicted in Fig. 15(b). Values of $`C_{\mathrm{max}}^2`$ close to $`1`$ imply that the number fluctuations at all lattice sites are small. Thus the Mott-insulator phase and the lobes of the quasi Bose-glass phase appear as regions of small number fluctuations (light shadings). In contrast, the rearrangement valleys in the quasi Bose-glass phase, which are associated with small $`C_{\mathrm{max}}^2`$, show increased number fluctuations. Number fluctuations of the order $`1`$ emerge only within a narrow band at small values of $`V/J`$. For large values of $`\mathrm{\Delta }/J`$ this corresponds to the localized phase discussed in Sec. V.2.

The behavior of the condensate fraction $`f_\text{c}`$ is shown in Fig. 15(c). In line with variations of the maximum number fluctuations, $`f_\text{c}`$ assumes very small values in the Mott-insulator phase. Within the quasi Bose-glass it decreases monotonically with increasing $`V/J`$, i.e., it does not show the lobe structure as the number fluctuations. Values of the condensate fraction above $`0.6`$ result only in a small area at small $`V/J`$ and $`\mathrm{\Delta }/J`$ and a very narrow stripe at $`V/J2`$. Thus only the localized phase and the expected superfluid phase exhibit sizeable condensate fractions.

So far, these quantities are determined by the properties of the ground state of Bose gas. Following the discussion in Sec. III.2 they therefore do not fully determine the superfluid properties of the system. Figure 16 shows the contour plots for two quantities which probe the excitation spectrum, these are the superfluid fraction and the energy gap. The plot for superfluid fraction $`f_\text{s}`$ in Fig. 16(a) pinpoints the boundaries of the superfluid phase. Large values of the superfluid fraction (dark shadings) result only in a restricted parameter region at small values of $`V`$ and $`\mathrm{\Delta }`$. Everywhere else—that is, in the Mott-insulator, the quasi Bose-glass, and the localized phase—the system is a perfect insulator. In comparing the boundaries of the superfluid phase defined in Fig. 16(a) with the structure of the other observables, like number fluctuations and condensate fraction in Fig. 15, it becomes evident that none of the ground state quantities can identify the superfluid region. Not only that the magnitude of the superfluid fraction cannot be determined, the whole structure of the plots is different.

This also applies to the energy gap $`E_{\text{gap}}`$ plotted in Fig. 16(b). Although the opening of a gap between the ground state and the first excited state is an indicator for the superfluid to Mott-insulator transition in the uniform lattice, there is no such correspondence for the irregular lattice. It is only the Mott-insulator phase that exhibits a large energy gap. The localized and the quasi Bose-glass phase, which are also insulating, have small energy gaps just like the superfluid phase.

The energy gap of the quasi Bose-glass phase deserves closer inspection. Figure 16(b) reveals that there is a small but finite energy gap within the stable lobes of the quasi Bose-glass phase. The size of the gaps, as well as the position and number of the lobes, depend on the differences between the various on-site energies. If the size of the unit cells is increased or if one goes to large random lattices then the differences between the $`ϵ_i`$ become smaller, the number of lobes increases, their size decreases, and the energy gap vanishes. Thus there is a continuous transition from the quasi Bose-glass to the strict Bose-glass phase with vanishing $`E_{\text{gap}}`$.

Another interesting observation results from the comparison of the superfluid fraction in Fig. 16(a) with the condensate fraction in Fig. 15(c). In an irregular lattices with weak interactions there is a region where the condensate fraction is larger than the superfluid fraction. An extreme case is the localized phase at $`V/J2`$, where $`f_\text{s}0`$ and $`f_\text{c}1`$. This means that the localized bosons in the different unit cells are still phase coherent but the superfluid flow is inhibited by the lattice topology. A similar phenomenon was recently observed in a Monte Carlo study of a continuous Bose gas with random impurities AsBo02  which act like the irregular lattice potential.

So far we have only considered the phase diagrams for a filling factor $`N/I=1`$. The qualitative structure of the phase diagram does not change if we go to larger integer filling factors. Quantitatively the superfluid region expands towards larger values of the interaction strength $`V/J`$ and the modulation amplitude $`\mathrm{\Delta }/J`$ if the filling factor is increased. For a filling factor $`N/I=4`$, for example, the contour corresponding to a superfluid fraction $`f_\text{s}=0.6`$ extends up to $`V/J15`$ and $`\mathrm{\Delta }/J40`$. The insulating phases are shifted accordingly.

### V.5 Non-commensurate filling factors

It is well-known that the Mott-insulator phase in a uniform lattice exists only for commensurate filling factors. If one adds a few particles on top of a commensurate configuration then these particles can be moved to any lattice site without changing the energy—even in the presence of strong repulsive interactions. In this way the additional particles gives rise to a finite superfluid fraction in the limit of large $`V/J`$.

This intuitive picture is confirmed by the exact calculations for the sinusoidal superlattice with $`I=10`$ sites and $`N=12`$ particles. Contour plots for the superfluid fraction and the energy gap in the $`V`$-$`\mathrm{\Delta }`$ plane are presented in Fig. 17. For the uniform lattice ($`\mathrm{\Delta }/J=0`$) the superfluid fraction shown in panel (a) remains finite even at large $`V/J`$.

If one starts to increase the amplitude $`\mathrm{\Delta }`$ of the sinusoidal modulation then the superfluid fraction decreases and eventually vanishes above $`\mathrm{\Delta }/J5`$. In this regime the growing lattice irregularity causes a pinning of the additional particles to the lattice sites with minimal on-site energies $`ϵ_i`$. The ground state has an increased occupation number at the deepest lattice well of each unit cell.

The behavior of the energy gap shown in Fig. 17(b) confirms this picture. For strong repulsive interactions the energy gap vanishes at small amplitudes $`\mathrm{\Delta }/J`$ in stark contrast to the commensurate system. With increasing amplitude of the sinusoidal modulation the energy gap starts to grow linearly with $`\mathrm{\Delta }/J`$. In this simple case $`E_{\text{gap}}`$ corresponds to the energy required to remove the additional particle from the lattice site it is affixed to. The slope of the increase is proportional to the energy difference between the two lowest $`ϵ_i`$. This behavior is equivalent to the structure of the lobes of the quasi Bose-glass phase.

Thus, in general, the change from a commensurate to a non-commensurate filling factor causes a vanishing of the Mott-insulator phase. It is replaced by a superfluid phase at small $`\mathrm{\Delta }/J`$, where the additional particles constitute a superfluid on top of an Mott-insulating layer. At larger $`\mathrm{\Delta }/J`$ the additional particles are pinned to the deepest lattice wells, superfluidity ceases, and an additional lobe of the quasi Bose-glass phase emerges. All other regions of the phase diagram are qualitatively unaffected by the change from a commensurate to a non-commensurate filling factor.

## VI Conclusions

We have shown that the superfluid fraction—the natural order parameter for superfluid to insulator phase transitions—is determined by the response of the system to an external perturbation, i.e., a twist in the boundary conditions, or equivalently an additional phase factor, in the tunneling term of the Bose-Hubbard Hamiltonian. This means that the superfluid fraction necessarily depends on the full excitation spectrum of the lattice system. The formal manifestation of this fact is given in Eq. (24), where the superfluid fraction is expressed as a sum of a first order term depending only of the ground state and a second order term which involves all excited states. The exact numerical solution of the Bose-Hubbard model for uniform one-dimensional lattices reveals that the second order term is crucial for the vanishing of the superfluid fraction in the Mott-insulator phase. The first order contribution gives only an upper bound for the superfluid fraction and can have a rather large value in the insulating phase.

This allows important conclusions for the prospects of measuring the superfluid fraction and the details of superfluid to insulator phase transitions. Present experimental observables are sensitive only to the ground state and therefore cannot probe the physics that governs the superfluid properties. The matter-wave interference pattern, for example, provides direct information on the quasimomentum distribution of the ground state but it does not measure the superfluid properties. In the case of the superfluid to Mott-insulator transition our numerical results reveal that the superfluid fraction vanishes much earlier than the interference fringes, i.e., an interference pattern is still visible in the Mott-insulator phase. This demonstrates that one has to clearly distinguish between superfluidity and coherence properties.

In the second part of the paper we have employed these tools to explore the zero-temperature phase diagram of Bose gases in non-uniform lattice potentials. We have demonstrated that even a simple superlattice potential which results from the superposition of two standing-wave lattices with different wavelengths gives rise to a very rich phase diagram. As function of the interaction strengths $`V/J`$ and the amplitude $`\mathrm{\Delta }/J`$ of the spatial variation of the on-site energies two additional insulating phases can be identified: a localized phase and a quasi Bose-glass phase. All insulating phases can be clearly distinguished through their characteristic signatures in matter-wave interference experiments or through the structure factor determined by Bragg diffraction of light.

Our results further support the view that Bose gases in optical lattices are a versatile tool for studying quantum mechanical many-body phenomena in strongly correlated systems. They allow us to address fundamental questions such as the connection between superfluidity and Bose-Einstein condensation. Moreover, they enable controlled studies of disorder-induced phenomena which result from a subtle competition between kinetic energy, two-body interaction, and lattice topology. To this end two-color lattices with relatively small unit cells seem a promising starting point because they facilitate a clear experimental distinction between the various phases.

## Acknowledgments

This work was supported by the DFG, the UK EPSRC, and the EU via the “Cold Quantum Gases” network. Keith Burnett thanks the Royal Society and Wolfson Foundation for support. We acknowledge useful discussions with S.A. Gardiner, J.A. Dunningham, T. Köhler, T. Gasenzer, and K. Braun-Munzinger.