Effective spin couplings in the Mott insulator of the honeycomb lattice Hubbard model

We have calculated the dependence of the magnetic exchange couplings on t/U using perturbative continuous unitary transformations (pCUTs) [22–24] along the lines of [10, 25]. The pCUT provides the magnetic couplings as series expansions in t/U. Note that the spectrum of the Hubbard model is symmetric at half-filling under the exchange $t\leftrightarrow -t$ and therefore only even order contributions are present [26]. We have determined all two-spin, four-spin and six-spin interactions up to order 14. In order 14 one has to calculate 345 topologically different graphs.

The bare series do not converge in the spin-liquid region of the Mott phase [20]. This is different from the recently analyzed case of the Hubbard model on the triangular lattice [10, 20] where the intermediate non-magnetic phase is already realized for smaller ratios of the bandwidth W to the interaction U. It is therefore mandatory to apply resummation schemes of the series which is very complicated to be performed reliably for the full set of exchange couplings. We have therefore applied the recently developed gCUTs [20] in order to derive non-perturbatively the effective spin model. In the following, we use the gCUT results to study the physics of the full low-energy spin model. The pCUT results are only used to analyze the most important terms in the effective model.

The general properties of the gCUT are discussed in [20]. Here we only mention the specific approach for the current problem. The basic idea is to use a CUT to map $\skew4\hat {H}$ unitarily to an effective Hamiltonian $\skew4\hat {H}_{\mathrm {eff}}$, which has the special property that the block without double occupancies representing the effective low-energy spin model is decoupled from the rest of the Hamiltonian. The convergence of the gCUT is not triggered by a small parameter but relies on the fact that the correlation length of charge fluctuations is finite in the Mott insulator.

In the gCUT one generates all topologically distinct connected graphs $\mathcal {G}_\nu$ of the lattice and one sorts them by their number of sites n. On each graph $\mathcal {G}_\nu$ a CUT is implemented by setting up the finite number of flow equations [27–29]. A continuous auxiliary variable l is introduced defining the l-dependent Hamiltonian $\skew4\hat {H}^{{\mathcal {G}}_\nu } (l):=U^\dagger (l)\skew4\hat {H}^{{\mathcal {G}}_\nu }U(l)$. Then the flow equation is given by

$$\begin{equation} \partial_l \skew4\hat{H}^{{\mathcal{G}}_\nu}(l) = \left[\skew4\hat{\eta} (l), \skew4\hat{H}^{{\mathcal{G}}_\nu}(l)\right], \end{equation} \tag{ 3 }$$

where $\hat {\eta } (l)$: = −U^†(l)(∂_lU(l)) is the anti-Hermitian generator of the unitary transformation. At the end of the flow l = ∞, one obtains an effective graph-dependent Hamiltonian $\skew4\hat {H}_{\mathrm {eff}}^{{\mathcal {G}}_\nu }$.

We have used the generator introduced by Wegner [27]

$$\begin{equation} \skew4\hat{\eta}^{\mathrm{Wegner}} (l) = \left[ \skew4\hat{H}^{\mathrm{d}} (l) , \skew4\hat{H}^{\mathrm{nd}} (l) \right], \end{equation} \tag{ 4 }$$

where the diagonal part $\skew4\hat {H}^{\mathrm {d}}$ is given by all matrix elements between states with the same number of double occupancies and $\skew4\hat {H}^{\mathrm {nd}}$ denotes all the remaining non-diagonal processes. We stop the flow for each graph once the residual off-diagonality (ROD) is below 10⁻⁹. Here the ROD is defined as the square root of the sum of all non-diagonal elements which connect to the low-energy subspace. Finally, one obtains the effective spin model $\skew4\hat {H}_{\mathrm {spin}}$ in the thermodynamic limit by subtraction of sub-cluster contributions and by embedding back the graphs on the honeycomb lattice [20].

The numerical effort scales exponentially with the number of sites n of a given graph. Here we have treated all graphs up to seven sites. This amounts to solving up to one million differential equations for the most demanding symmetry sector of a single graph for each value of t/U. Physically, one captures all charge fluctuations on this length scale. Let us remark that it is of course possible to only include a restricted set of graphs in the gCUT calculation. This will be used below for a direct comparison with CORE. Finally, one obtains for a given t/U the effective spin model non-perturbatively.

A view of the full effective spin model obtained by gCUTs for different values of t/U is shown in figure 1. Let us remark that for gCUT(7) the full low-energy model can be uniquely expressed in terms of spin operators, because multi-spin interactions involving more than six sites are not yet allowed for this cluster size.

Zoom In Zoom Out Reset image size

There are several implications which can be directly read off from figure 1. Firstly, only at rather small values of t/U ≲ 0.1 can a clear hierarchy in the amplitudes of the effective spin operators be seen. This corresponds to the purely perturbative regime where the importance of terms is fully given by the perturbative order where a spin coupling appears for the first time. In this regime the frustrated next-NN Heisenberg exchange J₂ (order 4 perturbation theory) is the leading correction to the dominant Heisenberg interaction J₁ (order 2 perturbation theory).

Secondly, the situation changes drastically in the intermediate coupling regime where the spin-liquid phase is realized. Here, no clear hierarchy in the amplitudes of the effective spin couplings can be detected. One is rather confronted with proliferating numbers of terms which is likely a consequence of the fact that the metal–insulator transition is second order. Consequently, the charge gap closes continuously and one expects an increasing number of effective spin couplings on increasing length scales when approaching the quantum critical point. Let us note that this is different for the Hubbard model on the triangular lattice [10]. In this case the metal–insulator transition is first order and one has no diverging length scale upon approaching the Mott transition from the insulating side.

Thirdly, and despite the complexity of the full effective spin model, one still expects that most of the properties of the Mott phase should be contained in a minimal magnetic model where one only considers the most important spin couplings having the largest amplitudes. Some evidence for this strategy is given in the next section by analyzing such a minimal magnetic model. Interestingly, we find that the perturbative hierarchy is not valid anymore for intermediate couplings: the most important correction to J₁ is not J₂ but rather six-spin interactions located on hexagons. These six-spin interactions arise in order 6 perturbation theory and originate dominantly from ring exchange processes around the shortest loop of even length on the honeycomb lattice corresponding to hexagons. Let us remark that this is very similar to the Hubbard model on the triangular lattice [10]. Here the most elementary loop of even length is a four-site plaquette. Consequently, four-spin interactions on the elementary plaquette are the dominant correction to the Heisenberg model.

In the remainder of this paper we are aiming at a minimal magnetic model which contains only the most important spin interactions in the regime of strong to intermediate coupling. In this parameter regime, the minimal model should display the low-energy properties of the Hubbard model on the honeycomb lattice. Additionally, we expect that going beyond the minimal model by considering the full graph decomposition used by gCUT (or alternatively by CORE) yields successfully even better results. The most important effective exchanges are depicted in figure 2: NN (J₁) and next-NN (J₂) two-spin interactions appearing in the J₁–J₂ model and the six-spin terms with couplings L₁,L₂,...,L₅. The corresponding gCUT and pCUT amplitudes for these couplings are displayed in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One clearly sees that the bare series of six-spin interactions are not converged in the interesting regime of intermediate t/U, confirming the conclusion obtained in [20]. In contrast, the results obtained by gCUT and self-similar extrapolation of the pCUT series are in fair agreement for all displayed couplings.

In the following, we want to further strengthen these findings by using CORE as an alternative tool to derive effective low-energy models.

It is also possible to derive non-perturbative effective theories by relying on the CORE technique [30]. CORE is a real-space renormalization procedure, applicable to generic lattice models, where effective models are derived by truncating local degrees-of-freedom. Such a construction requires the exact computation of low-lying eigenstates on finite connected graphs, similarly to what happens in the previously discussed gCUT procedure; implementation details are discussed in the literature [30, 31]. The resulting effective Hamiltonian is given as a cluster expansion

$$\begin{equation} {\cal H}_{\mathrm{eff}}= \sum_{g} h_g^{\mathrm{c}}, \end{equation} \tag{ 5 }$$

where the sum takes place over a set of graphs g, and h^c_g corresponds to the connected term that is obtained by subtracting contributions from embedded sub-clusters [30, 31]. Such an expansion must naturally be truncated at a certain maximum range [32], but its accuracy can in principle be controlled by analyzing the convergence of terms in the effective model with increasing range.

Since it is our goal here to obtain effective spin interactions for the Hubbard model (equation (1)) at half-filling, we select single sites as elementary blocks in the CORE expansion [30, 31] and retain only singly occupied states on all sites.

The clusters used in the CORE calculation are shown in figure 4. Here two setups are considered: (i) the first choice of clusters does not contain the four-site star graph displayed in figure 4(c). It therefore corresponds to the leading term of a full graph expansion in terms of hexagons. (ii) The second choice contains additionally the four-site star cluster (figure 4(c)). In the hexagon expansion mentioned before, such a graph would be included only on the level of three-hexagon clusters. Later we see that gCUT (restricted to exactly the same graphs) yields almost identical results for both choices of clusters. The second choice can then be regarded as an intermediate step between the minimal one-hexagon calculation and the full gCUT(7) calculation which includes all graphs shown in figure 4 as a subset.

Zoom In Zoom Out Reset image size

We start by discussing the first choice. The lowest-order contribution to the NN spin-exchange J₁ is simply computed by solving the Hubbard model on a two-site cluster (figure 4(a)): in this case only eigen-energies are required, since J₁ exactly corresponds to the singlet-triplet gap:

$$\begin{equation} J^{(2)}_1=\frac{U}{2} \left( \sqrt{1+\frac{16 t^2}{U^2}}-1 \right). \end{equation} \tag{ 6 }$$

Although trivially obtained, this result is non-perturbative, reducing to the well-known limit J₁ ≃ 4t²/U only when t/U ≪ 1 (see figure 5(b)). Let us mention that the same non-perturbative result is obtained when calculating J⁽²⁾₁ with gCUT(2).

Zoom In Zoom Out Reset image size

Longer-range couplings, along with corrections to shorter-ranged ones, are obtainable from the analysis of the larger clusters depicted in figures 4(b)–(d). In particular, in figure 5(a) we plot bare effective couplings, obtained by considering the six-site cluster in figure 4(d), as a function of t/U. Here the notion 'bare couplings' corresponds to the effective exchange couplings of a single cluster without considering subtractions and embeddings. We first focus on the next-NN coupling J₂ and note that, particularly in the range 0.233 ≲ t/U ≲ 0.286(4.3 ≳ U/t ≳ 3.9) where a QSL appears [11], J₂ ≲ 0.1J₁ and therefore this frustrating interaction is not large enough [15–19] to account for the absence of long-range Néel order before a semi-metal phase is stabilized [20]. However, we also note that, remarkably, much larger magnitudes are associated with the six-spin terms with couplings L₁,L₂,...,L₅ depicted in figure 2 (we also briefly remark that the signs for their amplitudes are in agreement with those obtained from a decomposition of the cyclic six-spin permutation, cf [33]). Furthermore, we observe that the amplitudes for longer-range two-spin, as well as four-spin, terms are found to be negligible in comparison to J₁ and L₁,L₂,...,L₅.⁷ One is thus inclined to conclude that such six-spin interactions are the main ingredients in stabilizing a QSL in the phase diagram for equation (1), a conjecture we aim at testing in section 4.

Let us mention that such six-spin interactions on local hexagons are generated in order 6 perturbation theory. In this order the couplings L₁–L₅ have the same absolute prefactor with an alternating sign which exactly corresponds to the six-spin interactions contained in the six-spin ring exchange permutation operator. But we stress that the latter operator also contains two-spin and four-spin interactions of the same magnitude which one has to contrast with our finding, both perturbatively and non-perturbatively, of suppressed two- and four-spin interactions. The full non-perturbative contribution of a single hexagon is therefore dominated by the six-interactions having the largest exchange couplings. Finally, we note that the non-perturbatively obtained couplings L₁–L₅ have different prefactors as can be seen in figure 5.

We proceed by analyzing the dependence of the effective couplings on the maximum range of the clusters. In order to do so, for each cluster depicted in figure 4 we subtract connected contributions from shorter-ranged, embedded, sub-clusters, according to the standard CORE recipe [30, 31]. We plot the results for the two-spin exchanges with couplings J₁ and J₂ as a function of t/U respectively in figures 5(b) and (c), and find a seemingly fast convergence with increasing cluster range up to the six-site cluster depicted in figure 4(d). Furthermore, the data in figures 5(b) and (c) are in good agreement with gCUT results, as we discuss below in section 3.3.

Finally, we comment on the fact that a few clusters larger than the ones in figure 4, such as the one comprising ten sites and forming an edge-sharing double hexagon, are amenable to numerical diagonalization, in principle implying that our CORE expansion could be extended to even longer ranges than considered here. To this end, a full graph decomposition must be implemented for the CORE approach; this is left as a task for future research. Here we would rather like to focus on the most minimal setup for describing the low-energy physics of the Hubbard model.

We compare now the effective models obtained from gCUT (section 3.1) and CORE (section 3.2). We first remark that both approaches yield qualitatively similar effective Hamiltonians, which display as dominant terms the two- and six-spin exchanges depicted in figure 2 (see figures 3 and 5). We therefore focus on such terms that may be regarded as defining a minimal effective model and plot their amplitudes as a function of t/U in figure 6. Let us mention that we do not show the constant term of the effective models which also displays a t/U dependence and which one therefore has to keep in the ED when comparing to the Hubbard model. Here three different sets of graphs are considered: (i) the leading graphs of a hexagon expansion not including the four-site star graph (CORE and gCUT), (ii) the leading graphs of a hexagon expansion including the four-site star graph (CORE and gCUT) and (iii) all graphs comprising up to seven sites (gCUT).

Zoom In Zoom Out Reset image size

We note that results of both methods are in good quantitative agreement up to couplings stabilizing the QSL phase [11] and beyond the perturbative regime [20] for all three choices. In particular, CORE and gCUT yield almost identical results for the restricted graph sets (i) and (ii). This is likely a consequence of the fact that the structure of the low-energy subspace is rather constrained due to the high symmetry of the considered clusters (the four-site graph as well as the six-site graph has a rotational symmetry) plus the number of different spin couplings defined on these graphs is rather small.

The most noticeable effect when going from the minimal choice to the full graph decomposition including all graphs up to seven sites is an increasing value for the next-NN two-spin exchange J₂. Nevertheless, we also remark that all three sets agree in that, unlike what was assumed in [12, 13], the effective values for the next-NN two-spin exchange J₂ are not large enough to destroy Néel order [15–19], suggesting that the QSL phase may instead be accounted for by the six-spin ring exchanges with amplitudes L₁,L₂,...,L₅ in figure 2. Furthermore, we show in the next section that the minimal effective spin model having the lowest J₂ gives the best agreement with the QMC results at intermediate couplings.

We start by discussing the properties of the minimal effective spin model obtained for the minimal set of graphs not including the four-site star graph. We stress again that CORE and gCUT give essentially identical low-energy spin models for this choice of graphs. In the following, we use the effective model obtained via the CORE algorithm (section 3.2).

We start by comparing the low-energy spectra for the Hubbard (equation (1)) and the effective (section 3) models on an N = 18 cluster for two values of t/U: t/U = 0.05 (deep into the Néel phase) and t/U = 0.25 (QSL phase, according to [11]). In figures 7(a) and (b), we plot the so-called Anderson tower of states for the Hubbard model. As we only make use of S_z symmetry (instead of the total spin), we plot the energy states versus S_z(S_z + 1), so that degenerate states with identical energies correspond to spin multiplets. We have also computed the one-particle gap, Δ_1P = (E₀(N/2 + 1) + E₀(N/2 − 1) − 2E₀(N/2))/2, where E₀(N_f) is the ground-state energy with N_f fermions. For t/U = 0.05 (figure 7(a)), the so-called Anderson tower of states is typical of a collinear Néel state ([34] and references therein) with the lowest excitation scaling as S(S + 1) (with a slope that should decrease with decreasing 1/N), followed by spin-wave excitations (with the energy scale $1/\sqrt {N}$). Moreover, on the scale of the figure, all states correspond to spin excitations since Δ_1P is quite large⁸. The same behavior is reproduced qualitatively and quantitatively by the effective spin model, as expected and shown in figure 7(c). Note that in the case of the effective spin model, we have been able to label each eigenstate with its total spin value S so that we plot the levels as a function of S(S + 1).

Zoom In Zoom Out Reset image size

When t/U is increased to 0.25, a distinct picture emerges (figure 7(b)). While the first excitations still correspond to S = 1 and 2 states at the Γ point, there is no longer a straight S(S + 1) line of excitations and, even more importantly, the gap Δ_1P to charge excitations is substantially lowered, beyond which a description solely based on spin degrees-of-freedom is no longer valid. However, according to the QMC results [11], for this coupling the one-particle gap should remain finite in the thermodynamic limit, so that our effective spin model could still describe the physics at the lowest energies: indeed, its low-energy spectrum, plotted in figure 7(d), reproduces the lowest spin excitations quite accurately. Furthermore, while the quantum numbers seen in the tower-of-states are consistent with standard collinear magnetic ordering, the lowest spin excitations do not display a clear S(S + 1) behavior, which might be an indication that there is no magnetic ordering.

It would also be desirable to compute correlators for the ground state of the effective model, but this is quite involved since one would also need to renormalize the operators when deriving an effective Hamiltonian [10, 20, 30, 31]. Nevertheless, some observables may straightforwardly be obtained from the ground-state energy per site, e₀, by relying on the Feynman–Hellmann theorem. For instance, the electronic double occupancy per site can be computed from ED of the effective model involving only spin degrees-of-freedom, as has been done in [10] using

$$\begin{equation} \langle n_{\uparrow} n_\downarrow \rangle = \frac{\mathrm{d} e_0}{\mathrm{d}U}. \end{equation} \tag{ 7 }$$

The results are shown in figure 8(a) and compared both to the exact result of the Hubbard model on N = 18 cluster and to QMC data⁹. Since this is a local quantity, its finite size effects are rather small as expected. We observe good agreement not only for small t/U (deep in the Néel phase), but also in the interesting regime t/U ≃ 0.25 where QSL is expected to emerge. As a side remark, let us emphasize that the discrepancy between the results on the Hubbard model and with the CORE effective one on N = 18 can be attributed to the role of short-length loops (of length 6) that exist on this small cluster, and are not captured by the effective model (see [10, 15] for a similar discussion). In fact, effective model results agree quite well with the data obtained on a much larger lattice with QMC.

Zoom In Zoom Out Reset image size

In [11], the analysis of the behavior of the kinetic energy density $E_{\mathrm {kin}}=\langle -t \sum _{\langle ij\rangle ,\sigma } (c^\dagger _{i\sigma } c_{j\sigma }+\mathrm {h.c.})\rangle /N$ versus U was also discussed in the context of the formation of local moments at strong U/t, in contrast with the itinerant regime for small U/t, whereas the QSL is stabilized in between. Using the Feynman–Hellmann theorem again, we can simply compute this quantity as

$$\begin{equation} \frac{\mathrm{d} E_{\mathrm{kin}}}{\mathrm{d}U} = - U \frac{\mathrm{d}^2 e_0}{\mathrm{d}U^2}. \end{equation} \tag{ 8 }$$

We have chosen a grid of U/t going from 3 to 8 by steps of 0.1 and computed the ground-state energy for various clusters. Using a finite-difference approximation, we obtain the derivative of the kinetic energy that we plot against t/U in figure 8(b). Our first comment is that we have some agreement with the exact QMC data [11] that show a maximum around 0.14 for U/t ≃ 5.0, although this quantity is expected to be quite sensitive to details since it is a second derivative of the ground-state energy. Note that QMC data are also quite noisy in this large t/U region, and we refer to the original publication [11] for additional data: the point is that a maximum of order 0.15 is reached for t/U ≃ 0.2. Our second remark is that we observe an anomaly around t/U = 0.3 which could signal that our CORE truncation is not safe beyond this value. This gives us confidence that at least for t/U ≃ 0.25, our approach could make sense, and this includes the spin-liquid phase. Lastly, since we observe this change of curvature around U/t = 4.5 in a spin-only model, it seems to us that the virtual charge fluctuations giving rise to our effective low-energy model are sufficient to account for this phenomenon. Additionally, it might be unsurprising that we underestimate this effect by our minimal magnetic model, because we expect that the large number of neglected additional spin operators as well as small renormalizations of the treated spin couplings could well lead to an enhanced itineracy of the ground state at intermediate t/U values.

As a conclusion on this part, while we are clearly not able to reach system lengths where QSL behavior is expected to occur—which would require several hundred sites—our simple effective model is able to reproduce both the low-energy properties and the ground-state local properties. It suggests that these properties may not be linked to the presence of a QSL and are robust short-distance features.

The results of the last subsections are very promising. CORE and gCUT give an almost identical minimal effective spin model which is obtained from the single hexagon graph and which compares well to the results obtained by QMC for the Hubbard model on the honeycomb lattice.

In the following, we want to go beyond this minimal choice of graphs and would like to see whether the above findings are confirmed and strengthened. We therefore compare the above results for the minimal model to EDs for the effective spin models obtained (i) by including the four-site star graph to the minimal set and (ii) by including all graphs up to seven sites. For case (i) we again use the results of the CORE algorithm, but we stress that the gCUT on the restricted set of graphs gives basically the same low-energy model. For the second case the results of gCUT(7) are used.

We start by comparing the ground-state energy per site as a function of t/U, which is displayed in figure 9. It can be clearly seen that all three different sets of graphs give a very similar energy per site in the full Mott phase, which also agrees well with the results of the Hubbard model. This is promising, but it is also a consequence of the fact that the ground-state energy per site is a rather insensitive quantity.

Zoom In Zoom Out Reset image size

We therefore turn to more sensitive quantities, namely the double occupancy per site and the derivative of the kinetic energy as discussed in the last subsection for the minimal effective model (figure 10). Consequently, details of the considered effective spin models matter, which is most clearly seen for the derivative of the kinetic energy. Interestingly and surprisingly, agreement between the results of the effective spin model and the results of the Hubbard model gets worse at intermediate couplings. Indeed, already the minimal set of graphs plus the four-site star graph yield the wrong tendency. Furthermore, ED of the full effective spin model from gCUT(7) does not even display the maximum in the derivative of the kinetic energy.

Zoom In Zoom Out Reset image size

Altogether, we find that the inclusion of more graphs in our calculation gives reduced agreement at intermediate couplings, in contrast to our expectation. This trend is most probably neither a numerical problem nor a difference between CORE and gCUT, because both techniques agree quantitatively for both restricted sets of graphs. In our opinion there are basically two possible scenarios. Firstly, an effective spin model can still be reliably derived at intermediate couplings using a hexagon expansion. Indeed, the minimal model only focusing on the leading one-hexagon graphs gives very nice results. It might then be possible that the full graph decomposition used for the gCUT by including all graphs up to seven sites does not give improved results, because the relevant low-energy physics is contained in graphs with longer loops (these are exactly the multi-hexagon graphs). Secondly, the trend seen in the full graph decomposition by gCUT is a real effect; that is, a controlled derivation of an effective low-energy spin model for the QSL at intermediate couplings remains a big challenge.