Designer Monte Carlo simulation for the Gross-Neveu-Yukawa transition

In this paper, we study the quantum criticality of Dirac fermions via large-scale numerical simulations, focusing on the Gross-Neveu-Yukawa chiral-Ising quantum critical point (QCP) with critical bosonic modes coupled with Dirac fermions. We show that finite-size effects at this QCP can be efficiently minimized via model design, which maximizes the ultraviolet cutoff and at the same time places the bare control parameters closer to the nontrivial fixed point to better expose the critical region. Combined with the efficient self-learning quantum Monte Carlo algorithm, which enables a nonlocal update of the bosonic field, we find that moderately large system size (up to $16\times 16$) is already sufficient to produce robust scaling behavior and critical exponents. The conductance of free Dirac fermions is also calculated, and its frequency dependence is found to be consistent with the scaling behavior predicted by the conformal field theory. The methods and model-design principles developed for this study can be generalized to other fermionic QCPs, and thus provide a promising direction for controlled studies of strongly correlated itinerant systems.

Figure 1

The fermion tight-binding models are shown in (a) and (b), with (a) only containing the nearest-neighbor hopping $t_1$ (Ref. [21]) and (b) containing both the nearest-neighbor hopping $t_1$ and the third-nearest-neighbor hopping between different sublattices $t_2$ (this work). The green dashed lines in the right panels mark the linear dispersion region, out of which the band structure starts to deviate from the Dirac linear dispersion. The $t_2$ hopping in (b) greatly increases the area of the linear dispersion region. The black arrows in the left panel represent the phase of the complex hopping strength $\theta =\frac{\pi }{4}$, and we set $t_1=2$ and $t_2=\frac{t_1}{27}$. The value of $t_2$ is obtained from the optimization process of trying to enlarge the linear dispersion area in the BZ. The lattice model is depicted in (c), where the coupling constants for bosons are $J_2=-\frac{J_1}{8},J_3=\frac{J_1}{63},J_4=-\frac{J_1}{896}$. The values of boson couplings are also obtained from the parameter optimization process. The bosons and fermions are coupled together via the ${\lambda }_{i,j}=\pm 1$ term shown in Eq. (2), where two next-nearest-neighbor fermions are coupled together with the boson between the two fermion sites. The sign of the coupling constant ${\lambda }_{i,j}$ alternates from one plaque to another, as marked in (c). Overall, one unit cell of this model contains two fermion sites and two boson sites. Parts (d) and (e) show dispersions of fermions and bosons, respectively (see Appendix pp3 for a detailed explanation of how to obtain these dispersions). Here, we presented both the bare dispersions (ignoring interactions) and the renormalized ones at the chiral Ising transition via QMC simulations (with $L=14$). In panel (d), it is easy to notice that the model shown in panel (b) greatly increases the linear-dispersion region.

Figure 2

(a) The Binder ratio $U(L,m)$ for different system sizes $L$ across the transition of $m$. Curves for two consecutive system sizes $L$ and $L+2$ cross at the finite-size transition point $m_c\left(L\right)$. (b) The correlation ratio $R(L,m)$ for different system sizes $L$ across the transition of $m$. Curves for two consecutive sizes $L$ and $L+2$ cross at the finite-size transition point $m_c\left(L\right)$. The finite-size analysis for these crossing points is shown in Fig. 3.

Figure 3

(a) Crossing-point analysis for the Binder and correlation ratios obtained in Fig. 2. The two $m_c\left(L\right)$ curves extrapolate to the same critical point $m_c^{*}=-0.062\left(3\right)$ at the thermodynamic limit (note that we set the two curves to meet at a single value of $m_c^{*}$ in the fitting), and the two powers in the extrapolation according to Eq. (6) consistently reveal $1/{\nu }_{\text{GN}}+\omega =1.8\left(1\right)$, with different nonuniversal coefficients $a$. (b) Crossing-point analysis of the Binder ratio $U_c\left(L\right)$ and the correlation ratio $R_c\left(L\right)$ according to Eqs. (7) and (8). The obtained correction exponent $\omega =0.8\left(1\right)$. Combining the exponents from (a) and (b), we find $1/{\nu }_{\text{GN}}=1.0\left(1\right)$. (c) Data collapsing for the order parameter measured at different $m$ and system sizes $L$. Here, we utilized the exponent ${\nu }_{\text{GN}}$ and the critical point $m_c^{*}$ obtained above, and the system sizes are $L=6$, 8, 10, 12, 14, and 16. A good collapse is achieved with dynamic critical exponent $z=1$ and the anomalous dimension ${\eta }_{\varphi }=0.59\left(2\right)$.

Figure 4

The imaginary-time decay of the fermion Green's function at the chiral Ising QCP. From the power-law decay of $G^f\left(\tau \right)\sim 1/{\tau }^{2+{\eta }_{\psi }}$ we obtain the anomalous dimension of the Dirac fermions ${\eta }_{\psi }=0.05\left(2\right)$.

Figure 5

The frequency dependence of the current-current correlation with $L=50$ and $T=1/200$ for free Dirac fermions, where the red (blue) data points are obtained from the designer model with $t_2=1/27$ (the conventional model with $t_2=0$). The black solid line is the fitting to the conformal theory prediction [Eq. (11)], which gives rise to ${\sigma }_{\infty }=\frac{1}{4}$ at the limit of ${\omega }_n\gg T$. The inset shows the fitted ${\sigma }_{\infty }$ as a function of the system size $L$. As $L$ increases, ${\sigma }_{\infty }$ from the designer model converges to the exact value (the dashed line) much faster than that of the conventional model.

Figure 6

Fermion and boson Green's functions and their fitting according to Eq. (C1) from which the finite-size gap $\mathrm{\Delta }\left(\mathbf{k}\right)$ at different momenta are obtained, These gaps, prepared on all the momenta along the high-symmetry-path, give rise to the dispersions $\omega \left(\mathbf{k}\right)$ shown in Figs. 1 and 1.