Kekulé valence bond order in the Hubbard model on the honeycomb lattice with possible lattice distortions for graphene

Graphene is often described by the Hubbard model on a honeycomb lattice. As a pioneering work on this model, Sorella and Tosatti elucidated a quantum phase transition from a semimetal (SM) to an antiferromagnetic Mott insulator (AFMI), which occurs at a finite strength of the interaction [1]. They further anticipated that their findings might be relevant for “the physics of strong correlations in the $\pi$-electron system in 2D graphite” [1]. Later, after the synthesis of graphene [2], not only the peculiar noninteracting band structure [3, 4], but also the many-body effects and their consequent quantum phase transitions in the Dirac electrons have been intensively studied [5, 6, 7, 8, 9]. Some of these studies were stimulated first by the possibility of a spin liquid phase [10, 11, 12, 13, 14, 15, 16] and later by an intriguing connection between the graphene physics and the celebrated Gross-Neveu model in high-energy physics [17, 18, 19, 20, 21, 22, 23].

While the effects of interaction on graphene have been to some extent understood at least based on lattice models, an experimental realization of AFMI in graphene, which would be highly promising for future device applications [24], has yet to be established. However, this does not necessarily prove that graphene is simply weakly correlated. Among many available estimations of the model parameters for graphene [25, 26, 27, 28, 29, 30, 31], adopting the partially screened on-site Coulomb interaction of $U_{00}=9.3\mathrm{eV}$ [29] and the widely accepted value of the hopping integral of $t\approx 2.7\mathrm{eV}$ [25, 26], we notice that their ratio is not far below the critical point of the Hubbard model on the honeycomb lattice $U_{\mathrm{c}}/t\simeq 3.8$ [13, 22, 23, 32, 33, 34]. This leads us to expect that AFMI is realized by applying strain, since strain would reduce $t$ with less effect on the Hubbard interaction $U$. It is noted that the application of strain has been explored as a way to introduce a gap in graphene not only by the Mott mechanism [35, 36, 37, 38, 39, 40].

Recently, an ab initio quantum Monte Carlo (QMC) calculation examined the effect of strain on graphene and found that a Kekulé-like dimerized state, rather than AFMI, is stabilized by increasing strain [41]. Being an insulator, this state should also be promising from an application point of view. Theoretically, while the density functional theory calculations suggest the AFMI state as the ground-state of the strained graphene [39, 42], it is fascinating that the more accurate description of the electron correlation, made possible by the QMC technique, predicts the Peierls-like ground state instead of the Mott insulator [41].

In this paper, as a complementary approach to the first-principles calculations, we investigate the Kekulé valence bond order based on the simple lattice model, that is, the Hubbard model on the honeycomb lattice with electron-lattice (e-l) coupling, using the auxiliary-field quantum Monte Carlo (AFQMC) method, which is numerically exact for the electron correlation. We consider lattice displacements which modulate the electron hoppings and treat them adiabatically. Even with this simplification, there may be a difficulty in determining the optimal bond-ordering pattern, as was the long-standing issue in similar models on a square lattice [43, 44, 45, 46, 47, 48]. However, this difficulty is largely circumvented on the honeycomb lattice. Frank and Lieb have shown that possible lattice distortions for graphene are periodic and reflection-symmetric, thus having at most six atoms per unit cell [49]. Hence, together with the constant-volume condition, it is sufficient to consider two different lattice lengths in this unit cell. We determine the optimal values of these lattice displacements by minimizing the free energy at finite temperatures under a constant volume and thus obtain the phase diagram for several strengths of the e-l coupling, which shows that the Kekulé valence-bond-solid (KVBS) phase is stabilized in the vicinity of $U_{\mathrm{c}}/t$. For a parameter corresponding to a rigid lattice as in the case of graphene, we find that the KVBS phase emerges as an intermediate phase between SM and AFMI, which suggests the correlation-driven SM-to-KVBS phase transition as predicted by the preceding ab initio QMC study [41].

The rest of the paper is organized as follows. In the next section, we introduce the model and explain how to obtain the most stable bond configuration. We then show the results of the QMC simulations, including the finite-temperature phase diagram, along with discussions in Sec. III. Finally, the summary is given in Sec. IV.

We study the e-l coupled Hubbard model described by the following Hamiltonian:

where $\hat{c}_{i\sigma}$ annihilates an electron with spin $\sigma$ ($=\uparrow,\downarrow$) at site $i$, $t$ denotes the transfer integral in the absence of the lattice distortion, $U$ is the Hubbard interaction, and $\hat{n}_{i\sigma}=\hat{c}^{\dagger}_{\sigma}\hat{c}_{\sigma}$ is a number operator. The sum $\langle i,j\rangle$ runs over all nearest neighbor sites of the honeycomb lattice. The lattice displacement denoted by $u_{ij}$ is defined as a normalized relative change in the bond length; the bond length between sites $i$ and $j$ is expressed as $a_{ij}=a(1+u_{ij})$, where $a$ is the lattice constant of the undistorted lattice. Thus, for example, if the bond shrinks, i,e, $u_{ij}<0$, the transfer integral increases by $-gu_{ij}$ at the cost of the elastic energy in the last term of Eq. (1) with $K$ being the elastic constant.

We consider the lattice degrees of freedom in the adiabatic limit, assuming that the lattice displacements $u_{ij}$ are frozen at the energy-minimizing configuration. On the other hand, we take into account both quantum and thermal fluctuations for the electron degrees of freedom. Specifically, $u_{ij}$ is determined so as to satisfy the condition,

where $F=-\beta^{-1}\ln Z$ is the free energy at an inverse temperature $\beta=1/T$, and $Z$ is the grand-canonical partition function, i.e., $Z=\textrm{Tr}\ e^{-\beta\left(\mathcal{\hat{H}}-\mu\hat{N}\right)}$, with $\hat{N}$ being the total number operator of electrons. By setting the chemical potential $\mu=0$, we study the model at half filling. The condition of Eq. (2) yields a self-consistent equation,

Here the bracket $\langle\cdot\rangle$ denotes calculating an expectation value, for which we employ the finite-temperature version of AFQMC method [50, 51, 52, 53]. Since the self-consistent equation is iteratively solved, the QMC simulation needs to be performed at each iterative step. Thus, in total, the required computational resources are rather large but feasible owing to the recent development of computational power. Several years ago, similar calculations were performed using the stochastic series-expansion QMC method [54, 55, 56], which is less computationally demanding than AFQMC, for the study of quasi-one dimensional molecular conductors [57, 58, 59, 60].

In addition to the computational cost, there is another issue to consider when using the QMC method in solving the self-consistent equation. Since the results of the QMC simulation necessarily involve statistical errors, the solution of the self-consistent equation also fluctuates. To improve the quality of the solution within the limited computational resources, we employ the Robbins-Monro algorithm, which is a simple and robust technique for root-finding problems from noisy observables [61, 62, 63]. We typically run about a hundred of the Robbins-Monro steps, in each of which a relatively short QMC simulation, e.g., 1000 Monte Carlo steps for equilibration followed by 1000 steps for measurements is performed on a finite-size cluster spanned by two lattice vectors of $L\bm{\tau_{1}}$ and $L\bm{\tau_{2}}$. Here $\bm{\tau_{1}}$ and $\bm{\tau_{2}}$ are primitive translation vectors (see Fig. 1) of the unit cell for the distorted honeycomb lattice as described below. The number of sites in the cluster is $N_{\mathrm{s}}=6L^{2}$, and we study the system with $L$ up to eight, which is moderate for the iterative calculations.

We take the unit cell which contains six sites and three different bonds as shown in Fig. 1. We emphasize that this is not an assumption of a supercell in a mean-field treatment. Instead, it is what Frank and Lieb have rigorously shown based on the symmetry argument [49]. Even if we take enlarged unit cells or consider all possible spatial configurations without the unit cell, the most stable configuration should converge to the above cell structure, possibly after a very long simulation time. In Fig. 1, the bond lengths are denoted by $a_{i}=a\ (1+u_{i})$ ($i$=1,2,3), and, without losing generality, we label these values as $u_{1}\leq u_{2}\leq u_{3}$. Then, apart from the trivial undistorted lattice, possible patterns of the lattice distortions fall into three categories: (1) maximum distortion of $u_{1}<u_{2}<u_{3}$ [Fig. 1(b)], (2) dimer-type distortion of $u_{1}<u_{2}=u_{3}$ [Fig. 2(a)], and (3) plaquette-type distortion of $u_{1}=u_{2}<u_{3}$ [Fig. 2(b)]. Here the state with the dimer-type distortion is a specific definition of our KVBS state, while the state with the plaquette-type distortion is referred to as the Kekulé resonating-valence-bond (KRVB) state. Furthermore, for simplicity, we study the system at a constant volume by imposing a constraint condition $\sum_{i}u_{i}=0$. This condition guarantees that the lengths of the primitive translation vectors are unchanged, i.e, $|\bm{\tau_{1(2)}}|=\sum_{i}a_{i}=3a$ (see Fig. 1), thus preserving the volume of the system. Subsequently, we can parameterize $u_{i}$ by two parameters, $u$ and $u^{\prime}$, as $u_{1}=-u-u^{\prime}/2$, $u_{2}=u/2-u^{\prime}/2$, and $u_{3}=u/2+u^{\prime}$. Solving the self-consistent equations, we obtain the optimized values of these parameters denoted as $\bar{u}$ and $\bar{u}^{\prime}$. In this way, the KVBS (KRVB) state is characterized by $\bar{u}>0$ and $\bar{u}^{\prime}=0$ ($\bar{u}=0$ and $\bar{u}^{\prime}>0$), and the state with the maximum distortion is considered as a mixture of these two Kekulé states, i.e., $\bar{u}>0$ and $\bar{u}^{\prime}>0$. We note that the constant volume condition is imposed just for simplicity. Since applying strain leads to an expansion of the volume, we could naively employ the lattice constant $a$ as a control parameter. In this case, the effect of the increasing $a$ is naturally incorporated in Eq. (1) by decreasing the transfer integral $t$, while the on-site Coulomb repulsion $U$ is expected to be less affected. Then, the effect of strain can be mimicked by increasing $U/t$, which is the conventional control parameter in the Hubbard model, instead of actually increasing $a$.

Let us begin by presenting the temperature dependence of the order parameters for the two Kekulé states as shown in Figs. 3 and 4 for $g/K=0.6$ and $g/K=0.4$, respectively. These values of the dimensionless parameter for the e-l coupling are representative for the soft and hard lattices in our model. In the first place, out of the three possible lattice distortion patterns, we observe that only the KVBS state is stabilized with $\bar{u}>0$ and $\bar{u}^{\prime}$=0 for all parameters we studied at low temperatures. We do not know why the KRVB state is not stabilized. Generally speaking, analytical methods would be eventually needed to understand the stability of different types of order [64, 65]. However, considering that the absence of the KRVB state is consistent with the ab initio QMC study [41], we suppose that it is not simply due to the model simplification or the technical treatment, such as the constant volume condition. For the soft lattice, as shown in Fig. 3, the phase transitions appear to be continuous for all $U/t$ and $L$, and the KVBS state is most stable with large $\bar{u}$ at $U/t=4\simeq U_{\mathrm{c}}/t$. On the other hand, in the case of the more rigid lattice (Fig. 4), the KVBS state is fragile. At $U/t=3$, the order parameter $\bar{u}$ depends on $L$, converging to the thermodynamic limit value only with $L\geq 8$. In addition, the converged value, i.e., $\bar{u}\simeq 0.055$, is quite small. At $U/t=5$, the temperature dependence of $\bar{u}$ for the smaller lattice ($L<6$) appears to signal that the transition is of first order. However, this apparent behavior can be ascribed to the finite-size effect. Since, for this large value of $U/t=5$, the antiferromagnetic (AF) correlation develops as $T/t$ decreases, its correlation length can reach the system size of small $L$ above the critical temperature of the Kekulé transition. In this case, the transition to the KVBS state which is triggered by further lowering $T/t$ can share features with the transition between KVBS and AFMI at the ground state, which is of first order as discussed later. However, the result of $L=6$ in Fig. 4(c) suggests that the Kekulé transition is more likely continuous.

The competition between KVBS and AFMI is also observed in the magnetic property. Figure 5 shows temperature dependence of the AF structure factor defined as

where $\bm{S}_{i}=\frac{1}{2}\sum_{s,s^{\prime}}c^{\dagger}_{is}(\bm{\sigma})_{ss^{\prime}}c_{is^{\prime}}$ is the spin operator at site $i$ with $\bm{\sigma}$ being the vector of Pauli matrices, and the sum $i\in A(B)$ runs over sites belonging to $A$ $(B)$ sublattices in the undistorted lattice [see Fig. 1 (a)]. Since in the usual Hubbard model without the e-l coupling the AF long-range order exists only in the ground state, $S_{\mathrm{AF}}$ for $g/K=0$ develops with decreasing $T/t$ and saturates at a low temperature where the AF correlation length exceeds the system size $L$. On the other hand, with the e-l coupling, $S_{\mathrm{AF}}$ shows a sudden drop at the KVBS transition temperature and does not develop thereafter, which indicates that the KVBS state has no AF long-range order.

We determine the critical temperatures for the KVBS transition from the temperature-dependence of $\bar{u}$ for several values of $U/t$ and $g/K$, and the results are summarized in the finite-temperature phase diagram of Fig. 6. The significant feature is that for small $g/K$, which corresponds to a more rigid lattice as in the case of graphene, the region of the KVBS phase is, although it is small, present in the vicinity of the critical point of the Hubbard model on the undistorted honeycomb lattice at $U_{\mathrm{c}}/t\simeq 3.8$ [13, 22, 23, 32, 34]. Thus, if the value of $U/t$ in real graphene is slightly smaller than $U_{\mathrm{c}}/t$, which is indeed anticipated from the first-principles calculations of $U=U_{00}=9.3\mathrm{eV}$ [29] and $t\approx 2.7\mathrm{eV}$ [25, 26], application of strain is expected to bring about the phase transition to the KVBS state.

To investigate more closely the vicinity of $U_{\mathrm{c}}/t$ where the order parameter $\bar{u}$ is small, we directly calculate the free energy as a function of $u$, i.e, $F(u)$ by numerically integrating its derivatives $\partial F/\partial u$, which is given by the left-hand side of Eq. (3), and determine $\bar{u}$ by locating the minimum of $F(u)$. This rather naive approach is computationally feasible since we have confirmed that the KRVB order does not appear and hence only one order parameter, $u$ for KVBS, is relevant in the free energy. We calculate the derivatives at about a hundred discrete $u$ points by the QMC simulations, of which the computational cost is comparable to that of the iterative method using the Robbins-Monro algorithm. The advantage of this method is that once the QMC simulations for the discrete $u$ points at a fixed value of $U/t$ are performed, we can determine $F(u)$ and $\bar{u}$ for any value of $g/K$ as shown in Fig. 7 for $U/t=3$, 4, and 5 at $T/t=0.05$. For $U/t=3<U_{\mathrm{c}}/t$, the finite-size effect is evident even up to $L=8$. This implies that a high resolution around the Dirac points in the momentum space is required to study the SM-KVBS transition at low temperatures. On the other hand, for $U/t>U_{\mathrm{c}}/t$, where the ground-state of the system with $\bar{u}=0$ is AFMI, the finite-size effect is small since the AF correlation length exceeds the system size $L$, and the transition from AFMI to KVBS is confirmed to be of first order.

We plot the KVBS order parameter $\bar{u}$ obtained by this approach at a fixed value of $g/K=0.4$ as a function of $U/t$ in Fig. 8(b), together with the critical temperature in Fig. 8(a) which is the same data used in Fig. 6. Now it is clear that the KVBS state emerges as an intermediate phase between the SM and AFMI phases, indicating that the Kekulé-like structural deformation is driven by the correlation before the AF long-range order sets in. This result is primarily consistent with the preceding ab initio QMC study [41].

Finally, let us explore the phase diagram at a fixed temperature $T/t=0.05$ shown in Fig. 9. As discussed above, the finite-size effect is strong in the phase boundary between the SM and KVBS phases. We calculate the critical value of $g/K$ for a large cluster $L=16$ in the noninteracting case and confirm that it is fairly close to that for $L=8$. Furthermore, there is less finite-size effect at $U/t\lesssim U_{\mathrm{c}}/t$. We thus expect that the phase boundary estimated for $L=8$ is sufficiently close to that in the thermodynamic limit. On the other hand, the critical value of $g/K$ separating the KVBS and AFMI phases hardly depends on the system size. This critical point is of first order and increases monotonically with $U$. These features are in accord with recent results of the Su-Schrieffer-Heeger Hubbard (SSHH) model on the square lattice, where quantum phonon dynamics is fully taken into account [66, 67, 68]. The first-order transition between the VBS and AF phases is also found in the spin-Peierls model on the honeycomb lattice [69], which corresponds to the SSHH model in the strong coupling limit ($U/t>>1$). Therefore, we expect that the first-order nature of the transition between these phases, being within the Landau-Ginzburg-Wilson paradigm, is ubiquitous especially near the adiabatic limit [69]. It is also pointed out that our phase diagram featuring a tricritical point with the SM, AFMI, and KVBS phases is similar to that obtained for a model of Dirac fermions devised from a designer Hamiltonian approach [70], except that the AFMI-KVBS transition is continuous because of the emergent SO(4) symmetry in the devised model.

To summarize, we have investigated the valence-bond-solid phase in the Hubbard model on the honeycomb lattice with the electron-lattice coupling, motivated by the recent ab initio calculation predicting the correlation-driven Kekulé valence-bond-solid phase in graphene. Considering the lattice displacements in the adiabatic limit, the model is readily trackable by the auxiliary-field quantum Monte Carlo method. Furthermore, based on the symmetry argument by Frank and Lieb, the possible distortions in the Hubbard model for graphene are restricted to have the unit cell containing at most six sites. Under this restriction as well as the constant-volume condition, we have obtained the finite-temperature phase diagram by determining the most stable displacement configuration. Out of the three possible bond-order patterns, only the Kekulé valence-bond-solid phase with the dimerized distortion is found to be stabilized. For a more rigid lattice, as in the case of graphene, the Kekulé valence-bond-solid phase exists in a narrow region near the critical point of the metal-insulator transition in the Hubbard model without the electron-lattice coupling. This implies that if the real graphene is located in the semimetallic side near this critical point in the phase diagram as a function of the interaction, applying strain, which should correspond to increasing the effective interaction, can bring about the valence-bond-order transition and hence introduce the gap in graphene. We have also discussed the phase diagram at a low temperature, where the semimetal, the antiferomagnetic Mott insulator, and the Kekulé valence-bond-solid phases are competing. The phase boundary between the antiferromagnetic Mott insulator and the Kekulé phase is confirmed to be of the first order, which falls in the Landau-Ginzburg-Wilson paradigm.