Temperature changes when adiabatically ramping up an optical lattice

When atoms are loaded into an optical lattice, the process of gradually turning on the lattice is almost adiabatic. In this paper we investigate how the temperature changes when going from the gapless superfluid phase to the gapped Mott phase along isentropic lines. To do so we calculate the entropy in the single-band Bose-Hubbard model for various densities, interaction strengths and temperatures in one and two dimensions for homogeneous and trapped systems. Our theory is able to reproduce the experimentally observed visibilities and therefore strongly supports that current experiments remain in the quantum regime for all considered lattice depths with low temperatures and minimal heating.


I. INTRODUCTION
During recent years substantial progress has been made in cold atom experiments. Among the greatest achievements are the experimental realization of the superfluid to Mott insulating transition with bosonic atoms in optical lattices [1] and the BEC to BCS crossover in fermionic quantum gases [2]. Due to the good tunability the excitement to use cold atomic gases as a so-called quantum simulator has been hard to temper. Such quantum simulators are systems which mimic -in a controllable and ultraclean way -simple strongly-interacting models such as the Hubbard model. With ultracold atom experiments the dream of using such quantum simulators to obtain new insight into the long-standing problems of other research areas seems 'almost' feasible. One of the notorious challenges of condensed matter physics that might be solved in such quantum simulations is the fermionic Hubbard model [3,4], believed to be relevant for high temperature superconductivity [5].
Beside the use as quantum simulators, the unprecedented tunability of cold atomic systems makes them very promising candidates for quantum computers. The experimental production of entanglement has been a large step forward into performing quantum computations [6].
The prerequisite of applying the quantum simulator to unknown physics, is a complete understanding of the present experiments. One would expect full agreement between the current generation of bosonic experiments and theory, since the properties of the underlying homogeneous Bose-Hubbard model are rather well studied. However, our understanding of the present experiments * Electronic address: pollet@itp.phys.ethz.ch is far from complete. The interpretation of the results is mainly complicated by two points: (i) The presence of an external trapping potential can cause the spatial coexistence of the superfluid and Mott phases. This replaces the quantum phase transition by a gradual crossover, which can obscure the characteristic signal of the quantum phases.
(ii) A change in temperature due to adiabatic or nonadiabatic origins can also hide the signature of the quantum phase transition, replacing it by a thermal transition. Considerable heating in current experiments would cast serious doubt [7,8,9] on former interpretations. Knowledge of the temperature in an optical lattice is therefore highly desirable, but whereas the temperature of a weakly-interacting bosonic gas in a parabolic trap can be measured accurately, no reliable temperature measurement exists in the presence of a deep optical lattice.
In this manuscript we reexamine the interpretation of recent experiments of bosons in optical lattices focusing on the effect of temperature in the presence of a trapping potential. We determine the lowest temperature that can be reached by loading the atoms adiabatically into the optical lattice. It is reasonable to assume that such processes are isolated and close to thermodynamic equilibrium at all times, since ramping up the lattice is typically slow (∼ 16 ms per increase E R in lattice laser intensity) compared to the tunneling rate of the bosons (∼ 1 kHz, except for the deepest lattices) [10]. This assumption gives a lower bound for the temperature which can be achieved in the experiments without applying further cooling techniques.
We compare the results of our approach with the measured experimental interference patterns over the whole parameter range. In contrast to previous work [7,8,12,13], we can take the full range of realistic heights of the lattice potential and the full trapping potential into account and our calculations are approximation-free and based on first principles.
After introducing the theoretical description of the quantum gas in section II, we discuss the Quantum Monte Carlo methods we developed in section III. This section can be skipped by the non-expert. In section IV we present results for the entropy in one dimension. We compare homogeneous systems with systems in the presence of realistic trapping potentials across the superfluid to Mott-insulating transition and we also study the Tonks-Giradeau gas. We make contact with analytic approximations where possible, and we study the validity of the local density approximation. Section V is dedicated to the entropy in two-dimensional homogeneous and trapped systems. We study the influence of the trapping potential, density, temperature, and of the actual experimental procedure on the visibility in section VI. We also compare the visibility obtained computationally with experimental data.

II. THEORETICAL DESCRIPTION
In order to prepare ultracold atoms into an optical lattice, the atoms are confined in an external trapping potential and cooled till they Bose-Einstein condense. Thereafter, additional lasers forming an optical lattice are turned on. One gradually increases the intensity of these lasers, slowly enough to remain in a low-energy state but also fast enough such that external influences are small. Adiabaticity can be checked by reversing the process, i.e. decreasing the intensity and comparing the state at a certain lattice height shows the same experimental signature as the corresponding state while switching on the lattice. Such experiments indicate that the loading of the atoms might to a good approximation be considered as adiabatic.
Atoms loaded in a sufficiently deep optical lattice are described by the Bose-Hubbard model [14], The notation i, j refers to the sum over nearestneighbor sites only. Bosons are created on site j by the operator b † j and the number of bosons on site j is counted by the number operator n j . The kinetic term describes hopping of the bosons with tunneling amplitude t, while the on-site repulsion has strength U . We will work in the canonical ensemble with a constant number of particles, N . The linear system size is L and we work in dimensions d = 1, 2. We restrict ourselves to a single-band Bose-Hubbard model, which is sufficient at low temperatures and lattice intensities around the superfluid-Mott transition, but it becomes approximative for very low repulsion or high temperatures.
The external trapping potential is included using ǫ i = v c r 2 , where v c describes the strength of the trapping and r the distance to the center of the lattice. If red detuned lasers are applied for creating the optical lattice potential the focusing of the lasers gives rise to an additional confinement. Under proper alignment of the two trapping potentials, the total trapping is given by ǫ(i) = mω 2 r 2 /2 [15] with where w is the waist of the lattice laser beam and V 0 its intensity expressed in single-photon recoil energies, Here λ is the wavelength of the lattice laser beam and m is the mass of the atoms. We took the waists isotropic and neglected corrections of the order of √ V 0 . The second term in Eq. (2) dominates already for moderately strong lattices.
The homogeneous Bose-Hubbard system shows a quantum phase transition from a superfluid to a Mott insulating state when the filling is commensurate [16]. The transition occurs at (U/t) c = 3.28(4) [17,18] in onedimensional systems, and at (U/t) c = 16.74(1) [19,20] in two-dimensional systems. Whereas in the superfluid state a continuous excitation spectrum exists at low energies, the Mott-insulator is characterized by a gap just above its ground state.
Analytically the entropy in the Bose-Hubbard model has been studied for non-interacting [12,21] and weakly interacting [12,13] bosons. The strongly interacting limit for the Mott-insulator and the Tonks regime have been considered for homogeneous [13,21] and trapped systems [13,22]. In three dimensions, the entropy deep in the superfluid and Mott insulating phases was calculated using effective masses [23]. Information on the trapped system was then obtained using the local density approximation (LDA). The Bose-Hubbard model can only be solved analytically in these limiting cases, where one is very deep in the superfluid or Mott insulating phase. Close to the phase transition numerical tools have to be employed, such as the QMC simulations performed here.

III. METHODS
In this chapter we discuss a number of methods to determine the entropy where E is the total (internal) energy of the system, β = 1/T the inverse temperature and F = − ln(Z)/β is the free energy. We use k B = 1. The main challenge is an accurate calculation of the partition function. It turns out that a combination of two methods discussed below gives the best results. Both are accurate, and do not involve any fitting nor noisy derivatives of numerical data, but they are efficient in a different temperature range. The results of both methods have been checked against each other for consistency. The flat histogram methods (Sec. III B) works best at high temperatures, while the thermodynamic integration method (Sec. III C) works better at low temperatures. This chapter is intended for the technically oriented readers and it is not necessary to read it in order to understand the discussion on the results in the next chapters.
A. The canonical worm algorithm All our simulations employ a canonical worm algorithm. A worm algorithm [24] in the path-integral representation is a quantum Monte Carlo algorithm where the decomposition of the partition function, is sampled indirectly by making local moves in the Green function sector, which is the extended configuration space of world lines with two open ends. Simulating the Bose-Hubbard model, we choose as diagonal part H 0 the potential energy, while the hopping terms are the perturbation V . In a canonical worm algorithm the operators of the equal-time Green function b i (τ )b † j (τ ) are propagated simultaneously. The extended partition function we sample reads where the terms W e (.) denote with E i = i 1 |H 0 |i 1 , and we have introduced sums over a complete basis set between any two off-diagonal operators. The terms W e are all positive and can thus be used as weights in a Monte Carlo sampling.
An efficient updating scheme has been presented in Ref. [25,26] and allows the straightforward computation of the kinetic and potential energy, density, compressibility, equal time Green function, and superfluid density. However, the partition function is not a thermodynamic average and is harder to compute.

B. Flat histogram methods
The goal of a flat histogram method is to obtain a density of states ρ(X) (where the coordinate X in classical simulations is usually the energy) directly by a random walk in X-space instead of performing a canonical simulation at fixed temperature. By sampling each value of X with a probability G(X) ∝ 1/ρ(X) we obtain a flat (constant) histogram H(X)ρ(X)G(X) = const.
In the Wang-Landau sampling scheme [27], a crude guess for G(X) is iteratively updated until it converges by a multiplicative factor f . During consecutive Wang-Landau iterations, f is reduced according to f → √ f when the current histogram H(X) is considered to be sufficiently flat. Then the histogram is reset, we have a more accurate estimator for the density of states, and the sampling restarts with the smaller f . The convergence of the scheme was proven by Zhou and Bhatt [28]. In particular, they showed that the minimum number of steps in each Wang-Landau iteration should scale as 1/ √ f . The generalization of the Wang-Landau scheme to quantum systems was discussed in Ref. [29]. Here, we generalize Eq. (4) of Ref. [29] to the path-integral formulation: The expansion order n corresponds to the number of kinks (particle hoppings) present in the path integral representation of a configuration. The original partition function (which is a function of the inverse temperature β) can be found back by setting λ = 1. The density of states g(n) corresponds to where W denotes the weight of a diagonal configuration, In such configurations all world-lines are continuous, and it occurs during the Monte Carlo run when the two open ends (worms) cancel each other on the same site and imaginary time.
Using the canonical worm algorithm [25,26], the density of states can be obtained as follows : A single Monte Carlo step is defined from a diagonal to a new diagonal configuration and has an acceptance factor q ′ . Taking the density of states into account, the acceptance factor should be modified to where g(x) is the density of states corresponding to the expansion order of the old configuration x. When the expansion order of the new configuration is larger than a predefined maximum expansion order, the update is rejected. When the Wang-Landau iteration is finished, we can obtain the partition function for all values of λ. For the Bose-Hubbard model that means that we obtain a whole set of partition functions through the scaling βt → βtλ, U/t → U/λt (thus βU is constant). A trap would also rescale as V → λV , which makes this scaling less useful in the trapped case and we use it just to obtain values at λ = 1. If we had worked in the grand-canonical ensemble, the chemical potential would scale analogously, and we lose all control over the particle number. Here we see a distinct advantage of the canonical ensemble over the grand-canonical ensemble.
The normalization of g(n) is fixed by calculating the partition function Z N in the canonical ensemble for the case without hopping (t = 0).

C. Thermodynamic integration
The second method we discuss is the thermodynamic integration method. We choose a set of inverse temperatures Then the partition function can be written as The partition function at infinitely hot temperatures Z 0 = Z β0=0 can be found by solving the combinatorial problem of placing N bosons on L lattice sites, The ratios in Eq. ( 12) can be estimated through the weights introduced in Eq. ( 5), where we sample all configurations σ at the temperature β j . In the canonical worm algorithm we have to measure with ∆β = β j − β j−1 and E d the time averaged potential energy of the configuration. We can thus compute the partition function at β j if we know the partition function at β j−1 . The accuracy of the scheme depends on the overlap between the system at β j and the one at β j−1 . If the overlap is small, the error on the partition function will increase rapidly and propagate systematically on to lower temperatures. In particular the first term β 1 should be chosen sufficiently close to zero, since there are only contributions if the expansion order is zero. The fluctuations in the (diagonal) energy in Eq. (15) are exponentially hard to control. We therefore choose our set of values of β such that ∆βE d < 1. At large values of U or large particle numbers, more β-points are required. When we are close to the ground state things get easier since there energy fluctuations are suppressed. In the limit that ∆β is infinitely small, the scheme reduces to an energy integration. Note however that the scheme remains exact when ∆β is finite.

D. Numerical strategy
Although we tried a number of alternatives, none of them were satisfactory. We briefly make some remarks about them. We tried to use and integrate the density fluctuation with respect to inverse temperature, but the scatter of the data was much bigger than the trend-line, which made this approach prohibitive without an adequate fitting procedure. Integrating the specific heat after differentiating the fitted curve through the energy was used in Ref. [23,30,31]. However, the division by the temperature is misbehaving at low T and the specific heat computed via the fluctuation formula c V = β 2 ( E 2 − E 2 ) is a quantity that converges slowly in the Monte Carlo simulation. We have also combined a grand-canonical directed loop algorithm in the path-integral representation [32]) with a quantum Wang-Landau reweighting scheme, though the fact that the Wang-Landau reweighting also changes the density, made this approach very cumbersome. A better attempt   was developing a canonical directed-loop algorithm in the stochastic series expansion [33] and combining it with a quantum Wang-Landau reweighting scheme. This approach has the advantage that one obtains all temperatures down to the one corresponding to the pre-chosen cut-off length at once. The drawback is that in a SSE representation the large values of U are also sampled, requiring extremely large orders even for moderate temperatures.
We obtained the highest accuracy by combining the two methods outlined in detail above. For high temperatures, up to βt ≈ 0.5, the combination of the canonical worm algorithm with the flat histogram (QWL) scheme is best and fast. For larger βt, and since we are interested in the entropy for a very large number of βt over a relatively small temperature range, thermodynamic integration is best. This method should only be used close to the ground state. Otherwise the fluctuations in Eq. (15) are hardly controllable which might lead to large systematic errors.
We compare the accuracy of both methods for a homogeneous one-dimensional system and for a trapped two-dimensional system in Table I and Table II,

IV. ISENTROPIC LINES IN ONE DIMENSIONAL SYSTEMS
We will start our discussion of the one-dimensional case with a homogeneous lattice and compare our results with analytic results in limiting cases. We will then gradually make the discussion more realistic (and more complicated) by including the parabolic confinement. The first step is discussing the entropy in a system of constant quadratic trapping, v c /t = const.. The second step is the case of an external parabolic trapping potential which is further strengthened by the focus of the lattice laser beams as is currently done in most experiments. We will compute the entropy with these two confining potentials starting in the superfluid and going to the Mott-insulator or the Tonks-Giradeau gas. We will also check the quality of a numerically based local density approximation and will find that in its regime of validity the speed-up in computing the entropy of trapped systems is considerable.

A. Homogeneous system
The dependence of the entropy on the temperature is closely related to the energy spectrum. Weakly interacting superfluid systems have an continuous energy spectrum. Therefore the entropy rises continuously with temperature.
Mott insulating states have an energy gap just above the ground state and entropy is exponentially suppressed up to temperatures of the order of k B T /U ≈ 0.1 [13,21]. Above the gap the bands of excited states lead to a finite entropy if temperature is high enough. Strongly interacting incommensurate systems also have gaps in the spectrum, but at substantially higher energies, signaled by a plateau in the The dependence of the entropy per site on the filling is shown in Fig. 2 for different interaction strengths at a moderately low temperature βt = 2. For weakly interacting systems (U/t = 2) the entropy reaches a maximum at half filling and decays monotonously for higher filling factors. Mott regions in strongly interacting systems (U/t = 12) appear as strong dips in Fig. 2, where the entropy is exponentially suppressed because the temperature is well below the gap. We also make a comparison with the atomic limit approximation, where only single particle-hole excitations are taken into account. The agreement with the Monte Carlo data is only qualitative, because the atomic limit misses higher order particlehole excitations which are important for U/t = 12. The atomic limit is expected to work well near commensurability, but already for densities 0.9 and 1.1 the deviation with Monte Carlo is considerable. The atomic limit approximation also incorrectly predicts a symmetric curve around n = 1. Figure 3 shows the entropy per site for a system of 40 particles on a lattice of 50 sites as a function of the interaction strength U/t and inverse temperature t/T = βt. Since we are away from integer filling we are always in the superfluid phase. At intermediate and high temperatures βt < 3, the entropy decreases when increasing the interaction strength. In contrast at lower temperature, βt 3, we see a non-monotonic behavior of the entropy with a minimum close to the critical interaction strength in a commensurate system. The same qualitative behavior was observed for a smaller lattice of L = 20 sites. The reason is the presence of the nearby quantum phase transition to a commensurate Mott state. The mass of the quasi-hole decreases when going away from the tip of the Mott lobe at (U/t) c to higher values of U/t. The quasi-hole absorbs more entropy and the entropy thus increases. Moving along isentropic lines we observe in units of t some mild heating when the initial temperature is reasonably high, βt < 3, but cooling if the initial temperature is sufficiently low, βt > 4.
The commensurate case in Fig. 4 shows the same behavior as Fig. 3 for low values of U/t < (U/t) c when we remain in the superfluid regime. However, when the Mott state develops, i.e. U/t > (U/t) c , the gap in the spectrum opens and the entropy at constant temperature decreases considerably. Along adiabatic trajectories the temperature in units of t shoots up near the transition point. This has been discussed for the homogeneous system in Ref. [21], from which the authors deduced that the temperature in present experiments must be of the order of U . However, Rey et al. [13] pointed out that in the presence of a parabolic confining potential less heating occurs in hard-core bosonic systems than in a homogeneous system. The next section addresses the same question for soft-core bosons.

B. Entropy distribution in a constant parabolic trapping potential
In the presence of an external trapping potential spatially separated quantum phases can coexist. In Fig. 5 we show two density profiles in the presence of a trapping potential. The first is for a superfluid state (Fig. 5a)) in The coexistence of spatially separated quantum phases can be understood in terms of a site-dependent chemical potential, the so-called local density approximation (LDA). Physical quantities are determined by using on each site the results obtained for a homogeneous system with the corresponding chemical potential. The sitedependent effective chemical potential provides a scan through the phase diagram. This approximation has been shown to work nicely for such quantities as the density or the variance in regions where the trapping potential varies slowly [34,35]. However, the LDA breaks down for steep trapping potentials and near the edges of Mott plateaus where numerical simulations are necessary to obtain reliable values [34]. To get a better understanding of the entropy distribution in an inhomogeneous system we developed a canonical and improved variant of the LDA, dubbed iLDA: we first calculate the exact density profile using a full QMC simulation for the trapped simulation. Then we take for every site the entropy from a homogeneous run corresponding to that density. This variant has the advantage that we start from the exact density profile, taking into account the rounding near the edges of the Mott plateaus due to the finite gradient of the trap. We tested the approach by comparing the total entropy of the trapped system calculated as the sum of the single sites to full numerical simulations. We found, as shown in Tab.III, that the iLDA can capture the qualitative trend of full calculations, but cannot reproduce the exact values. Using iLDA we obtained the entropy profiles for the density profiles of Fig. 5. In a system where superfluid and Mott-insulating regions coexist, we clearly see that the Mott-like regions are not able to accommodate entropy and the whole entropy is in the superfluid regions. If the whole system is superfluid (Fig. 5 (a)) the entropy varies only slowly from site to site and shows maxima at the filling close to n ≈ 1/2. Since the filling is larger than one in the center a dip in the entropy profile develops.
In Fig. 6 we show the entropy calculated in a full QMC calculation for a constant trapping potential. In Table IV we show the values of temperatures following an isentropic line as extracted from the data. The concept of 'adiabatic heating' is complicated by the different energy   Table IV for the corresponding temperatures.  Fig. 6. See Fig. 7 for the corresponding density profiles.
U/t βt βU 4 1.10(2) 4.4(1) 9 1.10(2) 9.9(2) 11 1.00 (5) (2) 18.2(7) scales (and units) used in the literature. It is important to define with respect to which energy scale the temperature is measured as illustrated in Table IV. We see that along isentropic lines in the superfluid phase (U/t < 10) temperature in units of the hopping t remains roughly constant. At higher temperatures there is some small heating in units of t, while at low temperatures we observe a little cooling.
This behavior can qualitatively be understood by looking at the density profiles and its variance along isentropic lines for different values of U/t, as shown in Fig. 7. Prior to the formation of a wide commensurate region, the temperature remains almost constant, for instance βt = 1.1 for U/t = 4 and U/t = 9 and βt = 1.0 for U/t = 11. Only when a considerable volume percentage of the system turns insulating the incommensurate edges cannot accommodate the entropy anymore which results in a rise of the temperature, i.e. βt = 0.8 for U/t = 15 to βt = 0.7 for U/t = 20. We checked that this effect is seen for a wide range of initial temperatures.
In units of the interaction strength U a cooling takes place. In particular, we see that for the chosen initial temperature the final temperature with respect to U stays below the temperature for which excitations in the Mott insulator are created in a homogeneous system. Therefore, the Mott insulator is stable up to the considered lattice height and the superfluid regions take most of the excitations.   Fig. 9. See Fig. 10 for the corresponding density profiles. 28.57 0.46(2) 13.1 (6) intensity). Table V gives the parameters for the temperature extracted along the isentropic line S = 13. We see that along this line the temperature in units of t increases while in units of U the temperature decreases. This means that the formation of the Mott-insulator is still possible starting at a low enough temperature, since the incommensurate regions take a lot of entropy. Figure 9 shows the entropy as a function of the temperature for different strengths of the optical lattice potential. Compared to Fig. 8 the number of atoms is increased to N = 140. This has the consequences that even at low temperature only small Mott-insulating regions can form and an incommensurate region survives in the center of the trap (Fig. 10). In Table VI we show the temperature along the isentropic line S = 20. In units of t we now get very strong heating, even in the superfluid phase, due to the high occupation which does not accommodate as much entropy as the low occupation region (cf. Fig 2).     Table VI for the corresponding temperatures.  However, in units of U we again find a temperature decrease.
One sees that the Mott-insulator is stable against the temperature change, since the temperature stays below 0.1U . Since the Mott-insulating regions are small and a broad incommensurate region survives in the center, the Mott transition does not play a central role in the behavior of the entropy curves.
Starting at a higher initial temperature βt = 0.8 and following the isentropic line S ≈ 50, the same qualitative effect of temperature increase in the units of t and decrease in the units of U can be seen in Tab. VII. However, at U/t ≈ 20, the temperature is still so high (k B T > 0.1U ) that no clean Mott-insulating region can be formed.

D. Tonks gas : one-dimensional trapped case
When the potential energy between the atoms of a onedimensional Bose-Hubbard model increases, the particles behave more and more like hard-core bosons. The limit of infinite repulsion and no multiple occupancies of a site, is called the Tonks-Giradeau gas. Quantities such as the en- ergy, average density, variance of the density can be computed accurately by assuming non-interacting fermions.
Other quantities such as the density matrix map to a noninteracting fermionic density matrix up to a phase factor coming from the Jordan-Wigner transformation. The experimental detection of the Tonks gas has demonstrated one of the fundamental concepts of quantum mechanics, namely the absence of a clear meaning of statistics in one-dimensional systems. The Tonks regime has been observed with [22] and without a lattice [36]. In the experiment with a lattice, the data was analyzed using fermionization, and good agreement with experiment was found in the region U/t > 5. The fermionization results were obtained at different temperatures along adiabatic lines. Temperature rose from βt = 2 for a lattice depth of V 0 = 4.6E R till βt = 0.77 for V 0 = 12E R as a consequence of the change in v c /t when ramping up the lattice. Consistently, it was argued in Ref. [7] that the temperature in the Tonks gas in a trap was of the order of the hopping t.
The Tonks problem was also studied by Monte Carlo simulations [37,38] at a low but constant temperature βt = 1, which is of the same order as the one used in the fermionization approach. The authors compared hardcore and soft-core bosons for a homogeneous lattice and for constant trapping. They found a gradual cross-over and found that the presence of a trap did not qualitatively change the Tonks onset.
The contrast in the experimental interference pattern was almost completely gone for the deepest lattices [22]. Although this is consistent with a strongly repulsive superfluid gas, one might fear that similar patterns are produced by either a Mott state or a thermal state due to a combination of soft-core bosons and an increased trapping depth.
Our analysis is again carried out in two steps. First, we calculated the entropy of a one-dimensional superfluid in the very low density limit trapped in a constant parabolic trap. From Fig. 11 we deduce that temperature remains remarkably constant. Thus, for a superfluid in the low density regime, adiabatic processes are (almost) isothermal when the external trapping is constant and weak. Second, we make the simulation of the experiment more realistic. We numerically evaluate the Bose-Hubbard parameters using the tight-binding approximation, and obtain the same parameters U and t as in Ref. [37]. In contrast to Ref. [37], we now also calculate the trapping parameter v c /t from the total axial trapping ω ax = 2π × 60 Hz [22] for all optical lattice depths. We find in Fig. 12 some heating, even in the low-density superfluid phase. Along an adiabatic line similar to the one taken in Ref. [22], we find that the temperature increases from βt = 2 at V 0 = 5E R (U/t = 6.9) to βt = 0.69 at V 0 = 12E R (U/t = 49.5). Thus, the temperature increase compares very well to the one calculated assuming hard-core bosons [22].
Summarizing, we confirm that the experiment has indeed reached the gradual cross-over toward the Tonks regime. The temperature remains of the order of the hopping t, even though a temperature increase of a factor of 3 is found due to the change in confinement strength v c /t when ramping up the lattice. It is exactly this increase in temperature that prevents the Mott domains from developing, since for V 0 = 12E R (U/t = 49.5, v c /t = 0.073 ) a broad Mott domain appears in the center of the trap for 15 particles and βt = 2. At a temperature βt = 0.69 the central density is 0.9 however.

A. Homogeneous two-dimensional superfluid
We again start our analysis in two dimensions with the homogeneous case. The phase transition from the superfluid to the Mott phase occurs for a n = 1 commensurate system at (U/t) c = 16.74(1) [19,20]. We calculate the entropy for a superfluid system close to commensurability, n = 0.8. In Fig. 13 the dependence of the entropy on the temperature and the interaction strength is shown. Its main behavior is similar to the one in a one-dimensional homogeneous case as reported in Fig. 1. At fixed interaction strength the entropy shows only a small increase for low temperatures. However, at a certain temperature it starts to increase strongly before it bends down again and a plateau is formed. This can as in the 1D case be related to the underlying band-structure in which first excitations in the lowest band can be created. Above a certain temperature the corresponding gap in the energy band-structure causes an intermediate saturation before at even higher temperatures further bands can be excited.
As a function of the interaction strength, the entropy shows -for constant but low temperature -a minimum  [19,20]. For low temperatures, we see an initial heating with increasing U/t, but around the transition point the presence of the Mott phase is felt and the system starts to cool, thanks to the lower effective mass. For larger values of U/t we see a further cooling for low temperatures since the Mott phase is far away and we go deeper inside the superfluid phase. The inset shows the entropy for low temperatures (same axes and symbols as in the main figure). close to the superfluid to Mott-insulator transition point of a commensurate system (inset of Fig. 13). As in the one-dimensional case this can be attributed to the effective mass change of the quasi-hole which has its maximum close to the phase transition point for the homogeneous system.
The density fluctuation shown in Fig. 14 is consistent with the behavior of the entropy shown in Fig. 13, using the relationship of Eq. (16). For infinitely hot tem- peratures the density fluctuation is independent of U/t (not shown). For low values of U/t density fluctuation goes down monotonously with βt. The normal-superfluid transition happens around βt ≈ 0.30(5) in our system of density 0.8 for a small lattice of 20×20 and U/t = 3. This transition belongs to the Kosterlitz-Thouless university class, and was studied in detail for the commensurate case in Ref. [19,20].
For large values of U/t we enter the quantum (thermal) critical regime determined by the quantum critical point U/t = 16.74 [19,20] and a minimum in the density fluctuations is reached. This is clearly seen for U/t = 15 around βt = 0.5. For larger values of U/t the minimum is reached at lower values of β. After this minimum is reached, the density fluctuations go slightly up with βt. The increase in the density fluctuations can be understood from the tendency of a dilute gas of vacancies (with respect to the n = 1 Mott state as a vacuum) trying to condense. In our canonical simulation, we will observe a tendency to increase the number of vacancies which will enhance pair formation. Thus theoretically the density fluctuations contain a lot of information about the system, but from the practical point of view the almost flat slopes in the quantum regime make this quantity a bad candidate for thermometry.

B. The superfluid to Mott insulating transition in a parabolic trap
In two dimensions we only consider trapping potentials taking into account the influence of the finite waist of the lattice laser on our sample. We use the parameters shown in Table. VIII. There are N = 200 bosons and the total system size is L = 20 × 20. These parameters are a compromise between increasing the trapping frequency while still obtaining meaningful results on a lattice of size L = 20 × 20. The filling in the system with N = 200 atoms for weak interactions is chosen such that it is close to n = 2 in the center for large U/t.
The entropy dependence on temperature is shown in Fig. 15 for different values of U/t. For low values of the temperature the slope for the curves is very flat. For higher temperature a clear increase in the entropy in the system can be seen. For low temperature the curves for different U/t almost coincide and only the curve for very   Table VIII for the corresponding values of  the trapping potential and Table IX for the corresponding  temperatures. strong interactions, i.e. U/t = 100, shows a considerably lower value of the entropy.
The behavior of the entropy can be explained by considering the density profiles of the system, shown in Fig. 16 along an adiabatic line S = 10. By increasing the optical lattice potential the central density goes down and we see the appearance of the n = 1 Mott-insulating region at U/t ≈ 25. These regions can be identified in the variance profile in Fig. 17. For stronger interactions, U/t ∼ 50, small Mott-insulating regions with n = 1 and n = 2 exist. For U/t = 70 the Mott-insulating regions cover already a large volume fraction of the system as clearly signaled in the variance profile. The formation of  Table VIII for the corresponding values of the trapping potential and Table IX for the corresponding temperatures. U/t βS=10t βS=10U βS=40t βS=40U βS=300t βS=300U the large Mott-insulating region causes the value of the entropy for the curve at U/t = 100 to lie below the others at low temperature. The temperature evolution along isentropic lines for different initial temperatures is shown more quantitatively in Table IX. When starting from a low temperature, we see almost no heating in units of t up to U/t = 70. The initially low entropy can be distributed over the remaining superfluid regions. Measuring the temperature in units of U leads to a cooling of the system below the value of U/t = 70.
In contrast, for interaction strengths larger than U/t = 70 almost the whole system is occupied by a commensurate region and the incommensurate regions have a negligible volume fraction. The energy cost to generate an excitation in this situation corresponds in the bulk to the large interaction energy and at the boundaries to the large potential energy cost resulting from the steep trapping potential. Hence the entropy cannot be well accommodated in the system and the temperature in units of t increases stronger than before when going to U/t = 100. In units of U the temperature first drops before it increases again or saturates for large U/t lattice potential.
Starting with a low temperature ensures that the temperature remains low enough for the presence of the Mott-insulating regions. In contrast when starting from a hot temperature (the right-most example in Table IX is already in the normal state), there is heating in units of t and only weak cooling in units of U . The quantum degeneracy regime is never reached. Note that the description by the one-band Hubbard model breaks down at such high temperatures.

VI. INTERPRETATION OF EXPERIMENTAL RESULTS
One of the standard experimental observation techniques consists of suddenly switching off the confinement and taking absorption images of the freely expanding gas after a finite flight time. The hereby resulting interference pattern is a reflection of the initial momentum distribution, where the factor w( k) is the Fourier transform of the Wannier function [14]. The use of the momentum profiles (and also of the visibilities) has been subject to debate, since in this quantity different effects as heating or the loss of coherence by stronger interactions have the same consequence. Furthermore, it is very difficult to extract information on the superfluid-Mott insulator transition point from these measurements in a trapping potential [17,34,39]. This is due to the spatially coexisting regions of superfluid and Mott-insulating character. A local measurement has to be implemented to obtain detailed information about the system. A first step into this direction has been taken by Fölling et al. [40] who measured the density in thin slices. Experimental progress was reported for a local measurement of the density in Ref.
[41], and another proposal was made in Ref. [42]. Such measurement techniques should be preferred over the visibility, which is only well suited for identifying the Mott and superfluid phase far away from the transition point. In Figs. 16 and 17 we show the typical evolution of the density profiles (with Mott plateaus for strongly repulsive systems) and variance profiles along isentropic lines.
A. Dependence of visibility on temperature and trapping potential Before we compare our results to the experimentally extracted quantities we would like to point out some features of the momentum distribution in a trapped systems at finite temperature. To do so we define the visibility V V = n max − n min n max + n min . where n min and n max are the values of the largest and smallest value of n( k).
Deep in the superfluid phase the darkest spot has almost n min ≈ 0, leading to a visibility close to unity, while in the Mott insulating phase the contrast between the brightest and the darkest spot is almost zero and one expects a very low visibility.
In Fig. 18 we show the visibilities for a trapped twodimensional Bose-Hubbard model at different densities. The calculation was done at a rather high temperature βt = 0.4. Even at this temperature regions with integer density and reduced compressibility exist. Precursors of the Mott-insulating regions can be seen in the density distribution and its variance around U/t ∼ 25. Looking first at the case of a constant trapping potential, we see that for the chosen number of particles a higher density leads to a higher visibility at low U/t. The visibility is lower at high U/t since the Mott region is larger. Taking the steepening of the trap into account we see that at low values of U/t the visibility is higher than the one obtained with the constant trapping for the considered particle number. The reason is the increased number of particles with density between one and two which form a superfluid edge between an n = 1 and an n = 2 Mott region. The n = 2 Mott region is absent in the calculation with the constant trapping potential. For the two curves labeled 'vc=mod' taking the change in the trapping potential into account, the visibility is well fitted by a logarithmic curve for large values of U/zt > 10 in agreement to the finding in Ref. [54]. The reason for this good agreement might be the suppression of the vol- ume fraction of the superfluid regions, such that only the Mott-insulating regions contribute to the decay of the visibility.
In Fig. 19 we show the visibility as a function of the inverse temperature for different strengths of the optical lattice potential. We see a strong decrease in the visibility if the temperature becomes of the order of βt ≈ 0.2 when the quantum regime is left. A comparable decrease in the visibility at low temperature by changing the interaction strength is only possible when it's very large, i.e. U/t = 200 leads to a visibility of the order of 0.2. However, we should note that at high temperatures the single-band approximation of Eq. (1) loses its validity.
We further show the visibility along different isentropic lines in Fig. 20 for the same parameters as taken in Table IX. The curves with low entropy (low initial temperature) both have a slope of approximately −1, but the onset value U/t is different. Even for the curve with high entropy, corresponding to an initial normal state, the data points for large value of U/t seem to be consistent with this slope, but differ for a small value of U/t (experiments show a constant visibility within error bars at low values of U/t).

B. Comparison to experimental results
We will now take the trap and the isentropic temperature change into account, and compare with experiments. This is the hitherto most realistic description done without numerical approximations, but for the 2D case we fail to take the same number of atoms and system sizes as in the experiment. Further we ignore time-of-flight collision effects.
We first compare our 2D calculations to the experimentally obtained visibilities. Up to now we have evaluated the visibility using n( k) at specific wave vectors k. In experiments instead an average over a square around the brightest and darkest spot is taken ( see Ref. [10] for the exact experimental procedure). The size of the area is a trade-off between signal and noise. The squares are chosen such that the contribution of the envelope of the Wannier functions (cf. Eq. (18)) is minimal, but not negligible. In Fig. 21  . This would correspond in the atomic limit to the presence of a small n = 2 plateau for the lower particle number and to a n = 2 plateau with almost half the number of particles for the higher number of particles [11].) and the visibility found in our model (labeled 'Th1' and 'Th2'). The curve 'Th1' corresponds to the same parameters as in Fig. 15 along the adiabatic line S = 10. The theoretical data is treated in the same way as explained in the text by taking the Wannier functions and averaging over a square of size 2π/5. The curve 'Th2' is a theoretical curve computed in the same way as 'Th1' and for the same system along an isentropic line S = 5 but for N = 50 particles such that the maximum occupation in the center of the trap does not exceed unity.
smaller value for the visibility. Fitting the slopes of the visibility with a power law gives roughly the same exponent for the theoretical and the experimental procedure, but the quality of the fit is better for the theoretical procedure.
Looking at three-dimensional experimental data, many groups find a visibility close to unity until U/zt ∼ 6.5, where z denotes the coordination number. For larger values of U/t the decrease in visibility is well approximated by V ∼ (U/zt) ν . The Zürich group [44,45] finds ν = −1.36 (5), while the Mainz group [46] finds ν ∼ −0.98 (7). The reason for the discrepancy is not fully understood. Similarly, the momentum distribution data on a log-log plot of the two-dimensional system studied in Ref. [47] were consistent with a behavior of (t/U ) and thus with a ν ≈ −1.
In Fig. 22 we show a comparison between our theoretical results for a two-dimensional system and an example of the experimental data [46] taken for a three dimensional system. We rescale the interaction strength by the coordination number taking the mean-field effect between the different dimensions into account. We compare two theoretical curves with N = 50 and N = 200 particles. In the strong interaction limit this corresponds to a density profile which has an n = 1 plateau for N = 50 and a profile with an n = 1 and an n = 2 plateau for N = 200.
The initial temperatures at U/t = 5 are βt = 1.7 and βt = 1.4, respectively. Note that the particle numbers and temperature do not directly correspond to the experimental ones. For a direct comparison of an experimental to a theoretical curve, three dimensional simulation with realistic particle numbers would be needed. However, we find a good agreement in Fig. 22 of the theoretical visibility treated in the way described above with the experimental data even up to relatively strong optical lattices. This shows that the experimental visibility can be explained by just taking the isentropic change of the temperature into account. In particular, the system does not have to leave the quantum regime to reach the low values of the visibility. The drop in the visibility can be explained by the formation of a broad Mott-insulating region (cf. Fig. 16 and 17).
In one-dimensional systems the procedure to get a quantity similar to the visibility has to be changed, since the interference pattern consists out of stripes. In Ref. [44] the superfluid to Mott-insulating transition in one-dimensional tubes was considered. The experimentalists extracted the coherent fraction of the atoms by taking the ratio of the content in the interference peaks and the background. We compared (not shown) the experimental data with the theoretical calculations where we treated the results for the momentum distribution in the same way. We found the same order of magnitude and the same qualitative trends ( e.g., the visibility for low lattice heights is between 0.5 and 0.6 in both cases). A detailed comparison however is hindered by the presence of many different parallel tubes in the experiments and by the large uncertainties stemming from the difficulty to fit the very sharp interference peaks.

VII. SUMMARY
An interpretation of experimental results on the superfluid to Mott insulator transition taking finite temperature effects into account has been given before by several authors, with differing conclusions. We have addressed the possibility of adiabatic temperature changes when ramping up the lattice for a Bose-Hubbard model using unbiased and first-principle quantum Monte Carlo simulations. We find that the expressions 'heating' and 'cooling' have to be taken with care, since they strongly depend on the unit in which temperature is measured.
For the homogeneous case, we found some small heating in units of t near commensurability in the superfluid phase for high temperatures βt < 0.5. For low temperatures, there is some small heating for low values of U/t < (U/t) c , but the system cools slightly for larger values of U/t when we are in the proximity of commensurability. At very low densities, temperature remains almost constant in the superfluid phase. When the density is commensurate, the temperature shoots up dramatically in units of t when the Mott gap opens in order to keep the accessible number of levels constant. This is in agree-ment with the findings of previous studies [7,13,21]. In a trapped system however the situation is different as already noted in [13] for hard-core bosons. We find that in a one-and a two-dimensional system in a trapping potential the main entropy contribution comes from incommensurate regions with low filling n ≈ 0.5. This is in good agreement with the finding for a three-dimensional system of the Amherst group where they find that the non-commensurate edges between the Mott plateaus become normal and accommodate the entropy [23]. We found that the temperature increased in units of t when the size of the commensurate regions is broader than several lattice sites in the case of constant trapping. For realistic trapping and densities around one or two, we found a temperature increase in units of t, even in the weakly interacting regime. However, in units of U temperature decreases or saturates, and its value lies deep in the quantum regime for the commensurate regions. Only if the incommensurate regions around filling n ≈ 0.5 are almost totally suppressed the entropy has to be taken by excitations in the commensurate regions or regions with higher filling, by which the Mott-insulating regions might be destroyed. Our result is in agreement with the results of F. Gerbier who finds that current experiments easily reach the thermal insulator regime (T < T * ≈ 0.2U , where Mott-like features persist but superfluidity is absent), and possibly the quantum region [43]. In contrast to the predictions of Ref. [8], we find that for realistic paramters no runaway temperature occurs.
Assuming adiabaticity in current experiments is in agreement with theory remaining in the quantum regime. In particular, the drop in the visibility of the interference pattern can be fully explained within this framework and no additional temperature rise has to be taken into account.
Future experiments using the spatially resolved measurements of the density and higher order correlations will be able to confirm the creation of Mott-insulating regions as the first evidence was reported in [48]. Unfortunately, we have seen that the integrated density fluctuations are not very sensitive to temperature changes when one is deep in the quantum regime.

VIII. OUTLOOK AND CONCLUSIONS
Our principal assumption that the loading of the lattice is approximately adiabatic needs to be verified considering the dynamical process at finite temperature. It will be equally important to extend our investigations to different and larger systems. A major goal will be to treat the full three-dimensional Bose-Hubbard system with a realistic number of particles. Systems with different types of particles will have to be addressed as well, since the development of a new energy scale that does not scale with U will have a negative influence on the possibility of reaching the ground state adiabatically. The visibility results of theory [49] at low temperature and experiment [45,50] are in disagreement for Bose-Fermi mixtures. It was believed that temperature is such that these systems are not in their ground state [49], and a study by Cramer et al. [51] hints at heating effects for a weak inter-species coupling. g the ground state. Antiferromagnetism and entropy were previously addressed for a homogeneous Fermi-Fermi system in Ref [31,52,53].
This work was motivated by the strongly different opinions that existed about temperature effects in the Bose-Hubbard model, ranging from a constant temperature to a runaway temperature. We have addressed all aspects of this discussion in homogeneous and trapped, one and two-dimensional Bose-Hubbard systems and the results of our work are uni-vocal: compatibility of experimental with numerical visibility curves (and density profiles) supports that the experimental initial temperature of the BEC is deeply in the superfluid phase, and that the quantum regime is not left even when adiabatically ramping up the optical lattice, even in the presence of a considerable Mott domain. The temperature scales almost linearly with U when the Mott domain is considerable in size, but temperature remains a very low fraction of the value of U . We have copied the experimental procedure of measuring the visibility in our simulations and found good agreement. It turns out that the Wannier function and the averaging over a small region of the interference pattern produce an almost constant shift of the visibility as a function of U/t. For deep lattices, the trapping potential becomes effectively steeper, and different Mott plateaus are forced to form. The main contributions to the visibility come from particle-hole excitations giving the visibility a slope of −1 as a function of U/t. The superfluid volume fraction is too low, but it still absorbs most of the entropy, as we could infer from the local density approximation. Our results strongly support that experimentalists have observed the superfluid and Mott phase and their crossover in a trapped system.