Fermi-liquid ground state of interacting Dirac fermions in two dimensions

TL;DR

This study uses large-scale quantum Monte Carlo simulations to demonstrate that the semimetallic phase of the 2D Hubbard model on a honeycomb lattice exhibits Fermi-liquid behavior, with a finite quasiparticle weight below a critical interaction.

Contribution

It provides the first large-scale numerical evidence confirming Fermi-liquid behavior in the 2D correlated semimetallic phase of the Hubbard model.

Findings

Quasiparticle weight remains finite below critical interaction.

Green's function exhibits algebraic decay consistent with Fermi-liquid theory.

Numerical simulations support Fermi-liquid description of 2D correlated metals.

Abstract

An unbiased zero-temperature auxiliary-field quantum Monte Carlo method is employed to analyze the nature of the semimetallic phase of the two-dimensional Hubbard model on the honeycomb lattice at half filling. It is shown that the quasiparticle weight $Z$ of the massless Dirac fermions at the Fermi level, which characterizes the coherence of zero-energy single-particle excitations, can be evaluated in terms of the long-distance equal-time single-particle Green's function. If this quantity remains finite in the thermodynamic limit, the low-energy single-particle excitations of the correlated semimetallic phase are described by a Fermi-liquid-type single-particle Green's function. Based on the unprecedentedly large-scale numerical simulations on finite-size clusters containing more than ten thousands sites, we show that the quasiparticle weight remains finite in the semimetallic phase below a critical interaction strength. This is also supported by the long-distance algebraic behavior ($\sim r^{-2}$, where $r$ is distance) of the equal-time single-particle Green's function that is expected for the Fermi liquid. Our result thus provides a numerical confirmation of Fermi-liquid theory in two-dimensional correlated metals.

Equations

$$\hat{H}={\hat{H}}_t+{\hat{H}}_U,$$

$${\hat{H}}_t$$

$${\hat{H}}_U=U\sum _i\sum _{\alpha =A,B}{\hat{n}}_{\alpha ,{\mathbf{r}}_i,\uparrow }{\hat{n}}_{\alpha ,{\mathbf{r}}_i,\downarrow }.$$

$$⟨\hat{\mathcal{O}}⟩=⟨{\Psi }_0\left|\hat{\mathcal{O}}\right|{\Psi }_0⟩=\underset{\tau \to \infty }{\lim }\frac{⟨{\Psi }_{\mathrm{T}}|e^{-\tau \hat{H}/2}\hat{\mathcal{O}}{e}^{-\tau \hat{H}/2}|{\Psi }_{\mathrm{T}}⟩}{⟨{\Psi }_{\mathrm{T}}|e^{-\tau \hat{H}}|{\Psi }_{\mathrm{T}}⟩},$$

$$e^{-{\Delta }_{\tau }{\hat{H}}_U}=C\sum _{s_{A,1}}\sum _{s_{B,1}}\cdots \sum _{s_{B,N_{\mathrm{cell}}}}\exp \left[\lambda \sum _{\alpha ,i}{s}_{\alpha ,i}\left({\hat{n}}_{\alpha ,{\mathbf{r}}_i,\uparrow }-{\hat{n}}_{\alpha ,{\mathbf{r}}_i,\downarrow }\right)\right]$$

$$e^{\pm {\Delta }_{\tau }{\mathbf{H}}_t}\approx {\mathcal{I}}_0(\rho {\Delta }_{\tau })\mathbf{I}+2\sum _{k=1}^M(\pm 1{)}^k{\mathcal{I}}_k(\rho {\Delta }_{\tau })T_k({\tilde{\mathbf{H}}}_t),$$

$${\left|\left|e^{\pm {\Delta }_{\tau }{\mathbf{H}}_t}-\left({\mathcal{I}}_0(\rho {\Delta }_{\tau })\mathbf{I}+2\sum _{k=1}^M(\pm 1{)}^k{\mathcal{I}}_k(\rho {\Delta }_{\tau })T_k({\tilde{\mathbf{H}}}_t)\right)\right|\right|}_{\max }<ϵ,$$

$$D_{AB,\sigma }(\mathbf{r})$$

$${\hat{H}}_t$$

$$D_{AB,\sigma }^{(0)}(\mathbf{r})=\frac{1}{N_{\mathrm{cell}}}\sum _{\mathbf{k}}⟨{\hat{c}}_{B,\mathbf{k},\sigma }^{†}{\hat{c}}_{A,\mathbf{k},\sigma }⟩e^{i\mathbf{k}\cdot \mathbf{r}}$$

$$\underset{T\to 0}{=}$$

$$D_{AA,\sigma }^{(0)}(\mathbf{r})=\frac{1}{N_{\mathrm{cell}}}\sum _{\mathbf{k}}⟨{\hat{c}}_{A,\mathbf{k},\sigma }^{†}{\hat{c}}_{A,\mathbf{k},\sigma }⟩e^{i\mathbf{k}\cdot \mathbf{r}}=\frac{1}{2}{\delta }_{\mathbf{r},0},$$

$$D_{AB,\sigma }(\mathbf{r})$$

$${\mathbf{G}}_{\sigma }(\mathbf{k},z)$$

$$=$$

$$D_{AB,\sigma }(\mathbf{r})$$

$$Z=\underset{|\mathbf{r}|\to \infty }{\lim }\frac{D_{AB,\sigma }(\mathbf{r})}{D_{AB,\sigma }^{(0)}(\mathbf{r})}.$$

$$D_{AB,\sigma }(r)\sim \frac{1}{r^2}$$

$$Z(L)=\frac{D_{AB,\sigma }({\mathbf{r}}_{\max })}{D_{AB,\sigma }^{(0)}({\mathbf{r}}_{\max })},$$

$$\frac{1}{N_{\mathrm{cell}}}\sum _{\mathbf{k}}\cdots \to \frac{S_{\mathrm{cell}}}{(2\pi {)}^2}\int d^2\mathbf{k}\cdots ,$$

$$D_{AB,\sigma }^{(0)}(\mathbf{r})$$

$$\mathbf{k}={\mathbf{k}}_v+\mathbf{q}.$$

$$h_{{\mathbf{k}}_K+\mathbf{q}}$$

$$-\frac{1}{2}\frac{S_{\mathrm{cell}}}{(2\pi {)}^2}{e}^{i{\mathbf{k}}_K\cdot \mathbf{r}}\int d^2\mathbf{q}\frac{i{q}_x+q_y}{q}{e}^{i\mathbf{q}\cdot \mathbf{r}}$$

$$\frac{1}{2\pi }{\int }_0^{\Lambda }dq{\int }_0^{2\pi }d\varphi e^{iqr\cos \varphi }=\frac{1}{r}{\int }_0^{r\Lambda }ds{J}_0(s),$$

$$\frac{1}{2\pi }\int d^2\mathbf{q}\frac{1}{q}{e}^{i\mathbf{q}\cdot \mathbf{r}}=\frac{1}{r}.$$

$$h_{{\mathbf{k}}_{K^{\prime }}+\mathbf{q}}$$

$$-\frac{1}{2}\frac{S_{\mathrm{cell}}}{(2\pi {)}^2}{e}^{i{\mathbf{k}}_{K^{\prime }}\cdot \mathbf{r}}\int d^2\mathbf{q}\frac{i{q}_x-q_y}{q}{e}^{i\mathbf{q}\cdot \mathbf{r}}$$

$$D_{AB,\sigma }^{(0)}(\mathbf{r})\simeq -\frac{S_{\mathrm{cell}}}{4\pi }\left(e^{i{\mathbf{k}}_K\cdot \mathbf{r}}\frac{r_x-i{r}_y}{r^3}+e^{i{\mathbf{k}}_{K^{\prime }}\cdot \mathbf{r}}\frac{r_x+i{r}_y}{r^3}\right).$$

$$D_{AB,\sigma }^{(0)}\left(n({\mathbf{\tau }}_1+{\mathbf{\tau }}_2)\right)\simeq -\frac{S_{\mathrm{cell}}}{2\pi }\frac{1}{r^2}$$

$$D_{AB,\sigma }\left(n({\mathbf{\tau }}_1+{\mathbf{\tau }}_2)\right)\simeq -Z\frac{S_{\mathrm{cell}}}{2\pi }\frac{1}{r^2}.$$

Full text

I Introduction

The characterization of different phases of matter is one of the essential issues in solid state physics. In the field of strongly correlated electrons, the correlation-induced metal-insulator transition Imada et al. (1998) is of particular importance since the itinerancy and localization of electrons Kohn (1964); Resta and Sorella (1999) can be regarded as a many-electron realization of the wave-particle duality, the fundamental concept of quantum mechanics.

The Hubbard model Gutzwiller (1963); Kanamori (1963); Hubbard (1963) is certainly one of the most important models in condensed matter physics since it has inspired many ideas and led to milestone achievements for understanding the fascinating properties of the metal-insulator transition. In particular, a semimetal-insulator transition occurs in the Hubbard model in a certain class of lattices where massless Dirac-like dispersion appears in the noninteracting limit, and has been therefore a subject of intense activity in recent years. Since such models can be constructed on bipartite lattices and thereby they are free from the negative-sign problem, the numerically exact auxiliary-field quantum Monte Carlo (AFQMC) method has played a major role in the study of this semimetal-insulator transition. In order to determine the ground-state phase diagram, most of the previous calculations have focused on the order parameters in the insulating phase, including the single-particle excitation gap and the antiferromagnetic spin-structure factor Meng et al. (2010); Sorella et al. (2012); Chang and Scalettar (2012); Assaad and Herbut (2013). Variants of such models have been further extended recently by coupling interacting Dirac fermions to Ising spins Sato et al. (2017) or by introducing disordered transfer integrals Ma et al. (2018).

On the theoretical side, the Green’s-function-based formalism Migdal (1957); Luttinger (1961); Nozières and Luttinger (1962); Luttinger and Nozières (1962); Abrikosov et al. (1975) of the Fermi-liquid theory Landau (1956) argues that one of the most important characteristics in a correlated metallic state is the quasiparticle weight $Z$ at the Fermi level, because finite $Z$ implies the existence of coherent zero-energy single-particle excitations. Although massless Dirac fermions exhibit only Fermi points instead of full Fermi surfaces, the quasiparticle weight $Z$ remains well defined Herbut et al. (2009), despite that the low-energy single-particle excitations and the electronic transport can be substantially different from those in simple metals Castro Neto et al. (2009); Das Sarma et al. (2011). In principle, $Z$ can be estimated from the imaginary-time-displaced single-particle Green’s function at the Dirac point with the AFQMC method Feldbacher and Assaad (2001); Assaad and Evertz (2008). However, the computation of imaginary-time-displaced quantities is considerably more expensive and suffers from much larger signal-to-noise ratio than the corresponding equal-time correlations. This is probably the main reason for preventing the calculation of $Z$ in the semimetallic phase with the AFQMC technique. In this regard, recently, three of us Otsuka et al. (2016) elucidated the quantum criticality emerging from the continuous semimetal-insulator transition, with large-scale zero-temperature AFQMC simulations Sorella et al. (1988, 1989); Sorella and Tosatti (1992); Sorella et al. (2012); Becca and Sorella (2017). However, no direct and systematic calculation of the quasiparticle weight for interacting Dirac fermions has been reported yet. It should also be noted that, in spite of the recent development of various numerical techniques and the continuous improvement of computer performances, a solid numerical evidence of the presence of quasiparticles and, by consequence, a clear validation of the Fermi-liquid theory, are still lacking for interacting fermions on any two-dimensional lattices.

In this paper, we first show that the quasiparticle weight $Z$ of the massless Dirac fermions at the Fermi level can be evaluated from the ratio of the interacting and noninteracting equal-time single-particle Green’s functions in the long-distance limit. The scheme is then demonstrated with the unbiased zero-temperature AFQMC simulation for the Hubbard model on unprecedentedly large finite-size clusters of the honeycomb lattice at half filling. Based on the numerical results for the quasiparticle weight, we address a fundamental and long-standing issue: whether the Fermi liquid can be realized in two spatial dimensions Varma et al. (1989); Anderson (1990, 1991). Our result implies that the Fermi-liquid picture is valid in the correlated semimetallic phase.

The rest of the paper is organized as follows. In Sec. II, we define the Hubbard model on the honeycomb lattice and describe the AFQMC method. In Sec. III, based on the Fermi-liquid theory, we show that the quasiparticle weight $Z$ of interacting massless Dirac fermions is calculated from the equal-time single-particle Green’s function. In Sec. IV, we provide the numerical results which strongly support the Fermi-liquid behavior in the semimetallic phase. In Sec. V, we summarize the paper and discuss the non-Fermi-liquid behavior in graphene. In Appendixes A and B, we analyze the long-distance behavior of the equal-time Green’s function in the semimetallic and insulating phases, respectively.

II Model and Method

II.1 Hubbard model on the honeycomb lattice

The Hamiltonian of the Hubbard model on the honeycomb lattice is given by

[TABLE]

where

[TABLE]

and

[TABLE]

Here, $\hat{c}_{\alpha,\boldsymbol{r}_{i},\sigma}^{\dagger}$ ($\hat{c}_{\alpha,\boldsymbol{r}_{i},\sigma}$) is a creation (annihilation) operator of a fermion at unit cell $i$, located at $\boldsymbol{r}_{i}=n_{i}^{(1)}\boldsymbol{\tau}_{1}+n_{i}^{(2)}\boldsymbol{\tau}_{2}$ (where $n_{i}^{(1)}$ and $n_{i}^{(2)}$ are integer), and sublattice $\alpha\,(=A,B)$ with spin $\sigma\,(=\uparrow,\downarrow)$, and $\hat{n}_{\alpha,\boldsymbol{r}_{i},\sigma}=\hat{c}_{\alpha,\boldsymbol{r}_{i},\sigma}^{\dagger}\hat{c}_{\alpha,\boldsymbol{r}_{i},\sigma}$ (see Fig. 1). $t$ is the hopping integral between the nearest-neighbor sites of the honeycomb lattice and $U$ is the strength of the on-site interaction. In this paper, we consider fermion density $n_{f}=1$, i.e., half filling, for which the Dirac points are located exactly at the Fermi level in the noninteracting limit.

Figure 1 shows the honeycomb lattice spanned by primitive translational vectors $\boldsymbol{\tau}_{1}=a(\frac{3}{2},\frac{\sqrt{3}}{2})$ and $\boldsymbol{\tau}_{2}=a(\frac{3}{2},-\frac{\sqrt{3}}{2})$ with $a$ being the lattice constant. A finite-size cluster of the linear dimension $L$ is defined by the two vectors $L\boldsymbol{\tau}_{1}$ and $L\boldsymbol{\tau}_{2}$, containing $N_{\rm cell}=L^{2}$ unit cells and hence $N_{\rm site}=2N_{\rm cell}=2L^{2}$ sites. We choose the clusters of $L=$ 8, 14, 20, 26, 32, 38, 44, 50, 62, and 74 under periodic boundary conditions, for which the closed-shell condition is satisfied Sorella (2015). The maximum size considered here thus contains 10952 sites, which is substantially (more than four times) larger than the previous largest AFQMC simulations of the two-dimensional Hubbard models Otsuka et al. (2016); Sorella et al. (2012).

II.2 Auxiliary-field quantum Monte Carlo method

We study the ground-state properties of the Hubbard model $\hat{H}$ with the zero-temperature AFQMC method Sugiyama and Koonin (1986); Sorella et al. (1988, 1989); Becca and Sorella (2017), where the ground-state expectation value of an operator $\hat{\mathcal{O}}$ is evaluated as

[TABLE]

where $|\Psi_{0}\rangle$ is the normalized ground state of $\hat{H}$, $\tau\geqslant 0$ is the projection time, and $|\Psi_{\rm T}\rangle$ is a trial wavefunction such that $\langle\Psi_{0}|\Psi_{\rm T}\rangle\not=0$. We choose as $|\Psi_{\rm T}\rangle$ the ground state of $\hat{H}_{t}$, i.e., the Fermi sea.

The imaginary-time evolution is performed with the second-order Trotter-Suzuki decomposition ${e}^{-\tau\hat{H}}=\prod_{l=1}^{N_{\tau}}({e}^{-\Delta_{\tau}\hat{H}_{t}/2}{e}^{-\Delta_{\tau}\hat{H}_{U}}{e}^{-\Delta_{\tau}\hat{H}_{t}/2})+O(\Delta_{\tau}^{2}),$ where $\tau$ is discretized into $N_{\tau}$ time slices with an interval $\Delta_{\tau}=\tau/N_{\tau}$ and $O(\Delta_{\tau}^{2})$ is the systematic error due to the imaginary-time discretization Trotter (1959); Suzuki (1976). At each time slice $l$, the discrete version of the Hubbard-Stratonovich transformation

[TABLE]

is applied, where $s_{\alpha,i}=\pm 1$ is the auxiliary field on sublattice $\alpha$ of the unit cell at $\boldsymbol{r}_{i}$, $\cosh(\lambda)={e}^{\Delta_{\tau}U/2}$, and $C=(e^{-\Delta_{\tau}U/4}/2)^{2L^{2}}$ Hirsch (1983); Hubbard (1959); Stratonovich (1957). When this equation is used to evaluate the full propagator $\prod\limits_{l}\exp(-\Delta_{\tau}\hat{H})$, an explicit imaginary-time ($l$) dependence of the field $s_{\alpha,i}=s_{\alpha,i}(l)$ appears in each time slice, according to Eq. (5). The multiple summation over $\{s_{\alpha,i}(l)\}$ is performed by the Monte Carlo method with the importance sampling. The negative-sign problem does not arise at half filling owing to the particle-hole symmetry Hirsch (1985). In this study, we set $\Delta_{\tau}t=0.1$ without attempting the extrapolation $\Delta_{\tau}\to 0$ because it already provides a satisfactory accuracy ($<2\%$) in all correlation functions studied. Large enough projection times $\tau t=50$ or equivalently $N_{\tau}=500$ ($\tau t=80$ or equivalently $N_{\tau}=800$) for clusters of $L\leqslant 20$ ($L\geqslant 26$) are used to obtain the converged $\tau\to\infty$ results in Eq. (4).

II.3 Sparse-matrix exponential

One of the most computationally expensive operations in the AFQMC method for large clusters is the multiplication of $e^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}$ (or $e^{\pm\Delta_{\tau}\boldsymbol{H}_{t}/2}$) to the wavefunction matrix or to the Green’s function matrix, where $\boldsymbol{H}_{t}$ is the (real-space) matrix representation of $\hat{H}_{t}$. Usually, $e^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}$ is treated as an $N_{\rm site}\times N_{\rm site}$ dense matrix with the spectral decomposition $e^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}=\boldsymbol{U}^{T}e^{\pm\Delta_{\tau}\boldsymbol{D}}\boldsymbol{U}$, where $\boldsymbol{U}$ is a $N_{\rm site}\times N_{\rm site}$ orthogonal matrix that diagonalizes $\boldsymbol{H}_{t}$, i.e., $\boldsymbol{H}_{t}\boldsymbol{U}=\boldsymbol{U}\boldsymbol{D}$. Although $e^{\pm\Delta_{\tau}\boldsymbol{D}}$ is diagonal, $\boldsymbol{U}$ is generally dense and thus $e^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}$ is dense. Therefore, the computational cost of the matrix-matrix multiplication scales as $O(N_{\rm site}^{3})$. Here, we describe an alternative multiplication scheme of $e^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}$ which is efficient for large clusters by taking full advantage of the sparseness of $\boldsymbol{H}_{t}$.

In this scheme, we expand the matrix exponential as a polynomial of degree $M$, i.e.,

[TABLE]

where $\boldsymbol{I}$ is the identity matrix, ${\mathcal{I}}_{k}(\rho\Delta_{\tau})$ is the $k$th order modified Bessel function of the first kind, $\rho$ is the spectral radius of $\boldsymbol{H}_{t}$ ($\rho=3|t|$ in the present case), and $\tilde{\boldsymbol{H}}_{t}=\boldsymbol{H}_{t}/\rho$. $T_{k}(\tilde{\boldsymbol{H}}_{t})$ is the $k$th order Chebyshev polynomial of the first kind, which can be obtained iteratively as $T_{0}(\tilde{\boldsymbol{H}}_{t})=\boldsymbol{I}$, $T_{1}(\tilde{\boldsymbol{H}}_{t})=\tilde{\boldsymbol{H}}_{t}$, and $T_{k}(\tilde{\boldsymbol{H}}_{t})=2\tilde{\boldsymbol{H}}_{t}T_{k-1}(\tilde{\boldsymbol{H}}_{t})-T_{k-2}(\tilde{\boldsymbol{H}}_{t})$ for $k\geqslant 2$ Tal-Ezer and Kosloff (1984); Vijay and Metiu (2002); Iitaka and Ebisuzaki (2003); Weiße and Fehske (2008). A similar orthogonal-polynomial expansion of the Boltzmann factor with the Legendre polynomial has been employed in a finite-temperature dynamical density-matrix-renormalization-group method Sota and Tohyama (2008). As shown below, we find that, for large $N_{\rm site}$, the multiplication of ${e}^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}$ with manipulating $\boldsymbol{H}_{t}$ as a sparse matrix in the right-hand side of (6) is faster than the direct multiplication of the dense matrix ${e}^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}$, even when machine accuracy is reached with large enough $M$.

Figure 2(a) shows the computational time of one space-time Monte Carlo sweep with the two multiplication schemes for fixed $N_{\rm\tau}=100$, $\Delta_{\tau}t=0.1$, and $U/t=3.5$. The same initial auxiliary field configuration $\{s_{\alpha,i}(l)\}$ with the same random seed for the same random number generator is used for both schemes. The stabilization (i.e., orthonormalization) of the wavefunction Sorella et al. (1989); Imada and Hatsugai (1989); White et al. (1989) is made every 10 time slices. $M=8$ ($M=7$) is used for the expansion with $\pm\Delta_{\tau}t$ ($\Delta_{\tau}t/2$) to achieve an accuracy of $<10^{-13}$ (see below). Since $\boldsymbol{H}_{t}$ has only $z_{\rm c}$ ($z_{\rm c}$: the coordination number, i.e., $z_{\rm c}=3$ for the honeycomb lattice) nonzero matrix elements in each column and row (thus, totally $z_{\rm c}N_{\rm site}$ nonzero elements), the computational cost of the multiplication of $e^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}$ to an $N_{\rm site}\times N_{\rm site}$ dense matrix scales as $O(z_{\rm c}MN_{\rm site}^{2})$ when the polynomial expansion scheme in Eq. (6) is employed. A convenient speedup larger than one is achieved for $N_{\rm site}\gtrsim 2000$ in our computing environment and increases with $N_{\rm site}$. Modern processors have the possibility to perform several independent tasks, called ”threads” within the same computational unit. As shown in Fig. 2(b), the speedup with threading is also as effective as that in the dense-matrix case. Here, the compressed-row-storage (CRS) format (see for example Ref. Hager and Wellein (2008)) is used to store the nonzero matrix elements of $\boldsymbol{H}_{t}$ within the polynomial expansion scheme.

In analogy with the high-temperature series expansion Imada and Takahashi (1986); Jaklič and Prelovšek (2000), the convergence of the polynomial expansion in Eq. (6) with relatively small $M$ is evident because usually $\Delta_{\tau}$ is taken small ($\Delta_{\tau}t\ll 1$) in the AFQMC simulation. Given a desired accuracy $\epsilon$ for the polynomial expansion, $M$ can be determined to satisfy

[TABLE]

where $||\boldsymbol{A}||_{\rm max}=\max_{ij}|A_{ij}|$ is the maximum norm of $\boldsymbol{A}$ and $A_{ij}=[\boldsymbol{A}]_{ij}$. We set $\epsilon=10^{-13}$ and find that $M=8$ ($M=7$) is the minimum value that satisfies the inequality in Eq. (7) for $\Delta_{\tau}t=0.1$ ($\Delta_{\tau}t/2=0.05$ when $e^{\pm\Delta_{\tau}\boldsymbol{H}_{t}/2}$ is expanded), irrespectively of the system size. This implies that the polynomial expansion is well controlled and even does not introduce the additional systematic error by terminating the expansion at finite $M$ as $\epsilon$ is negligibly smaller than the statistical error. Finally, we note that, if the degree $M$ is the same, the Chebyshev polynomial expansion in Eq. (6) gives better accuracy than the Taylor expansion ${e}^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}\approx\sum_{k=0}^{M}\frac{(\pm\Delta_{\tau})^{k}}{k!}\boldsymbol{H}_{t}^{k}$, in the sense that the matrix norm of the difference from the exact ${e}^{\pm\Delta_{\tau}\boldsymbol{H}_{t}}$ is smaller, for the model studied here.

Recently, a different approach to reduce the computational effort of fermionic quantum Monte Carlo (QMC) simulations, dubbed as effective momentum ultra-size QMC, has been proposed and successfully used in some model systems Liu et al. (a, b). This approach is designed to capture the low-energy physics for original lattice models of interest.

III Quasiparticle weight

The main quantity considered here is the equal-time single-particle Green’s function

[TABLE]

where $\boldsymbol{r}$ denotes a relative spatial position of two unit cells at $\boldsymbol{r}^{\prime}$ and $\boldsymbol{r}^{\prime}+\boldsymbol{r}$, and the average $\langle\cdots\rangle={{\rm Tr}[{e}^{-\hat{H}/T}\cdots]}/{{\rm Tr}[{e}^{-\hat{H}/T}]}$ is defined at a finite temperature $T$ for the clarity of the following formulation. The zero-temperature limit will be taken only at the end of the calculation. Since $D_{AB,\sigma}(\boldsymbol{r})$ represents the probability amplitude that a hole created on sublattice $A$ in the unit cell at $\boldsymbol{r}^{\prime}$ propagates to sublattice $B$ in the unit cell at $\boldsymbol{r}+\boldsymbol{r}^{\prime}$, the long-distance behavior of $D_{AB,\sigma}(\boldsymbol{r})$ should enable us to distinguish whether the system is semimetallic or insulating. Indeed, as shown in Appendixes A and B, $D_{AB,\sigma}(\boldsymbol{r})$ decays algebraically, with a prefactor proportional to $Z$, in the semimetallic phase, while it decays exponentially in the insulating phase.

III.1 Noninteracting limit

First, we analyze $D_{AB,\sigma}(\boldsymbol{r})$ in the noninteracting limit. For this purpose, we diagonalize $\hat{H}_{t}$ as

[TABLE]

where $h_{\boldsymbol{k}}=t\left(1+{e}^{-{i}\boldsymbol{k}\cdot\boldsymbol{\tau}_{1}}+{e}^{-{i}\boldsymbol{k}\cdot\boldsymbol{\tau}_{2}}\right)$, $\hat{\psi}_{+,\boldsymbol{k},\sigma}=\frac{1}{\sqrt{2}}\left(\hat{c}_{A,\boldsymbol{k},\sigma}+{e}^{{i}\theta_{\boldsymbol{k}}}\hat{c}_{B,\boldsymbol{k},\sigma}\right)$, $\hat{\psi}_{-,\boldsymbol{k},\sigma}=\frac{1}{\sqrt{2}}\left(\hat{c}_{A,\boldsymbol{k},\sigma}-{e}^{{i}\theta_{\boldsymbol{k}}}\hat{c}_{B,\boldsymbol{k},\sigma}\right)$, ${e}^{{i}\theta_{\boldsymbol{k}}}=h_{\boldsymbol{k}}/|h_{\boldsymbol{k}}|$, and $\hat{c}_{\alpha,\boldsymbol{k},\sigma}=N_{\rm cell}^{-1/2}\sum_{i}\hat{c}_{\alpha,\boldsymbol{r}_{i},\sigma}{e}^{-{i}\boldsymbol{k}\cdot\boldsymbol{r}_{i}}$. The bonding- and antibonding-band energies are $-|h_{\boldsymbol{k}}|$ and $|h_{\boldsymbol{k}}|$, respectively. The zero-energy modes protected by the chiral symmetry Semenoff (2012); Hatsugai et al. (2013) appear at two inequivalent momenta, $K$ and $K^{\prime}$ points, which are specified by the vectors $\boldsymbol{k}_{K}=\frac{1}{a}(\frac{2\pi}{3},\frac{2\pi}{3\sqrt{3}})$ and $\boldsymbol{k}_{K^{\prime}}=\frac{1}{a}(\frac{2\pi}{3},-\frac{2\pi}{3\sqrt{3}})$, respectively. $D_{AB,\sigma}(\boldsymbol{r})$ in the noninteracting limit is now evaluated as

[TABLE]

where the superscript “$(0)$” denotes that the quantity is in the noninteracting limit. $n_{\rm F}(E)=1/({e}^{E/T}+1)$ is the Fermi distribution function, which arises from the occupation of the fermions $\left\langle\hat{\psi}^{\dagger}_{\pm,\boldsymbol{k},\sigma}\hat{\psi}_{\pm,\boldsymbol{k},\sigma}\right\rangle=n_{\rm F}(\pm|h_{\boldsymbol{k}}|)$. The summand in Eq. (10) exactly at the $K$ and $K^{\prime}$ points is zero because $n_{\rm F}(|h_{\boldsymbol{k}_{K(K^{\prime})}}|)-n_{\rm F}(-|h_{\boldsymbol{k}_{K(K^{\prime})}}|)=0$ and thereby these two momenta are excluded from the summation in Eq. (11).

We should note that, on the contrary to $D_{AB,\sigma}^{(0)}(\boldsymbol{r})$, $D_{AA,\sigma}^{(0)}(\boldsymbol{r})$ at half filling gives merely a trivial $\boldsymbol{r}$ dependence, i.e.,

[TABLE]

because $\left\langle\hat{c}_{A,\boldsymbol{k},\sigma}^{\dagger}\hat{c}_{A,\boldsymbol{k},\sigma}\right\rangle=1/2$. Here, $\delta_{\boldsymbol{r},\boldsymbol{0}}=1$ when $\boldsymbol{r}=\boldsymbol{0}$ and zero otherwise. This is also the case when the interaction $U$ is finite because $\langle\hat{c}_{A\boldsymbol{k}\sigma}^{\dagger}\hat{c}_{A\boldsymbol{k}\sigma}\rangle=1/2$ as long as the particle-hole symmetry is preserved. Therefore, $D_{AA,\sigma}(\boldsymbol{r})$ and similarly $D_{BB,\sigma}(\boldsymbol{r})$ do not show any long-distance propagation of a hole that can discriminate the nature of the different ground states.

III.2 Interacting case

In order to analyze $D_{AB,\sigma}(\boldsymbol{r})$ in an interacting system, we now express this quantity with the single-particle Green’s function $G_{AB,\sigma}(\boldsymbol{r},{i}\omega_{\nu})$ in the Matsubara-frequency representation Ezawa et al. (1957); Matsubara (1955), i.e.,

[TABLE]

where ${i}\omega_{\nu}=(2\nu+1)\pi{i}T$ with $\nu$ integer is the fermionic Matsubara frequency, $G_{AB,\sigma}(\boldsymbol{r},{i}\omega_{\nu})=N_{\rm cell}^{-1}\sum_{\boldsymbol{k}}G_{AB,\sigma}(\boldsymbol{k},{i}\omega_{\nu}){e}^{{i}\boldsymbol{k}\cdot{\boldsymbol{r}}}$, and the frequency sum is converted to the contour integral. The contour $\mathcal{C}$ is chosen so as to include all the singularities of $G_{AB,\sigma}(\boldsymbol{k},z)$, which lie on the real axis, and therefore does not enclose the Matsubara frequencies.

We now assume that the single-particle Green’s function near the Fermi level has a Fermi-liquid-type pole Abrikosov et al. (1975), which should be consistent with the particle-hole symmetry of the model, i.e.,

[TABLE]

where $\tilde{h}_{\boldsymbol{k}}=(v_{\rm F}/v_{\rm F}^{(0)})h_{\boldsymbol{k}}$ with $v_{\rm F}$ and $v_{\rm F}^{(0)}(=3|t|a/2)$ being the Fermi velocity of the interacting and noninteracting systems, respectively, and $Z$ is the quasiparticle weight at the nodal Dirac point. The incoherent part is a function of $z$ and the singularities lie well away from the Fermi level.

By substituting $G_{AB,\sigma}(\boldsymbol{k},z)$ of Eq. (19) into Eq. (13) and performing the contour integral, we obtain in the large-distance limit ($|\boldsymbol{r}|/a\gg 1$) that

[TABLE]

Here, the incoherent part does not contribute to $D_{AB,\sigma}(\boldsymbol{r})$ in the long-distance limit. This is because the singularities of the incoherent part appear away from the Fermi level and thus the contribution of the incoherent part to $D_{AB,\sigma}(\boldsymbol{r})$ decays exponentially in $|\boldsymbol{r}|$ (see Appendix B). Note that the $K$ and $K^{\prime}$ points are excluded from the summation in Eq. (20), as in the noninteracting case. This justifies the use of finite-size clusters with $L=3n+2$ (or $L=3n+1$, where $n$ is integer) for our AFQMC simulations, where the closed-shell condition in the noninteracting limit is convenient for accurate simulations Sorella (2015).

The form of $D_{AB,\sigma}(\boldsymbol{r})$ in Eq. (20) is quite natural as it matches the simple substitution of the quasiparticle operators $\hat{c}_{B,\boldsymbol{r},\sigma}^{\dagger}\mapsto\hat{q}_{B,\boldsymbol{r},\sigma}^{\dagger}=\sqrt{Z}\hat{c}_{B,\boldsymbol{r},\sigma}^{\dagger}$ and $\hat{c}_{A,\boldsymbol{r},\sigma}\mapsto\hat{q}_{A,\boldsymbol{r},\sigma}=\sqrt{Z}\hat{c}_{A,\boldsymbol{r},\sigma}$ into $D_{AB,\sigma}^{(0)}(\boldsymbol{r})$ Fabrizio (2007). The quasiparticle weight $Z$ at the Fermi point in the thermodynamic limit is now simply evaluated via the ratio of the equal-time single-particle Green’s functions in the long-distance limit, i.e.,

[TABLE]

Since the Fermi velocity $v_{\rm F}$, another unknown quantity, does not appear here, $Z$ can be estimated independently of $v_{\rm F}$.

IV Numerical results

Employing the AFQMC method, we now examine numerically the long-distance behavior of $D_{AB,\sigma}(\boldsymbol{r})$. As shown in Appendixes A and B, $D_{AB,\sigma}(\boldsymbol{r})$ decays in $r=|\boldsymbol{r}|$ as

[TABLE]

in the Fermi liquid, while $D_{AB,\sigma}(r)$ decays exponentially in the insulating state. Figure 3 shows the cluster-size ($L$) dependence of $D_{AB,\sigma}(r_{\rm max})$ for $U/t=$3.5, 3.7, and 4, where $r_{\rm max}=|\boldsymbol{r}_{\rm max}|$ is the maximum distance available in a given finite-size cluster of linear dimension $L$ (see Fig. 1). We take $\boldsymbol{r}_{\rm max}$ in the $x$ direction to remove the phase factors in $D_{AB,\sigma}(r)$ (for details, see Appendix A). The lines are linear fits to the data of the form $-\alpha(L_{\rm min})\ln L+b(L_{\rm min})$, where $\alpha(L_{\rm min})$ and $b(L_{\rm min})$ are fitting parameters with $L_{\rm min}$ being the minimum $L$ used for the fit. As summarized in Fig. 4, $\alpha(L_{\rm min})$ approaches to $2$ for $U/t=3.5$ and $3.6$, as expected for the Fermi liquid, while $\alpha(L_{\rm min})$ increases with $L_{\rm min}$ for $U/t=3.8$, $3.9$ and $4$, indicating the insulating behavior. Only in the vicinity of the quantum critical point $U_{\mathrm{c}}/t\simeq 3.7$ separating the semimetal and the antiferromagnetic insulator Sorella et al. (2012); Otsuka et al. (2016), we observe the non-Fermi-liquid behavior characterized by the non-trivial exponent of $2+\eta_{\psi}$, where $\eta_{\psi}\simeq 0.2$ Otsuka et al. (2016) is the fermion anomalous dimension 111 Note that the value of the fermion anomalous dimension $\eta_{\psi}$ (as well as other critical exponents) is still controversial as it ranges from $0.071(2)$ to $0.242$, depending on numerical and analytical techniques used Janssen and Herbut (2014); Otsuka et al. (2016); Zerf et al. (2017); Knorr (2018) . Therefore, these results already imply that the semimetallic phase is the Fermi liquid.

Next, we evaluate the quasiparticle weight on finite-size clusters,

[TABLE]

as recently applied by the authors to identify the semimetallic state on a triangular lattice Otsuka et al. (2018). For the Fermi-liquid ground state, the quasiparticle weight in the thermodynamic limit, i.e., $Z=\lim_{L\to\infty}Z(L)$, is finite. Figure 5 shows $Z(L)$ as a function of $1/L$ and lines are second-order polynomial fits of the form $\sum_{n=0}^{2}c_{n}L^{-n}$ to the data with $\{c_{n}\}$ being fitting parameters determined by the least-squares method. The extrapolated values of $c_{0}=Z$ and their error bars in the thermodynamic limit are also shown at $1/L=0$ for the semimetallic phase where the Fermi-liquid-like asymptotic behavior is observed in $D_{AB,\sigma}(r)$ (see Fig. 4). We find that these extrapolated values are consistent, within two standard deviations, with our previous results Otsuka et al. (2016) which are estimated from the jump of the momentum distribution function and indicated by stars in Fig. 5. Our new calculations with Eq. (23) performed on the larger clusters are however more accurate as the error bars are more than six-times smaller, supporting the validity of the Fermi-liquid theory in the semimetallic phase of the Honeycomb lattice.

V Conclusions and discussions

In conclusion, we have shown by the AFQMC method that a Fermi-liquid ground state is realized in the semimetallic phase of the Hubbard model on the honeycomb lattice at half filling. This conclusion is obtained by studying the asymptotic behavior of the equal-time single-particle Green’s function $D_{AB,\sigma}(r)\sim 1/r^{2}$ and by providing firm numerical indication of a finite quasiparticle weight $Z$ in the semimetallic phase. The finite $Z$ immediately implies the presence of the quasiparticles, each of which carries a spin $1\over 2$ and a charge $-e$ (for many electron systems) with the Fermi surface unaltered from the noninteracting one, due to the particle-hole symmetry Luttinger (1960); Seki and Yunoki (2017). In the vicinity of the quantum critical point, the non Fermi liquid behavior characterized with a non-trivial exponent is also probed directly by the asymptotic behavior of $D_{AB,\sigma}(r)$.

Considering the Hubbard model as the minimal model for graphene Schüler et al. (2013), our results imply a realization of Fermi liquid in graphene, which has been often assumed, for example, in Ref. Katsnelson (2008). However, because of the vanishing density of states at half filling, the unscreened long-range Coulomb interactions are certainly important for a more realistic modeling of graphene to examine a possible non-Fermi-liquid behavior accompanied with the diverging Fermi velocity González et al. (1999); Kotov et al. (2012); Ulybyshev et al. (2013); Wu and Tremblay (2014); Tang et al. (2015); Tupitsyn and Prokof’ev (2017); Tang et al. (2018); Buividovich et al. (2018). Indeed, an anomalous increase of the Fermi velocity in graphene has been reported experimentally Elias et al. (2011). The Hubbard-type models with long-range Coulomb interaction Hohenadler et al. (2014) on the honeycomb lattice might be promising to investigate the non-Fermi-liquid state in graphene and also other possible many-body electronic states in carbon-based low-dimensional materials such as condensed excitonic states Phan and Fehske (2012); Varsano et al. (2017).

Acknowledgements.

We acknowledge Tomonori Shirakawa for useful discussions. This work has been supported in part by Grant-in-Aid for Scientific Research from MEXT Japan (under Grant Nos. 26400413 and 18K03475), RIKEN iTHES Project, and the Simons Collaboration on the Many Electron Problem. The numerical simulations have been performed on the HOKUSAI supercomputer at RIKEN (Projects No. G17030, No. G17032, No. G18007, and No. G18025) and on the K computer at RIKEN Center for Computational Science (R-CCS) through the HPCI System Research Project (Projects No. hp160159, No. hp170079, No. hp170162, No. hp170308, No. hp170328, and No. hp180098). K. S. acknowledges support from the JSPS Overseas Research Fellowships.

Appendix A

$D_{AB,\sigma}(\boldsymbol{r})$ in the semimetallic phase

In this Appendix, we show that $D_{AB}(\boldsymbol{r})$ decays algebraically in $|\boldsymbol{r}|$ for $|\boldsymbol{r}|/a\gg 1$ in the semimetallic phase. First, we consider the noninteracting limit. To examine the asymptotic form of $D^{(0)}_{AB,\sigma}(\boldsymbol{r})$, we replace the sum over discrete $\boldsymbol{k}$ in Eq. (11) by the integral over continuous $\boldsymbol{k}$ in the whole first Brillouin zone, i.e.,

[TABLE]

where $S_{\rm cell}=3\sqrt{3}a^{2}/2$ is the area of the unit cell. This is justified in the thermodynamic limit and useful for analyzing the low-energy and long-distance behavior. In the thermodynamic limit, Eq. (11) now reduces to

[TABLE]

Since the long-distance behavior of $D_{AB,\sigma}^{(0)}(\boldsymbol{r})$ is dominated by the low-energy spectrum around the Dirac ($K$ and $K^{\prime}$) points, we measure momentum $\boldsymbol{k}$ from the Dirac points ($v=K,K^{\prime}$) as

[TABLE]

Expanding $h_{\boldsymbol{k}}=h_{\boldsymbol{k}_{K}+\boldsymbol{q}}$ around the $K$ point with respect to $\boldsymbol{q}=(q_{x},q_{y})$ and taking up to the linear term in $\boldsymbol{q}$ yield

[TABLE]

The contribution to $D_{AB,\sigma}^{(0)}(\boldsymbol{r})$ from the momentum around the $K$ point is thus evaluated as

[TABLE]

where $q=|\boldsymbol{q}|$, $\boldsymbol{r}=(r_{x},r_{y})$, and $r=|\boldsymbol{r}|$. Here, the integral in the second line is treated as

[TABLE]

where $s=qr$, $J_{0}(s)$ is the zeroth-order Bessel function of the first kind, and $\Lambda$ is a cutoff momentum of order $\Lambda\sim 1/a$. The upper bound of the integral satisfies $r\Lambda\gg 1$ because our interest is in the long-distance ($r/a\gg 1$) behavior. Since the long-distance behavior of the hole propagation should not be affected by the cutoff momentum $\Lambda$, it is possible to set $r\Lambda\to\infty$. Then, the integral of the Bessel function can be performed as $\int_{0}^{\infty}{d}sJ_{0}(s)=1$ and Eq. (29) results in $1/r$, as in the Fourier transform (or the Hankel transform) of the Coulomb potential in two dimensions

[TABLE]

Therefore, the propagation of a hole is long ranged.

Similarly, around the $K^{\prime}$ point, $h_{\boldsymbol{k}_{K}+\boldsymbol{q}}$ can be expanded as

[TABLE]

The contribution to $D_{AB,\sigma}^{(0)}(\boldsymbol{r})$ from the momentum around the $K^{\prime}$ point is thus evaluated as

[TABLE]

The asymptotic form of $D^{(0)}_{AB,\sigma}(\boldsymbol{r})$ for $r/a\gg 1$ is given by the sum of (28) and (32), i.e.,

[TABLE]

Since the contributions from the $K$ and $K^{\prime}$ points interfere with each other, $\boldsymbol{r}$ dependence of $D_{AB,\sigma}^{(0)}(\boldsymbol{r})$ is in general complicated. Nevertheless, among several directions of $\boldsymbol{r}$, one can find that $\boldsymbol{r}$ in the $x$ direction, i.e., $\boldsymbol{r}=n(\boldsymbol{\tau}_{1}+\boldsymbol{\tau}_{2})=(3na,0)$ with $n$ integer, gives a simple asymptotic form

[TABLE]

for $r/a\gg 1$. Figure 6 shows $D^{(0)}_{AB,\sigma}\left(n(\boldsymbol{\tau}_{1}+\boldsymbol{\tau}_{2})\right)$ calculated directly on an $L=1000$ cluster using Eq. (11), which is compared with its asymptotic form in Eq. (34). The agreement of the two results for $r/a\gg 1$ verifies the algebraic decay of $D^{(0)}_{AB,\sigma}(\boldsymbol{r})$, including the coefficient $S_{\rm cell}/2\pi$.

In the case of an interacting system, it is apparent from Eq. (20) that the asymptotic form of $D_{AB,\sigma}(\boldsymbol{r})$ for $r/a\gg 1$ under the assumption of Eq. (19) is given as

[TABLE]

Therefore, in principle, the quasiparticle weight $Z$ can be estimated from the asymptotic behavior of the equal-time single-particle Green’s function itself, without referring to the noninteracting Green’s function.

As shown in Fig. 6, $D_{AB,\sigma}(\boldsymbol{r})$ of the noninteracting system approaches its asymptotic form only at a very long distance in a large cluster. This might also be the case for the interacting systems. Therefore, the direct observation of the asymptotic behavior of $D_{AB,\sigma}(\boldsymbol{r})$ is difficult within the cluster sizes affordable at present within the AFQMC method. Nevertheless, with an appropriate finite-size-scaling analysis, we can obtain useful and reliable predictions on the asymptotic behavior, within the available cluster studied by AFQMC. Indeed, we have found that the quasiparticle weight can be estimated more accurately from the finite-size scaling of the ratio of $D_{AB,\sigma}(\boldsymbol{r})$ between the interacting and noninteracting systems as in Eq. (23), instead of directly fitting the asymptotic behavior of $D_{AB,\sigma}(\boldsymbol{r})$. On the other hand, the exponent characterizing the asymptotic behavior of $D_{AB,\sigma}(\boldsymbol{r})$ in the semimetallic phase can be estimated with reasonable accuracy, also for the noninteracting system, in the way shown in Figs. 3 and  4.

Appendix B $D_{AB,\sigma}(\boldsymbol{r})$ in the insulating phase

In this Appendix, we show that $D_{AB,\sigma}(\boldsymbol{r})$ decays exponentially in $r$ for $r/a\gg 1$ in the insulating phase. The derivation is essentially the same as that in Appendix A. The main difference due to the finite single-particle excitation gap is that the integral over $\boldsymbol{q}$ (the momentum measured from the Dirac point), which yields a massless (Coulomb-potential-like) form for the semimetallic phase as in Eq. (30), now yields a massive (Yukawa-potential-like) form for the insulating phase as in Eq. (41)

To examine the asymptotic form of $D_{AB,\sigma}(\boldsymbol{r})$ in the insulating phase, we model the single-particle Green’s function $\boldsymbol{G}_{\sigma}(\boldsymbol{k},z)$ with the same analytical form of an antiferromagnetically ordered state, i.e.,

[TABLE]

where $(-1)^{\sigma}\equiv+1\,(-1)$ for $\sigma=\uparrow(\downarrow)$ and $\Delta$ is the gap function corresponding to the staggered magnetization that breaks the chiral symmetry Semenoff (2012); Hatsugai et al. (2013). Here, we assume that the magnetization is along the $z$ spin-quantization axis with real $\Delta\,(>0)$, for simplicity. The energy dispersion is obtained by solving $\det\boldsymbol{G}_{\sigma}^{-1}(\boldsymbol{k},z)=0$ with respect to the frequency $z$, i.e., $\pm\sqrt{|\tilde{h}_{\boldsymbol{k}}|^{2}+\Delta^{2}}$, and thus it is massive. In particular, the single-particle excitation gap at the $K$ and $K^{\prime}$ points is $2\Delta$.

Inserting the model single-particle Green’s function into Eq. (13) and taking the zero temperature limit, we can obtain the equal-time single-particle Green’s function for the insulating phase, i.e.,

[TABLE]

By expanding $\tilde{h}_{\boldsymbol{k}}$ around the $K$ point as in Eq. (27), we find that the contribution to $D_{AB,\sigma}(\boldsymbol{r})$ from the momenta around the $K$ point is given as

[TABLE]

where, with the same argument for Eq. (29), the integral over $\boldsymbol{q}$ is performed, as in the Fourier transform (or the Hankel transform) of the Yukawa potential in two dimensions, i.e.,

[TABLE]

with

[TABLE]

The propagation of a hole is thus short ranged in the insulating phase due to the finite single-particle excitation gap $\Delta$.

With the propagation range $\xi$ of a hole in the insulating phase, Eq. (40) can be written as

[TABLE]

Similarly, the contribution to $D_{AB,\sigma}(\boldsymbol{r})$ from the momentum around the $K^{\prime}$ point is evaluated as

[TABLE]

Adding (43) and (44) yields the asymptotic form

[TABLE]

In the limit of $\xi\to\infty$, i.e., $\Delta\to 0$, Eq. (45) reduces to the noninteracting limit in Eq. (33). In conclusion, the equal-time single-particle Green’s function $D_{AB,\sigma}(\boldsymbol{r})$ decays exponentially in $r$ in the single-particle-gapful system with a characteristic length scale $\xi$ given in Eq. (42).