Revealing Fermionic Quantum Criticality from New Monte Carlo Techniques

TL;DR

This review discusses recent advances in numerical methods, especially quantum Monte Carlo techniques, that have enabled the study of fermionic quantum criticality in complex condensed matter systems, providing new insights and potential experimental guidance.

Contribution

It highlights novel model designs and algorithm improvements that allow unbiased large-scale simulations near quantum critical points, advancing understanding of fermionic quantum critical phenomena.

Findings

Development of improved Monte Carlo algorithms for fermionic systems

Large-scale numerical simulations revealing properties of quantum critical points

Potential experimental implications for superconductors and cold atom systems

Abstract

This review summarizes recent developments in the study of fermionic quantum criticality, focusing on new progress in numerical methodologies, especially quantum Monte Carlo methods, and insights that emerged from recently large-scale numerical simulations. Quantum critical phenomena in fermionic systems have attracted decades of extensive research efforts, partially lured by their exotic properties and potential technology applications and partially awaked by the profound and universal fundamental principles that govern these quantum critical systems. Due to the complex and non-perturbative nature, these systems belong to the most difficult and challenging problems in the study of modern condensed matter physics, and many important fundamental problems remain open. Recently, new developments in model design and algorithm improvements enabled unbiased large-scale numerical solutions to be achieved in the close vicinity of these quantum critical points, which paves a new pathway towards achieving controlled conclusions through combined efforts of theoretical and numerical studies, as well as possible theoretical guidance for experiments in heavy-fermion compounds, Cu-based and Fe-based superconductors, ultra-cold fermionic atomic gas, twisted graphene layers, etc., where signatures of fermionic quantum criticality exist.

Figures

Equations

$$S=S_{\text{fermion}}+S_{\text{boson}}+S_{\text{coupling}}$$

$$S_{\text{fermion}}=\int d\mathbf{r}dt\sum _{\alpha }\left(c_{\alpha }^{†}i{\partial }_t{c}_{\alpha }-\frac{{\hslash }^2}{2m_e}\nabla c_{\alpha }^{†}\nabla c_{\alpha }\right),$$

$$S_{\text{boson}}=\frac{1}{2}\int d\mathbf{r}dt\left[({\partial }_t\varphi {)}^2-(\nabla \varphi {)}^2-m{\varphi }^2-\frac{u}{2}{\varphi }^4\right].$$

$$S_{\text{coupling}}=g\int d\mathbf{r}dt\varphi s_z,$$

$${\chi }_0(\mathbf{q},\omega )=\frac{1}{{\omega }^2-q^2-{\xi }_c^{-2}},$$

$$\chi (\mathbf{q},\omega )=\frac{{\chi }_0}{a_{0}i|\omega |/q+a_1{\omega }^2-a_2(q^2+{\xi }_c^{-2})},$$

$$H=H_{\text{fermion}}+H_{\text{boson}}+H_{\text{coupling}},$$

$$H_{\text{fermion}}$$

$$H_{\text{boson}}$$

$$H_{\text{coupling}}$$

$$H_{\text{fermion}}$$

$$H_{\text{boson}}$$

$$H_{\text{coupling}}$$

$$e^{\Delta \tau h{s}_i^x}=\cosh (\Delta \tau h)1+\sinh (\Delta \tau h)s_i^x,$$

$$⟨S_z^{\prime }|e^{\Delta \tau h{s}_i^x}|S_z⟩=\Lambda e^{\gamma s_z^{\prime }{s}_z}.$$

$$⟨S_z|e^{\Delta \tau h{s}_i^x}|S_z⟩$$

$$⟨-S_z|e^{\Delta \tau h{s}_i^x}|S_z⟩$$

$$\gamma$$

$${\Lambda }^2$$

$$Z$$

$${{\mathcal{W}}_{\mathcal{C}}}^{\text{boson}}=$$

$$(\prod _i\prod _{⟨\tau ,{\tau }^{\prime }⟩}\Lambda e^{\gamma {\sum }_i{s}_{i,\tau }^z{s}_{i,{\tau }^{\prime }}^z})$$

$${{\mathcal{W}}_{\mathcal{C}}}^{\text{fermion}}=\det [1+\prod _{\tau }(e^{-\Delta \tau T}\prod _i{e}^{V\{s_{i,\tau }^z\}})]$$

$$H_f=\sum _k[ϵ(k)-\mu ]c_{ka}^{†}{c}_{ka};$$

$$H_{fb}=\lambda \sum _{k{k}^{\prime }}{c}_{ka}^{†}{M}_{ab}{c}_{k^{\prime }b}{\varphi }_{k-k^{\prime }}.$$

$${\varphi }_k=\frac{1}{N}\sum _i{\varphi }_i{e}^{-i\mathbf{k}\cdot {\mathbf{r}}_i}.$$

$$H_{\text{fermion}}=$$

$$-t_3\sum _{⟨⟨⟨ij⟩⟩⟩,\lambda ,\sigma }{c}_{i,\lambda ,\sigma }^{†}{c}_{j,\lambda ,\sigma }+h.c.-\mu \sum _{i,\lambda ,\sigma }{n}_{i,\lambda ,\sigma },$$

$$Z_{{\mathbf{k}}_F}\approx \frac{1}{1-\frac{\text{Im}\Sigma ({\mathbf{k}}_F,i{\omega }_0)}{{\omega }_0}},$$

$$\chi (h,T,\mathbf{q},{\omega }_n)=\frac{1}{L^2}{\int }_0^{\beta }d\tau \sum _{ij}{e}^{i{\omega }_n\tau -i{\mathbf{qr}}_{ij}}⟨s_i^z(\tau )s_j^z(0)⟩.$$

$$\chi (h,T,\mathbf{q},{\omega }_n)=\frac{1}{c_t{T}^2+c_h\left|h-h_c\right|+c_q{q}^2+c_{\omega }{\omega }^2+\Delta (\mathbf{q},{\omega }_n)}.$$

$$\Delta (\mathbf{q},{\omega }_n)=c_{\mathrm{HM}}\frac{|{\omega }_n|}{\sqrt{{\omega }_n^2+(v_{f}q{)}^2}}$$

Full text

1. I Introduction

2. II Theoretical setting

3. II.1 Model

4. II.2 the Hertz-Millis-Moriya theory and beyond

5. III Numerical methodologies

6. III.1 Designer-Hamiltonians

7. III.2 Determinantal quantum Monte Carlo

8. III.3 Self-learning Monte Carlo, EMUS, etc

9. IV Results

10. IV.1 Phase Diagrams

11. IV.1.1 FM-QCP

12. IV.1.2 AFM-QCP with $3\mathbf{Q}=\Gamma$

13. IV.1.3 AFM-QCP with $2\mathbf{Q}=\Gamma$

14. IV.2 Non-Fermi liquid behavior

15. IV.3 Quantum critical scaling analysis

16. IV.3.1 FM-QCP

17. IV.3.2 AFM-QCP with $3\mathbf{Q}=\Gamma$

18. IV.3.3 AFM-QCP with $2\mathbf{Q}=\Gamma$

19. IV.4 Superconductivity

20. V Quantum criticalities of fermions coupled to phonons, $\mathbb{Z}_{2}$ or $U(1)$ gauge fields

21. V.1 Dirac fermions coupled with phonons

22. V.2 Dirac fermions coupled with $\mathbb{Z}_{2}$ or $U(1)$ gauge fileds

23. VI Discussion and prospectives

I Introduction

In the past few decades, the discovery of novel quantum materials and many-body states whose unconventional properties do not fitted into traditional quantum many-body paradigms such as the Fermi liquid theory of metals and the Landau-Ginzburg-Wilson framework for phases and phase transitions, urged for new theoretical principles and insights to better understand these new quantum states and to guide the search for novel quantum materials.

In conventional materials, such as a noble metal, the net effect of interactions is limited to modify certain quantitative properties of electrons, such as the effective mass, which results in the formation of the so-called “quasiparticle". In conventional materials, quasiplarticles at the Fermi surface has infinite lifetime at zero temperature, which implies that these quasiparticles are effectively free from scatterings and thus behave as non-interacting particles. In the field theory description, the infinite lifetime leads to a finite quasiparticle residue, and as long as the quasiparticle residue remains finite, the non-interacting picture remains asymptotically accurate at low-energy and thus a complicated many-body quantum system is now transformed into a non-interacting system, whose physical properties can be easily understood and predicted. This approach, known as the Fermi liquid (FL) theory, has been proved to be highly effective in the study of conventional materials and it is the foundation for our understanding about many-body electronic systems.

However, in the past few decades, it is found that the Fermi-liquid theory fails in many quantum materials, such as the Cu- and Fe-based superconductors, heavy-fermion metal, and transition-metal alloys on the brink of magnetic orders Stewart (2001); Custers et al. (2003); Löhneysen et al. (2007). These materials support the so-called ‘non-Fermi-liquid (nFL)" or “strange metal" phase, in which the quasiparticle residue vanishes as the temperature is reduced to zero, instead of remaining finite as in a FL. The vanishing quasiparticle residue violates the key assumption of the FL theory, and thus results in novel quantum properties in direct contrast to FL predictions. In the studies of nFLs, one of the central focuses is systems with itinerant quantum critical points (IQCPs), which have attracted extensive research efforts for nearly half century dating back to the celebrated Hertz-Millis-Moriya (HMM) framework Hertz (1976); Millis (1993); Moriya (1985). In the study of correlated quantum materials, quantum criticality in itinerant electron systems is of great importance and interests Hertz (1976); Millis (1993); Moriya (1985); Stewart (2001); Chubukov et al. (2004); Belitz et al. (2005); Löhneysen et al. (2007); Chubukov and Maslov (2009), and it plays an important role in the understanding of anomalous transport, strange metal and nFL behaviors Metzner et al. (2003); Senthil (2008); Holder and Metzner (2015); Metlitski et al. (2015); Xu et al. (2017a) in various quantum materials, such as heavy-fermion materials Custers et al. (2003); Steppke et al. (2013), Cu- and Fe-based high-temperature superconductors Zhang et al. (2016); Liu et al. (2016); Gu et al. (2017); Wang et al. (2018a), the recently discovered pressure-driven quantum critical point (QCP) between magnetic order and superconductivity in transition-metal monopnictides, CrAs Wu et al. (2014), MnP Cheng et al. (2015), CrAs1-xPx Matsuda et al. (2018) and other Cr/Mn-3d electron systems Cheng and Luo (2017) and the more recent discoveries in twisted angle graphene heterostructures Cao et al. (2018, 2019); Shen et al. (2019). However, after decades of extensive efforts Hertz (1976); Millis (1993); Moriya (1985); Stewart (2001); Metzner et al. (2003); Abanov et al. (2003); Abanov and Chubukov (2004); Chubukov et al. (2004); Belitz et al. (2005); Löhneysen et al. (2007); Chubukov and Maslov (2009); Metlitski and Sachdev (2010a, b); Sur and Lee (2015, 2016); Schlief et al. (2017); Lunts et al. (2017); Lee (2018); Schlief et al. (2018), itinerant quantum criticality still remains among the most challenging subjects in condensed matter physics due to its nonperturbative nature, and many fundamental questions and puzzles still remain open.

It is worthwhile to point out that similar to classical critical phenomena and thermal phase transitions, the theoretical understanding of quantum criticality and electronic quantum criticality requires combined efforts from both analytic theory and large-scale numerical simulations. In the past few decades, a lot of great progress have been made on the theory side Hertz (1976); Millis (1993); Moriya (1985); Stewart (2001); Metzner et al. (2003); Abanov et al. (2003); Abanov and Chubukov (2004); Chubukov et al. (2004); Belitz et al. (2005); Löhneysen et al. (2007); Chubukov and Maslov (2009); Metlitski and Sachdev (2010a, b); Sur and Lee (2015, 2016); Schlief et al. (2017); Lunts et al. (2017); Lee (2018); Schlief et al. (2018). However, numerically testing and verifying these predictions remain a highly challenging task for IQCPs in 2D or higher dimensions, mainly due to limitations in all available numerical techniques. The absence of reliable and unbiased numerical verification makes it difficulty to verify key assumptions adopted in theoretical studies and thus making it highly challenging to refine our theoretical knowledge about IQCPs.

Recently, thanks to major new progress in quantum Monte Carlo (QMC) techniques, accurate and unbiased numerical solutions becomes achievable for various IQCPs. It is found that designer models of fermions coupled to critical bosonic fields offer a pathway to access a wide variety of IQPCs, and in the same time avoid the notorious sign-problem large-scale QMC simulations possible Berg et al. (2012). Utilizing this approach, a variety of fermionic quantum critical points are studied, with critical/soft fluctuations such as antiferromagnetic, nematic or phononic fluctuations Berg et al. (2012); Schattner et al. (2016a, b); Li et al. (2016). In the study of unconventional superconductors, such as Cu- or Fe- based high temperature superconductors, these soft low-energy fluctuations have been experimentally observed. However, whether these fluctuations are the main driving force for the formation of high temperature superconductivity remains an important open questions, and these QMC simulations are expected to shed new lights on this long-standing question. Remarkably, all these QMC simulations shares one common observation: as one turns on the coupling between fermions and critical fluctuations, a superconducting dome always arises and covers the QCP regardless of the symmetry breaking patterns, which suggests that there might exist a deep and fundamental connection between high temperature superconductivity and these quantum critical fluctuations.

In this review, we focus on a different aspect of IQCPs. Instead of exploring superconducting domes induced by IQCPs, our objective is to understand the IQCPs themselves, i.e. to identify the critical theory and its universal scaling behavior and to explore the non-Fermi liquids induced by the critical flucuations. For this purpose, we need to access the close vicinity of a QCP, which requires suppressing the superconducting dome to expose the QCP. Experimentally, superconductivity can be in principle suppressed by a strong magnetic field, but for QMC studies, this approach is not the optimum pathway, because handling gauge fields is numerical challenging. Instead, in QMC, we found that a more convenient approach is to use designer Hamiltonians, especially using bilayer models that contain two identical copies of fermions Xu et al. (2017a). Although the mechanism is not yet fully understand, in all bilayer models that we have examined, this setup suppresses the superconducting dome down to extremely low temperature unaccessible to QMC simulations. In addition to suppressing superconductivity, to obtain accurate critical phenomena and scaling exponents, large system sizes, in both spatial and temporal directions, are necessary, which requires new numerical techniques to access larger systems and lower temperature.

Along this line of efforts, new designer models and QMC simulation techniques are developed to examine in detail the properties of IQCPs. For example, to study strange metal behaviors near the ferromagnetic quantum phase transition, we developed the self-learning quantum Monte Carlo method Liu et al. (2017a); Xu et al. (2017b) to solve this problem and observed clear signatures of non-Fermi-liquid behaviors near the IQCP Xu et al. (2017a). More remarkably, quantum critical analysis indicates it is a new critical point with an anomalous scaling exponent Xu et al. (2017a). These results not only sharpen the theoretical understanding about IQCP, but also offer important theoretical guidance to experimental investigations on various strongly correlated metals. The similar technique is also applicable for antiferromagnetic QCPs, via coupling antiferromagnetic fluctuations with itinerant fermions Berg et al. (2012); Liu et al. (2018, 2019, 2018), which can be used to investigate the fundamental phenomena such as antiferromagnetic fluctuations mediated charge density wave, psesudo-gap phase and unconventional superconductivity. All these questions can now be hopefully addressed with unbiased numerical simulations. Moreover, with this technique, one can also study the effect of compact gauge fields coupled to a Fermi surface, such that the interdisciplinary research between condensed matter and high energy physics on the subject of quantum spin liquid with emergent gauge field and the confined and deconfined state of matter can be explored Xu et al. (2019). All these fast developments, by and large, are calling for a topical review such as this in front of the readers, to summarize the key developments and point out future directions.

The content of the review is based on several recent numerical works. The discussion about ferromagnetic itinerant QCP and the associated nFL quantum critical region (QCR), is mainly based on Ref. Xu et al. (2017a). The discussion on antiferromagnetic IQCP, in realization of 2d triangle and square lattices and their similarity and difference, are mainly based on Refs. Liu et al. (2018, 2018). The developments in the QMC methodologies, including basic outline of the determinantal QMC (DQMC) for fermion coupled to critical boson problems are based on Refs. Xu (2017); Liu (2018), and the framework of self-learning Monte Carlo (SLMC in short) methods are based on the Refs. Liu et al. (2017a); Xu et al. (2017b); Chen et al. (2018) as well as the momentum-space based elective momentum ultra-size quantum Monte Carlo method (EMUS in short) are based on the Refs. Liu et al. (2019, 2018). The other recent developments, including the application of SLMC on electron-phonon coupled Holstein-type of problems Chen et al. (2018, 2019) and the question of U(1) gauge fields coupled to matter fields in the form of Dirac fermions and the deconfinement-to-confinement transition discoveried in the such model simulations Xu et al. (2019), are discussed as well.

Because this review mainly focuses on QMC studies, we will also provide a pedagoical QMC code package developed over the past few years by some of us, which contains the lattice construction, local update of conventional DQMC and SLMC steps with learning of effective model and the cumulative update scheme afterwards, as well as the EMUS construction of the momentum patches. The package can be found in the GitHub repository Xu et al. (2018) and the corresponding manual can be found in the Ref. Liu (2018). Readers who are interested in understanding the numerical details of these new Monte Carlo techniques and would like to reproduce our results and carry out further investigations, are strongly encouraged to download, install and run the package and modify it to his/her purposes. And the developers are more than willing to offer help if questions/problems/comments arise.

II Theoretical setting

The description of quantum phase transitions and quantum criticality in metallic systems was pioneered by the Hertz-Millis-Moriya theory Hertz (1976); Millis (1993); Moriya (1985). However, it was soon found that Hertz-Millis-Moriya (HMM) formalism may not be adequate, and for a number of important situations, higher order terms ignored in the Hertz-Millis-Moriya are found to result in unstable RG flows away from the Hertz-Millis-Moriya fixed point, such as the ferromagentic, antiferromagnetic and nematic quantum critical points in (2+1)D. And more advanced theoretical frameworks, aiming at capturing these quantum critical points and more importantly the non-Fermi-liquid behavior associated with them, are developed Metzner et al. (2003); Senthil (2008); Metzner et al. (2003); Abanov et al. (2003); Abanov and Chubukov (2004); Chubukov et al. (2004); Belitz et al. (2005); Chubukov and Maslov (2009); Metlitski and Sachdev (2010a, b); Dalidovich and Lee (2013); Lee (2009a); Sur and Lee (2016); Schlief et al. (2017); Lee (2018); Schlief et al. (2018) over the years. Although this review mainly focuses on the numerical end, we will briefly introduce the theoretical setting and describes its fundamental challenges in simple terms.

II.1 Model

For strongly-correlated systems, although theoretical descriptions often start from fermions with direct fermion-fermion interactions (e.g. four-fermion interactions), to describe a quantum phase transition as often used in the Hertz-Millis-Moriya approach Hertz (1976); Millis (1993); Moriya (1985), the most efficient way is to use an effective model where directly fermion-fermion interactions are replaced by certain boson-mediated interactions. These bosons can be thought of as the auxiliary field in the Hubbard-Stratonovich transformation. Here, we choose the auxiliary filed to share the same symmetry as the order parameter of the phase transition, i.e. an order-parameter field. Thus, the expectation value of the boson operator serves as the order parameter and these bosonic fields describes soft critical fluctuations in the vicinity of the QCP. It is worthwhile to emphasize that this effective model don’t fully replicate the details of the original model with fermion-fermion interactions. In particular, contributions from other quantum fluctuations beyond the order parameter field is ignored in the effective theory. However, near the QCP, universal properties, which are independent of microscopic details, are expected to survive.

Interestingly, this idea of using an effective model also plays an important role in QMC studies. For QMC simulations, directly simulating four-fermion interactions often turns to be challenging due to the arising of the sign problem. However, the effective model used in the Hertz-Millis-Moriya approach (with boson-mediated interactions) turns out to be much more QMC friendly and the sign problem can often be avoided. Thus, although a generic and efficient way to avoid the sign problem is still absent for the Hubbard or Hubbard-like models, the sign problem can be easily avoided in many of these boson-mediated effective model. As mentioned above, these effective models are precisely the setup utilized in most of the theoretical studies and they are expected to capture universal properties at IQCPs. Thus, these effective models with boson-mediated interactions opens a pathway to study these IQCPs via large scale numerical simulations.

To demonstrate the effective theory used in the study of fermionic quantum criticality, here we show one example based on an Ising ferromagnetic transition in continuous space, but the same construction generically applies to other quantum phase transitions as well as lattice systems. Consider the following action.

[TABLE]

Here $S_{\text{fermion}}$ is the action of a non-interacting fermion gas

[TABLE]

where $c_{\alpha}^{\dagger}$ and $c_{\alpha}$ are the fermion creation and annihilation operators with spin index $\alpha=\uparrow$ or $\downarrow$. $m_{e}$ and $\hbar$ are the effective mass of the fermion and the Planck constant respectively. The second term $S_{\text{boson}}$ is the action of a bosonic field. This boson field respects the same symmetry as the order parameter of quantum phase transition and this action includes all symmetry allowed terms of the boson field, as well as symmetry allowed quantum dynamics. For an Ising transition, we can use the $\phi^{4}$-theory here

[TABLE]

Without fermions, this action describe a bosonic quantum phase transition with the same symmetry breaking pattern. In the example here, a bosonic QCP arises as we decrease the value $m$ to negative. The last term $S_{\text{coupling}}$ is the coupling term between fermions and bosons. The boson field is coupled with a fermion bilinear, which respects the same symmetry as the order parameter. In the example here, the fermion bilinear is z-component of the fermion spin $s_{z}=c_{\alpha}^{\dagger}\sigma^{z}_{\alpha,\beta}c_{\beta}$ where $\sigma^{z}$ is the Pauli matrix and the action is

[TABLE]

where $g$ is the coupling constant.

In this model, the interactions between fermions are mediated by the bosonic field $\phi$. And in return, fermionic fluctuations also renormalize the bosonic part of the action. At $g=0$, when the bosons and fermions decouple, the susceptibility of the bosonic order parameter field (at the tree level) is

[TABLE]

where $\xi_{c}$ is correlation length of bosoinc field, which diverges at the critical point. As we turn on the boson-fermion coupling ($g\neq 0$), a phase transition is still expected as we reduce the value of $m$, although the critical threshold $m_{c}$ is now renormalized. This construction demonstrate a generic approach to realize and study fermionc quantum phase transitions.

II.2 the Hertz-Millis-Moriya theory and beyond

As boson mediates interactions between fermions, the fermion also mediates nontrivial couplings among the bosons. This can be seen by integrating out the fermionic degrees of freedome in the action of Eq. (1), which leads to a bosonic effective theory for a fermionic QCP. In the Hertz-Millis-Moriya theory, this is done within the framework of the random phase approximation (RPA). Namely, only the so-called bubble diagram is included while higher order contributions are ignored.

With this assumption, the boson susceptibility receives significant renormalization

[TABLE]

where $\chi_{0}$, $a_{1}$ and $a_{2}$ are renormalized parameters, whose values goes to unity in the $g\rightarrow 0$ limit where the susceptibility without fermions in the form of Eq. (5) is recovered. Aside from renormalizing parameters at the quantitative level, the most significant contribution from the fermions lies in the new quantum dynamics it introduced to $\phi$, i.e. the term $i|\omega|/q$, known as Landau damping. The coefficient of this term, $a_{0}$, depends on the coupling strength $g$. It is finite for $g\neq 0$ and vanishes in the $g\to 0$ limit. For finite $g$, in the low-energy limit (small $\omega$), this Landau-damping term dominates over the old $\omega^{2}$ term. In the language of renormalization group, this term is more relevant and thus making the old $\omega^{2}$ term subleading and negligible. After dropping the sub-leading $\omega^{2}$ term, a dimension counting indicates that this quantum critical point shall have the dynamic critical exponent $z=3$ Hertz (1976). This Landau damping term and the dynamic critical exponent $z=3$ are generically expected for any phase transition that preserves the (lattice) translational symmetry (e.g. Ising feromagnetic or nematic). For QCPs at which the (lattice) translational symmetry is spontaneously broken (e.g. charge- or spin- density wave states), the Landu damping term takes a slightly different form $\propto i|\omega|$, which often results in $z=2$ Hertz (1976).

As a result, with in the Hertz-Millis-Moriya approximation, the fermions push the effective dimension of a fermionic QCP to $d+z$. For $d\geq 2$, no matter the value of $z$ is $2$ or $3$, $d+z$ will be above or equal to the upper critical dimension $4$, and thus mean-field exponents is expected regardless of microscopic details, which is a key prediction of the Hertz-Millis-Moriya theory. This prediction is highly generic. Although RG flows may lead to nontrivial temperature dependence beyond the RPA approximation Millis (1993), at $T=0$, the RPA prediction and mean-field exponent is always expected within this framework.

However, it is worthwhile to highlight that this strong and universal prediction is based on one important assumption: as we integrate out the fermions, the HMM approach assumes that only leading order one-loop diagram (i.e. the bubble diagram) needs to be included. This assumption can be justified in a Fermi liquid phase. However, at a fermionic QCP, non-Fermi-liquid behaviors often emerge and it is highly unclear whether this one-loop approximation remains valid or not. In particular, as higher order diagrams are taken into account, latter studies indicate that certain fermionic QCPs exhibit a highly nontrivial scaling behavior, with scaling exponent significantly different from mean-field expectations Abanov et al. (2003); Metlitski and Sachdev (2010a, b); Lee (2009a, 2018). In particular, for antiferromagnetic critical point, these anomalous dimensions are computed within the framework of $1/N$ expansion, where $N$ is the number of the fermion hotspots, which will be defined below. Anomalous dimensions are observed even in the first order correction ($\propto 1/N$) Abanov et al. (2003), and later, more careful analysis indicates that the dynamic critical exponent also receives an anomalous dimensions Metlitski and Sachdev (2010b). And self-consistent theory is developed to predict the RG flow of the dynamic critical exponent Lee (2018). Similar breakdown of HMM scaling was also predicted for some other fermionic QCPs Metlitski and Sachdev (2010a).

Although this non-mean-field exponents has been predicted for years, in contrast to the long efforts and rapid progress in analytic theory, progress on the numerical front had a slower pace and no signature of anomalous dimension has been achieved due to system size and other limitations, until very recently. Below, we will summarize some of the latest progress along this line and in particular, we focuses on two questions: (1) whether it is possible for numerical techniques to achieve high enough resolution to distinguish mean-field and non-mean-field exponents and (2) whether we can determine numerically which fermionic QCPs have mean-field exponents (i.e. the Hertz-Millis-Moriya theory remains valid) and which fermionic QCPs exhibit the breakdown of the Hertz-Millis-Moriya theory.

III Numerical methodologies

Before showing the numerical results that answer the second question above, in this section we address the first question on the development of numerical methodologies.

III.1 Designer-Hamiltonians

To put the problem in Eq. (1) on a lattice and solve with quantum Monte Carlo, one would write down the following microscopic Hamiltonian,

[TABLE]

with

[TABLE]

where $H_{\text{fermion}}$ is the free fermion Hamiltonian term; $H_{\text{boson}}$ is the bosonic Hamiltonian realized in the form of a transverse-field quantum Ising model and $H_{\text{coupling}}$ is the couple interaction term between the Ising spin and the fermion. $i$, $j$ and $k$ are the lattice site indices and $\alpha$ contains the set of the orbital and spin and other internal freedom. The free fermion Hamiltonian $H_{\text{fermion}}$ is the hopping term in real space, after Fourier transformation, we can obtain the energy band in momentum space

[TABLE]

The Ising spin term $H_{\text{boson}}$ makes use of the tuning parameter $h_{x}/J$ to obtian an quantum phase transition between the magnetic ordered phase and the disordered paramagnetic phase at zero temperature. The coupling term $H_{\text{coupling}}$ plays a vital role in introducing the critical fluctuations into the fermions and vice versa. The fermion-spin coupling will move the location of the quantum critical point. More importantly, because the fermion-spin coupling is a relevant perturbation in RG, it will change the critical exponents and transfers the quantum critical points into a different universality class.

To be more specific, let us consider a 2D system with three layers of square lattices, including two fermion layers with index $\lambda$ and spin index $\sigma$, and one Ising-spin layer lattice inserted in between the two fermion lattice layers, as schematically shown in Fig. 1. At each lattice site, the Ising spin couples with fermion spins. Such type of specific designer-model, as studied in Refs. Xu et al. (2017a); Liu et al. (2018, 2019, 2018), can be written as

[TABLE]

In Eq. (LABEL:eq:specificformula) the basis of the fermion Hilbert space is $\tau\otimes\sigma\otimes i$. The coupling between fermion and boson is onsite in $H_{\text{coupling}}$. $c_{i}^{\dagger}$ is a four component spinor, $c_{i}^{\dagger}=\begin{bmatrix}c_{i1\uparrow}^{\dagger}&c_{i1\downarrow}^{\dagger}&c_{i2\uparrow}^{\dagger}&c_{i2\downarrow}^{\dagger}\end{bmatrix}$, and $\tau_{0}$ is a $2\times 2$ identity matrix in layer/orbital space, $\sigma_{z}$ is the third Pauli matrix in spin space, respectively. $\sigma_{i\lambda}^{z}=\frac{1}{2}(n_{i\lambda\uparrow}-n_{i\lambda\downarrow})$ is the $z$ component of the fermion spin at orbit $\lambda$ on site $i$. This system preserves a $Z_{2}$ Ising symmetry corresponding to the Ising spin flipping $s_{i}^{z}\rightarrow-s_{i}^{z}$, which will be spontaneously broken in the ordered phase.

As shown in the Fig. 1 (b), the Ising spins develop long-range Ising order in the ordered phase. Deep in the ordered phase, (small $h$ and low $T$), quantum fluctuations of the Ising fields diminish and thus, the Ising spins can be treated as a static and average background, i.e. a mean field. Through the coupling term $H_{\text{coupling}}$, this mean field formed by Ising spins induce the same symmetry breaking pattern to the fermions, and this is how a fermionic phase transition is achieved in these type of models. This construction is highly flexible and generic. Any type of symmetry breaking phase be achieved, as long as proper boson (Ising) fields that break the desired symmetry is utilized. In systems with itinerant fermions, quantum phase transition can be largely classified into two categories depending on whether the ordered phase preserves (lattice) translational symmetry or spontaneously breaks it. Both these two types of quantum phase transition has been studied recently using the QMC approach described here. Examples in the first category includes ferromagnetic Xu et al. (2017a) and nematic Schattner et al. (2016a) phase transitions in metals, while antiferromagnetic quantum phase transition (with itinerant fermions) Gerlach et al. (2017); Liu et al. (2018, 2018) belongs to the second category. As will be shown in latter sections, these symmetry-breaking metallic phases are still FLs with well-defined quasiparticles, after certain Fermi surface reconstructions, e.g. splitted, distorted or folded Fermi surfaces due to the symmetry-breaking and/or the induced Brillouin zone folding.

The real surprise of these designer models [Eq. (LABEL:eq:specificformula)] arises near the quantum critical point in the quantum critical region (QCR), which is also the most challenging regime for theoretical treatments. Inside QCR, the critical bosons cannot be viewed as a mean field any more, and their strong quantum critical fluctuations mediate intensive effective interactions between fermions, which pushes the system into the strongly-correlated regime and is also the origin of non-Fermi-liquid behaviors in these systems. These strong interactions make theoretical description of such a quantum critical regime highly challenging, and exactly solution is not expected. In fact, almost the entire efforts in the past four or five decades, starting from HMM as explained in Sec. II.2, focuses on employing or inventing new and advanced theoretical (renormalization group) machinery, to extract the effective low-energy interactions in these systems. Such effective fermionic and bosonic interactions are not only non-local in space and also non-local in time, rendering the entire problem unperturbative and consequently prohibiting controlled analytical solutions. It is such kind of theoretical challenge, motivated the development of unbiased quantum Monte Carlo methods, to tackle these problems and to bring new insights into our understanding of fermionic quantum criticality. Along the way, better numerical methodologies are inspired and invented with guidance from field theory analysis, from which new results that are bringing us coherent understanding between both numerical and theoretical perspectives, are revealed. In the two subsections followed, we will discuss the quantum Monte Carlo methodologies used and developed during this interesting process.

III.2 Determinantal quantum Monte Carlo

To solve the model in Eq. (LABEL:eq:specificformula), one makes use of the determinant quantum Monte Carlo method. For the purpose of this review, the basic literatures that introduce the method are Refs. Assaad and Evertz (2008); Xu (2017); Liu (2018). DQMC plays a vital role in the study of correlated electron models. As an unbiased method, DQMC can solve sign-problem-free strongly-correlated fermion models with large-scale parallelization to access the thermodynamic limit. Historically, it was first introduced to condensed matter physics systems from high-energy physics community by Blankenbecler, Scalapino and Sugar, so the method is also referred to as BSS method Blankenbecler et al. (1981). DQMC is based on the canonical ensemble within a path integral formalism. Starting from $d$ dimensional quantum grand canonical ensemble $Z=\mathrm{Tr}\{e^{-\beta H}\}$, and use the path integral formalism to tranform the partition function into the $(d+1)$ dimensional classical partition. In the $(d+1)$ dimensional classical partition, the partition function $Z$ is represented as the sum over weights of configuration, which can be seen as the probability of configurations.

In the past, DQMC is mostly used to solve Hubbard or Hubbard-like models with explicit four-fermion interactions Assaad and Evertz (2008); Varney et al. (2009); Meng et al. (2010); He et al. (2016). And to be able to construct the path-integral formalism of the partition function, one needs to introduce auxiliary field to decouple the four-fermion interaction (usually refered as discrete Hubbard-Stratonovich (HS) transformation) Hirsch (1983, 1985). After the decoupling, bosonic HS fields will emerge along with fermion bilinears and it is the space-time configuration of the HS fields that one could use to compute the fermion determinant and the determinant plays the role to the configurational weight, upon which the Markov-chain of Monte Carlo is built.

However, in the new types of problem as shown in Eq. (LABEL:eq:specificformula), the spin-fermion model already acquires the decoupled form, or in another word, the boson fields in Eq. (LABEL:eq:specificformula) already couples to the fermion bilinears, so we do not need to perform the HS transformation. The configurational space in such cases are spanned by the space-time configurations of the boson fields, and for each boson configuration, one can trace out the fermions by evaluating the fermion determinant. Below we first discuss how to construct the bosonic weight of the bare transverse field Ising model $H_{\text{boson}}$.

It is easily to see that from

[TABLE]

one obtains

[TABLE]

For the Ising field in $H_{\text{boson}}$ $S_{z}=\pm 1$, so we can get

[TABLE]

then $\Lambda$ and $\gamma$ can be evaluated as

[TABLE]

Defining $|S_{z}^{N}\rangle$ as the product state of Ising spin at each site in time slice $\tau=N\Delta\tau$, $\beta=M\Delta\tau$ and impose the periodic condition in time direction $|S_{z}^{M}\rangle=|S_{z}^{1}\rangle$. We can rewrite the bare transverse field Ising model partition function

[TABLE]

where $S_{i,\tau}^{z}=\pm 1$ is the Ising field lives in $d+1$ dimensional space after path integral. With the partition function of $H_{\text{boson}}$ obtained, one can readily write down the partition function for the entire system of Eq. (LABEL:eq:specificformula) since here there is no four-fermion interaction terms, as

[TABLE]

where

[TABLE]

is the weight of bare bosonic part in Eq. (17), and

[TABLE]

is the Femion weight which depends on the space and imaginary time distribution of the Ising field $\{S^{z}_{i,\tau}\}$. Eq. (18) can be sampled with Markov process according to the local update of $\{S_{i,\tau}^{z}\}$ spins. And with the DQMC update scheme developed over the years, such evaluation of the ratio of the fermion determinant, with propagation of the $N\times N$ matrices from [math] to $\beta$ Assaad and Evertz (2008), acquires the computational complexity scaling with the system parameters as $O(\beta N^{3})$ in the ideal case, if the autocorrelation between different configurations were not a problem.

III.3 Self-learning Monte Carlo, EMUS, etc

From the above discussion it is obvious that the evaluation of determinant inevitably renders the $O(\beta N^{3})$ scales of operation, which is the computational bottleneck that governs whether DQMC could access larger system sizes and lower temperatures, and what make things worse is the critical slowing down inside the QCR. While physically critical regions are difficult as the fluctuations are the strongest here, but intelligent Monte Carlo update scheme can always help to greatly reduce the autocorrelation among the Monte Carlo configurations and effectively reduces the physical time spent. Successful examples include the Wolff and Swenden-Wang algorithms of the Ising model at classical critical points Swendsen and Wang (1987); Wolff (1989) and loop-updates and stochastic series expansions in the quantum bosonic and spin systems Evertz (2003); Sandvik and Syljuåsen (2003). The successful cluster or non-local Monte Carlo updates all share the merit that the Monte Carlo move actually captures the correct modes in the low-energy landscape of the problem at hand, for example the Wolff cluster resembles the scale-invariance at the critical point and the size of the cluster varies in all length scales. It is even more so in the quantum Monte Carlo simulations in loop-updates or stochastic series expansion. Hence, for the spin-fermion model of our interest, it will be of crucial importance that a non-local update can be designed to respect the effective low-energy physics of the QCR and to overcome the autocorrelation and reduce the physical computational time. This goal, as will be explained here, is partially achieved by the self-learning Monte Carlo method and its subsequential developments such as elective momentum ultra-size Monte Carlo (EMUS) method.

The self-learning Monte Carlo update scheme (SLMC) Liu et al. (2017a, b); Xu et al. (2017b); Nagai et al. (2017); Shen et al. (2018); Chen et al. (2018) is designed to speed up the Monte Carlo simulations, both in classical Liu et al. (2017a, b); Bojesen (2018), quantum few-body Nagai et al. (2017); Shen et al. (2018) and quantum many-body Xu et al. (2017b); Chen et al. (2018, 2019) systems. As shown in Fig. 2, in the context of strongly-correlated fermion systems, the SLMC is implemented in the following manner. One first performs the standard DQMC simulation on the model in Eq. (LABEL:eq:specificformula), and then train an effective boson Hamiltonian that contains long-range two-body interactions both in spatial and temporal directions. The effective Hamiltonian serves as the proper low-energy description of the problem at hand with the fermion degree of freedom integrated out. We then use the effective Hamiltonian to guide the Monte Carlo simulations, i.e., we will perform many sweeps of the effective bosonic model (as the computational cost of updating the boson model is $O(\beta N)$, drastically lower than the update of fermion determinant which scales as $O(\beta N^{3})$), and then evaluate the fermion determinant of the original model in Eq. (18) such that the detailed balance of the global update is satisfied. As shown in our previous works Liu et al. (2017a); Xu et al. (2017b); Liu et al. (2018, 2019, 2018); Chen et al. (2018, 2019), the SLMC can greatly reduce the autocorrelation time in the conventional DQMC simulation and make the larger systems and lower temperature accessible.

As will be shown in the next section, SLMC worked very well for the spin-fermion models, and as we started to explore more related models, further numerical improvements, inspired by the theoretical understanding of the fermionic quantum critical points, were also made possible. Among them, elective momentum ultra-size quantum Monte Carlo method (EMUS) is of particular advantage in dealing with fermionic quantum critical point subject to antiferromagnetic or other finite momentum bosonic fluctuations. EMUS is developed by us in Ref. Liu et al. (2019), it is inspired by the awareness that in a finite $\mathbf{Q}$ quantum critical point, the critical bosonic fluctuations mainly couple to fermions near the hot spots, which are the points on the FS connected by the finite $\mathbf{Q}$ bosonic fluctuations as shown in Fig. 3 in the case of square lattice band structure subject to antiferromagnetic $\mathbf{Q}=(\pi,\pi)$ fluctuations. Thus, instead of including all the fermion degrees of freedom, we ignore fermions far away from the hot spots and focus only on momentum points near the hot spots in the simulation. This approximation will produce different results for non-universal quantities compared with the original model, such as critical field $h_{c}$ or critical temperature $T_{c}$, as one actually performed a hard UV cutoff to the problem at hand. However, for universal quantities at the IR, such as scaling exponents, which are independent of microscopic details and the high energy cutoff, EMUS has been shown to generate consistent values with those obtained from standard DQMC Liu et al. (2019, 2018), both for the itinerant QCPs on triangle and square lattices.

In DQMC, configurations of the bosonic modes $\phi_{i}$ are stochastically sampled, with weights obtained through integrating out fermion modes. Conventionally, real-space single particle modes, $c_{ia}$, are chosen as the basis on which the calculation is performed. While in EMUS, we instead use the momentum space. To this end, we rewrite Eqs. (LABEL:eq:specificformula) in momentum space,

[TABLE]

and

[TABLE]

Here, $\phi_{k}$ denotes the $k$-component of the Fourier transform of the bosonic field $\phi_{i}$:

[TABLE]

Rewriting the problem in the momentum space gives us a freedom of choosing arbitrary $k$ points in the summation in Eqs. (21) and (22). When to study the low-energy and long-wavelength physics, we choose IR fermion modes that are particularly relevant, and throw away other modes without worrying a proper UV completion of a lattice model. Particularly, for studying AFM-QCP as shown in Fig. 3 (a), only fermion modes near hot spots, where two patches of FSs are connected by the ordering wave vector $\mathbf{Q}_{i}$ [see Fig. 3 (b)], are relevant to universalities. As usually found in the analytic calculation of such problems Metlitski and Sachdev (2010b); Abanov and Chubukov (2004), one keeps modes in patches around these hot spots, and neglect other modes in the BZ, as shown in the patches in Fig. 3 (b). In this way, the number of fermion modes used in computing the effective weights is greatly reduced, from the total size (or total volume) $N$ to the patch size (or patch volume) $N_{f}$. Therefore, while retaining the same IR physics of what a lattice DQMC simulation of system size $N$ can achieve, EMUS has drastically lifted the computational burden.

In EMUS, as local couplings in real space become non-local in the momentum basis [Eq. (22)], making a simple local update in the standard DQMC costs $\beta N\cdot O(\beta N_{f}^{3})$ computational complexity, one can no longer use such kind of local update. Fortunately, the cumulative update scheme in the SLMC developed Liu et al. (2017b, 2018, 2019, 2018); Chen et al. (2018, 2019) come to help. Such a cumulative update is a global move of the Ising spins and gives rise to the computational complexity at most $O(\beta N_{f}^{3})$ for computing the fermion determinant. Since $N_{f}$ can be much smaller than $N$, a speedup of the order $(\frac{N}{N_{f}})^{3}\sim 10^{3}$ of EQMC over DQMC, with $\frac{N}{N_{f}}\sim 10$, can be comfortably achieved.

IV Results

In this section, we present several of our recent works on the fermionic quantum critical points, revealed with the quantum Monte Carlo methods discussed in the previous section. Most of these results can be found in the several recent publications, including the Ising ferromagnetic quantum critcal point (FM-QCP) studied in Ref Xu et al. (2017a), the Ising antiferromagneitc quantum critical point (AFM-QCP) on the triangular lattice (with 3$\mathbf{Q}=\Gamma$) studied in Ref Liu et al. (2018) and the Ising antiferromagnetic quantum critical point on the square lattice (with $2\mathbf{Q}=\Gamma$) studied in Ref. Liu et al. (2018). From the point of view of a pedagogical narrative, here we took the content of each paper and restructured them in a coherent manner, such that the overall logical flow and the internal relation of the physical results beyond each individual papers manifest.

IV.1 Phase Diagrams

IV.1.1 FM-QCP

We start from the ferromagnetic quantum critical point (FM-QCP). In the general designer Hamiltonian illustrated in Eq. (LABEL:eq:specificformula), we take $J<0$ to realize a FM-QCP and particularly, we are interested in the case where the Fermi surface away from perfect nesting which can be realized by tuning chemical potential. We draw the microscropic lattice model again in Fig. 4 (a).

In the simulation, we set $t=1$, $J=-1$ and chemical potential $\mu=-0.5$, resulting in a fermion density $\langle n_{i\lambda}\rangle\approx 0.8$, with FS away from perfect nesting. The coupling strength between the fermion and bosons is set to $\xi=1$. By tuning the transverse field $h$, we realized an itinerant PM-FM transition as shown in the phase diagram of Fig. 5(a). In the temperature-transverse filed ($T-h$) phase diagram, the black dots are the phase boundary of bare transverse Ising model without coupling to fermions ($\xi=0$ case) where the finite temperature paramagnetic (PM) to ferromagnetic (FM) phase transition belongs to 2D Ising universality class, and the zero temperature PM-FM transition (large black dot) belongs to 3D Ising universality class. When turning on the couplings between the Ising spins and fermions ($\xi=1$ here), the phase boundary shift a little bit to larger $h$ direction (red points), and the finite temperature PM-FM phase transition is still 2D Ising as at finite temperature the Fermion can be trivially integrated out without influencing the criticality, while at zero temperature (large red point), there exist no controllable analytical tools to take care of critical fluctuations of these strongly coupled fermions and bosons, and the criticality here is believed to be different from 3D Ising universality and in fact different from any known ones. Our QMC numerical results reveal the properties of such itinerant 3D FM-QCP and indeed discovered the scaling behavior near the QCP with new critical exponents significantly deviate from 3D Ising with critical FS of non-FL and large anomalous dimensions, as will be discussed in Secs. IV.2 and IV.3.

IV.1.2 AFM-QCP with $3\mathbf{Q}=\Gamma$

We choose a triangular lattice to realized a AFM-QCP with $3\mathbf{Q}=\Gamma$-by setting $J>0$ in the general designer Hamiltonian in Eq. (LABEL:eq:specificformula), as shown in Fig. 4 (b). The notation of $3\mathbf{Q}=\Gamma$ means that three times length of the antiferromagnetic fluctuations $\mathbf{Q}=(\frac{2\pi}{3},\frac{2\pi}{\sqrt{3}})$ will bring the wavevector back to the $\Gamma$ point in the extended BZ scheme. Once we apply the same mode on square lattice, it is obvious that twice length of the $\mathbf{Q}=(\pi,\pi)$ wavevector will bring the fluctuations back to $\Gamma$ in the extended zone, therefore dubbed AFM-QCP with $2\mathbf{Q}=\Gamma$. As will become clear later in Sec. IV.3, these two types of AFM-QCPs actually have different scaling behaviors, consistently revealed both from QMC simulation and theoretical argument.

We set $t=1$, $J=1$, $\mu=-0.5$ (electron density $\langle n_{i,\lambda}\rangle\sim 0.8$ and explore the $T$-$h$ phase diagram as shown in Fig. 5(b). When the Ising spin and fermion is not coupled ($\xi=0$ case), the phase diagram is already well studied Moessner and Sondhi (2001); Isakov and Moessner (2003); Wang et al. (2017), and the phase boundary is denoted as dashed lines Fig. 5(b). At low temperature and small $h$, a clock phase is formed, which is characterized by a complex order parameter $me^{i\theta}=m_{1}+m_{2}e^{i4\pi/3}+m_{3}e^{-i4\pi/3}$ where $m_{\alpha}$, with $\alpha=1,2,3$, equals to $\frac{1}{3N}\sum^{N/3}_{i=1}s^{z}_{i,\alpha}$, representing magnetization of the three sublattices. In clock phase, the mentioned order parameter has a finite momentum $\mathbf{Q}=(\frac{2\pi}{3},\frac{2\pi}{\sqrt{3}})$. At $T=0$, increasing $h$ destroys clock phase through a continuous quantum phase transition at $h_{c}=1.63(1)$ and a U(1) symmetry is emerged at QCP, making the quantum phase transition guided by the (2+1)D XY universality class Isakov and Moessner (2003); Wang et al. (2017). When increasing temperature, the clock phase finally melts and thermal melting goes through an intermediate BKT phase. When the Ising spin fermion is coupled ($\xi=1$ case), the phase boundary shifted to a larger ordered region and the QCP is now at $h_{c}=1.83(1)$ as shown in Fig. 5(b). As the Ising spin clock phase generates a SDW order (with wavevector $\mathbf{Q}$) in the fermionic sector, the fermion band is folded with six FS pockets and six pairs of hot spots connected by the $\mathbf{Q}$ wavevector (the detailed FS structure can be seen in Fig.1 of Ref. Liu et al. (2018)). In contrast to FM QCP case with a $\mathbf{Q}=0$ where the entire FS is critical, here with a finite momentum $\mathbf{Q}$ QCP, only the hot spots are critical, and we saw non-Fermi liquid behavior at hot spots, as will be discussed in Sec. IV.2 and Sec. IV.3. We have performed complementary DQMC and EMUS simulations on the model, as shown in Ref. Liu et al. (2018) and Ref. Liu et al. (2019) respectively.

IV.1.3 AFM-QCP with $2\mathbf{Q}=\Gamma$

Antiferromagnetic quantum critical point can be realized in many forms. Here we implement it in the form of square lattice with $\mathbf{Q}=(\pi,\pi)$ fluctuations. We realized a $2\mathbf{Q}=\Gamma$ AFM-QCP by setting $J>0$ in the general designer Hamiltonian in Eq. (LABEL:eq:specificformula) and the microscopic lattice model is depicted in Fig. 4(a).

To tune the FS, in particularly to make it represent that of the high temperature superconductor cuprates, we design fermion hopping part as

[TABLE]

and choose $t_{1}=1.0$, $t_{2}=-0.32$, $t_{3}=0.128$, $\mu=-1.11856$ (electron density $\langle n_{i,\lambda}\rangle\sim 0.8$), according to Ref. Chowdhury and Sachdev (2014).

We set $J=1$ and implement both DQMC and EMUS here to explore the $T-h$ phase diagram as shown in Fig. 5(c). The AFM-QCP with $2\mathbf{Q}=\Gamma$ is realized by turning on Ising spin fermion coupling (we set $\xi=1$ here). The standard DQMC with SLMC update scheme greatly reduce the autocorrelation time and make the larger systems and lower temperature accessible. We use it provides unbiased results as a benchmark. On the other hand, EMUS concentrates on fermions with momentum near the hot spots, therefore it enjoys a speed up of the order $(\frac{N}{N_{f}})^{3}$ and makes even larger system size accessible. Fig. 6 shows FS obtained by $G(\mathbf{k},\beta/2)\sim A(\mathbf{k},\omega=0)$ with both DQMC and EMUS. Higher momentum resolution near hot spots is obtained in EMUS, and it plays a vital role in unveiling critical behavior of this QCP.

As will be discussed in detail in Sec. IV.2 and Sec. IV.3, although the AFM-QCP with $3\mathbf{Q}=\Gamma$ and $2\mathbf{Q}=\Gamma$ share similarities, such as the formation of FS pockets and hot spots. They give rise to very different quantum criticalities. It turns out the length of the antiferromagnetic wavevector actually matters in terms of the effects of the critical fluctuations. In the triangular lattice case, the QCP is found to be consistent with the HMM prediction, within our current energy and momentum resolution. However, in the square lattice case, quantum criticality with large anomalous dimension is discovered with scrutiny. The discovery of anomalous dimension is of great significance in the theoretical understanding of the fermionic quantum critical points. As many field theoretical developments beyond the HMM, such as those in Refs. Abanov et al. (2003); Abanov and Chubukov (2004); Metlitski and Sachdev (2010b); Metzner et al. (2012); Schlief et al. (2017); Lee (2018), the anomalous dimension at the AFM-QCP (actually also the FM-QCP but it is even more difficult to handle analytically) have been suggested to have various forms, such as a function of the number of hot spots $N_{\text{h.s.}}$, but the validity of such renormalization group (RG) calculations are not known a priori, since it is still the fundamental question that a proper small parameter for the RG flow for fermionic quantum critical point is in absence. Therefore, our quantum Monte Carlo results actually provide the first set of unbiased attempts to provide concrete results for the analytical approaches to compare with, and from such comparison, further adjustments from both theoretical and numerical sides, can be made to reach the eventually solution of the problem. In this regard, our results shown in the review, in particular in Sec. IV.3, will be of great importance for the further developments of the community.

IV.2 Non-Fermi liquid behavior

Now let’s focus on how the fermionic degrees of freedom react at the fermion QCPs. In particular, we would like to find out whether nFL behavior with the low-temperature fermionic quasi-particle weight vanishes at the Fermi surface or hot spots Oganesyan et al. (2001); Senthil (2008); Lee (2009b); Metlitski and Sachdev (2010a, b); Dalidovich and Lee (2013); Schlief et al. ; Lee (2017) due to strong damping of the fermions induced by critical bosonic fluctuations, can be seen from our simulations. Following Chen et al. (2012), the first Matsubara frequency ($\omega_{0}=\pi T$) of self-energy gives a good estimation of quasi-particle weight

[TABLE]

as shown in Fig. 7(a) for $h=h_{c}$ (squares) and $h>h_{c}$ (circles) for the FM-QCP case. The behavior of quasi-particle weight with temperature in the PM phase ($h>h_{c}$) and at QCP ($h=h_{c}$) is totally different, that it remains close to unity at low temperature in the PM phase, indicating well-defined quasi-particles, while it is strongly suppressed and extrapolates to zero as $T\rightarrow 0$ at QCP, which is a key signature of a non-Fermi liquid. The different quasi-particle behaviors in the PM phase and at QCP is also revealed in the fermions self-energy and Green’s function data as shown in Fig. 7(b) and (c). The increasing of both $-\text{Im}(\Sigma)$ and $-\text{Im}(G)$ with decreasing $\omega_{n}$ at QCP indicates strong damping of the fermions at low frequencies but without any signature of gap opening. Such enhancement of the self-energy clearly suggests the strong correlation at the QCP, which means either this is a new fermionic QCP without previously knowledge or that when the temperature/energy scales go lower than what have accessed by now the system might develop a first order transition Belitz et al. (2005). But either the case is interesting and have not been addressed with any unbiased approaches before.

In the FM-QCP, since the bosonic fluctuations are at $\mathbf{Q}=0$, the entire FS become critical and as shown in Fig. 7 the quasiparticles indeed lose coherence at the FM-QCP. The AFM-QCPs, in this regard, are different and our results show that it is only at the hot spots points, i.e., the momentum points connected by the antiferromagnetic wavevectors, the nFL behavior can be observed at the AFM-QCP, as shown in Fig.2 in Ref. Liu et al. (2018) for triangular lattice case and in Fig.4 in Ref. Liu et al. (2018) in the square lattice case. But at the momenta away from the hot spots, the FL behavior of the quasiparticle keep intact as their quasiparticle fraction measured in Eq. (25) remains close to 1. From these results, one could also tell that actually the FM-QCP is harder to handle than the AFM-QCP in that at the FM-QCP all the fermionic modes on the FS are strongly coupled to the critical bosonic fluctuations and it is therefore difficult to develop controlled RG calculation scheme. At the AFM-QCP, only finite number of hot spots fermionic modes are critically coupled with finite $\mathbf{Q}$ bosonic fluctuations (in the case in Sec. IV.1.3, $N_{\text{h.s.}}=16$). As will be seen in Sec. IV.3, a finite $\mathbf{Q}$ QCP is theoretically better controlled in comparison to QCPs with $\mathbf{Q}=0$, and recent numerical results show good agreement with theory predictions at the qualitative level. However, at the quantitative level, numerical results still differ from theory predictions by a numerical factor, which requires further investigations and efforts from both numerical and analytical sides.

IV.3 Quantum critical scaling analysis

After showing the fermionic degrees of freedom at the quantum critical points, we now move on to the most difficult and yet the most direct observable that reveal the nature of the QCP – the dynamic magnetic susceptibilities. In fact, all the previous analytical works, are in one way or the other, trying to find controlled scheme to integrate out the fermionic degrees freedom at the quantum critical point and arrive at a bosonic susceptibility that captures the unique divergence at the QCP. Therefore, our results in this section provide solid examples that analytical approaches can be used to test with and might inspire further analytical developments.

IV.3.1 FM-QCP

We again start with the FM-QCP. Our results in Ref. Xu et al. (2017a) show that the critical behavior of itinerant FM-QCP is strongly deviates from 3D Ising universality class. The first evidence comes from the PM-FM transition temperature, as shown in Fig. 5(a), the phase boundary of the FM onset temperature scales as $T_{N}(h)\sim|h-h_{c}|^{c}$. We have $c=0.77(4)$ for the FM-QCP case while that for 3D Ising is $c=\nu z\approx 0.63$, also note that due to itinerant nature of the QCP, the exponent $c$ is no longer expected to obey the relation $c=\nu z$ Löhneysen et al. (2007); Belitz et al. (2005).

A more rigorous analysis of scaling behavior is based on the scaling analysis of the data inside the QCR, in which the bosonic critical modes become strongly renormalized due to the coupling to gapless fermions. Our QMC results indicate that the scaling behavior deviates from the HMM predictions significantly and we propose a modified Hertz-Millis scaling formula which fits our QMC data for the Ising spin susceptibility at all the momenta and frequencies simulated near QCP in the long wavelength and low frequencies limit.

As discussed in Sec. II.2, even at the level of random phase approximation (RPA), non-trivial quantum dynamics has already entered the bosonic susceptibility Hertz (1976); Millis (1993); Moriya (1985). We, however, will measure the quantity unbiasedly in the quantum Monte Carlo simulations, the Ising spin susceptibility is given by,

[TABLE]

At the RPA level, $\chi(h,T,\mathbf{q},\omega_{n})$ has the following hypothesized form near the QCP,

[TABLE]

where $c_{t}$, $c_{h}$, $c_{q}$, $c_{\omega}$ are constants. Here, the $c_{q}q^{2}+c_{\omega}\omega^{2}$ comes from the bare action of the Ising degrees of freedom, and the $\Delta(\mathbf{q},\omega_{n})$ term is the contribution of the fermionic fluctuations, and has following form

[TABLE]

in the isotropic and low momentum/energy limit. The $c_{\mathrm{HM}}$ is a constant and $v_{f}$ is the Fermi velocity. The form in Eq. (28) has the following meaning that due to the coupling between critical bosons with fermions, the limits of $q=0$, $\omega\rightarrow 0$ and $\omega=0$, $q\rightarrow 0$ no longer commute with each other, which also points out the $z=1$ Lorentz symmetry present at the bare 3D Ising boson critical point no longer hold in the coupled case, that the FM-QCP is indeed a new quantum criticality. Such statement will be verified with our numerical data.

Our QMC results share some characteristics with the above form, such as the singular behavior that $\lim_{q\to 0}\lim_{\omega_{n}\to 0}\chi^{-1}(\mathbf{q},\omega_{n})$ differs from $\lim_{\omega_{n}\to 0}\lim_{q\to 0}\chi^{-1}(\mathbf{q},\omega_{n})$ by a constant $c_{\mathrm{HM}}=0.20(4)$, but importantly, we found that the RPA form shall be adjusted to a form with finite anomalous scaling dimensions,

[TABLE]

with $a_{q}=1.85(3)$ ($\eta=0.15(3)$ from $a_{q}=2-\eta$) . This is different both from the exponent for an Ising transition in $(2+1)$ dimensions, $a_{q}=1.96$ ($a_{q}=2-\eta$ and $\eta=0.04$ for 3D Ising), and from the HMM RPA value of $a_{q}=2$. The presence of an anomalous exponent is another key finding of the study in Ref. Xu et al. (2017a). The singular behavior that in the limits of $q=0$, $\omega\rightarrow 0$ and $\omega=0$, $q\rightarrow 0$ and that they do not commute are also found in a closely related quantum Monte Carlo investigation of nematic QCP Schattner et al. (2016a), but due to the limitation of system size and temperature, the errorbar in the susceptibility in Ref. Schattner et al. (2016a) is too large (at the level of 0.3) to make any concrete statement on the absence and existence of the anomalous dimension we find here.

We fit all the constants and exponents in Eq. (29) based on our QMC data. For example, when setting $q=0$ and $\omega_{n}=0$ and follows the temperature and magnetic field dependence, we can obtain $c_{t}=0.13(1)$, $a_{t}=1.48(4)$, $c_{h}=0.7(1)$ and $\gamma=1.18(4)$, these results can be seen in Fig.8 and 9 in Ref. Xu et al. (2017a). Then one can set $\omega=0$, $h=h_{c}$ and at a very low temperature, to monitor the momentum dependence of $\chi(q)$ in Eq. (29), the results are shown in Figs. 8 (a). From two different system sizes, we can obtain $a_{q}=1.85(3)$ and $c_{q}=1.00(2)$. The finite $a_{q}$ suggests a finite anomalous dimension $\eta=0.15(3)$. In the same vein, one can set $q=0$, $h=h_{c}$ and at a very low temperature, to monitor the frequency dependence of the $\chi(\omega)$ in Eq. (29), the results are shown in Fig. 8 (b), from the fit, we obtain $c_{\omega}=0.10(2)$ and $c_{\text{HM}}=0.20(4)$. The finite $c_{\text{HM}}$ indicate the two limits $q=0$, $\omega\rightarrow 0$ and $\omega=0$, $q\rightarrow 0$ do not commute, at least upto the system sizes and temperature accessed within our QMC simulation.

One way to show that the above fitting is robust is to use the obtained exponents and coefficients to collapse all $\chi$ data for all $q$, $\omega_{n}$, $h$ and $T$ simulated according to Eq. (29). Such a collapse is shown in Fig. 8 (c), where all the data collapse onto the diagonal line, especially for small $q$, $\omega_{n}$, low temperature $T$ and $h\sim h_{c}$. Such a collapse strongly suggests that the FM-QCP revealed with our finite size QMC simulations is a novel quantum critical point and the results obtained will be useful for future analytical works to test with.

One more thing we want to remark here is about the exact form of $\Delta(\mathbf{q},\omega_{n})$, we found that as long as $\lim_{\omega_{n}\to 0}\Delta(\mathbf{q},\omega_{n})=0$ and $\lim_{q\to 0}\Delta(\mathbf{q},\omega_{n})=c_{\mathrm{HM}}=0.20(4)$ is kept, the impact of specific functional form of $\Delta(\mathbf{q},\omega_{n})$ is little within numerical error bars of the QMC data, and this uncertainty in $\Delta(\mathbf{q},\omega_{n})$ implies that the exact form of $\Delta(\mathbf{q},\omega_{n})$ can be much more complicated than the simple RPA approximation and the current system sizes and temperature are still not sufficient to distinguish their differences (e.g., the precise value of the dynamical critical exponent $z$). To further explore the nature of the FM-QCP, future simulations with larger system sizes and lower temperatures are still necessary and will be performed with even better Monte Carlo techniques, such as updated version of self-learning Monte Carlo. Fortunately, this does not affect our analysis at the static limit ($\omega_{n}=0$), and all the conclusions in this limit, including the anomalous dimensions $\eta=2-a_{q}=0.15(3)$, stand without any dependence on the details of $\Delta(\mathbf{q},\omega_{n})$.

IV.3.2 AFM-QCP with $3\mathbf{Q}=\Gamma$

For AFM-QCP with finite $\mathbf{Q}$, there is no singularity behavior under different orders of $\omega\rightarrow 0$ and $q\rightarrow 0$ limits, such equality applies for both $3\mathbf{Q}=\Gamma$ and $2\mathbf{Q}=\Gamma$ cases. For the former, we found the spin susceptibility can be well captured by the form

[TABLE]

At low temperature and frequency, the linear $T$ and $\omega$ terms dominant, this is in consistent with HMM prediction for antiferromagnetic cases. At higher temperature and frequency, the quadratic $T^{2}$ and $\omega^{2}$ terms dominant, in consistent with bare boson (2+1)D XY universality class for the antiferromagnetic Ising transverse field model on triangular lattice Moessner and Sondhi (2001); Wang et al. (2017); Isakov and Moessner (2003). Such crossover behavior from low to high energy scale is well established in Eq. (30) and in our QMC data. One thing we want to highlight here is we do not find anomalous dimension (see Fig. 9) in the $3\mathbf{Q}=\Gamma$ AFM-QCP case. As discussed in Sec. II.2 and will be discussed in Fig. 11 in Sec. IV.3.3, this is consistent with the understanding that an important scattering channel is missed in this case, whose effect will be seen in the $2\mathbf{Q}=\Gamma$ AFM-QCP case. Our Monte Carlo simulation results in this and next sections are the first set of unbiased numerical data to reveal such interesting and subtle difference among the AFM-QCPs.

IV.3.3 AFM-QCP with $2\mathbf{Q}=\Gamma$

The scaling analysis of AFM-QCP with $2\mathbf{Q}=\Gamma$ is presented in this section. Comparing with the AFM-QCP with $3\mathbf{Q}=\Gamma$ in Sec. IV.3.2, AFM-QCP with $2\mathbf{Q}=\Gamma$ shows interesting features deviated further away from HMM. We found the magnetic susceptibility with the following form with an anomalous dimension can guide our fitting of QMC data very well.

[TABLE]

From the momentum dependence of the magnetic susceptibility, as shown in Fig. 10, we obtain a finite anomalous dimension $\eta=0.125$. For the frequency dependence, a crossover behavior is observed, and for quantity $\chi^{-1}\equiv\chi^{-1}(T,h_{c},0,\omega)-\chi^{-1}(T,h_{c},0,0)$, it is well described by

[TABLE]

At low frequency, the first term with power law form $\omega^{0.875}$ dominant, consistent with the anomalous dimension $\eta=0.125$ found in momentum dependence. At higher frequency, the $\omega^{2}$-term dominants, back to the bare bosonic susceptibility of (2+1)D Ising universality class. We want to remark here that although our QMC data is consistent with dynamical exponent $z=2$ of HMM, significant derivations in the form of anomalous dimension in the dynamical exponent as predicted in Ref. Metlitski and Sachdev (2010b) cannot be excluded. Interestingly, the $O(2)$ itinerant AFM-QCP case is considered in another recent numeric work Gerlach et al. (2017), where the dynamical exponent $z=2$ is found as well, consistent with the one-loop RG results, but higher order correction indeed exists in the latest RG calculation Schlief et al. (2017). However, for the Ising case we considered here, such similar calculations are still absent. It will be of great theoretical and numerical interests, that such a higher order RG calculation, can be performed in the AFM-QCP with Ising symmetry and $2\mathbf{Q}=\Gamma$ case as well.

The anomalous dimension $\eta$ found in AFM-QCP with $2\mathbf{Q}=\Gamma$ case but not in AFM-QCP with $3\mathbf{Q}=\Gamma$ case is highly nontrivial and worth to be further discussed. Below we try to understand such difference and the implication of the presence of $\eta$ in the $2\mathbf{Q}=\Gamma$ case in more detail.

First, between $2\mathbf{Q}=\Gamma$ and $3\mathbf{Q}=\Gamma$, the constraints dictated by the momentum conservation law are different. The QCPs with $2\mathbf{Q}=\Gamma$ deviates from the HMM already at the level of four-boson vertex correction Abanov et al. (2003), as shown in Fig. 11 (a). This four-boson vertex gives two topologically different bosonic self-energy diagrams as shown in Fig. 11 (b) and (c) when $2\mathbf{Q}=\Gamma$, and particularly the diagram shown in Fig. 11 (c) results in logarithmic corrections and is responsible for the breakdown of the HMM scaling. However for the $3\mathbf{Q}=\Gamma$ case in the triangular lattice model, this crucial diagram is prohibited by the momentum conservation law, hence the deviations from the HMM is not expected, at least within the same level of approximation.

Second, let us discuss the meaning of the particular value of the anomalous dimension we found. Comparing to existence theories, the Heisenberg AFM-QCPs with $SU(2)$ symmetry has been considered in the literature Abanov et al. (2003); Metlitski and Sachdev (2010b), while the Ising spin case is not yet carefully analyzed. We compare our numerical results with existing theoretical predictions from Heisenberg AFM-QCPs, qualitative but not quantitative agreement is expected because of the symmetry difference. In the calculations in Refs Abanov et al. (2003); Metlitski and Sachdev (2010b), the anomalous dimension depends on the angle between Fermi velocity at the hot spots and the AF wavevector $\mathbf{Q}$ of order parameter. In our model, we have the angle close to $45^{\circ}$. But the RG calculations in Refs. Abanov et al. (2003); Metlitski and Sachdev (2010b) actually assume the angle equals zero degree at the fixed point, and only when that happens, the RG predicted $\eta=1/N_{h.s.}$ where $N_{h.s.}$ is the number of hot spots. In our model, we have $N_{h.s.}=16$, thus RG predicted $\eta=1/16$, while our QMC results find qualitatively consistency but with a value close to $2/N_{h.s.}$. The reason of the difference, possible comes from the symmetry difference or any other contributions, is an open issue, and further RG calculations on the AFM-QCP with Ising symmetry is highly desirable.

Another prediction from RG calculation is that near the QCP, the Fermi surface at hot spots will rotate towards nesting Abanov et al. (2003), and it further increases the anomalous dimension and can even renormalize dynamical exponent $z$ Abanov et al. (2003); Metlitski and Sachdev (2010b). Such a rotation of Fermi surface is indeed observed in very recent QMC studies Liu et al. (2018). However, the rotation we observed is small due to very slow RG flow and thus the resulting changes in anomalous dimension and dynamical exponent is beyond numerical resolution at the current stage. Even larger system sizes and lower temperatures, or better designer Hamiltonians, are needed to be accessed in future QMC simulations.

IV.4 Superconductivity

In the vicinity of QCP, other instabilities like superconductivity may emerge, which usually may bury the QCP, pollute QCR and prevent us to do a faithful scaling analysis down to very low temperature. The construction we make seems to enjoy a very pristine QCP that does not develop any superconductivity down to lowest temperature we studied with coupling constant $\xi=1$. Further increasing $\xi$, it may finally drives a superconductivity instabilities, which is also find in other related numerical studies, such as in Refs. Schattner et al. (2016b); Lederer et al. (2017); Gerlach et al. (2017).

Fig. 12 shows the superconductivity correlations for different coupling constant in FM-QCP case. We do find as the coupling constant $\xi$ increases, the pairing correlations start to increase with system size near QCP, which is a signature of pairing instability.

V Quantum criticalities of fermions coupled to phonons, $\mathbb{Z}_{2}$ or $U(1)$ gauge fields

The philosophy of designer Hamiltonian in Eq. (7) is generic, besides couple fermion to critical bosonic modes in the forms of ferromagnetic, antiferromagnetic and nematic fluctuations, the boson modes can be extended to other situations such as Holstein photons Johnston et al. (2013); Chen et al. (2018, 2019) to study the electron-phonon interaction mediated metal to charge density wave (CDW) transition Johnston et al. (2013); Chen et al. (2018, 2019) and metal to superconductivity transition Esterlis et al. (2018); Costa et al. (2018); Li et al. (2018), and these interaction-driven phases. Along this line of research, new quantum Monte scheme in the form of SLMC which introduces symmetry-enforced global update of phonon field to reduce the autocorrelation time Chen et al. (2018) and Langevin dynamics Batrouni (1987); Batrouni and Scalettar (2019) which performs efficient global update to reduce the computational complexity have been successfully introduced. The situation, that simulating Holstein model at 2D or higher was extremely difficult compared with other DQMC simulations for, say, Hubbard model Hohenadler and Lang (2008), has been profoundly changed in the last few years.

More interestingly, the exact form of critical bosons can even go beyond the conventional condensed matter wisdom and acquire incarnations with high-energy flavor, such that they could play the role of gauge degree of freedom and carry gauge symmetry and topological orders. In this way, the situation of matter fields couples to gauge fields can be realized in our setting and the rich physics of fractionalization of electrons to anyons coupled to the emergent gauge fields of $\mathbb{Z}_{2}$ Assaad and Grover (2016); Gazit et al. (2017, 2018); Chen et al. (2019); Gazit et al. (2019) and $U(1)$ Xu et al. (2019); Wang et al. (2019) symmetries in many strongly correlated systems such as $\mathbb{Z}_{2}$ and $U(1)$ topological orders Wen (1991); Senthil and Fisher (2000); Rantner and Wen (2002); Hermele et al. (2005) and their material and model realizations in quantum spin liquids Fu et al. (2015); Feng et al. (2017); Hermele et al. (2004); Sun et al. (2018); Wang et al. (2018b); Huang et al. (2018) and deconfined quantum criticalities Senthil et al. (2004); Sandvik (2007); Qin et al. (2017); Ma et al. (2018), and the fundamental question of the existence of deconfinement at (2+1)D quantum electrodynamics (QED3) with fermionic matter in high-energy physics Polyakov (1977); Burkitt and Irving (1988); Fiebig and Woloshyn (1990); Armour et al. (2011), can now be addressed with unbiased quantum Monte Carlo simulations. The development in this direction is fast and profound.

In this section, we try to briefly summarize the activities in the aforementioned areas in recent years, based primarily upon the works of our own.

V.1 Dirac fermions coupled with phonons

As mentioned above, although the electron-phonon coupled systems, realized in the Holstein model or Holstein-Hubbard model, play important role in the understanding of metal to charge density wave (CDW) transition Johnston et al. (2013); Chen et al. (2018, 2019); Zhang et al. (2019) and metal to superconductor transition Esterlis et al. (2018); Costa et al. (2018); Li et al. (2018) and these interaction-driven phases themselves, the quantum Monte Carlo simulations of Holstein model in 2D or higher dimensions are notoriously difficult due to strong autocorrelation times Hohenadler and Lang (2008). The continuous phonon fields render the local update of conventional DQMC inefficient and the growth of autocorrelation time in the DQMC of Holstein model can be shown much faster than that of the Hubbard model, see Fig.2 (a) in Ref. Chen et al. (2018) of 2D Holstein model at metal-to-CDW transition temperature for example, the measured autocorrelation time $\tau\sim L^{5}$, much larger than the typical $\tau\sim L^{2}$ for the local update of Metropolis algorithms. Recently, we applied SLMC method on 2D Holstein model and found an efficient way of designing the effective bosonic Hamiltonian that enforces the $Z_{2}$ symmetry of the on-site energy potential of the phonons, in this way, the SLMC of Holstein model can reduce the autocorrelation time to $\tau\sim L^{2}$ Chen et al. (2018) and makes the 2D lattice systems simulatable.

With such development, we investigated the electron-phonon coupled problem on 2D honeycomb lattice Chen et al. (2019). This model is relevant to the electron-phonon coupling in graphene and the twisted angle graphene multilayers with superlattices Cao et al. (2019, 2018); Shen et al. (2019). The Dirac fermion in the graphene honeycomb lattice is robust against the weak coupling with phonons. But when the coupling is strong enough, the fluctuations of the phonons modes would generate effective attraction between electrons and eventually lead to a CDW insulator. In the work of Ref. Chen et al. (2019), such a metal-to-CDW transition is studied with unprecedentedly large system sizes (qauntum phase transition at $T=0$ with electron number of $2\times L\times L$ upto $2\times 15\times 15$), and the transition is revealed with dynamically generated mass of Dirac fermion and belongs to $N=2$ flavors (2+1)D chiral Gross-Neveu universality He et al. (2018); Chandrasekharan and Li (2013); Zerf et al. (2017). The critical exponents are determined and consistent with previous quantum Monte Carlo results with the same universality but physically different systems He et al. (2018); Chandrasekharan and Li (2013). The comparison of the exponents with high order perturbative RG calculation Zerf et al. (2017) is still not completely satisfactory and this would require future works with even larger system sizes, also more controlled calculation from the theoretical end. Judging from the developments and achievements summarized in this review, we are confident that eventually with the mutual stimulation and developments from numerical and analytical communities, the inconsistency in the chiral Gross-Neveu universities of the gapped out Dirac fermions with spontaneous fermion bilinear of symmetry-breaking, would be cured.

V.2 Dirac fermions coupled with $\mathbb{Z}_{2}$ or $U(1)$ gauge fileds

Also as aforementioned, in the other direction where fermions are coupled to $\mathbb{Z}_{2}$ or $U(1)$ gauge fields, the system becomes directly related to both high energy and condensed matter physics. This is because in the high energy physics, there exists a fundamental question of whether the deconfined phase exists or not in the setting of fermion coupled to the gauge field in QED3, where it is known that the pure gauge field in QED3 will always be confined due to monopole proliferation Polyakov (1977), but with the help of the matter fields, the situation might be different and this question relies on the unbiased quantum Monte Carlo simulation to verify. And in the condensed matter physics setting, the situation of the fermions couple to gauge fields is believed to provide the low energy descriptions of many condensed matter systems ranging from high-temperature superconductors Lee et al. (1998, 2006) to $\mathbb{Z}_{2}$ or $U(1)$ quantum spin liquids Fu et al. (2015); Feng et al. (2017); Hermele et al. (2004); Sun et al. (2018); Wang et al. (2018b); Huang et al. (2018) and deconfined quantum critical points Senthil et al. (2004); Sandvik (2007); Qin et al. (2017); Ma et al. (2018), etc. Due to these reasons, the quantum Monte studies of matter fields coupled to gauge fields are of great interests to very broad audience.

The Dirac fermions coupled with $\mathbb{Z}_{2}$ gauge fields in 2D systems are studied in Refs. Gazit et al. (2017); Assaad and Grover (2016); Gazit et al. (2018, 2019). In Refs. Gazit et al. (2017, 2018), the gauge constrain is enforced in the Monte Carlo simulation and in Ref. Assaad and Grover (2016), the gauge constrain is dynamically generated at low temperature. But overall, similarly rich phase diagrams are established at different number of fermion flavors and the strength of the coupling of gauge field with matter field (Dirac fermions). There always exists a deconfined phase where the Dirac cones in the matter fields persist and coexist with an emergent $\mathbb{Z}_{2}$ topological order. As the coupling becomes weak, the systems enter a confined phase via deconfinement-to-confinement phase transition, and depending on the number of fermion flavors, the confined phases develop various symmetry-breaking in the form of CDW, valence bond solid (VBS) and AFM phases, in all these phases, the Dirac fermion is gapped out. In the work of Ref. Gazit et al. (2018), a on-site Hubbard $U$ term is added into the Hamiltonian, and the large $U$ limit is also a AFM phase but with $\mathbb{Z}_{2}$ topological order. Therefore, there is another deconfinement-to-confined phase transition between two different AFM phases. The transition from deconfined phase to AFM with confinement is found to possibly host an emergent $SO(5)$ symmetry Gazit et al. (2018). All these results are of great interests. And the question of Fermi surface coupled to $\mathbb{Z}_{2}$ topological order via the help of $\mathbb{Z}_{2}$ matter field has recently been address by some of us Chen et al. (2019).

More recently, Dirac fermions coupled $U(1)$ gauge fields problem is studied by some of us Xu et al. (2019); Wang et al. (2019). Our designer Hamiltonian, in the spirit of Eq. (7), is a 2D quantum rotor model which couples Dirac fermions with compact $U(1)$ gauge field on 2D square lattice,

[TABLE]

where $N_{f}$ is the number of fermion flavors, $\hat{L}_{ij}$ are canonical angular momentum and $\hat{\theta}_{ij}$ are coordinate operators, and they satisfy relation $[\hat{L}_{ij},e^{\pm i\hat{\theta}_{ij}}]=\pm e^{\pm i\hat{\theta}_{ij}}$. The Hamiltonian defined in Eq. (33) can be formulated in a path integral with action

[TABLE]

where

[TABLE]

This model is sign problem free when number of fermion flavors $N_{f}$ is even due to pseudo-unitary symmetry which guarantee each determinant is a real number (see Appendix A. in Ref. Xu et al. (2019) for the rigorous proof). We focus on $K=t=1$ and explore the phase diagram when turning $J$ for different $N_{f}$ as showed in Fig. 13. We studied cases when $N_{f}=2$, 4, 6 and 8 respectively, and discovered the $U(1)$ deconfined phase (U1D) and confinement transition in each case. The properties of U1D phase are consistent with the proposal of algebraic spin liquid, where various competing orders (such as AFM order, VBS order etc.) all have identical power-laws algebraic correlation in real space. Importantly, we found the decaying power match perfectly with the larger-$N_{f}$ perturbative renormalization group expression Rantner and Wen (2002); Hermele et al. (2005, 2007); Xu and Sachdev (2008). The continuous confined transitions from U1D to AFM or VBS we found should be described by QED3-Gross-Neveu O(2) or O(3) universality, depending on the symmetry group that the fermion bilinears break in the confined phase, and further carefully study of the critical properties of these transitions via QMC simulations and analytical calculations is certainly worthwhile Wang et al. (2019). Recently perturbative renormalization group calculations to higher orders have been carried out in attempt to accquire the critical properties of the deconfinement to confinement transition in form of QED3-Gross-Neveu universality classes Janssen and He (2017); Ihrig et al. (2018); Zerf et al. (2018); Gracey (2018); Boyack et al. (2019); Zerf et al. (2019) and to address the relevance or irrelevance of the monopole operators to the stability of the U1D phase Song et al. (2018a, b); Dupuis et al. (2019), show substantial interests and great ongoing efforts along this direction.

The discovery of the stable U1D phase at low fermion flavor numbers (possibly starting from $N_{f}=2$) is of vital importance as the situation corresponds to the experimental relevant case. Further more, the phase diagram is of importance, and we believe that it will trigger a number of future investigations on the numerical and analytical fronts, both in the solid state and the high energy communities.

VI Discussion and prospectives

In this review, we summarize the recent progress in developing and employing unbiased quantum Monte Carlo techniques, to investigate the fermion quantum criticalities that can be realized in the designer models such as Eqs. (1),  (7) and  (LABEL:eq:specificformula). In earlier parts of this review, we provided a highly condensed summary about the long history and the intensive theoretical efforts that have been devoted to understand these highly interesting quantum systems, as well as key challenges and difficulties that hold us from obtaining full and final solutions. Nevertheless, great theoretical progress has been made since the HMM framework, and new concepts in the RG and in nFL, hot spots, anomalous transport and quantum critical scaling behaviors prevail.

On the other hand, it took the numerical community quite some time to understand and appreciate these modern theoretical developments and finally is able to incorporate them into the developments of numerical methodologies, exemplified here in the large-scale quantum Monte Carlo simulations. Instead of the treating explicit four-fermion interactions like the Hubbard interaction, which was the focus of the numerical community since the early date of DQMC, in the new designer models discussed in this review, one can directly couple the FS to various critical bosonic modes driven by the consideration of the physical question one would like address, ferromagnetic, antiferromagnetic, nematic critical bosonic fluctuations and QCPs, just to name a few. In short, one can design Hamiltonians and to directly simulate the situations that have been asked in the field theoretical analyses/hypotheses mentioned above. Such awareness has greatly liberated the mindset of the numerical community and new models and results from QMC simulations thus flourishes. Along this process, better numerical methods are also invented from the mutual inspiration and dialogues between numerical and theoretical communities, SLMC, EMUS are the successfully examples.

The works presented here, contain the fermionic QCP in 2D systems on square, triangle lattices, with the bosonic fluctuations in ferromagnetic, antiferromagnetic forms acquired Ising symmetry. We found that in the FM-QCP case, the entire FS becomes critical and the quasiparticle vanishes at the critical FS, the quantum critical scaling clearly acquires finite anomalous dimension different from the 3D Ising universality and that of the HMM prediction, the Lorentz symmetry between the space and time is also lost at the FM-QCP. All this evidence is pointing to a new and unknown universality. Our results provide the valuable references for the further development in the field theoretical treatment to this difficult problem.

As for the AFM-QCPs, we found an interesting difference between the $2\mathbf{Q}=\Gamma$ and $3\mathbf{Q}=\Gamma$ situations, where the former is also clearly deviated from the HMM and acquires the anomalous dimension, the latter is actually consistent with the HMM prediction which enhances the dynamical exponent from $z=1$ in the bare boson case to the $z=2$ in the coupled case. The $2\mathbf{Q}=\Gamma$ is unique in the sense that this is for the first time, an finite anomalous dimension with $\eta=0.125$ is discovered in the unbiased numerical simulations and shining light on the correctness in the previous RG calculation that the $\eta$ might be obtained from perturbative calculation as a function of $N_{\text{h.s.}}$. However, the discrepancy still exists and this might due to the fact that the existing theoretical calculations are performed with $O(3)$ bosons whereas our results are for $Z_{2}$ bosons. Nevertheless, it is clear that from here, both theoretical and numerical communities can adjust their models and try to perform the first exact comparison with exactly the same low-energy physical condition and to see whether the same fixed point can be reached.

Above are the main content of the review, and in Sec. V, we also summarized the recent activities in other related directions such as the fermions coupled to phonons and fermions coupled to gauge degree of freedoms. The progresses in these directions, especially the later ones, go even beyond the boundary of the condensed matter physics and are bridging out to the fundamental question such as the existence of deconfinement in the QED3 in high energy physics. The existence of the $\mathbb{Z}_{2}$ and $U(1)$ gauge field theory with fermion modes provide the concrete examples of the unconventional metallic phases which might stem from the fractionalization of electrons into anyons with associated emergent gauge fields in many strongly correlated systems that are the breeding ground of the quantum state of matter beyond conventional paradigms of FL and spontaneous symmetry breaking. Examples include topological ordered states and quantum spin liquids, Cuprate high-temperature superconductors, etc.

In addition to results summarized above, this progress in model design and numerical techniques offers a platform to explore other challenging problems in strongly correlated systems. In particular, this platform finally enables direct comparison between unbiased numerical simulations and theoretical predictions, and thus a closed feedback loop between theory and numerical studies can be formed. We can now utilize unbiased numerical results to verify and refine our theoretical knowledge, and the progress on the theory side will in return provide more guidance for numerical exploration. For classical critical phenomena, such a feedback loop has played a vital role since decades ago and has resulted in highly fruitful results and achievements. In quantum systems, we have good reasons to believe that through the same type of theory-numerical iterations, we are getting closer than ever before to the fully understanding of strongly correlated phenomena, such as fermionic quantum criticality and beyond. From a humble beginning, an odyssey of discovery has lead us to this stage, that based on the efforts mentioned in this review, the numerical community and theoretical one can finally sit down and start to compare each others results. From this point, one could expect the solution of the fermionic QCPs, or, put in the better terms, the concrete framework of the fermionic QCPs beyond that of the HMM, as concrete as that of the Landau-Ginzburg or Wilson-Fisher, are expected to be reached in the near future.

Acknowledgements.

We would like to express our gratitude towards Fakher Assaad, Yoni Schattner, Erez Berg, Richard Scalettar, George Batrouni, Martin Hohenadler, Simon Trebst and Snir Gazit upon the illustrative discussions and kind support from the numeric perspectives throughout the simulations of the works presented in this review. Likewise, we would like to acknowledge Andrey Chubukov, Subir Sachdev, Steven Kivelson, Max Metlitski, Sung-Sik Lee and Michael Scherer on the very helpful discussions and education on theoretical understanding of modern developments of the itinerant quantum critcalities. XYX, ZHL, GPP and ZYM acknowledge the supports from the Ministry of Science and Technology of China through the National Key Research and Development Program (2016YFA0300502), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB28000000), and the National Science Foundation of China (11421092,11574359,11674370). XYX is also thankful for the support of HKRGC through grant C6026-16W. YQ acknowledges support from Minstry of Science and Technology of China under grant numbers 2015CB921700, and from National Science Foundation of China under grant number 11874115. K.S. acknowledges support from the National Science Foundation under Grant No. EFRI-1741618 and the Alfred P. Sloan Foundation. We thank the Center for Quantum Simulation Sciences in the Institute of Physics, Chinese Academy of Sciences, the Tianhe-1A platform at the National Supercomputer Center in Tianjin and Tianhe-2 platform at the National Supercomputer Center in Guangzhou for their technical support and generous allocation of CPU time.