Critical points of the three-dimensional Bose-Hubbard model from on-site atom number fluctuations

We discuss how positions of critical points of the three-dimensional Bose-Hubbard model can be accurately obtained from variance of the on-site atom number operator, which can be experimentally measured. The idea that we explore is that the derivative of the variance, with respect to the parameter driving the transition, has a pronounced maximum close to critical points. We show that Quantum Monte Carlo studies of this maximum lead to precise determination of critical points for the superfluid-Mott insulator transition in systems with mean number of atoms per lattice site equal to one, two, and three. We also extract from such data the correlation-length critical exponent through the finite-size scaling analysis and discuss how the derivative of the variance can be reliably computed from numerical data for the variance. The same conclusions apply to the derivative of the nearest-neighbor correlation function, which can be obtained from routinely measured time-of-flight images.


Results
He near the lambda transition 9 . The curve is given by 0 0127, B ≈ 460 J/(mol K) ⋅ , and T c ≈ 2.17 K. The A coefficient for < T T c and T T c > is approximately given by −447 J/(mol K) ⋅ and −471 J/(mol K) ⋅ , respectively. Its asymmetry produces the lambda shape of the specific heat curve. Specific heat is not divergent at the critical point because the exponent α < 0. Measurements of specific heat as large as 120 J/ ⋅ (mol K) were reported. A similar shape is observed in our results for Var ∂ η .
boundary conditions, would require placing cold atoms in a three-dimensional optical lattice enclosed in an optical box trap. The tools needed for experimental creation of such a trap have been recently developed [17][18][19][20][21][22] . A comprehensive review of the properties of this model can be found in ref. 23 . In short, its physics depends on the filling factor n, i.e. the mean number of atoms per lattice site, and the ratio of the tunneling coupling J to the interaction energy U η = . J U / (2) We are interested in integer filling factors, for which there is a quantum phase transition 15 between the Mott insulator and the superfluid phase at the critical point η = J c c /U c . The system is in the Mott insulator phase for c η η < and in the superfluid phase for η η > c . The critical point at unit filling factor was theoretically studied via perturbative expansions 24,25 , QMC simulations 26 , the non-perturbative renormalization group approach 27 , and the projection operator approach 28 . The critical points at higher filling factors were studied perturbatively in ref. 24 for = n 2, 3 and in ref. 25 for an arbitrary filling factor. The results of all above-mentioned papers can be summarized as follows We will now summarize experimental results on the critical points of the 3D Bose-Hubbard model 5,29,30 . The experiment 5,31 is done in an optical lattice of wavelength λ = 852 nm with 87 Rb atoms in the | = = 〉 F m 2, 2 F state. Their s-wave scattering length 32 is 5.45 (26) nm. The position of the critical point for the unit filling factor was estimated to correspond to lattice heights between 10E R and 13E R , where the recoil energy E R is defined as k 2 2  /2m with π = k 2 /λ and m being the mass of the considered atom. This may be written as 11.5(9)E R , where the standard deviation has been estimated 33 by dividing maximum uncertainty of 1.5E R by 3. Using the formulas for J and U coefficients from 23 , we find c η to be 0.04 (1). The number reported in the bracket provides one standard deviation due to uncertainties in the lattice height and scattering length. It has been obtained through the uncertainty propagation formula. Nearly identical experimental setup was used in 29 . It was found there that the lattice heights corresponding to critical points for double and triple filling factors are 14.1(8)E R and 16.6(9)E R , respectively. Applying the same procedure as above, albeit with λ = 850 nm, we get c η for the double and triple filling factors 0.021(5) and 0.011(2), respectively. Finally, we come to ref. 30 , where again 87 Rb atoms, but in a different hyperfine state, are studied. It was found there that the superfluid-Mott insulator transition takes place at η = .
0 029(2) c for the unit filling factor. To put these results in perspective, we can compare them to the mean-field predictions, which for our system are 34 This yields η c equal to 0.029, 0.017, and 0.012 for n = 1, 2, and 3, respectively. Therefore, we see that more accurate experimental results are needed for characterization of beyond mean-field effects in the position of the critical points. It should be also said that in all above-mentioned experiments external harmonic trapping is imposed on the system. At the very least, it enhances finite-size effects making detailed comparision between experiments and the theory based on Hamiltonian (1) difficult. Such comparision is additionally complicated by the fact that the 3D Bose-Hubbard model captures only the leading-order behavior of cold atoms in optical lattices 35 . As a result, more precise experimental results on the critical points would presumably ask for a bit more advanced theoretical description of the system. Our method for locating the critical points should be immediately applicable to such non-standard versions of the 3D Bose-Hubbard model.
Besides critical points, quantum phase transitions are also characterized by critical exponents, which are supposed to be the same within a given universality class. The quantum 3D Bose-Hubbard model belongs to the universality class of the classical 4D XY model 15 . To the best of our knowledge, however, detailed studies of the critical properties of the latter model have not been presented in the literature so far. This is in sharp contrast to the properties of the lower-dimensional XY models, which have been studied in great detail 36 . The difficulties presumably arise here from the complexity of numerical studies of such a high-dimensional model. Furthermore, we note that the upper critical dimension of the XY model is four. This means that the mean-field theory, whose dynamical z and correlation-length ν critical exponents are given by should provide the lowest-order approximation to behavior of the 4D XY model. As it will turn out below, our relatively small system-size simulations are unable to capture corrections to the mean-field values of the critical exponents.
the observable of interest. We will study here variance of the on-site atom number distribution where the site index i can be chosen arbitrarily due to the translational invariance of our model. Such an observable can be conveniently computed with QMC algorithms. It can be also experimentally measured 37-39 . www.nature.com/scientificreports www.nature.com/scientificreports/ Alternatively, since we are actually interested in the derivative of the variance, one may focus on the derivative of the nearest-neighbor correlation function and use the mappinĝˆ † † η η η The higher-order zero-temperature perturbative calculations of the variance in the Mott insulator phase of the 3D Bose-Hubbard model were numerically performed in comprehensive work of Teichmann et al. 25 .
Quantum Monte Carlo simulations. We perform QMC simulations, which we briefly describe in the Methods section (see 40 for a cold-atom-oriented review of this subject). This allows us to study physics on the superfluid side of the transition, where dependence of the variance on η is most interesting for our purposes. Additionally, this approach allows us to get nonzero-temperature results, which is of interest from the experimental perspective.
We perform our studies in lattices of size 3 where the linear system size L 4 16 ≤ ≤ for = n 1, 2 and L 4 12 ≤ ≤ for n 3 = . Most of the time, we investigate systems at temperature k T B / = . U 0 02, where k B is the Boltzmann constant. Such temperature can be converted to Kelvins by considering typical experimental conditions. To do so, one may assume that the lattice of wavelength λ = 532 nm is populated by either 87 Rb or 174 Yb atoms having the s-wave scattering lengths of 5.45 nm and 5.56 nm, respectively 32,41 . Using then formulas from 23 , we find that U c /E R for n 1, 2, 3 = equals 0.47, 0.54, 0.59 for 87 Rb and 0.48, 0.55, 0.60 for 174 Yb. These numbers have been obtained by assuming that critical points are given by (3). Since we work near critical points, we may take U c /k B as the unit of temperature. In nano Kelvins, it equals 184, 211, 230 (93, 107, 117) for = n 1, 2, 3 in 87 Rb ( 174 Yb). We have checked that the same results can be obtained by generating Wannier functions and then integrating them over in order to get J and U coefficients 16 . www.nature.com/scientificreports www.nature.com/scientificreports/ www.nature.com/scientificreports www.nature.com/scientificreports/ We see from these calculations that temperature k T B /U 0 02 = . corresponds to a few nano Kelvins in typical rubidium and ytterbium systems. While such low temperatures are certainly experimentally challenging, it does not mean that our studies are completely free from nonzero-temperature corrections. For example, the 3D Bose-Hubbard model was studied through QMC simulations in ref. 26 at temperature about twenty times smaller than our k T B /U 0 02 = . . These studies were done at the unit filling factor in systems, whose sizes were similar to the ones used by us. The critical point was extracted from finite-size scaling of the excitation gap. The relative difference between the position of the critical point found in our work and the one from ref. 26 is about 0.5%. We will thus skip systematic discussion of nonzero-temperature effects in our computations. It should be stressed, however, that our approach to finding critical points can be applied to "warmer" systems as well, where nonzero-temperature scaling analysis can be deployed 3 .
The variance for the filling factors n = 1 and n = 2, 3 is presented in Figs 2 and 3, respectively. We see there its steep increase around the expected positions of critical points (3). To locate the points, where changes of the variance proceed most rapidly, we compute the first derivative of the variance. Such a procedure encounters a technical problem: the derivative is sensitive to fluctuations of data that is being differentiated (Fig. 4). Therefore, such data has to be smoothed first, which we do by fitting the Padé approximant of order (m, s) 42 Usefulness of this procedure is illustrated in Fig. 4, where we see that the derivative of the Padé approximant provides a smooth curve that can be easily subjected to detailed analysis. We mention in passing that the very same procedure could be applied to data for the variance coming from experimental measurements.
Fitting noisy data with (10) requires the approximation order to be adapted to the problem. While small orders may result in a bad fit due to insufficient number of fitting parameters, choosing too large orders causes problems as well. In the latter case, such extra flexibility leads to reproduction of noise-induced fluctuations of data points instead of averaging the fluctuations out.
Choosing the optimal order of the Padé approximant is not difficult in our computations. Indeed, we have found that for every combination of (L, n, T) parameters, there exists a stability island in the set of all reasonable approximation orders. By taking the order of approximation within the island, stable results for the variance and its derivative are obtained. We have found that typically our QMC data sets can be reasonably fitted with m 8 ≤ , ≤ s 9. We have also found that by considering a denser numerical grid, or by reducing the QMC noise by increasing the sample size, stable results can be obtained. We have applied this strategy for creation of Figs 2-8, where Padé approximants of the fixed (8,8) order are employed.
The system-size and temperature dependence of Var ∂ η is presented in Figs 5 and 6 for the filling factors n 1 = and = n 2, 3, respectively. The first thing we notice there is the lambda-shape of the plots reminiscent of the specific-heat plot in liquid 4 He (Fig. 1).
Then, we find that the first derivative of the variance has a maximum near the critical point on the superfluid side of the transition, say at η ⁎ . We see that ⁎ η shifts towards the critical point as the system size increases. The same is observed when temperature decreases. Moreover, Var( ) η ∂ η ⁎ grows with the system size and inverse tem- www.nature.com/scientificreports www.nature.com/scientificreports/ perature. All these observations suggest that the maximum of the derivative of the variance encodes the position of the critical point. This is not the first time when the derivative of an experimentally-accessible quantity is used for finding the critical point of the 3D Bose-Hubbard model. Indeed, the derivative of experimentally-measured visibility of the time-of-flight interference pattern was used for such a purpose as well 29 .
More quantitatively, we study the position of the maximum of Var ∂ η by fitting C to numerical results. The idea here is that the parameter a estimates the position of the maximum in the thermodynamically-large system (c 0 > ). The typical fits that we perform are shown in Fig. 7a-c, where the positions of the maxima have been extracted from Padé approximants of order (8,8). To check sensitivity of these results to the order of approximants (10) We have obtained ± . = .
± . = a n n n 0 03430 0 00008 for 1 0 02020 0 00006 for 2 0 01447 0 00003 for 3 , where the error bars are chosen to capture all the results in the parameter range given by (12). Standard deviations of fitted coefficients for = = m s 8 are typically a bit smaller (Fig. 7). A quick look at (3) reveals that these results provide positions of critical points, which makes us confident that ⁎ η converges to c η in the thermodynamic limit.  Table 1. Coefficients obtained by fitting (18) to ∂ η Var(η*). Error bars capture spread of the results due to the varying order of Padé approximation (12). QMC results for k B T/U = 0.02 have been used to prepare this table, the same system sizes as in Fig. 7 have been employed in the fitting.
www.nature.com/scientificreports www.nature.com/scientificreports/ which can be explained through the standard finite-size scaling argument 8,43 . Namely, we assume that near the critical point the system properties depend on the linear system size and the ratio between the linear system size and the correlation length ξ, say c where f is a non-singular scaling function. According to this ansatz, the position of the extremum of ∂ η Var scales and we see that the function h(L) captures system-size dependence of Var ∂ η at ⁎ η . The exponent in this equation matches the c parameter if we use the mean-field value of the critical exponent ν (5). More accurate studies are needed for checking if there are beyond-mean-field corrections to the critical exponent ν in the 3D Bose-Hubbard model. We believe that the key limitations here come from the relatively small system sizes that can be numerically handled.
Finally, for the sake of completeness, we provide the results for the fitting parameter b again under the variation of the order of Padé approximation (12) ± . = .
± . = . b n n n 0 059 0 005 for 1 0 029 0 002 for 2 0 028 0 002 for 3 (17) Next, we will discuss scaling of ⁎ Var( ) η ∂ η with the linear system size. We fit to numerics. Typical results that we obtain are presented in Fig. 7d-f, where again Padé approximants of order (8,8) have been employed. Next, we quantify influence of the order of Padé approximation on these results. Proceeding similarly as with (13), (14) and (17), we get the results summarized in Table 1. All of them suggest that ⁎ η ∂ η Var( ) slowly increases with the linear system size reaching a finite value in the thermodynamic limit. This has an interesting consequence that can be readily spotted in Figs 5 and 6.
Namely, we see that the curves showing ∂ η Var for different system sizes at constant temperature cross near the critical point. This can be explained by the finite-size scaling ansatz (15) and take into account that h(L) weakly depends on the linear system size reaching a finite value in the thermodynamic www.nature.com/scientificreports www.nature.com/scientificreports/ limit. The latter remark follows from the fact that h(L) is proportional to , which we have just discussed. We mention in passing that similar-looking crossing of curves near the critical point was used for finding the position of the critical point from QMC data for the excitation gap 26 .
Further insights into Var ∂ η can be obtained by setting k T B /U 0 = and using the Feynman-Hellmann theorem to arrive at 14 where  is the ground-state energy per lattice site. This expression is closely linked to the one for specific heat oftentimes studied in the context of classical phase transitions (see e.g. Fig. 1 and the discussion around it). Indeed, specific heat per lattice site can be written as 44 where  is the free energy per lattice site. Its singular part is typically assumed to scale as where α is the specific heat critical exponent. A quick look at (19) and (20) reveals the mapping between the two expressions. It is then unsurprising that the singular part of  is usually assumed to scale as The exponent α is linked to the z and ν critical exponents through the quantum hyperscaling relation where d is the system's dimensionality 3 .
Var c in the infinite system, which would imply that h L L ( ) / ∼ α ν . As a result, 0 α = would be compatible with our numerics in the large-L limit. Such a value can be obtained by putting mean-field critical exponents (5) into (22). There are, however, at least two reasons to be cautious here.
First, the upper critical dimension of the Bose-Hubbard model is three and so it is expected that there will be corrections to the mean-field scaling laws. As a result, it is unclear to us what is the actual value of α.
Second, even if α would be zero, the presence of logarithmic singularities in the derivatives of the ground-state energy could not be ruled out without detailed analysis. For example, such a situation takes place in the one-dimensional quantum Ising model, where 0 α = due to z d 1 ν = = = 45,46 . The singular part of its ground-state energy per lattice site Ising  turns out to be proportional to g g g g ( ) ln www.nature.com/scientificreports www.nature.com/scientificreports/ thermodynamically-large chain, where g is the magnetic field driving the transition and g c is the critical point. As a result, ∂ g 2 Ising  diverges logarithmically with | − | g g c . In the finite chain, there is the extremum of  . The latter property differs from what we seem to observe in the 3D Bose-Hubbard model.
One possible explanation of our puzzling observation that Var( ) ⁎ η ∂ η reaches a finite value in the thermodynamic limit might be that the system sizes that we consider are much too limited rendering our extrapolations unreliable. Still, the use of our fitting results for interpolation purposes should be very well justified and useful.
Finally, we would like to briefly discuss dependence of our results on the filling factor n. The idea that we explore here comes from Teichmann et al. 25 , where it was found through perturbative studies that deeply in the Mott insulator phase the variance of the on-site atom number operator is a function of This implies that Var  ∂ η should be a function of (23) as well. Our QMC numerics, which we present in Fig. 8, perfectly follows these predictions away from the critical point on the Mott insulator side of the transition. We also see in this figure that the mapping η η fails a bit in the superfluid phase for low-n data that we explore. Nonetheless, judging from quite good overlap between the n = 2 and 3 results, it is reasonable to expect that the mapping will be accurately supported by numerical simulations in both phases in the limit of n 1  . Further studies are needed for establishing this observation.

Discussion
We have studied equilibrium properties of the 3D Bose-Hubbard model focusing our attention on the variance of the on-site atom number operator and its derivative with respect to the parameter driving the superfluid-Mott insulator transition. Our results have been obtained in systems with the mean number of atoms per lattice site equal to one, two, and three. They come from Quantum Monte Carlo simulations.
The key finding of this work is that the derivative of the variance has a pronounced maximum close to critical points. For example, in a very small lattice of linear size 4, when the number of atoms equals the number of lattice sites, the position of the maximum estimates the position of the critical point with 10% relative accuracy (Fig. 4). The mismatch between the two decreases quadratically with the linear system size and it can be further suppressed by simple extrapolation to the thermodynamic limit.
Besides discussing the position of the critical point, which is an interesting albeit non-universal feature, we have found that even in small systems the critical exponent ν can be extracted from the finite-size shift of the maximum of the derivative of the variance. This is interesting because knowledge of this exponent can provide important information about the universality class of the 4D XY model. This is the least-studied universality class of the XY model. Limited knowledge of its properties stems from numerical shortcomings, clearly seen in our work, and difficulties in finding physical systems, where it can be experimentally approached. The latter can be done in condensed matter and atomic physics setups. In the condensed matter context, it was proposed that some properties of either strongly underdoped cuprate superconductors or 4 He in nanoporous media can be captured by the 4D XY universality class [47][48][49] . In the atomic physics context, cold atoms in a three-dimensional optical lattice are the best example of a system whose scaling properties should mirror those of the 4D XY model.
We view cold atom setups, simulating the 3D Bose-Hubbard model, as the cleanest and most promising platform for future quantitative studies of the 4D XY universality class. In fact, measurements of the on-site atom number fluctuations in the 3D Bose-Hubbard systems have been recently reported 38,39 . Direct comparision of these results to our findings is difficult because setups studied in 38,39 are non-uniform due to the external trapping potential adding a local chemical-potential-like term to Hamiltonian (1). We are hopeful, however, that blending of techniques presented in these references with the recent optical box trapping advances can lead to successful creation of the homogeneous 3D Bose-Hubbard quantum simulator. Such a system could be large-enough to overcome small-size limitations plaguing numerical simulations. As a matter of fact, quantum simulators, at the very least, are supposed to do just that.

Methods
We use the Directed Worm Algorithm from the ALPS software package 50,51 . This algorithm samples the path-integral representation of a density matrix of a grand canonical ensemble (GCE) with configurations called worldlines. Since we work with systems having fixed filling factor n, computation of any average in a lattice of the linear size L requires rejecting those worldlines, where the total numer of particles differs from nL 3 . To improve the sample count of the remaining fraction, the chemical potential is adjusted to set the expected GCE density of the system to n particles per site. The statistical error of the determined variance is significantly reduced by adopting periodic boundary conditions, where observables do not depend on the lattice site. As a result, the variance of the on-site atom number operator can be averaged over the resulting ensemble and over all L 3 lattice sites, which is exactly what we do.
Due to the amount of computational power needed, our QMC simulations are limited to system sizes and temperatures discussed below equation (9). Early symptoms of these limitations can be spotted in Figs 5 and 6. We see there that for the largest systems considered, there is a small warp in the derivative of the variance slightly to the left of the dotted lines marking positions of critical points. Therefore, it is important to check that positions of the maxima, which we extensively study, are stable under change of parameters of our QMC simulations. Several tests are thus performed. First, we vary the total number of later averaged worldlines reaching typically the level of 10 7 to 10 8 . Second, when generating the worldlines, only every m-th trajectory is included into the final ensemble (if it additionally contains nL 3 particles), to ensure that subsequent worldlines in the ensemble are independent