Quantum fluctuations beyond the Gutzwiller approximation in the Bose-Hubbard model

We develop a quantum many-body theory of the Bose-Hubbard model based on the canonical quantization of the action derived from a Gutzwiller mean-field ansatz. Capitalizing on the ability of this latter to describe both the Mott insulator and the superfluid phases, our theory is a systematic generalization of the Bogoliubov theory of weakly interacting gases and provides accurate results across the whole phase diagram. We characterize the superfluid-insulator phase transition studying the two-point correlation functions, the local number fluctuations, and the superfluid stiffness across the whole phase diagram. Two different universality classes are recovered at integer and non-integer filling and the density fluctuations are successfully compared to accurate quantum Monte Carlo data. To conclude we highlight the potential of our theory in view of including interactions between collective modes to describe their finite lifetime and their quantum optical properties.

The Hubbard model is one of the most celebrated models of quantum condensed matter theory. The main reason is probably the widespread belief that its twodimensional fermionic version holds the key to understand how high-temperature superconductivity emerges upon doping a Mott insulator [1,2]. Its central feature is the archetypal competition between the kinetic energy term, which favors delocalized states, and the local Coulomb repulsion which favors localization [3,4].
In the two-dimensional fermionic model this competition is somewhat hidden by the presence of other phases bridging between the Mott insulator, the superconducting and the metallic states, including the celebrated pseudogap [5] and charge-ordered [6] phases. While a quantitative description of these many-body phases remains challenging even to modern numerical approaches, the competition between the kinetic and interaction energies directly manifests itself as a quantum phase transition between a superfluid and a Mott insulator in the bosonic version of the model, the so-called Bose-Hubbard (BH) model [7]. As a consequence of its paradigmatic nature, this transition has attracted an enormous experimental interest in the last years using cold atoms trapped in optical lattices [8][9][10][11][12][13][14][15] and, more recently, using interacting photons in the novel context of non-equilibrium phase transitions [16][17][18].
On the theoretical side one of the easiest approaches to the BH model is based on the Gutzwiller ansatz [19][20][21], where the variational wave function is site-factorized and the (complex) weight c n (r) is a variational parameter. While some important features both in the superfluid and in the insulating phase are qualitatively captured by this wave function, sophisticated techniques are typically necessary to obtain quantitative predictions for the observables, in particular those involving quantum correlations between distant sites.
In the weakly-interacting superfluid regime, the Bogoliubov approach -based on the quantization of the fluctuations around the mean-field Gross-Pitaevskii equation (GPE) of pure Bose-Einstein condensates [15,22] -provides an accurate description of the properties of the gas, including non-trivial quantum correlations [23,24] and interactions between quasi-particles [25,26].
In this paper we combine the successful features of the Gutzwiller and Bogoliubov approaches, introducing a description of the BH model in terms of quantized fluctuations around the Gutzwiller solution. The collective phonon excitations of the Bogoliubov theory are replaced here by the more complex excitations of the BH model, in particular the Goldstone and Higgs modes on the superfluid side and particle/hole excitations on the insulating one [27][28][29][30]. In spite of the local nature of the Gutzwiller ansatz (1), the quantum description of the excitations and, in particular, of their zero-point fluctuations turns out to correctly reproduce the non-local many-body correlations of the different phases as well as the different critical behaviours of the commensurate and incommensurate phase transition [7,27].
As compared to analogous Gutzwiller approaches recently developed for Fermi systems [31], the advantage of our formalism is that it directly includes quantum fluctuations of the collective modes and can naturally incorporate effects beyond linearized fluctuations. This is essential to successfully tackle problems such as the finite quasi-particles lifetime via their nonlinear interactions and the collective quantum correlations between the products of their decay.
We consider the BH model in a three dimensional (d = arXiv:1908.03470v1 [cond-mat.quant-gas] 9 Aug 2019 3) lattice composed by N sites: where J is the hopping amplitude, U the on-site interaction, and µ the chemical potential, while r, s labels all pairs of nearest-neighboring sites.
A variational approach based on the Gutzwiller ansatz (1) leads to the following Lagrangian functional for the Bose-Hubbard model [21]: where the dot indicates temporal derivative, ψ (r) = ψ(r) = n √ n c * n−1 (r) c n (r) is the mean-field order parameter and H n = (U/2) n (n − 1) − µ n. This reformulation of the Bose-Hubbard many-body problem allows to adopt the complex-valued Gutzwiller parameters c n (r) as the coordinates of the problem. Noting that c * n (r) = ∂L/∂ċ n (r) are the canonical momenta corresponding to c n (r), canonical quantization of the theory [32,33] requires promoting the coordinates c n (r) to operatorsĉ n (r) and imposing equal-time canonical commutations ĉ m (r),ĉ † n (s) = δ r,s δ m,n . In analogy with the Bogoliubov approximation for dilute Bose-Einstein condensates [15,22], we use the Gutzwiller ground state c 0 n obtained by minimizing the energy Ψ G |Ĥ|Ψ G as our starting point and writê c n (r) =Â r c 0 n + δĉ n (r). By inserting the expansion of c n (r) into Eq. (3) and keeping only terms up to the quadratic order in the fluctuations we obtain a pseudo-Hermitian quadratic form, i.e. Hermitian according to the η = diag[Id, −Id] product, as in standard Bogoliubov theory [22]. One can restrict to local fluctuations orthogonal to the ground state n c 0 n * δĉ n (r) = 0 and then [34] for the operator ensuring the proper normalization of the Gutzwiller ansatz.
Thanks to homogeneity, fluctuations can be Fourier transformed as δĉ n (r) ≡ N −1/2 k∈BZ e ik·r δĈ n (k), leading to the quadratic Hamiltonian: where the different components δĈ n (k) are gathered into the vector δĈ(k), E 0 is the mean-field ground state energy, and L k is the matrix of the pseudo-Hermitian quadratic form describing the dynamics of small fluctuations of wave vector k around the Gutzwiller mean-field solution [21,30].
Even though this paper is focused on the quadratic expansion (4), inclusion of the higher-order terms in δĈ n arising from the hopping term in (3) does not pose fundamental conceptual difficulties [35] and will be the subject of forthcoming publication.
A suitable rotation of the Gutzwiller operators can then be used to recast the quadratic form (4) into a diagonal formĤ corresponding to the different many-body excitation modes with frequencies ω α,k , labeled by their momentum k and the branch index α. Bosonic commutation rules [b α,k ,b † α ,k ] = δ k,k δ α,α between the annihilation b α,k and creationb † α,k operators are enforced by choosing the the usual Bogoliubov normalization condition As usual [22], the sum over the α modes is restricted to positive-norm modes [36]. In the MI phase the two lowest branches are the gapped particle and hole excitations. As the SF phase is approached, one of them transforms into the gapless Goldstone mode of the broken U(1) symmetry in the SF, while the other (gapped) one, often referred to as the Higgs mode, is related to the fluctuations of the amplitude of the order parameter.
In the same way as the quantum Bogoliubov theory goes beyond a linearization of the mean-field GPE [22], the operator-valuedb α,k are the main novelty of our quantum Gutzwiller theory as compared to previous mean-field Gutzwiller approaches based on C-numbervalued excitation amplitudes [21,30].
In the following we show the power of our theory by investigating relevant observables. Notably, we will display accurate results in regimes where non-local quantum correlations (which were not included in the initial meanfield Gutzwiller ansatz) play a crucial role.
We write the Bose field operatorψ(r) by replacing in the expression for the order parameter used in (3) the corresponding operatorsĉ n (r). To the lowest order in the fluctuations δĉ n (r) one has: where ψ 0 is the order parameter in the ground state and the particle (hole) amplitudes U α,k (V α,k ) involve a combination of the u α,k,n (v α,k,n ) [30]. At the quadratic order in the δĉ's, the normalized zero-temperature one- body coherence function reads and typical curves in different regions of the phase diagram [ Figure 1(a)] are plotted in Figure 1(b-d).
In the deep SF phase [panel (b)], the spectral weight is saturated by the Goldstone mode alone so that our prediction for g (1) (r) reduces to the Bogoliubov one. This is no longer true in the intermediate region 2 dJ/U ∼ O(1) of a strongly interacting superfluid, where the Gutzwiller approach includes also the contribution of other excitation modes [29].
In the MI phase [panels (c,d)], in spite of the locality of the Gutzwiller ansatz (1), our approach is able to capture the finite exponentially decreasing coherence, g (1) (r) ∼ exp (−r/ξ). Such result is in stark contrast with the r = 0 on-site coherence obtained within the mean-field Gutzwiller theory.
These same panels illustrate how the critical behaviour at the superfluid-insulator transition depends on whether we cross it at integer or non-integer filling. If we cross the transition outside the tip of the Mott lobe, the correlation length ξ of the MI grows but remains bounded [panel (c)]. Once on the SF side, long-range order suddenly appears. On the other hand, ξ diverges if we approach the SF phase through the tip of the Mott lobes and a powerlaw g (1) is found at the critical point [panel (d)]. This remarkably different behaviour is related to the universality class of the transition [27,37]. The critical point of the integer-filling transition at the tip of the Mott lobe is a O(2) Lorentz invariant point where both the particle and hole modes become gapless (before turning into the Goldstone and Higgs modes on the SF side), which explains the divergent coherence length. In any other critical point on the edge of the lobe, only one excitation becomes gapless, namely the hole (particle) excitation below (above) the tip. Since the non-trivial short-distance coherence of the MI is due to virtual particle-hole excitations, the exponential decay of g (1) (r) is dominated by the gap of the particle (hole) excitation, which remains finite. Beyond the critical point, the transition to the superfluid phase occurs via the proliferation of a (initially dilute) gas of holes (particles) on top of the MI, which automatically condense at zero temperature and form a hole (particle) superfluid [30].
Superfluid density -The proper definition [38,39] of the superfluid density fraction F s = n s /n involves the response of the system to a transverse phase twist [40] and the current-current response function are calculated by expanding the current operatorĴ x (r) = i J[ψ † (r + e x )ψ(r) − h.c.] in powers of δĉ's. The only finite contribution to Λ xx J is the 4 th -order one which reveals the crucial role played by the coupling between different collective modes in suppressing the superfluid stiffness and creating a sort of normal component [41]. The result of the quantum Gutzwiller method is illustrated as a red thick line in Figure 2: the superfluid fraction is unity in the deep SF, approaches zero at the critical point and then remains zero in the MI regime. Interestingly, the superfluid fraction remains always larger than the condensed fraction (green line) defined as the r → ∞ limit of (8). Further insight is obtained by isolating the two contributions in (9): the normal current correlation Λ defined in (11) (black thick line) exhibits a non-monotonic behaviour as a function of J/U , since it tends to zero at the extrema of the phase diagram and reaches its maximum at the MI-SF transition. In the strongly-interacting SF regime, the Goldstone-Higgs vertex almost saturates the sum (11) and leads to a complete suppression of F s . In the MI phase, the vanishing F s results from the perfect cancellation of the short-range correlations in the virtual particle/hole by the residual kinetic energy.
For completeness we also report the Bogoliubov prediction (dashed blue line in Figure 2), to which our results approach in the weakly-interacting limit 2 dJ/U 1 and, as noted in [42], leads to a vanishing normal current response and therefore to a large superfluid density even in the Mott phase.
On the SF side, the antibunching g (2) (0) < 1 due to the repulsive on-site interactions grows stronger for decreasing 2dJ/U and well matches the Bogoliubov prediction (gray dot-dashed line) in the deep SF. On the MI side, while the mean-field Gutzwiller theory (dashed blue line) predicts a perfect antibunching g (2) (0) = 0, our quantum theory is able to account for the virtual excitation of doublon-hole pairs leading to a g (2) ∝ J 2 behaviour at low 2dJ/U , in excellent agreement with perturbative calculations (green short-dashed line) [20]. Even more remarkably close and across the critical point the theory is in very good agreement with low-temperature quantum Monte Carlo predictions [43] (thin black dotted line), where no other semi-analytical theory is available. The role of the quantum fluctuations is further illustrated by the off-site density correlations function shown in the inset: while mean-field theory predicts a trivial g (2) (r = 0) = 1, negative non-trivial correlations due to virtual particle-hole pairs (on the MI side) and to quantum density fluctuations on the SF side are found, these latter in agreement with Bogoliubov theory in the weakly interacting regime.
In this paper we have introduced a simple and powerful analytical tool to study the many-body physics of quantum interacting particles on a lattice based on a canonical quantization of the fluctuations around a Gutzwiller ansatz. The power of the method has been validated on the archetypal case of the Bose-Hubbard model. In spite of the locality of the initial mean-field Gutzwiller ansatz, our quantized treatment correctly includes very non-local features such as the superfluid stiffness of the superfluid phase, the spatial structure of the virtual particle-hole pair excitations on top of a Mott insulator, and the different behaviour of the correlation functions at various critical points.
Moreover, our formalism can be straightforwardly extended to treat inhomogeneous configurations, more complex forms of hopping and interaction terms, as well as to finite temperature and/or time-dependent problems. Going beyond the quadratic expansion around the meanfield, our method can in fact naturally incorporate additional nonlinear terms describing interactions between quasi-particles [35] so to describe, e.g., their temporal decay into entangled pairs via multi-branch Beliaev decay processes [44]. An exciting long term perspective will be to extend the formalism to those driven-dissipative models that can now be realized in photonic systems [16][17][18].