Phase diagram of the square lattice bilayer Hubbard model: a variational Monte Carlo study

We investigate the phase diagram of the square lattice bilayer Hubbard model at half filling with the variational Monte Carlo method for both the magnetic and the paramagnetic case as a function of interlayer hopping t_perp and on-site Coulomb repulsion U. With this study we resolve some discrepancies in previous calculations based on the dynamical mean field theory, and we are able to determine the nature of the phase transitions between metal, Mott insulator and band insulator. In the magnetic case we find only two phases: An antiferromagnetic Mott insulator at small t_perp for any value of U and a band insulator at large t_perp. At large U values we approach the Heisenberg limit. The paramagnetic phase diagram shows at small t_perp a metal to Mott insulator transition at moderate U values and a Mott to band insulator transition at larger U values. We also observe a reentrant Mott insulator to metal transition and metal to band insulator transition for increasing t_perp in the range of 5.5t<U<7.5t. Finally, we discuss the obtained phase diagrams in relation to previous studies based on different many-body approaches.


Introduction
Understanding the origin of transitions from a metal to a Mott or a band insulator in correlated systems has been a topic of intensive debate in the past few years. Various generalizations of the Hubbard model have been investigated for this purpose, like the extended Hubbard model [1], the ionic Hubbard model in one and two dimensions [2][3][4][5][6], the two-band Hubbard model [7] and the bilayer Hubbard model [8][9][10][11][12][13][14][15]. The latter has been especially important in relation to the bilayer high-T c cuprates [16,17]. Previous investigations of the above models have been carried out primarily by employing dynamical mean-field theory (DMFT) [18] and its cluster extensions [19][20][21][22][23]. While DMFT already captures a significant amount of the key properties in correlated systems, it is extremely important to analyse these models with completely unrelated many-body methods in order to get a deeper understanding of the underlying physics. In this work we investigate the phase diagram of the square lattice bilayer Hubbard model at halffilling with the variational Monte Carlo (VMC) method. With this study (i) we are able to resolve some discrepancies between previous DMFT and cluster DMFT studies and (ii) we find new aspects of the Mott to band transition not captured in previous studies.
The bilayer Hubbard model on the square lattice is given by the following Hamiltonian: , , is the electron density. t is the nearest neighbour hopping parameter in the plane, ⊥ t denotes the hopping between planes and U is the on-site Coulomb repulsion; see figure 1. Note that in the following we measure all quantities in units of t. The U = 0 dispersion, x y k for the bonding/antibonding (±) bands is at half-filling perfectly nested, i.e. ε ε = + ± ∓ k Q k with π π = Q ( , ), as illustrated in figure 1. As a consequence, the ground state [24] is an ordered antiferromagnet for any > U 0, as long as a Fermi surface exists, which is the case for < ⊥ t t / 4 in the limit → U 0. In this work we analyse the magnetic phase diagram at T = 0 as well as the paramagnetic case, which we investigate by suppressing the long-range magnetic ordering. This latter investigation, while it is done at T = 0, is relevant for predictions at finite and small temperatures where long-range magnetic order is absent as a consequence of the Mermin-Wagner theorem [25].
In the limit of large interaction strength, → ∞ U , the Hubbard Hamiltonian in equation (1) reduces to the Heisenberg Hamiltonian,

∑∑
. The bilayer Heisenberg Hamiltonian has itself been a subject of extensive research, and it has been found to undergo an order-disorder transition [26] at = ⊥ J J / 2.552, which corresponds to = ⊥ t t / 1.588 in terms of the Hubbard Hamiltonianʼs hopping parameters [27,28].
Studies of the complete phase diagram of the bilayer Hubbard model have been conducted by Fuhrmann et al using DMFT [8] and by Kancharla and Okamoto with cluster DMFT [9]. The authors of reference [8] concentrated on the paramagnetic phase at finite temperature and found a metallic phase at small U and small ⊥ t values as well as an insulating phase as ⊥ t increases (the band insulator) or U increases (the Mott insulator), but no clear separation between the Mott and band insulating phases was found.
The authors of reference [9] considered clusters of sizes × 2 2 (i.e. two sites per plane) within cluster DMFT and performed exact diagonalization calculations to solve the impurity problem, which allowed them to investigate T = 0 and to obtain a magnetic phase diagram as well as a T = 0 phase diagram where magnetism is suppressed. These results show several important differences with respect to the earlier phase diagram given by Fuhrmann et al [8]. First of all, Kancharla et al [9] distinguish between a Mott and a band insulator by looking at the behaviour of the charge gap as a function of ⊥ t . Also, they find a Mott insulator at = ⊥ t 0 for any value of > U 0, even if magnetism is suppressed. They mention the cluster DMFTʼs intralayer spatial correlation as the reason for this being found with cluster DMFT but not in the DMFT results of Fuhrmann et al [8]. However, at = ⊥ t 0 the planes are completely decoupled and should have the same properties as a single plane. For the single plane, different methods predict that the system becomes a Mott insulator at some finite U, as shown for instance in reference [29], where the critical U was calculated using the dynamical cluster approximation with a variety of cluster sizes. Also, Tocchio et al [30] and Capello et al [31] showed with the VMC method (which in the first reference includes Fermi surface renormalization effects) the appearance of a Mott insulator at a finite U. A possible reason for the occurrence of a Mott insulator for any finite U in reference [9] might be the small cluster size of × 2 2 used in cluster DMFT. Having two sites in each plane breaks the fourfold rotational symmetry of the square lattice and results in an artificially enhanced local pair within each plane for any > U 0. This can notably affect the phase diagram of a system, as reported by Lee et al for the two-orbital Hubbard model [32]. A particularly interesting feature of the phase diagram in reference [9] is that for a certain range of U values the system goes through two phase transitions as ⊥ t is increased, first from a Mott insulator to a metal and then from a metal to a band insulator, while at large U the system exhibits a direct transition from a Mott to a band insulator. Analogous features have also been proposed in reference [3] for the ionic Hubbard model. If magnetic order is allowed, it is noteworthy that no magnetic ordering is found in the cluster DMFT [9] for small U and ≲ < ⊥ t 2 4, even though there is a perfect nesting between the Fermi surfaces of the bonding and antibonding bands with π π = Q ( , ). A metallic phase at small U in the magnetic phase diagram has also been obtained in a determinant quantum Monte Carlo (DQMC) study [10]. This result could be a consequence of the finite temperatures that have been used in DQMC studies, since they could be large enough to destroy the tiny magnetic order at small U, where the gap is exponentially small.
Overall, DMFT and previous DQMC results do not give a conclusive picture yet, making it worthwhile to also employ other methods in order to gain a deeper understanding of the physics in the bilayer Hubbard model.

Methods
The variational Monte Carlo method was introduced by McMillan in 1965 to calculate the ground state of liquid He 4 [33], and in 1977 applied to a fermionic system for the first time [34]. Its basic idea is to use the Rayleigh-Ritz principle [35] to approximate the ground state through a variational many-body wavefunction. It is a Monte Carlo method because a stochastic sampling is used to evaluate the sum over a high dimensional configuration space. A detailed description of how the variational Monte Carlo method can be applied to the Hubbard model may be found for instance in reference [36]. The VMC approach has also played a central role when examining the limit of large U of the Hubbard model [37][38][39], in the context of high temperature superconductivity.
The choice of the variational many-body wavefunction is crucial in order to obtain reliable results. Here, we define a variational state Ψ | 〉 that consists of two parts: a Slater determinant Φ | 〉 and a Jastrow correlator J acting on Φ | 〉: Here the Slater determinant Φ | 〉 ensures the antisymmetry of the wavefunction while the Jastrow factor J modifies its amplitude to take into account electronic correlations.
The state Φ | 〉 is the ground state of a variational mean-field HamiltonianĤ var which may include up to five different terms: nearest neighbour hopping within the planes, hopping between the planes, superconducting pairing in the planes with d-wave symmetry, pairing between the planes and an antiferromagnetic term, according to the following [40][41][42]: , , Here r labels the sites within the planes, l is the plane index, e x y ( ) is the unity vector along the x y ( ) direction and S i z indicates the z component of the spin operator on site i. Note that the square lattice bilayer model is a bipartite lattice, with τ ∈ i ( ) {1, 2} labelling the sublattice of site i, so a different spin orientation is preferred for each of the two sublattices when μ > 0 m . We would like to mention that the following particle-hole transformation has been used in order to diagonalize the variational Hamiltonian: This is possible because we chose the spins to align along the z direction in the antiferromagnetic term of equation (13). The alternative choice of aligning the spins along the x direction-see for instance [43]does not allow one to study magnetism and pairing together as a single Slater determinant, but it is often preferred because the variational state is improved by the application of a spin Jastrow that couples spins in a direction orthogonal to the ordering one; see [44].
The Jastrow factor J implements a long-range density-density correlation which has been shown to be essential in the variational description of Mott insulators [45]. In order to account for the bilayer nature of the system we used a modified version of the Jastrow factor with a different set of variational parameters for intraplane (v) and interplane ( ⊥ v ) correlations:  with the v(ij) and the ⊥ v ij ( ) being optimized independently for every distance | − | r r i j . The Metropolis algorithm [46] with single-particle updates has been used to generate the electronic configurations, while the optimization of the variational wavefunction was done using the stochastic reconfiguration method [47,48], which allows us to independently optimize every variational parameter in J , as well as ⊥ t (var) , Δ, Δ ⊥ and μ m in the mean-field state. The inplane hopping parameter t is kept fixed to set the energy scale. A lattice size of 10 × 10 sites per plane was used, unless stated otherwise.

Results
The non-magnetic and the magnetic phase diagrams obtained using VMC simulations are presented in figure 2.
The non-magnetic phase diagram (figure 2, left panel) shows a metallic phase at small U and ⊥ t , which goes, as expected, into a band insulating phase at = ⊥ t 4. The critical value of ⊥ t that is needed to make the system band insulating decreases as U is increased. At small ⊥ t the system undergoes a metal to Mott insulator transition as U is increased, with the critical U ranging from about 5.5 at = ⊥ t 0 to 7.5 at = ⊥ t 0.7. For large U the system undergoes a Mott to band insulator transition at ≈ ⊥ t 0.7, where the critical value of ⊥ t is independent of U and only weakly dependent on the systemʼs size. Hysteresis in the variational parameters when going through the Mott to band insulator transition suggests the transition to be of first order. An interesting feature is that for ≲ ≲ U 5.5 7.5 the system first undergoes a Mott insulator to metal transition and then a metal to band insulator transition as ⊥ t is increased. In contrast, only two phases are found in the magnetic phase diagram (figure 2, right panel). The system is a Néel ordered antiferromagnetic Mott insulator as long as ⊥ t is smaller than some critical value, and a paramagnetic band insulator for larger ⊥ t . The critical interlayer hopping is = ⊥ t 4 at U = 0 and decreases as U is increased. At U = 10 it reaches a value of about = ⊥ t 1.9, which suggests that the critical interlayer hopping approaches the Heisenberg limit of = ⊥ t 1.588 for → ∞ U . We now discuss the details of the calculations undertaken in order to obtain the above presented phase diagrams.

The non-magnetic phase diagram
In order to obtain the non-magnetic phase diagram it is useful to first analytically diagonalize the variational single-particle Hamiltonian whose eigenstates are used to build the determinantal part of the wavefunction. The argument here is that if the Slater determinant is already gapped, no correlator can make the wavefunction conducting again. Therefore, one can identify a band insulator purely by looking for a gap in the band structure of the determinantal part. Note also that a superconductor has a gap in the mean-field state, without being an insulator. However, the only parameter in our variational wavefunction that could induce superconductivity is the inplane pairing Δ that is always gapless in the nodal direction, since it has d-wave symmetry. Therefore, we do not risk accidentally classifying a superconductor as an insulator, and the existence of a gap in the mean-field state is indeed equivalent to the system being insulating.
In the non-magnetic case the variational Hamiltonian consists of four terms for the intraplane/interplane hopping and pairing:

t t var (var)
Diagonalizing it analytically using the particle-hole transformation gives the following band structure: x y x y k (2,4) (var) 2 2 There are two bands with only positive/negative energies. At half-filling there are enough electrons to populate exactly two bands, so the two negative bands are always completely filled and the two positive bands entirely empty. Consequently, the only way not to have a band gap is by having the bands touch each other at zero energy. The easiest way to understand the influence of the different parameters is to look at figure 3, where the bands are plotted for different values of the variational parameters. It can easily be seen that there are two ways of opening a gap in the mean-field part of our variational wavefunction, namely having a > ⊥ t 4 (var) , or a non-zero Δ ⊥ which corresponds to the formation of singlets between the planes. We now discuss the variational Monte Carlo simulation results used to draw the non-magnetic phase diagram in figure 2. Figure 4 shows the optimized interplane pairing Δ ⊥ as a function of the interplane hopping ⊥ t . Starting with = ⊥ t 4 at U = 0 the region with a non-zero Δ ⊥ extends to lower ⊥ t as U is increased. For ⩾ U 8 the jump is at a constant ≈ ⊥ t 0.7. As any non-zero Δ ⊥ opens a gap in the mean-field part of the wavefunction, the region of a non-zero Δ ⊥ maps out the band insulating part of the phase diagram in figure 2.
At variance with the band insulator case, a Mott insulating region is characterized by a gapless mean-field state, while the insulating nature is driven by the electronic correlations that are included in the Jastrow factor J . In order to discriminate between a Mott insulator and a metal, we use the following single-mode ansatz for the wavefunction of the excited state, which goes back to Richard Feynmanʼs work on excitations in liquid helium [49], and was later successfully applied to fermionic systems [50,51]: wheren q is the Fourier transform of the particle density with = q q q ( , , 0) x y . By calculating the energy of the excited state, one can derive a formula for an upper bound of the charge gap E g , that relates it to the static structure factor = 〈ˆˆ〉 − N n n q ( ) q q [52]: g q 0 2 Figure 5 shows the charge gap E g as a function of the Coulomb repulsion U for different interplane hoppings ⊥ t . As expected from the known monolayer results, we find the system to be a Mott insulator for large enough values of U. Note that the critical U needed to make the system Mott insulating, i.e. when the charge gap E g starts to grow as a function of U, increases from about U = 5.5 at = ⊥ t 0 to U = 7.5 at = ⊥ t 0.6. Note that a finite charge gap for ≲ U 5.5 and = ⊥ t 0 is just an artefact of the limited number of q points that are available for the extrapolation to = q 0 and indeed it decreases as the system size increases. We point out that a sizable in-plane Δ can be found only in the Mott insulating region, indicating that the pairing within the planes is to be understood in terms of the resonating valence bond theory [53,54], in which d-wave pairs are formed, but not phase coherent. Indeed, it is the presence of the Jastrow factor of equation (16), that allows the Φ | 〉  J wavefunction to describe a Mott insulator through the opening of a charge gap without any symmetry breaking. Figure 6 shows the renormalization of the variational interplane hopping ⊥ t (var) . In general the variational ⊥ t (var) is not too different from the ⊥ t in the original Hubbard Hamiltonian, but there are two exceptions to this rule. The first one is that in the Mott insulating phase (at large U and small ⊥ t ) the variational ⊥ t (var) is renormalized to quite small values. Together with the lack of a Δ ⊥ in this region, this indicates that the mean-field part of the  (21)) as a function of the Coulomb repulsion U.
wavefunction is almost that of two decoupled Hubbard planes. The other exception is that for small U and large ⊥ t the variational ⊥ t (var) is renormalized to values larger than the original ⊥ t . This is due to the fact that in the limit → U 0 the ground-state wavefunction is independent of the value of ⊥ t with the bonding band filled as long as > ⊥ t 4. The band gap and the Mott gap opening via Δ ⊥ , presented in figures 4 and 5 respectively, suffice for mapping out the regions of band and Mott insulators in the phase diagram. The physics of the respective transitions can be confirmed by looking at the density of double occupancies presented in figure 7, as a function of the interplane hopping ⊥ t . For intermediate values of U, the density of double occupancies first rises sharply and then drops off abruptly again as the ⊥ t is increased, signalling that there is a metallic phase in between the Mott and band insulating phases for ≲ ≲ U 5.5 7.5. The double-occupation density was also calculated within 2× 2 cluster DMFT [9], but, contrary to our results of figure 7, it was found that it decreases in the metallic phase. We believe that our results correspond to the physics of the Mott insulating phase, which suppresses double occupancies.
At large ≳ U 8 the Mott insulator goes directly into the band insulator at ≈ ⊥ t 0.7. The Mott insulating wavefunction is characterized by an in-plane Δ > 0 and a Δ = ⊥ 0, while the band insulating wavefunction has a sizable Δ ⊥ but no in-plane Δ. A strong hysteresis in these variational parameters was observed when going through the transition, so both a Mott and a band insulator could be obtained for ⊥ t close to the transition. The optimal wavefunction is the one with the lowest energy, as plotted in figure 8. Finding hysteresis means that both the Mott and the band insulating wavefunctions are local energy minima in our variational parameter space, thus suggesting that the Mott to band insulator transition is of first order.

The magnetic phase diagram
The simulations through which the magnetic phase diagram was obtained were performed in the same way as those for the non-magnetic case, except that the site and spin dependent chemical potential μ m of equation (13) was no longer fixed to zero during the optimization. Figure 9 shows the magnetic potential μ m as a function of the interplane hopping ⊥ t . For small U, an antiferromagnetic order is found for <  antiferromagnetic region to smaller ⊥ t , indicating that the critical interplane hopping goes to the Heisenberg limit of = ⊥ t 1.588 for → ∞ U . One can show, by analogy to the non-magnetic case in the previous section, that any nonzero μ m makes the system insulating, by analytically diagonalizing the variational single-particle Hamiltonian:ˆ=ˆ+ˆ+⊥ The two negative bands are filled and the system can only be conducting if these bands touch the empty bands at zero energy. Looking at the equations for the energy bands, one can see that this cannot happen for μ ≠ 0 m , and hence one can use a non-zero μ m as a criterion for an insulating state. Note that we classify the ordered state in the phase diagram (figure 2, right panel) as a Mott insulator, even though we have a gap in the mean-field state. This is due to the fact that the antiferromagnetic ordering is correlation induced.
Comparing the interplane Δ ⊥ in figure 10 with the plot of the magnetic potential μ m in figure 9, we see that a non-zero Δ ⊥ is found in the entire paramagnetic region. Both the μ m and the Δ ⊥ open a gap in the mean-field part of the wavefunction, so the square lattice bilayer Hubbard model is always an insulator at > U 0 if magnetic order is allowed. It is interesting to note that at large ≳ U 8 and < ≲ ⊥ t 1 2 there is a small region with both a non-zero magnetic potential μ m and an interplane pairing Δ ⊥ . This, together with the fact that no hysteresis in the variational parameters was observed at the order-disorder transition, suggests that the transition is indeed continuous.

Conclusion
In summary we have calculated the magnetic and non-magnetic phase diagram of the square lattice bilayer Hubbard model using the variational Monte Carlo method, as summarized in figure 2. Moreover, our results suggest that the Mott insulator to band insulator transition is of first order in the non-magnetic phase diagram, while it becomes continuous when magnetic order is allowed. Comparison of our results to the ones obtained with DMFT [8] and cluster DMFT [9] reveals that our non-magnetic phase diagram includes some features observed in DMFT and 2 × 2 cluster DMFT but also new distinct properties, as follows. While, in agreement with 2 × 2 cluster DMFT, there is a region in which the system first goes from a Mott insulator to a metal and then to a band insulator as ⊥ t is increased, we do not find that this region extends down to U = 0. Instead, for the decoupled planes we find a metal to Mott insulator transition at a critical ≈ U 5.5, which agrees with the DMFT results of Fuhrmann et al, and the dynamical cluster approximation results for the single layer of Gull et al [29]. For large U our results agree with those of 2 × 2 cluster DMFT in that there is a direct transition from a Mott to a band insulator, but our critical ⊥ t is smaller by about a factor of 3. The reason for this might be the cluster with two sites in each plane used by the authors of reference [9] that breaks the fourfold rotational symmetry of the square lattice and creates an artificially enhanced local pair within each plane, ultimately stabilizing the in-plane Mott phase against the interplane dimers of the band insulating phase.
The magnetic phase diagram, however, shows clearly a behaviour different from the one predicted by 2 × 2 cluster DMFT calculations. The most obvious difference is that we no longer find a metallic phase if magnetic ordering is allowed. Instead we find a Néel ordered Mott insulator, which we attribute to the perfect nesting between the bonding and antibonding bandʼs Fermi surfaces. The reason for the variational Monte Carlo approach stabilizing magnetic ordering as compared to cluster DMFT with two sites per plane might be the much more explicit treatment of long-range correlations.
The variational Monte Carlo results improve our understanding of the bilayer Hubbard model, but further investigation may be necessary to clarify the origin of the sizable differences between VMC and DMFT results.