Slater and Mott insulating states in the SU(6) Hubbard model

TL;DR

This study uses large-scale quantum Monte Carlo simulations to explore the transition from Slater to Mott insulating states in the half-filled SU(6) Hubbard model, revealing a crossover from antiferromagnetic to valence bond solid order as interaction strength increases.

Contribution

It provides the first detailed numerical analysis of the SU(6) Hubbard model, identifying the nature of phase transitions and critical exponents in this high-symmetry system.

Findings

Antiferromagnetism persists in SU(2) and SU(4) but transitions to valence bond solid in SU(6).

Single-particle gap scales linearly with U in the Mott regime.

Critical exponents for antiferromagnetic transition are estimated as ν=0.60 and η=0.44.

Abstract

We perform large scale projector determinant quantum Monte-Carlo simulations to study the insulating states of the half-filled SU(6) Hubbard model on the square lattice. The transition from the antiferromagnetic state to the valence bond solid state occurs as increasing the Hubbard $U$. In contrast, in the SU(2) and SU(4) cases antiferromagnetism persists throughout the entire interaction range. In the SU(6) case, antiferromagnetism starts to develop in the weak interacting regime based on the Slater mechanism of Fermi surface nesting. As $U$ passes a crossover value $U^*/t\approx 9$, the single-particle gap scales linearly with $U$, marking the onset of Mott physics. In the Mott regime, antiferromagnetism becomes to be suppressed as $U$ increases, and vanishes after $U$ passes the critical value $U_{\rm AF,c}/t=13.3\pm 0.05$. The critical exponents are obtained via critical scalings as $\nu_{\rm AF}=0.60\pm 0.02$ and $\eta_{\rm AF}=0.44\pm 0.03$. As $U$ further increases, the valence bond solid ordering appears exhibiting the anomalous dimension $\eta_{\rm VBS}=0.98\pm 0.01$.

Equations

$$H=-t\sum _{⟨i,j⟩,\alpha }\left(c_{i\alpha }^{†}{c}_{j\alpha }+h.c.\right)+\frac{U}{2}\sum _i{\left(n_i-\frac{N}{2}\right)}^2,$$

$$m_r=\frac{1}{6}\left(\sum _{\alpha =1}^3{n}_{r\alpha }-\sum _{\alpha =4}^6{n}_{r\alpha }\right),$$

$$S_{mQ}=\frac{1}{L^2}\sum _{r{r}^{\prime }}⟨m_r{m}_{r^{\prime }}⟩(-1{)}^{r-r^{\prime }}.$$

$$d_{r\hat{e}}=\sum _{\alpha =1}^6\left(c_{r\alpha }^{†}{c}_{r+\hat{e}\alpha }+c_{r+\hat{e}}^{†}{c}_{r\alpha }\right),$$

$$d_{\hat{x}(\hat{y})}=\frac{1}{L^2}\sum _r{d}_{r\hat{x}(\hat{y})}(-1{)}^{r_x(r_y)}.$$

$$S_{vbs}=\frac{1}{L^2}\sum _{r{r}^{\prime },\hat{e}}⟨d_{r\hat{e}}{d}_{r^{\prime }\hat{e}}⟩(-1{)}^{r_{\hat{e}}-r_{\hat{e}}^{\prime }}.$$

$$d_{r\hat{e}}\propto \frac{t}{U}\sum _{\alpha \beta }{c}_{r\alpha }^{†}{c}_{r\beta }{c}_{r+\hat{e}\beta }^{†}{c}_{r+\hat{e}\alpha },$$

$$m_Q^2(L)=S_{mQ}/L^2,O_{vbs}^2(L)=S_{vbs}/L^2.$$

$$f(L)=a+\frac{b}{L^c},$$

$${\xi }_{\mathrm{AF}}=\frac{1}{q}\sqrt{\frac{m_Q^2}{m_{Q+q}^2}-1},$$

$${\xi }_{\mathrm{AF}}(U,L)=Lf\left[|U-U_{\mathrm{AF},\mathrm{c}}|L^{1/{\nu }_{\mathrm{AF}}}\right],$$

$$⟨(-{)}^r{m}_r{m}_0⟩$$

$$m_Q^2\sim \frac{1}{L^{d+z-2+{\eta }_{\mathrm{AF}}}},$$

$$U_{\mathrm{AF},\mathrm{c}}$$

$${\nu }_{\mathrm{AF}}$$

$${\eta }_{\mathrm{VBS}}=0.98\pm \mathrm{0.01.}$$

Full text

I introduction

How repulsive interactions turn a partially filled electron band into the insulating state is an important question of strong correlation physics. In the presence of nested Fermi surfaces, the antiferromagnetic (AF) order appears at infinitesimal interactions based on Fermi surface nesting. The single-particle gap is at the same order of the AF gap function Slater (1951). Such a state is known as the Slater insulator exhibiting strong charge fluctuations. On the other hand, charge fluctuations are frozen in the strong interaction regime due to the large charge gap linearly scaling with the repulsive interaction Mott (1949), and such a state is Mott insulating. The low energy physics lies in the spin channel arising from the superexchange among local spin moments. Mott insulators can even exhibit no symmetry breaking, for example, the 1D Hubbard model at half-filling exhibits power-law AF correlations and charge gap, but without long-range AF ordering.

However, in strongly correlated electron systems, the above two pictures of insulators are often blended together Imada et al. (1998); Lee et al. (2006). For example, in the half-filled SU(2) Hubbard model on the square lattice, Fermi surface nesting leads to the Slater AF state at weak $U$, while the strong $U$ side is effectively described by the Heisenberg model and attributed to the Mott AF insulator. Both regimes exhibit the commensurate Neel ordering smoothly connected by a crossover Hirsch (1985); *Hirsch_PRL_1989; Pruschke and Zitzler (2003).

In recent years, the rapid development of ultra-cold atom physics provides a new route to investigate strong correlation physics. It was proposed that cold fermions with multiple spin components are ideal systems to study high symmetries that typically are met in high energy physics Wu et al. (2003); *Wu_MPLB_2006; Honerkamp and Hofstetter (2004); Wu (2012). For example, the spin-3/2 fermion systems possess a generic Sp(4) symmetry without fine-tuning, which is further enlarged to SU(4) when the interaction is spin-independent. These symmetries play an important role to study novel quantum magnetism beyond SU(2) Wu et al. (2003); *Wu_MPLB_2006; Wu (2005); Hung et al. (2011); Xu and Wu (2008); Wang et al. (2014); Chen et al. (2005). The study of high symmetric ultra-cold fermions has been attracting considerable interests both experimental and theoretical recently Wu et al. (2003); *Wu_MPLB_2006; Honerkamp and Hofstetter (2004); DeSalvo et al. (2010); *Taie_PRL_2010; Krauser et al. (2012); *Taie_NP_2012; *Zhang_S_2014; Cazalilla and Rey (2014); *Laflamme_AP_2016.

As for the SU(4) Hubbard model at half-filling on the square lattice, i.e., two fermions per site, a pervious determinant quantum Monte Carlo shows that the AF order is non-monotonic as $U$ increases: After reaching a maximal at $U/t\approx 8$, the AF order starts to decrease but remains finite throughout the interaction range simulated $U/t\leq 20$ Wang et al. (2014). Meanwhile, the system exhibits no valence bond solid (VBS) ordering. A recent study directly on the SU(4) Heisenberg model with the one column self-conjugate representation also shows the survival of the AF order Kim et al. (2019), hence, the AF order should also persist through the entire interaction range. The SU(4) and SU(6) Hubbard models of Dirac fermions in the honeycomb lattice and the $\pi$-flux square lattice exhibit the transition from the Dirac semi-metal phase to VBS state and show the absence of the AF order Zhou et al. (2016, 2018). In contrast, the half-filled SU(6) Hubbard model in the square lattice behaves very differently. The QMC simulations show a transition from the AF state at weak $U$ to the VBS state at strong $U$ Wang et al. (2014); Lang (2015). In the weak $U$ limit, the AF is a direct consequence of the Fermi surface nesting and the Van Hove singularity, while the VBS state is a manifestation of the Mott physics. How such a Slater to Mott transition occurs is an interesting and open question, which is the main aim of the present work.

On the other hand, the quantum phase transition from the AF state to the VBS one is argued to be continuous as a result of the deconfined criticality beyond the Landau-Ginzburg-Wilson paradigm Senthil (2004); *Senthil_PRB_2004; *Levin_PRB_2004. Such a prediction has been supported by numerical simulations in recent years. Sandvik (2007); Melko and Kaul (2008); *Sandvik_PRL_2010; *Kaul_PRL_2012; *Pujari_PRL_2013; *Shao_S_2016; *Nahum_PRX_2015; *Wang_PRB_2016; *Assaad_PRX_2016. But there are also works claiming the first order phase transition Kragset et al. (2006); *Kuklov_AP_2006; *Kuklov_PRL_2008; *Sen_PRB_2010; *Papanikolaou_PRL_2010; *Block_PRL_2013; *DEmidio_PRB_2016; *DEmidio_PRL_2017. Nearly all these models studied so far are based on quantum spin models in which charge fluctuations are frozen. It would be interesting to directly investigate the transition between the AF and VBS states based on the fermionic Hubbard model, which takes into account both charge and spin fluctuations.

An additional motivation of this work is the AF order studied below belongs to the self-conjugate representation of SU(N), which could be described by the $U(N)/[U(N/2)\otimes U(N/2)]$ nonlinear $\sigma$-model MacFarlane (1979); Duerksen (1981); Read and Sachdev (1989a); Read_NPB_1989; Read_PRB_1990. The symmetry class is different from the widely studied CPN-1 model which respects the fundamental representation corresponding to the case of $U(N)/[U(1)\otimes U(N-1)]$. Therefore, the AF phase transition here (if continuous) belongs to a different universality class, whose critical exponents would be desired to calculate to characterize such a university class.

In this work, we apply the projector determinant QMC free of the sign problem to study the half-filled SU(6) Hubbard model on the square lattice. Our main results are shown in Fig. 1. From the slope of the single-particle gap, a crossover from the Slater to Mott insulating regime around $U^{*}/t\approx 9$ accompanied by the obvious enhancement of the AF order. The AF order reaches the maximum around $U/t\approx 10$ and then starts to drop as $U$ further increases. In the Mott insulator side, the vanishing of the AF order occurs at $U_{\rm AF,c}/t=13.3\pm 0.05$, and simulations show a continuous transition with the critical exponents $\nu_{\rm AF}=0.60\pm 0.02$ and $\eta_{\rm AF}=0.44\pm 0.03$. As for the appearance of the VBS order, $\eta_{\rm VBS}=0.98\pm 0.01$, and more nature of this transition will be deferred for a future study.

II The model definition and QMC parameters

The SU($N$) Hubbard model on the square lattice at half-filling is defined as

[TABLE]

where $c_{i\alpha}$ is a fermion annihilation operator with $i$ the site index and $\alpha$ the flavor index satisfying $1\leq\alpha\leq N$. The $t$-term represents hoppings between the nearest neighbour sites, and $t$ is the hopping integral. The $U$-term describes the on-site Hubbard interaction as usual and the onsite particle number $n_{i}=\sum_{\alpha=1}^{N}c_{i\alpha}^{\dagger}c_{i\alpha}$. Eq. 1 satisfies the particle-hole symmetry, hence, the average particle number per site is fixed at $N/2$, and the chemical potential $\mu$ is not shown explicitly. When $N=2$, Eq. 1 goes back to the usual spin-$\frac{1}{2}$ Hubbard model. In this article, we focus on $N=6$.

The half-filled SU($N$) (with even $N$) Hubbard model on a bipartite lattice is free of sign problem in auxiliary field QMC simulations as a result of the particle-hole symmetryCai et al. (2013a), which enables us to perform large scale simulations. Details of the algorithm can be found elsewhere Assaad and Evertz (2008); Wang et al. (2014); Zhou et al. (2016), and will not be repeated here. In our simulations, the projection time $\beta=2L$ is used, which is long enough to achieve convergence for a given linear lattice size $L$ up to $24$. The discrete time slice $\Delta\tau=0.05$ is chosen. For each group of parameters, the simulation is performed on 24 cores with 1000 Monte Carlo steps for warming up and no less than 1000 steps for measurements on each core. The exception is the case of $L=24$ which is performed on 48 cores with 500 Monte-Carlo steps for warming up and no less then 500 steps for measurements.

For later convenience, we define the following correlation functions. For the AF order, due to the SU(6) symmetry, we take the diagonal component of the spin-moment operator

[TABLE]

whose largest eigenvalue is normalized to 1/2 as in the case of SU(2). The Fourier component of spin-moment at $Q=(\pi,\pi)$ corresponds to the AF order parameter $m=\frac{1}{L^{2}}\sum_{r}m_{r}(-1)^{r}$. Due to the finite sizes of QMC simulations, there is no spontaneous symmetry breaking, and we measure spin structural factor defined as the equal-time spin-spin correlation,

[TABLE]

To describe the VBS order, we define the kinetic dimer operator,

[TABLE]

where $\hat{e}=\hat{x},\hat{y}$. Then the VBS order parameters $d_{\hat{x}}$ and $d_{\hat{y}}$ are defined as the Fourier components at $(\pi,0)$ and $(0,\pi)$, respectively,

[TABLE]

Again, we directly measure the structure factor of dimer-dimer correlation

[TABLE]

In large $U$ limit (Heisenberg limit), the kinetic dimer order is equivalent to the spin-Peierls VBS defined as

[TABLE]

for the SU($N$) Heisenberg models Read and Sachdev (1989a); *Read_NPB_1989; *Read_PRB_1990; Marston and Affleck (1989) through the 2nd order perturbation theory. (For finite $U$, charge fluctuations may cause the inequivalence of these two kinds of VBS definitions.) Based on the AF and VBS structure factors $S_{mQ}$ and $S_{vbs}$, we further denote

[TABLE]

In the presence of long-range ordering of AF and VBS, $m^{2}_{Q}$ and $O_{vbs}^{2}$ exhibit non-vanishing values in the thermodynamic limit $L\to\infty$, respectively.

III QMC simulation results

We first present the single particle gap $\Delta_{c}$ extracted from the slope of $\ln G(\tau,k)$, where $G(\tau,k)$ is the single particle Green’s function defined as $G(\tau,k)=-\frac{1}{L^{2}}\sum_{rr^{\prime}}\langle\mathcal{T}_{\tau}c_{r\alpha}(\tau)c_{r^{\prime}\alpha}^{\dagger}\rangle\mathrm{e}^{ik\cdot(r-r^{\prime})}$. The momentum $k$ is taken at $(\pi,0)$ on the Fermi surface. The results for $L=14$ are shown in Fig. 1. When $\Delta_{c}\gtrsim 1$, it shows very little size dependence, because it describes the local charge fluctuations with a very short charge coherence length estimated as $\xi_{c}\sim t/\Delta_{c}\lesssim 1$. Hence, the results at $L=14$ already can be taken as the thermodynamic limit. An interesting observation is the nearly linear dependence of $\Delta_{c}$ on $U$ at $U>U^{*}\approx 9t$, whose slope is very close to 1/2, indicating the characteristic feature of the Mott insulator. This is consistent with that in the atomic limit, i.e., $t/U\rightarrow 0$, which is simply $U/2$, the energy cost by adding or removing an electron on the half-filled Mott insulating background. On the other hand, in the regime $U<U^{*}$, $\Delta_{c}$ keeps at very small values, which is consistent with the AF insulators based on the Fermi surface nesting as in a Slater insulator. Therefore, we take $U^{*}$ as a crossover from the Slater to Mott regimes since no symmetry breaking occurs.

Next we consider the AF and VBS orderings near the AF-VBS transition by performing the finite size extrapolation of $m^{2}_{Q}$ and $O_{vbs}^{2}$ at $L\to\infty$. Without a precise knowledge of finite size effects in prior, we have tried different fitting functions. It turns out that the usual polynomial (neither square nor cubic) functions of $1/L$ used in Ref. [Neuberger and Ziman, 1989; Sandvik, 1997] fail to fit the data. Instead, a simple (non-integer) power law function

[TABLE]

works pretty well, where $a$, $b$, and $c$ are fitting parameters. We suspect that this is due to strong quantum fluctuations in the interaction parameter regime (from $U/t=10$ to $16$) near the quantum phase transition. The complex excitations would significantly change the finite size effect. The finite size scalings are shown in Fig.  2. The extrapolated values of the AF and VBS order parameters are plotted in Fig. 1. Both of them drop to zero as the interaction parameter approaches a small regime of $13<U/t<13.5$ from the opposite directions. We have also checked the evolutions of total and kinetic energies. Neither of them shows an obvious discontinuous behavior, which suggests continuous transitions.

Certainly, there exist a few possibilities regarding to the nature of these two ordering transitions: $a$) There is only a direct 2nd order phase transition as in the framework of the deconfined criticality, i.e., the two orderings share the same critical value of $U$; $b$) they exhibit two separate but very close 2nd order phase transitions with a quantum disordered phase in between; $c$) the same as in $b$) but with a small coexistence regime of both orders; $d$) a weak 1st order transition between them. To further address the nature of these transitions, we perform the following scaling analysis.

We first consider the scalings from the AF side. The following definition of the AF correlation length is employed based on spin-spin correlations Sandvik (2010b)

[TABLE]

where $Q=(\pi,\pi)$ is the ordering wavevector, and $q$ is a small deviation from $Q$ chosen as $q=(2\pi/L,0)$. In Fig. 3($a$), we plot $\xi_{\rm AF}/L$ versus $U$ at different values of $L$, and find that they cross at $U_{\rm AF,c}/t=13.3\pm 0.05$, which is taken as the transition point for the AF order.

Based on the critical value of $U_{\rm AF,c}$, we further perform the data collapse as plotted in Fig. 3($b$) according to the scaling function of $\xi_{\rm AF}$,

[TABLE]

where the exponent of the divergence of correlation length is determined to be $\nu_{\rm AF}=0.60\pm 0.02$. Such an exquisite scaling behavior is a strong hint to a continuous phase transition.

At a quantum critical point, the two-point correlation function in $d+1$ dimensions is expected to be algebraic decay as

[TABLE]

where $z$ is the dynamic critical exponent and $\eta_{\rm AF}$ is the anomalous dimension Sondhi et al. (1997). After the Fourier transformation, the structure factor at finite size $L$ scales as

[TABLE]

at large enough values of $L$. In Fig. 3 ($c$), $m^{2}_{Q}$ is plotted versus $L$ on a log-log coordinate around $U_{\rm AF,c}$, which exhibits a good linear behavior up to $L=24$. From their slopes, $z+\eta_{\rm AF}=1.44\pm 0.03$ is found. In our simulations, $z$ is difficult to determine accurately since it requires the time evolutions of two particle Green’s functions which, however, tend to be gapless at the critical point. If we adopt the $z=1$ directly following the prediction of the deconfined critical theorySenthil (2004); *Senthil_PRB_2004; *Levin_PRB_2004, we arrive at the anomalous dimension $\eta_{\rm AF}=0.44\pm 0.03$.

Based on the above analysis, the AF transition exhibits quite clear evidence of a 2nd order phase transition. Here we summarize the critical value of $U_{\rm AF,c}$ and the two critical exponents for the AF transition as

[TABLE]

As for the correlations for the VBS orderings, unfortunately, it is difficult to obtain high quality data for the scaling of the VBS correlation length to determine the critical value of $U_{\rm VBS,c}$ and compare it with $U_{\rm AF,c}$. Nevertheless, since the VBS transition is very close to the AF one, we present the log-log plot of $O^{2}_{vbs}(L)\sim(1/L)^{d+z-2+\eta_{\rm VBS}}$ for $U$ in the vicinity of $U_{\rm AF,c}$. The fitting of the slopes yields $z+\eta_{\rm VBS}=1.98\pm 0.01$, which corresponds to

[TABLE]

The above anomalous dimensions $\eta_{\rm AF}$ and $\eta_{\rm VBS}$ are different from those obtained from the SU(6) $J$-$Q$ model Kaul and Sandvik (2012) indicating they indeed belong to different universality classes. Our case is based on the fermionic SU(6) Hubbard model, and in its Mott insulating state each site is in the self-conjugate representation, i.e., 3 fermions per site, while the $J$-$Q$ model is equivalent to the non-compact CP5 model in which neighboring states belong to the SU(6) fundamental and anti-fundamental representations.

In order to identify the type of the VBS order, we plot the histograms of the VBS configurations during QMC simulations, as shown in Fig. 4(b). Deep inside the VBS state (for large $U$ and large $L$), the histogram shows larger weights at $\text{arg}(d_{x}+id_{y})=0$ than $\pi/4$, indicating that the VBS belongs to the columnar type. However, near the phase boundary, the histograms are difficult to tell which type of the VBS is, which in fact is consistent with an emergent $U(1)$ symmetry in the framework of deconfined criticality. Sandvik (2007)

IV summary and discussion

In summary, we have performed a large scale projector QMC simulations on the half-filled SU(6) Hubbard model in the square lattice. As $U$ increases, we have found a crossover at $U^{*}/t\approx 9$ from the Slater-AF to the Mott-AF insulators. As $U$ further increases, a (signature of) continuous phase transition at $U_{c}/t=13.3\pm 0.05$ from the Mott-AF to Mott-VBS states.

Several remarks of these numerical observations are given as follows: (1)The finite size extrapolations in this work are based on the power law fitting in Eq. 9, which is different from most studies, especially the SU(2) Heisenberg model where the cubic-order polynomial works very well Neuberger and Ziman (1989); Sandvik_PRB_1997. The difference may be rooted in the stronger quantum fluctuations from the higher symmetry group $SU(N)$. A full understanding requires more sophisticated knowledge of the excitation properties of the SU(N) Hubbard which is left in future studies. (2) We have also tried to obtain the universal plots of $\xi_{\rm VBS}/L$ and $L^{z+\eta_{\rm VBS}}O_{vbs}^{2}$ versus $|U-U_{c}|L^{1/\nu_{\rm VBS}}$ but failed. Of course, this may be caused by insufficient lattice sizes up to $L=24$ in our simulations. However, another possible reason may be the recently proposed two-length scaling hypothesis Shao et al. (2016), which requires very large lattice sizes difficult to reach for the determinant QMC simulations for fermions. (3) By symmetry analysis, the observed AF-VBS phase transition belongs to a broader universality class governed by a $U(N)/[U(m)\otimes U(N-m)]$ nonlinear sigma model beyond the CPN-1* model corresponding to $m=1$. The full understanding, e.g. critical exponents, of the $m>1$ models calls for more elaborated theoretical efforts via e.g. $1/N$ expansionWang et al. (2019), or, renormalization group analysisDas (2018) in the future.

V acknowledgement

D. W. is supported by NSFC under grant No. 11874205. L. W. is supported by supported by the Ministry of Science and Technology of China under the Grant No. 2016YFA0302400 and the Strategic Priority Research Program of Chinese Academy of Sciences Grant No. XDB28000000. C. W. is supported by AFOSR FA9550-14-1-0168. The numerical calculations were performed on Tianhe-II platform at the National Supercomputer Center in Guangzhou.