Ground-state phase diagram of the half-filled bilayer Hubbard model

Employing a combination of functional renormalization group calculations and projective determinantal quantum Monte Carlo simulations, we examine the Hubbard model on the square lattice bilayer at half filling. From this combined analysis, we obtain a comprehensive account of the ground-state phase diagram with respect to the extent of the system's metallic and (antiferromagnetically ordered) Mott-insulating as well as band-insulating regions. By means of an unbiased functional renormalization group approach, we exhibit the antiferromagnetic Mott-insulating state as the relevant instability of the free metallic state, induced by any weak finite on-site repulsion. Upon performing a careful analysis of the quantum Monte Carlo data, we resolve the difficulty of identifying this antiferromagnetic ground state for finite interlayer hopping in the weak-coupling regime, where nonmonotonic finite-size corrections are shown to relate to the two-sheeted Fermi surface structure of the metallic phase. On the other hand, quantum Monte Carlo simulations are well suited to identify the transition between the Mott-insulating phase and the band insulator in the intermediate-to-strong coupling regime. Here, we compare our numerical findings to indications for the transition region obtained from the functional renormalization group procedure.

Figure 1

Phase diagram of the Hubbard model on the square lattice bilayer (inset) at half filling. Red dots display the phase boundary obtained from projective DQMC (the red solid line is a guide to the eye). The blue dashed line indicates the reduction of antiferromagnetic instabilities upon increasing $t_{\perp }$ at intermediate values of $U/t$ in the fRG approach (cf. Sec. 4b for details).

Figure 2

Fermi surfaces of the square lattice bilayer tight-binding model ($U=0$) at half filling for different values of $t_{\perp }/t=0,1,2,3$ (from left to right). Every point in one band is connected to the other band by the nesting vector $\mathit{Q}=(\pi ,\pi )$. The grid points denote the momenta included in a finite system of linear length $L=6$ with periodic boundary conditions. For the system sizes accessible in the DQMC simulations, these represent the only values of $t_{\perp }/t$ with discrete momenta located on the Fermi surface (larger dots).

Figure 3

The $U=0$ bilayer density of states (DOS) $\rho (\epsilon ,t_{\perp })$ for $t_{\perp }=2t$ (solid red line) in comparison with the single-layer DOS (dashed black line). The dotted vertical line marks the position of the Fermi level.

Figure 4

Left panel on top: Interaction vertex and spin convention. Below: Loop contributions from the particle-particle channel (a), the crossed particle-hole channel (b), and the direct particle-hole channel (c). Right panel: Patching scheme of the Brillouin zone for a total of $N_p=32$ patches. Depending on the energy band the coupling function is evaluated on a wave vector at the Fermi level indicated by the dots (red: upper band, blue: lower band).

Figure 5

Effective interaction vertex near the critical scale for the AF-SDW in units of $t$, exhibiting a sharply pronounced momentum structure. The axes are numbered according to the number of the patch, cf. Fig. 4, where wave vectors $k_1$ are depicted vertically and $k_2$ are depicted horizontally. We fix $k_3$ to be on patch 1. Left panel: Effective vertex for ${\lambda }_1={\lambda }_2={\lambda }_3={\lambda }_4$. Middle panel: ${\lambda }_1={\lambda }_3,{\lambda }_2={\lambda }_4$. Right panel: ${\lambda }_1={\lambda }_4,{\lambda }_2={\lambda }_3$.

Figure 6

Left panel: fRG critical scale ${\Lambda }_c$ as a function of $t/U$ for the single-layer case (solid line) and the bilayer case with $t_{\perp }=1t,2t,3t$ (from top to bottom) using $N_p=48$ patches, exhibiting exponential decrease of ${\Lambda }_c\sim e^{-\mathrm{const}.\times t/U}$. This functional dependence suggests an instability for any value of $U>0$ in accord with HFMFT. Right panel: The fRG critical scale ${\Lambda }_c$ decreases by several orders of magnitude when $t_{\perp }/t$ is increased from 0 to 4. We show this behavior for three choices of $U/t=3,U/t=2$, and $U/t=1$ by the solid, dashed, and dotted lines, respectively.

Figure 7

Static spin-spin correlations along a triangular path (inset) for different coupling strengths $U/t$. The calculations were performed for a bilayer system with interlayer hopping $t_{\perp }=2t$ and linear length $L=24$.

Figure 8

Finite-size scaling of the staggered magnetization for $t_{\perp }=t$ (top) and $t_{\perp }=2t$ (bottom). For couplings $U/t\ge 4$, linear fits were performed (solid lines).

Figure 9

Comparison of the antiferromagnetic structure factor of the weakly interacting system and the noninteracting static susceptibility ${\chi }_0^{+-}\left(\mathit{Q}\right)$ at finite temperature $T={10}^{-4}t$ (blue dots) for different linear lengths $L$. A quadratic function may be fitted to the susceptibility data for lengths $L$ that are multiples of 12 (blue solid line).

Figure 10

Comparison of the single-particle gaps from HFMFT (solid and dashed lines), QMC (red dots and black squares), and the critical scale from fRG (blue triangles). Left: Single-layer case. Right: Bilayer system with $t_{\perp }=2t$.

Figure 11

Finite-size scaling analysis of the rescaled antiferromagnetic structure factor around the AF-BI transition, assuming $\beta =0.369$ and $\nu =0.711$. The crossing points indicate the critical coupling $U_c/t$ of the transition.

Figure 12

Data collapse analysis of the rescaled antiferromagnetic structure factor around the AF-BI transition, assuming $\beta =0.369$ and $\nu =0.711$.

Figure 13

Apparent crossings of the rescaled structure factor data for $t_{\perp }=0$ (left panel) and $t_{\perp }=2t$ (right panel) in the lower-$U$ region, indicative of the crossover in the nature of the AF state.

Figure 14

Left panel: Effective low-energy interaction vertex in units of $t$ for the orbital combination ${\lambda }_1={\lambda }_2={\lambda }_3={\lambda }_4$ for three choices of $t_{\perp }\in \{2t,3.5t,3.9t\}$ (from top to bottom) at fixed on-site interaction $U=6t$. The conventions are identical to the ones in Fig. 5. Right panel: Interaction vertex profile for fixed $k_1$ at patch 1 for the $t_{\perp }$ from the left panel. This corresponds to a horizontal cut of the interaction vertices shown in the left panel. The flat line in the bottom (green) shows the value of the initial interaction.

Figure 15

Flow of the width $C_{\Lambda }$ for three choices of $t_{\perp }\in \{3t,3.5t,3.9t\}$ (solid, dashed, and dotted line, respectively) at fixed on-site interaction $U=6t$.