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

TL;DR

This paper demonstrates that the critical points of the 3D Bose-Hubbard model can be precisely identified using the derivative of on-site atom number fluctuations, which are experimentally accessible, and extracts critical exponents through finite-size scaling.

Contribution

It introduces a method to accurately locate quantum critical points using atom number fluctuation derivatives and validates it with Quantum Monte Carlo simulations.

Findings

Precise determination of critical points for superfluid-Mott insulator transition.

Extraction of the correlation-length critical exponent from fluctuation data.

Applicability of the method to experimentally measured quantities like time-of-flight images.

Abstract

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.

Tables

Table 1. Table 1: Coefficients obtained by fitting ( 18 ) to ∂ η Var ​ ( η ∗ ) subscript 𝜂 Var superscript 𝜂 \partial_{\eta}\!{\rm Var}(\eta^{*}) . 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 subscript 𝑘 𝐵 𝑇 𝑈 0.02 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.

	A	B	C
$n=1$	$23.3\pm 0.6$	$-44\pm 4$	$1.2\pm 0.1$
$n=2$	$41\pm 1$	$-64\pm 4$	$1.1\pm 0.1$
$n=3$	$54\pm 1$	$-101\pm 7$	$1.3\pm 0.1$

Equations

$$\hat{H}=-J\sum _{⟨\mathbf{i},\mathbf{j}⟩}({\hat{a}}_{\mathbf{i}}^{†}{\hat{a}}_{\mathbf{j}}+{\hat{a}}_{\mathbf{j}}^{†}{\hat{a}}_{\mathbf{i}})+\frac{U}{2}\sum _{\mathbf{i}}{\hat{n}}_{\mathbf{i}}({\hat{n}}_{\mathbf{i}}-1),$$

$$[{\hat{a}}_{\mathbf{i}},{\hat{a}}_{\mathbf{j}}^{†}]={\delta }_{\mathbf{ij}},[{\hat{a}}_{\mathbf{i}},{\hat{a}}_{\mathbf{j}}]=0,{\hat{n}}_{\mathbf{i}}={\hat{a}}_{\mathbf{i}}^{†}{\hat{a}}_{\mathbf{i}},$$

$$\eta =J/U.$$

$$\eta_{c}\approx\left\{\begin{array}[]{ll}0.034&\ \text{for}\ n=1\\ 0.02&\ \text{for}\ n=2\\ 0.014&\ \text{for}\ n=3\end{array}\right..$$

$${\eta }_c=\frac{2n+1-2\sqrt{n(n+1)}}{6}.$$

$$z=1,\nu =1/2,$$

$$\mathrm{Var}=⟨{\hat{n}}_{\mathbf{i}}^2⟩-⟨{\hat{n}}_{\mathbf{i}}{⟩}^2,$$

$$\frac{\partial }{\partial \eta }\mathrm{Var}=6\eta \frac{\partial }{\partial \eta }C,C={⟨{\hat{a}}_{\mathbf{i}}^{†}{\hat{a}}_{\mathbf{j}}+a_{\mathbf{j}}^{†}{\hat{a}}_{\mathbf{i}}⟩|}_{⟨\mathbf{i},\mathbf{j}⟩},$$

$$\langle\hat{n}^{2}_{\bf i}\rangle-\langle\hat{n}_{\bf i}\rangle^{2}\approx\left\{\begin{array}[]{ll}24\eta^{2}+3960\eta^{4}&\ \text{for}\ n=1\\ 72\eta^{2}+33408\eta^{4}&\ \text{for}\ n=2\\ 144\eta^{2}+131400\eta^{4}&\ \text{for}\ n=3\end{array}\right..$$

$$L^3,$$

$$\mathrm{Var}(\eta )=\frac{{\sum }_{i=0}^m{A}_i{\eta }^i}{1+{\sum }_{i=1}^s{B}_i{\eta }^i}.$$

$$a+b{L}^{-c}$$

$$7\le m,s\le 9.$$

$$a=\left\{\begin{array}[]{ll}0.03430\pm 0.00008&\ \text{for}\ n=1\\ 0.02020\pm 0.00006&\ \text{for}\ n=2\\ 0.01447\pm 0.00003&\ \text{for}\ n=3\end{array}\right.,$$

$$c=\left\{\begin{array}[]{ll}2.11\pm 0.08&\ \text{for}\ n=1\\ 1.97\pm 0.08&\ \text{for}\ n=2\\ 2.23\pm 0.05&\ \text{for}\ n=3\end{array}\right.,$$

$${\partial }_{\eta }\mathrm{Var}=h(L)f\left(\frac{L}{\xi }\right),\xi \sim |\eta -{\eta }_c{|}^{-\nu },$$

$$|{\eta }^{\ast }-{\eta }_c|\sim L^{-1/\nu }$$

$$b=\left\{\begin{array}[]{ll}0.059\pm 0.005&\ \text{for}\ n=1\\ 0.029\pm 0.002&\ \text{for}\ n=2\\ 0.028\pm 0.002&\ \text{for}\ n=3\end{array}\right..$$

$$A+B{L}^{-C}$$

$$\frac{\partial }{\partial \eta }\mathrm{Var}=-2\eta \frac{{\partial }^2}{\partial {\eta }^2}\left(\frac{\mathcal{E}}{U}\right),$$

$$-T\frac{{\partial }^2}{\partial T^2}\mathcal{F},$$

$$|\eta -{\eta }_c{|}^{2-\alpha }.$$

$$\alpha =2-\nu (d+z),$$

$$\tilde{\eta }=\sqrt{n(n+1)}\eta .$$

Full text

The field of quantum phase transitions is significant for both fundamental and practical reasons Sac (a, b); Con . On the fundamental side, it provides insights into rich physics of the most complicated many-body quantum systems whose description challenges the most advanced theoretical and numerical studies. On the practical side, understanding of quantum phase transitions is critical for explanation of properties of numerous materials of potential technological importance such as high-temperature superconductors. One of the key obstacles to progress in the field of quantum phase transitions is our inability to efficiently solve the models describing strongly-correlated systems.

The recent progress in cold atomic physics suggests that it could be possible to approach this problem from another angle. Namely, one can use cold atomic setups as hardware for simulation of condensed matter models. The idea to use an easy-to-control quantum system to simulate another quantum system that is hard-to-study in conventional setups dates back to Feynman Fey . Its proof of principle experimental demonstration was achieved about two decades ago in the cold atom experiment Gre (a). Ever since this publication, the field of quantum simulation has attracted a lot of attention by providing both a fertile ground for interdisciplinary research and a plausible promise of delivering flexible experimental hardware for solving various challenging models.

At the center of these studies is the implicit assumption that such models can be quantitatively studied in cold atom simulators. Despite fascinating experimental efforts that have been reported so far Lew ; Blo , the current characterization of quantum phase transitions in best cold atom simulators is nowhere near the level of accuracy that is found in top condensed matter experiments.

To justify this statement and to place the research discussed in this paper in a broader context, we compare measurements of the lambda transition in liquid 4He to measurements of the superfluid–Mott insulator transition in the two-dimensional (2D) cold atom cloud. Such a comparison is meaningful because the two systems are expected to belong to the same universality class, the one of the 3D $XY$ model.

The experimental studies of the lambda transition in 4He have started about a century ago and continue to this date. The good reason for this long-term interest is that they provide the most stringent tests of the renormalization-group (RG) theory, which is the cornerstone of our understanding of phase transitions. The RG predictions are best established very close to the critical point Car (a), where it is difficult to do measurements. One of the key current problems in accurate testing of the RG theory at the lambda transition is that gravity broadens the transition (the transition temperature depends on pressure, which varies across the sample due to the gravitational field). To remove this obstacle, the measurements were taken in a zero-gravity environment in a Space Shuttle. This allowed for checking the RG predictions up to the relative distances $|1-T/T_{c}|\approx 10^{-9}$ Lip , where $T_{c}\approx 2.17$K marks the position of the critical point ($T_{c}$ can be extracted out of the specific heat data from ref. He (4) with relative accuracy of about $10^{-8}$; ref. Lip does not provide its value). These measurements, whose main outcome we briefly summarize in Fig. 1, provided the specific heat critical exponent $\alpha=-0.0127(3)$. Such a value is close to the exponent that one can theoretically obtain from the studies of the 3D $XY$ model Vic (a). It provides significant evidence about the universality class of the lambda transition. On the less optimistic side, these experiments are quite complicated in execution as they are performed over $200$ km above the Earth. It is thus hard to expect that the progress in these studies can be sustained.

These results can be compared to the state of the art experiments in the 2D cold atom clouds placed in a periodic potential formed by standing laser beams (an optical lattice) Por (a, b). These experiments aimed at determination of the position of the critical point from the disappearance of the condensate fraction in the Mott insulator phase. They located the critical point with $10$%–$15$% relative accuracy, which precluded experimental studies of the cold atom system as near the critical point as in the above-reported measurements in liquid 4He. Moreover, the condensate fraction in these experiments was estimated by the relative size of the main peak in time-of-flight images. Such an estimation is heuristic, as the condensate fraction is rather defined as the largest eigenvalue of the single-particle density matrix, which was not studied in refs Por (a, b). Finally, these references did not report any measurements of the critical exponents, which excluded experimental identification of the universality class of the studied system.

On the more positive side, recent developments in the experimental manipulations of cold atomic samples suggest that substantial progress in the studies of phase transitions in these systems could be expected. For that purpose, we need efficient tools that can be used to extract quantitative rather than qualitative features of phase transitions out of experimental data. This observation motivates the work presented here.

The purpose of this work is to discuss a practical scheme leading to precise determination of the critical points of the 3D Bose-Hubbard model. In addition to that, the method that we explore provides the correlation-length critical exponent, albeit with lower accuracy. The idea behind these studies was proposed in ref. BDS , where we tested it in the two-dimensional Bose-Hubbard model. The 3D Bose-Hubbard model has been intensively discussed over the past decade because it provides basic description of the properties of cold atoms in optical lattices. Despite all this interest, precision experimental results on its critical properties are hard-to-find in the literature.

The outline of this paper is the following. We will briefly summarize below basic features of the 3D Bose-Hubbard model, state of the art experimental and theoretical results on positions of its critical points, and our scheme for finding them. Then, we will discuss our results coming from Quantum Monte Carlo (QMC) simulations.

Results

The model. The Hamiltonian of the 3D Bose-Hubbard model reads Fis ; Jak

[TABLE]

where $\langle{\bf i},{\bf j}\rangle$ describes nearest-neighbor sites ${\bf i}$ and ${\bf j}$ in the cubic lattice and periodic boundary conditions are imposed on the creation and annihilation operators. A physical realization of such a model, albeit with open 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 Rai ; Had ; Hom (a, b, c, d).

A comprehensive review of the properties of this model can be found in ref. Kru . 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$

[TABLE]

We are interested in integer filling factors, for which there is a quantum phase transition between the Mott insulator and the superfluid phase at the critical point $\eta_{c}=J_{c}/U_{c}$. The system is in the Mott insulator phase for $\eta<\eta_{c}$ and in the superfluid phase for $\eta>\eta_{c}$ Fis .

The critical point at unit filling factor was theoretically studied via perturbative expansions Mon ; Eck , QMC simulations San , the non-perturbative renormalization group approach Dup , and the projection operator approach Sen . The critical points at higher filling factors were studied perturbatively in ref. Mon for $n=2,3$ and in ref. Eck for an arbitrary filling factor. The results of all above-mentioned papers can be summarized as follows

[TABLE]

We will now summarize experimental results on the critical points of the 3D Bose-Hubbard model Gre (a); Ger ; Ket . The experiment Gre (a, b) is done in an optical lattice of wavelength $\lambda=852$ nm with 87Rb atoms in the $|F=2,m_{F}=2\rangle$ state. Their s-wave scattering length is $5.45(26)$ nm Jul . 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 $\hbar^{2}k^{2}/2m$ with $k=2\pi/\lambda$ 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 by dividing maximum uncertainty of $1.5E_{R}$ by $\sqrt{3}$ sqr . Using the formulas for $J$ and $U$ coefficients from Kru , we find $\eta_{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 Ger . 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 $\lambda=850$ nm, we get $\eta_{c}$ for the double and triple filling factors $0.021(5)$ and $0.011(2)$, respectively. Finally, we come to ref. Ket , where again 87Rb atoms, but in a different hyperfine state, are studied. It was found there that the superfluid–Mott insulator transition takes place at $\eta_{c}=0.029(2)$ 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 Sto

[TABLE]

This yields $\eta_{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 (LABEL:H) 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 Kub . 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 Fis . 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 Vic (b). 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 $\nu$ critical exponents are given by

[TABLE]

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

[TABLE]

where the site index ${\bf 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 Tak ; 3DB ; exp . 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 mapping

[TABLE]

which can be obtained through the Feynman-Hellmann theorem (it is strictly valid in systems that are either thermodynamically-large or periodic). The nearest-neighbor correlation function can be extracted out of the routinely measured time-of-flight images in cold atom simulators Blo ; 3DB .

The variance can be perturbatively calculated deeply in the Mott insulator phase at temperature of absolute zero Kru

[TABLE]

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. Eck .

Quantum Monte Carlo simulations. We perform QMC simulations, which we briefly describe in the Methods section (see Pol 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 $\eta$ 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

[TABLE]

where the linear system size $4\leq L\leq 16$ for $n=1,2$ and $4\leq L\leq 12$ for $n=3$. Most of the time, we investigate systems at temperature $k_{B}T/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 $\lambda=532$ nm is populated by either 87Rb or 174Yb atoms having the s-wave scattering lengths of $5.45$ nm and $5.56$ nm, respectively Jul ; Yba . Using then formulas from Kru , we find that $U_{c}/E_{R}$ for $n=1,2,3$ equals $0.47$, $0.54$, $0.59$ for 87Rb and $0.48$, $0.55$, $0.60$ for 174Yb. 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 87Rb (174Yb). 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 Jak .

We see from these calculations that temperature $k_{B}T/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. San at temperature about twenty times smaller than our $k_{B}T/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. San 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 Con .

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)$ Pad

[TABLE]

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 $8\leq m,s\leq 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 $\partial_{\eta}\!{\rm 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 4He (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 $\eta^{*}$. We see that $\eta^{*}$ shifts towards the critical point as the system size increases. The same is observed when temperature decreases. Moreover, $\partial_{\eta}\!{\rm Var}(\eta^{*})$ grows with the system size and inverse temperature. 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 Ger .

More quantitatively, we study the position of the maximum of $\partial_{\eta}\!{\rm Var}$ by fitting

[TABLE]

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 Figs 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 done calculations for

[TABLE]

We have obtained

[TABLE]

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 $\eta^{*}$ converges to $\eta_{c}$ in the thermodynamic limit.

For the $c$ coefficient under constraint (12), we have got

[TABLE]

which can be explained through the standard finite-size scaling argument Car (b, a). 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 $\xi$, say

[TABLE]

where $f$ is a non-singular scaling function. According to this ansatz, the position of the extremum of $\partial_{\eta}\!{\rm Var}$ scales

[TABLE]

and we see that the function $h(L)$ captures system-size dependence of $\partial_{\eta}\!{\rm Var}$ at $\eta^{*}$. The exponent in this equation matches the $c$ parameter if we use the mean-field value of the critical exponent $\nu$ (5). More accurate studies are needed for checking if there are beyond-mean-field corrections to the critical exponent $\nu$ 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)

[TABLE]

Next, we will discuss scaling of $\partial_{\eta}\!{\rm Var}(\eta^{*})$ with the linear system size. We fit

[TABLE]

to numerics. Typical results that we obtain are presented in Figs 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 $\partial_{\eta}\!{\rm Var}(\eta^{*})$ 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 $\partial_{\eta}\!{\rm Var}$ for different system sizes at constant temperature cross near the critical point. This can be explained by the finite-size scaling ansatz (15), if we note that $\partial_{\eta}\!{\rm Var}(\eta_{c})\sim h(L)$ and take into account that $h(L)$ weakly depends on the linear system size reaching a finite value in the thermodynamic limit. The latter remark follows from the fact that $h(L)$ is proportional to $\partial_{\eta}\!{\rm Var}(\eta^{*})$, 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 San .

Further insights into $\partial_{\eta}\!{\rm Var}$ can be obtained by setting $k_{B}T/U=0$ and using the Feynman-Hellmann theorem to arrive at BDS

[TABLE]

where $\cal E$ 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 Bax

[TABLE]

where $\cal F$ is the free energy per lattice site. Its singular part is typically assumed to scale as $|T-T_{c}|^{2-\alpha}$, where $\alpha$ 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 $\cal E$ is usually assumed to scale as

[TABLE]

The exponent $\alpha$ is linked to the $z$ and $\nu$ critical exponents through the quantum hyperscaling relation

[TABLE]

where $d$ is the system’s dimensionality Con .

Combining (19) with (21), one gets $\partial_{\eta}\!{\rm Var}\sim|\eta-\eta_{c}|^{-\alpha}$ in the infinite system, which would imply that $h(L)\sim L^{\alpha/\nu}$. As a result, $\alpha=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 $\alpha$.

Second, even if $\alpha$ 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 $\alpha=0$ due to $z=\nu=d=1$ Lie ; Pfe . The singular part of its ground-state energy per lattice site ${\cal E}_{\text{Ising}}$ turns out to be proportional to $(g-g_{c})^{2}\ln|g-g_{c}|$ in the thermodynamically-large chain, where $g$ is the magnetic field driving the transition and $g_{c}$ is the critical point. As a result, $\partial^{2}_{g}{\cal E}_{\text{Ising}}$ diverges logarithmically with $|g-g_{c}|$. In the finite chain, there is the extremum of $\partial^{2}_{g}{\cal E}_{\text{Ising}}$ at say $g^{*}$. When the system size increases, $g^{*}\to g_{c}$ and $|\partial^{2}_{g}{\cal E}_{\text{Ising}}(g^{*})|\to\infty$. The latter property differs from what we seem to observe in the 3D Bose-Hubbard model.

One possible explanation of our puzzling observation that $\partial_{\eta}\!{\rm Var}(\eta^{*})$ 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. Eck , 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

[TABLE]

This implies that $\partial_{\tilde{\eta}}\!{\rm 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 $\eta\mapsto\tilde{\eta}$ 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\gg 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 $\nu$ 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 4He in nanoporous media can be captured by the 4D $XY$ universality class Fra ; Bro ; Egg . 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 3DB ; exp . Direct comparision of these results to our findings is difficult because setups studied in 3DB ; exp are non-uniform due to the external trapping potential adding a local chemical-potential-like term to Hamiltonian (LABEL:H). 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 alp (a, b). 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 from each other. We check sensitivity of our results to this so-called skip $m$ parameter by changing it from $64$ to $256$. Third, we fit Padé approximants to data sets that differ from each other by range and grid density of the parameter $\eta$. All these tests make us confident that the results that we present for the maxima of the derivative of the variance are well-converged.

Acknowledgments

We thank Yoshiro Takahashi for a conversation about ref. 3DB . We also thank Tama Ma for discussions about the details of the ALPS implementation of the Directed Worm Algorithm. OAP and BD were supported by the Polish National Science Centre (NCN) grant DEC-2016/23/B/ST3/01152. MŁ was supported by the Polish National Science Centre (NCN) project 2016/23/D/ST2/00721. Numerical computations in this work were supported in part by PL-Grid Infrastructure.