# Mesoscopic Effects in Quantum Phases of Ultracold Quantum Gases in Optical Lattices

###### Abstract

We present a wide array of quantum measures on numerical solutions of 1D Bose- and Fermi-Hubbard Hamiltonians for finite-size systems with open boundary conditions. Finite size effects are highly relevant to ultracold quantum gases in optical lattices, where an external trap creates smaller effective regions in the form of the celebrated “wedding cake” structure and the local density approximation is often not applicable. Specifically, for the Bose-Hubbard Hamiltonian we calculate number, quantum depletion, local von-Neumann entropy, generalized entanglement or Q-measure, fidelity, and fidelity susceptibility; for the Fermi-Hubbard Hamiltonian we also calculate the pairing correlations, magnetization, charge-density correlations, and antiferromagnetic structure factor. Our numerical method is imaginary time propagation via time-evolving block decimation. As part of our study we provide a careful comparison of canonical vs. grand canonical ensembles and Gutzwiller vs. entangled simulations. The most striking effect of finite size occurs for bosons: we observe a strong blurring of the tips of the Mott lobes accompanied by higher depletion, and show how the location of the first Mott lobe tip approaches the thermodynamic value as a function of system size.

## I Introduction

Quantum phase transitions (QPTs) are often treated in the thermodynamic limit, where the concepts of a critical point, critical exponents, and scaling relations are very clear sachdev1999 (). However, many physical many-body systems for practical applications these days are finite. For example, nano-devices can be smaller than 100 nm. Then a crystal with a lattice constant of a fraction of a nanometer has just 2 or 3 orders of magnitude to work with in scale. In nuclear physics one has QPTs in finite nuclei which are very far from being in the thermodynamic limit caprio2008 (). Experiments on trapped ions, recently used to simulate many body Hamiltonians, normally work with small systems kimK2009 (). Ultracold atoms in optical lattices provide a fourth example, where the number of lattice sites is typically 50 to 100 in each direction, so that in one dimension (1D) they are highly mesoscopic. These systems have been proposed as quantum simulators to solve hard quantum many body problems, including the positive- Fermi-Hubbard Hamiltonian for general filling hofstetter2002 (); jaksch2004 (); lewensteinM2007 (); bloch2008 (). Progress towards this goal requires understanding of how finite-size effects change the quantum phase diagram. In this paper, we explore such mesoscopic or finite size effects for both Bose- and Fermi-Hubbard Hamiltonians.

In addition to finite-size effects, ultracold atomic quantum simulators are non-uniform, as ordinarily there is a weak overall harmonic trap keeping atoms from spilling off the edge of the lattice. This non-uniformity gives rise to a “wedding-cake” or tiered structure of Mott plateaus separated by superfluid layers for bosons batrouniGG2002 (); gemelke2009 (), and similar layered structure for fermions schneiderU2008 (); jordens2008 (). The non-uniformity creates even smaller regions of Mott insulator or superfluid, often just 10 or 20 sites in length. For these small regions the quantum phase diagram is substantially different than in the thermodynamic limit: for example, the quantum depletion, (defined precisely below in Sec. II) shows that the Mott lobe is not closed but instead goes through a very broad crossover into the superfluid regime, as we show in Sec. III.1. Figure 1 shows the sharp phase boundary of a Bose Hubbard system in the thermodynamic limit (white circles kuhner2000 ()) and the continuously variable quantum depletion (defined in Sec. II) obtained in our calculations for 51 sites. As can be seen in the figure, the tip of the Mott lobe is significantly blurred even for this relatively large lattice.

Approximating the effects of a slowly varying trap by using a local chemical potential in the local density approximation, , is not very useful for cold atoms in 1D. The breakdown of the local density approximation has been explored with quantum Monte Carlo and exact diagonalization elsewhere (see rigol2009 () and references therein); our work, which isolates finite size effects in small uniform regions, is complementary to previous studies which have explored the harmonic trap structure as a whole batrouniGG2002 (); rigol2004 (); rigol2009 (). The success of quantum simulators relies in part on understanding all possible sources of error trotzky2009 (), in the same sense as one might do so for a high precision measurement of time-reversal symmetry, for example. Therefore going beyond the local density approximation rigol2009 () is key to the eventual success of quantum simulators, especially in lower dimensions. Our use of TEBD offers a wide array of quantum measures, many of which are not directly accessible to quantum Monte Carlo kuhner1998 (); sengupta2005 (); fubini2006 () or exact 1D methods Lieb1968 (); gu2004 (), including various entanglement entropies and fidelity. Multiple measures provide maximal information, useful for a thorough analysis; in particular, fidelity and fidelity susceptibility make no assumptions about states or symmetries, but simply ask about many body wavefunction overlap in a model-independent way venuti2008 ().

It has already been shown that the mean-field Gutzwiller approximation, which assumes a product state between sites of the lattice, is not very useful in 1D for Hubbard Hamiltonians lewensteinM2007 (). But how far does one have to move away from the Gutzwiller approximation in order to get correct results? A high level of effort has been put into alternative methods, such as generalized dynamical mean field theory in two dimensions for Bose-Fermi mixtures titvinidze2009 (). The convergence parameter in TEBD, which can be couched in terms of the degree of spatial entanglement, allows one to go beyond the Gutzwiller approximation in a methodical manner, and we present a thorough study of this issue. A second question concerning numerical methods is the use of the canonical ensemble vs. the grand canonical ensemble. The two only formally give the same result in the thermodynamic limit. The canonical ensemble is numerically preferred for TEBD and related methods daley2004 () because the number of Fock states is much smaller and therefore calculations are more efficient. We explicitly explore the effects of ensemble dependence Hubbard Hamiltonian phase diagrams.

This article is outlined as follows. In Sec. II we present the Bose- and Fermi-Hubbard Hamiltonians, the many quantum measures we use, and a brief description of and motivation for our implementation of TEBD. In Sec. III we present and discuss our numerical results for mesoscopic effects in ground states of the Bose- and Fermi-Hubbard Hamiltonians, respectively. In Sec. IV we compare canonical to grand canonical ensembles and Gutzwiller to entangled TEBD simulations; these comparisons also serve as a convergence study. Finally, in Sec. V, we conclude.

## Ii Models, Measures, and Methods

We treat two fundamental Hamiltonians for ultracold quantum gases in optical lattices, the Bose- and Fermi-Hubbard Hamiltonians. In both cases we take the tight-binding lowest-band approximation. The Bose-Hubbard Hamiltonian is

(1) | |||||

where is the hopping strength, is the interaction strength, is the chemical potential, , are bosonic destruction/creation operators satisfying bosonic commutation relations on-site and commuting between sites, is the bosonic number operator, and denotes a sum over nearest neighbors. The corresponding Fermi-Hubbard Hamiltonian in particle-hole symmetric form is

(2) | |||||

where are fermionic destruction and creation operators satisfying fermionic anticommutation relations; and have the same meaning as in the bosonic case; and we have assumed -wave interactions and spin 1/2, so that . We note that there is an exact mapping under the canonical transformation from the negative- to positive- Fermi-Hubbard Hamiltonian at half filling (). This will appear visually in our phase diagrams in Sec. III.2.

To find the ground states of these models, we use TEBD tebdOpenSource () in 1D and propagate in imaginary time, starting with a state with uniform weights. TEBD has one main convergence parameter, the number of eigenvalues retained in the reduced density matrix. In all cases we have checked convergence by at least doubling , and seeing if we obtain the same results in the plots; all plots have at least a of 16, and we have spot-checked a of 32. We emphasize that as we have a very high resolution of () pixels in our diagrams, we cannot go to the preferred typical values of in the hundreds, in common use for TEBD calculations; also, grand canonical calculations are inherently more restrictive, due to the much larger state space. We use both canonical and grand canonical ensembles depending on circumstances. The latter is useful to understand the proper interpolation method to be used for the canonical ensemble for smaller numbers of sites. In the canonical ensemble, chemical potential is derived by finite differencing, ; for just a few sites, this would lead to only a few points of resolution on the chemical potential axis without interpolation. That is, a finite difference approximation for the chemical potential is only directly useful for large ; interpolation can only be justified by comparing to the grand canonical results. We illustrate this point in Sec. IV.1. For fermions, there is an additional advantage of grand canonical calculations in that, in practice, we find that imaginary time propagation is less likely to get stuck in a local minimum.

For fermions, the local, or on-site dimension of the Hilbert space is fixed at since we consider only a single band. For bosons provides a second convergence parameter. In practice we truncate at , meaning 0 to 5 bosons per site, and we have checked for specific points in the phase diagrams and found no visible changes. We implement imaginary time propagation in the second order Suzuki-Trotter scheme to obtain the ground state. Finally, we parallelize over the parameter space via MPI for bosons, and for fermions. We typically perform our calculations on between and cores. We note that TEBD becomes exact, aside from the error of the Suzuki-Trotter decomposition, for .

Our quantum measures can be divided into four classes: moments, correlation, entropies, and fidelity. Although these classes of measures can have significant overlap in the information they provide, nevertheless there are always cases where they are distinct. For example, the spin singlet is maximally entangled but is not correlated. Thus it is useful to consider all classes.

First, for moments, the energy, , and filling factor, , are defined as

(3) | |||||

(4) |

where is the number of lattice sites. To calculate and for fermions, one must also sum over the spin index .

Second, for correlations in bosons we calculated the quantum depletion, , which is a measure of how far away the system is from the Penrose-Onsager definition of superfluidity penrose1956 ():

(5) |

where are the eigenvalues of the single particle density matrix ordered so that . For correlations in fermions we use three standard measures scalettar1989 (); varney2009 (). The s-wave pairing correlation, , is given by

(6) | |||||

(7) |

The charge-density wave correlation, , is given by

(8) | |||||

(9) |

The antiferromagnetic structure factor is given by

(10) | |||||

(11) |

Third, for entropies, the average local von Neumann entropy kitaev2006 (), or average local entropy of entanglement, , is

(12) |

where the local reduced density matrix, , is defined as

(13) |

and is a pure-state density matrix for the whole system, as we work at zero temperature. Note that the base of the log normalizes the entropy to lie between zero and unity. We also studied the Q-measure meyer2002 (); brennenGK2003 () or average local impurity, , a special case of generalized entanglement barnumH2003 (); barnumH2004 (), which gives similar results to the von Neumann entropy (see Sec. III):

(14) |

Fourth, we calculated both fidelity, , and fidelity susceptibility buonsante2007 (); venuti2008 (); rigol2009b (), . The fidelity is an overlap measure on the states, and is given by

(15) |

where is a control parameter in a quantum phase diagram, in our case either chemical potential or hopping. The fidelity susceptibility is a second derivative of the fidelity in the control parameter, and we normalize it to the unit length:

(16) |

This proves to be a more useful quantity than the fidelity, as it is independent of the control parameter step size , and a divergence of in the thermodynamic limit signals a quantum phase transition. Note that bears no relation to the convergence parameter for TEBD, . In the figures below we take to be for the Bose-Hubbard Hamiltonian and for the Fermi-Hubbard Hamiltonian.

## Iii Quantum Measures on Ground States of Hubbard Hamiltonians

### iii.1 Bose-Hubbard Hamiltonian

The Bose-Hubbard Hamiltonian phase diagram has been calculated in the mean field approximation in higher dimensions fisher1989 () and with both the DMRG method monien1998 (); kuhner2000 () and Monte Carlo in one dimension. Trapped Bose systems with an external harmonic rigol2009 () and linear potential carr2009j () have been treated with Monte Carlo and TEBD, respectively. However, these studies did not focus on finite size effects without aiming at the thermodynamic limit.

In Fig. 2 are illustrated a suite of quantum measures for and 40 sites, from left to right in column form. The rows represent the different quantum measures as labeled on the right hand side: from top to bottom these are density, entropy, Q-measure, depletion, fidelity, and fidelity susceptibility. The calculation of these measures and their definitions are given in Sec. II. We emphasize that quantum Monte Carlo cannot be used to obtain entanglement. Only matrix-product-state-based techniques can do so for general local Hamiltonians, including TEBD and DMRG schollwock2005 (), which lead to an equivalent result, although with different numerical implementations. Bethe-ansatz-based techniques have been used to calculate entanglement in certain limiting cases (near ) in the Fermi-Hubbard Hamiltonian gu2004 (), but not generally. The case of is exact diagonalization for and , but at the result is already approximate, since states would be required for exact diagonalization in matrix product state form (see Sec. II). Each panel of Fig. 2 is of resolution in . The largest system of took several weeks to complete in the grand canonical ensemble on nodes with separate memory, each node consisting of 8 cores with shared memory, all on the Golden Energy Computing Organization geco () high performance computing cluster “Ra”. For this reason we do not go to for the whole phase diagram, but only make spot checks, since simulations scale as . We calculated both canonical and grand canonical results, and show the canonical ones in Fig. 2.

Turning to the physical interpretation of Fig. 2, we begin with number. In the canonical ensemble a definite total number is guaranteed. This can be seen first in the average filling, which proceeds in multiples of , as illustrated in the top row of the figure, starting from the vacuum and a filling of zero. The vacuum region is evident in the lower portion of the panels in the top row of where it is shaded a deep blue. The density increases in discrete jumps as the chemical potential is increased; interpolation is used between curves of definite number difference in order to fill in the picture, as described in Sec. II and Sec. IV.1.

In the second row we show the quantum depletion. One clearly sees the highly depleted Mott lobes and the lowly depleted outer superfluid region. Note that there are only two Mott lobes displayed in these figures; the apparent third, lowest such lobe is in fact the vacuum, , for which quantum depletion is ill-defined. In keeping with the 1D nature of the problem we consider, we observe that even the superfluid region is still depleted at a level of about 30 % far away from the Mott lobes. Again, the smooth nature of the Mott-insulator-superfluid transition near the tip is in strong contrast to the sharper transition in other regions.

In the third and fourth rows of Fig. 2 we show entropies. For completeness we compare the von Neumann entropy or entropy of entanglement, which is now a standard measure used in quantum information nielsenMA2000 () and also proposed for the study of quantum phase transitions osterlohA2002 (); kitaev2006 (); dengS2008 (), to the Q-measure, which has been more closely connected to local observables barnumH2003 (); sommaR2004 (). We find both measures useful and nearly but not exactly the same; we choose the Q-measure for the rest of our study. Note that the local von Neumann entropy of the Mott insulating regions is low by this measure, which reflects the nearly pure number state on each site in the Mott phase. On the other hand, the superfluid region displays high entropy, because the superfluid state is delocalized across sites, leading to mixed states on each site when measured separately. This contrasts with the usual association of low thermal entropy with a superfluid phase.

In the second through fourth rows one can clearly pick out the Mott-insulating region. It has a sharp phase boundary for small . Kühner, White and Monien kuhner2000 () first provided an accurate evaluation of the critical value, , for the Kosterlitz-Thouless (KT) transition at the tip of the Mott lobe in the infinite system limit; in our notation this corresponds to a value of . Near the entropy changes more continuously for 3, 5, and 10 sites, while for 20 sites the lobe begins to close. Finally at 40 sites, on the far right of the figure, the region of continuous change is very narrow. We clarify that by “sharper boundary” as compared to “more continuous” we are speaking in relative terms, specifically of the difference in entropies between states with commensurate and incommensurate filling. There is no truly sharp transition since this is a finite system; all “quantum phase transitions” seen here are in fact formally crossovers between regions with different characteristics.

Finally, in the fifth and sixth rows we show the fidelity and fidelity susceptibility . In our present implementation we compare the neighboring pixels of the figures in for curves of fixed and , as described in Sec. II. This excludes a comparison between changes in total number; such a comparison, which we also performed in the grand canonical ensemble, is not shown, because features in and are dominated by changes in the conserved quantity. In the thermodynamic limit the superfluid region of the phase diagram forms a continuous tensor product space and the fidelity is everywhere zero even for arbitrarily small . This is the well-known Anderson orthogonality catastrophe anderson1967 (). Thus, the fidelity, through its scaling and qualitative behavior, reflects for a finite system the critical behavior of an infinite system. Indeed, the Mott lobes appear clearly in this set of measures, as can be seen in the figure.

In order to provide a more quantitative estimate of the tip location for finite systems than a by-eye estimation from Fig. 2, we study quantum measures and derivatives thereof for unit filling in the canonical ensemble. Three figures illustrate our results. In Fig. 3 are shown derivatives of quantum measures which could be used to define a crossover; in rows from top to bottom these are derivatives of average local von Neumann entropy, Q-measure, and quantum depletion; and lastly, the fidelity susceptibility. The fidelity susceptibility provides a particularly sharp signature of the Mott phase, a result that was noted previously in the context of perodic boundary conditions and the completely connected graph buonsante2007 (). Within the Mott lobe, assumes values close to those given by second order perturbation theory rigol2009b () in the limit that :

(17) |

for open boundary conditions. For periodic boundary conditions the system size cancels and the expression is . These perturbations represent particle-hole fluctuations, which are weighted more heavily in higher Mott lobes due to Bose stimulation. This causes to have a higher value in the upper Mott lobe in Fig. 2 as compared to the lower lobe.

As the number of sites increases from 3 to 40 the extrema move towards the thermodynamic limit . Figure 4 summarizes the most useful extrema of these curves in an easy to read single plot showing that quantum depletion and fidelity susceptibility appear to converge most rapidly to their asymptotic values in the thermodynamic limit. However, is affected by the boundary conditions, as can be seen in Eq. 17 as well as in Fig. 5, where we compare periodic and open boundary conditions using both exact diagonalization and TEBD; this was also a significant issue in DMRG and perturbation theory studies, where earlier calculations monien1998 () obtained a less accurate estimate of the critical point due to the choice of boundary conditions kuhner2000 (). This figure also serves as a check on our code. We emphasize that open boundary conditions are closer to the conditions of ultracold atom experiments at the present time.

### iii.2 Fermi-Hubbard Hamiltonian

The Fermi-Hubbard Hamiltonian in 1D has an exact solution Lieb1968 (); Lieb2003 () via the Bethe ansatz; in principle, one can extract any information from such a solution, including all measures of interest, and the Bethe ansatz has been applied to finite as well as infinite domains. Because of their computational complexity, Bethe ansatz results are not widely used for performing partial traces of the type that we consider here gu2004 (). Therefore it is useful to consider finite size effects for a suite of measures simultaneously, including measures not previously accessed through exact solution techniques. TEBD or DMRG methods also enable one to “turn on” non-mean-field effects, i.e. entanglement, in discrete steps in the form of the parameter , as we show in Sec. IV.2.

Figure 6 displays ground state solutions of the Fermi-Hubbard Hamiltonian [Eq. (2)] in the grand canonical ensemble for 3, 5, 9, 10, and 20 sites from left to right in columns, following the same pattern as Fig. 2 but with measures appropriate to fermions. Since we take , for 3 and 5 sites TEBD is equivalent to exact diagonalization, as ; in fact, we explicitly checked our results with a separate exact diagonalization routine. The first row shows the density. For positive we observe that all number states are allowed; as in Sec. III.1, the grand canonical ensemble settles on a total-number-conserving state for these small systems due to the discreteness of the spectrum. For example, for there are seven colored regions for increasing and greater than 0, of which all regions are visible up to ; these correspond to 0,1,,6 fermions. For there are eleven colored regions, of which 2 are barely visible due to our choice of scale, and for higher we observe the same pattern. For odd numbers of sites the exact mapping from positive to negative in Eq. (2) at zero chemical potential does not hold, and therefore there is no symmetry around the vertical axis at for , while for we observe this symmetry. For negative , particles enter only in pairs, excepting at for even, as can be seen from the figure.

In the second row of Fig. 6 we show the entanglement entropy. One clearly observes the particle-hole symmetry between negative and positive chemical potential. For odd numbers of sites, the most highly entangled regions for positive are those in which there is one extra particle or hole above half filling; for even systems it is two extra particles or holes, as can be seen in the stripe-like pattern for , for example. States nearer the band insulator, with just a few particles (holes) compared to the vacuum (maximal filling) are less entangled. Entanglement also decreases in the Mott region for increasing since the system moves closer to a product state. There is a discontinuity in the -measure between regions of half filling and slightly away from half filling which vanishes at =0. This can be understood as the requirement that the derivatives approaching from above and below are required to satisfy by particle-hole symmetry. For negative the most highly entangled regions are near zero chemical potential and near half-filling. Therefore, in general the pattern for both positive and negative is that the entanglement is largest near half filling. We show only the Q-measure, having established in Sec. III.1 that the average local von Neumann entropy does not show substantial additional information.

In the third row of Fig. 6 we show the charge-density correlation (CDC). This measure shows how close the system is to a “checkerboard” occupying every other site. For an odd number of sites, as visible in , and negative , the CDC is higher for two pairs than for one pair, by the nature of the measure [see Eq. (9)] – this is due to the fact that there is one extra pair which cannot be placed. For even numbers of sites, as seen in it is particle-hole symmetric, as expected, yielding the same result upon flipping all plots around the horizontal axis. CDCs do not follow the positive- to negative- mapping of the Fermi-Hubbard Hamiltonian. CDCs decrease as is increased, since the system is closer to a Mott product state. In the fourth row of Fig. 6 we show the pairing correlations, which are of interest only in the negative- region where pairing naturally occurs due to the negative interaction.

In the fifth row of Fig. 6, we show the antiferromagnetic structure factor. This measure is the spatial fourier transform of the spin-spin correlation function evaluated at , and is motivated by the fact that a canonical transformation maps the half-filled positive- Hubbard model in the strong coupling limit onto a spin-1/2 antiferromagnetic Heisenberg model with in an external magnetic field

(18) |

This measure jumps sharply at half filling where the antiferromagnetic order is expected to exist, grows positively and reaches a maximum for some value of , and then decreases to 0 for large . The symmetry about is preserved because the operators involved are which vanish for on-site singlet pairs.

In the sixth row of Fig. 6 we show the fidelity susceptibility. Due to our use of the grand canonical ensemble, only different total number states have an overlap significantly different from unity for small systems. For larger systems, fidelity susceptibility becomes less dominated by differences in number states as the number of jumps in number approaches the resolution in , as was also the case for the Bose-Hubbard Hamiltonian in Fig. 2.

## Iv Key Numerical Comparisons

### iv.1 Grand Canonical vs. Canonical Ensembles

The canonical ensemble is much preferable for large systems because it is more efficient. However, in order to obtain a chemical potential one must approximate the derivative in . Our volume is automatically fixed by performing the derivative for a given system size . The derivative is then evaluated by finite differences, either a simple forward difference of , or a higher order finite difference. At any level of finite difference approximation, for smaller systems there is only a very low resolution for the chemical potential axis. This leads to a lack of information about quantum measures. This issue is highlighted in Fig. 7 by showing what the actual data is for the canonical ensemble (right panel) as compared to the grand canonical ensemble (left panel) for 20 sites. In the canonical calculation, each horizontal curve represents an allowed value of . The grand canonical ensemble allows us to motivate the interpolation of the canonical ensemble result to go past the finite-difference result and obtain a continuous plot. Our interpolation algorithm is to take the value from the nearest curve according to the standard distance measure . Figure 8 shows the actual canonical ensemble data used to produce Fig. 2.

The initial state chosen for our imaginary time propagation has the same projection onto every Fock state in the Hilbert space allowed by our numerical truncation. In our grand canonical calculations, this includes states which do not have the same total number, . However, for small systems, number states are sufficiently different in energy that imaginary time propagation converges on a state with a definite . On the border between number states one requires a very high to deal with the near degeneracy of states with different . This leads to the regular pattern of bright spots seen in the left panel of Fig. 7. The sole effect of higher is to remove bright spots from the entropy and Q-measure. The density of such low-entropy spikes increases with increasing because a higher is needed for a larger system – more sites means the possibility for more spatial entanglement. This is another reason why canonical ensemble calculations are preferred for bosons. We have verified that all panels in Fig. 2 match the grand canonical picture, with the exception of the bright spots, which do not occur in canonical ensemble calculations.

Although we do not illustrate it in Fig. 7, the variance in the total number, , where , is zero to 5 or 6 digits of precision. This is in fact a measure of convergence for our grand canonical simulations; we require all Schmidt values (eigenvalues of the reduced density matrix) be converged to 10 digits in imaginary time propagation, and obtain this total number variance, having started with an initial state with weights on all total numbers allowed by the cut-off , i.e., from 0 atoms to atoms.

### iv.2 Gutzwiller Mean Field vs. Entangled Time-Evolving Block Decimation

Since mean field methods are usually the most tractable approaches to solution of many body problems, identifying their domain of validity is important. For instance, for the three-dimensional Bose-Hubbard Hamiltonian one can obtain the phase diagram via the Hubbard-Stratonovich transformation fisher1989 (), which is a mean field method, or alternately, via the Gutzwiller mean field approximation, as discussed in Sec. I. On the other hand, the Gutzwiller method fails for the one-dimensional Bose-Hubbard Hamiltonian lewensteinM2007 (). It is natural to ask, how far beyond mean field does one need to go in order to obtain correct results for essential quantum measures?

To address this question, we repeat the calculations of Sec. III.2 for the Fermi-Hubbard Hamiltonian for 10 sites only, starting with the mean field calculation of and then doubling successively. Figure 9 shows the results: each column corresponds to while the rows correspond to the same quantum measures as Fig. 6. Note that can be found in the third column of Fig. 6, so we do not repeat it here. We make several remarks concerning this figure. The density is a first moment and only requires lower values of , as can be seen in the first row. For in the first column, the entanglement shown in the second row is zero to 16 digits of accuracy, reflecting our use of double precision floating point arithmetic in our implementation of TEBD. The entanglement does not converge until , as can be determined from a comparison of Fig. 9 and Fig. 6.

In the third row, for large positive we obtain the correct result for charge density correlations at half-filling even for very low because the Mott phase is nearly a product state. The pairing phase for negative , on the other hand, requires higher , as can be seen in the fourth row. This can be understood from the fact that superfluids and superconductors are highly entangled in position space; similarly, in the Bose Hubbard Hamiltonian one must go to high in order to resolve the superfluid region far away from the Mott lobes. By “spatial entanglement” we refer to exactly the Schmidt number used in TEBD. The fifth row shows anti-ferromagnetic order.

Finally, in the sixth row we observe that the main features of the fidelity susceptibility appear already for , because fidelity measures in the grand canonical ensemble are mainly sensitive to total number states in our system; correct total number states are obtained at very low . We observe that particle-hole symmetry is maintained even at the mean field level of , while the positive- to negative- mapping at half-filling and zero chemical potential is only good for at least .

Although Fig. 9 shows that is sufficient to obtain by-eye converged values of all quantum measures considered in Fig. 6, a more careful quantitative study is desirable. In Fig. 10 we show this study for both fermions and bosons for 10 sites. Both the average error (left panels) and maximum error (right panels) over all pixels in for fermions and for bosons are evaluated. The error for any given pixel is defined as

(19) |

where is the quantum measure of interest evaluated at Schmidt number and is an integer, in our case taken to be 1, 2, 3, 4, 5. That is, we evaluate measure at of 16, 8, 4, 2 and compare it to the evaluation of the same measure at the previous smaller value of . The maximum and average errors are given by

(20) | |||||

(21) |

For fermions we calculate these two errors for density, Q-measure, and pairing in the grand canonical ensemble for 10 sites, corresponding to Fig. 9. For brevity, we do not show an analogous figure to Fig. 9 for bosons, but calculate maximum and average error for sites and the grand canonical ensemble shown in Fig. 7 in the third column. As can be seen from the figures, the average error is about 1 % for for bosons, and the maximum error is about 10 %. For fermions it somewhat larger for the same . The maximum error appears large, but upon a detailed investigation we find that this error is localized in the vacuum or band insulator regions; this means that the max error is dominated by imaginary time propagation effects. Other regions have an error corresponding to the average error, including the important half-filling Mott region for positive . In general, local measures like density or entanglement converge faster than non-local measures like quantum depletion, as can be seen in the figure. Also, the bosonic measures converge more quickly than do the fermionic measures.

## V Conclusions

We explored finite size effects in the Bose- and Fermi-Hubbard Hamiltonians with TEBD and imaginary time propagation, studying a range of quantum measures including moments, correlations, entanglement, and fidelity. We found that entanglement nielsenMA2000 (); weiTC2005 () followed the same pattern as more traditional measures for quantum phase transitions, such as depletion for the Bose-Hubbard Hamiltonian and the superfluid-Mott-insulator transition, and pairing, charge-density, and antiferromagnetic structure factor for the Fermi-Hubbard Hamiltonian in the Mott transition; the exception was the tip of the Mott lobe, which was best described by the depletion, in terms of converging to the thermodynamic result more quickly as a function of system size. Fidelity susceptibility provides a useful alternative measure for characterizing the phase transition, and is a non-local measure.

We compared results computed within the grand canonical ensemble and canonical ensembles: the chemical potential can be approximated in the canonical ensemble by a finite difference, a method only useful for large . However, by first performing a grand canonical calculation and then using it to motivate an interpolation of the canonical ensemble we could subsequently use the more efficient canonical calculations to obtain more highly converged results. The canonical ensemble is also more useful for fidelity susceptibility calculations, which are otherwise dominated by changes in the conserved quantity of total number.

We performed quantitative convergence studies and showed that particle-hole symmetry is maintained in the particle-hole symmetric form of the Fermi-Hubbard Hamiltonian for all values of , but the positive- to negative- mapping at half-filling and zero chemical potential is correct only for at least , meaning it is not satisfied in the mean field approximation. We found that for bosons the tip of the Mott lobe is wide-open for smaller system sizes of less than about 20 sites, so that exploring quantum phase transitions by changing chemical potential finds sharp effects while changing hopping find a very broad crossover.

We acknowledge useful discussions with Ippei Danshita, Ludwig Mathey, Ryan Mishmash, Alejandro Muramatsu, and Marcos Rigol. This work was supported by the National Science Foundation under Grant PHY-0547845 as part of the NSF CAREER program (LDC, MLW) and by the Aspen Center for Physics (LDC), by the Summer Undergraduate Research Fellow (SURF) program at the National Institute of Standards and Technology (DGS, RCB), and by the NSF under Physics Frontiers Center grant PHY-0822671 (RCB, CWC).

## References

- (1) S. Sachdev, Quantum Phase Transitions (Cambridge University Press, New York, 1999).
- (2) M. A. Caprio, P. Cejnar, and F. Iachello, Ann. Phys. 323, 1106 (2008).
- (3) K. Kim, M.-S. Chang, R. Islam, S. Korenblit, L.-M. Duan, and C. Monroe, arXiv:0905.0225 (2009).
- (4) W. Hofstetter et al., Phys. Rev. Lett. 89, 220407 (2002).
- (5) D. Jaksch and P. Zoller, Annals of Physics (2004).
- (6) A. Lewenstein, M. Sanpera et al., Adv. Phys. 56, 243 (2007).
- (7) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
- (8) T. D. Kühner, S. R. White, and H. Monien, Phys. Rev. B 61, 12474 (2000).
- (9) G. G. Batrouni et al., Phys. Rev. Lett. 89, 117203 (2002).
- (10) N. Gemelke, X. Zhang, C.-L. Hung, and C. Chin, arXiv:0904.1532 (2009).
- (11) U. Schneider et al., Science 322, 1520 (2008).
- (12) R. Jördens et al., Nature (London) 455, 204 (2008).
- (13) M. Rigol, G. G. Batrouni, V. G. Rousseau, and R. T. Scalettar, Phys. Rev. A 79, 053605 (2009).
- (14) M. Rigol and A. Muramatsu, Phys. Rev. A 69, 053612 (2004).
- (15) S. Trotzky et al., arXiv:0905.4882 (2009).
- (16) T. D. Kühner and H. Monien, Phys. Rev. B 58, R14741 (1998).
- (17) K. Sengupta and N. Dupuis, Phys. Rev. A 71, 033629 (2005).
- (18) A. Fubini et al., e-print cond-mat/0605602 (2006).
- (19) E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).
- (20) S.-J. Gu, S.-S. Deng, Y.-Q. Li, and H.-Q. Lin, Phys. Rev. Lett. 93, 086402 (2004).
- (21) L. C. Venuti et al., Phys. Rev. B 78, 115410 (2008).
- (22) I. Titvinidze, M. Snoek, and W. Hofstetter, Phys. Rev. B 79, 144506 (2009).
- (23) A. J. Daley, C. Kollath, U. Schollwock, and G. Vidal, J. Stat. Mech. - Theor. Exp. 2004, P04005 (2004).
- (24) http://physics.mines.edu/downloads/software/tebd/.
- (25) O. Penrose and L. Onsager, Phys. Rev. 104, 576 (1956).
- (26) R. T. Scalettar et al., Phys. Rev. Lett. 62, 1407 (1989).
- (27) C. N. Varney et al., arXiv:0903.2519 (2009).
- (28) A. Kitaev and J. Preskill, Phys. Rev. Lett. 96, 110404 (2006).
- (29) D. A. Meyer and N. R. Wallach, Journal of Mathematical Physics 43, 4273 (2002).
- (30) G. K. Brennen, Quant. Inf. Comp. 3, 619 (2003).
- (31) H. Barnum, E. Knill, G. Ortiz, and L. Viola, Phys. Rev. A 68, 032308 (2003).
- (32) H. Barnum et al., Phys. Rev. Lett. 92, 107902 (2004).
- (33) P. Buonsante and A. Vezzani, Phys. Rev. Lett. 98, 110601 (2007).
- (34) M. Rigol, B. S. Shastry, and S. Haas, arXiv:0906.5347 (2009).
- (35) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
- (36) H. Monien, M. Linn, and N. Elstner, Phys. Rev. A 58, R3395 (1998).
- (37) L. D. Carr and M. K. Oberthaler, Phys. Rev. Lett., under review (2009).
- (38) U. Schollwock, Rev. Mod. Phys. 77, 259 (2005).
- (39) http://geco.mines.edu/.
- (40) M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, UK, 2000).
- (41) A. Osterloh, L. Amico, G. Falci, and R. Fazio, Nature 416, 608 (2002).
- (42) G. O. S. Deng and L. Viola, e-print http://arxiv.org/abs/0802.3941 (2008).
- (43) R. Somma et al., Phys. Rev. A 70, 042311 (2004).
- (44) P. W. Anderson, Phys. Rev. Lett. 18, 1049 (1967).
- (45) E. H. Lieb and F. Y. Wu, Physica A 321, .
- (46) T.-C. Wei et al., Phys. Rev. A 71, 060305 (2005).