Solving fermion sign problem in quantum Monte Carlo by Majorana representation
In this paper, we discover a new quantum Monte Carlo (QMC) method to solve the fermion sign problem in interacting fermion models by employing Majorana representation of complex fermions. We call it “Majorana QMC” (MQMC). MQMC simulations can be performed efficiently both at finite and zero temperatures. Especially, MQMC is fermion sign free in simulating a class of spinless fermion models on bipartite lattices at half filling and with arbitrary range of (unfrustrated) interactions. Moreover, we find a class of fermionic models with odd , which are sign-free in MQMC but whose sign problem cannot be in solved in other QMC methods such as continuous-time QMC. To the best of our knowledge, MQMC is the first auxiliary field QMC method to solve fermion sign problem in spinless (more generally, odd number of species) fermion models. We conjecture that MQMC could be applied to solve fermion sign problem in more generic fermionic models.
Introduction: Interacting fermionic quantum systems with strong correlations and/or topological properties have attracted increasing attentionsXGWenbook (); Fradkinbook (). Nonetheless, in two and higher spatial dimensions, strongly interacting quantum systems are generically beyond the reach of analytical methods in the sense of solving those quantum models in an unbiased way. As an intrinsically-unbiased numerical method, quantum Monte Carlo simulation plays a key role in understanding physics of strongly correlated many-body systemsScalapino-81 (); Blankenbecler-81 (); Rebbi-81 (); Hirsch-81 (); Hirsch-85 (). Unfortunately, in simulating fermionic many-body systems, QMC often encounters the notorious fermion minus-sign problemLoh-90 (); Troyer-05 (), which arises as a consequence of Fermi statisticsfootnoteFermi (). Undoubtedly, generic solutions of fermion sign problems would lead to a great leap forward in understanding correlated electronic systemsTroyer-05 ().
Many QMC algorithms are based on converting an interacting fermion model into a problem of free fermions interacting with background auxiliary classical fields; the Boltzmann weight is the determinant of free fermion matrix which is a function of auxiliary fields and which can be positive, negative, or even complex. In such determinant QMC (DQMC), when the determinants are rendered to be positive definite, we say a solution to the fermion sign problem is found. For spinful electrons, conventional strategy of solving fermion sign problem is to find a symmetric treatment of both spin components of electrons such that the Boltzmann weight can be written as the product of two real determinants with the same sign and is then positive definiteHirsch-86 (); Hands-00 (); Wu-05 (); Berg-12 (); Assaad-12 (); Wu-12 (). For spinless or spin-polarized fermion models, it is usually much more difficult to solve fermion sign problem because the Boltzmann weight contains only a single determinant and the usual strategy used for even species of fermions cannot be directly applied here.
In this paper, based on Majorana representation of fermions, we propose a genuinely new auxiliary field QMC approach to solve fermion sign problem in spinless fermion models. We observe that each complex fermion can be represented as two Majorana fermions. Consequently, we can express spinless fermion Hamiltonians in Majorana representation and then perform Hubbard-Stratonovich (HS) transformations to decouple interactions by introducing background auxiliary fields. Under certain conditions such as particle-hole symmetry, we can find a symmetric treatment of two species of Majorana fermions, namely the free Majorana fermion Hamiltonian obtained after HS transformations is a sum of two symmetric parts each involving only one species of Majorana fermions, such that the Boltzmann weight is a product of two identical real quantities and is then positive definite. This is the basic idea of the Majorana approach to solve fermion sign problem in spinless or spin-polarized fermion models which we call “Majorana QMC” (MQMC). Note that the MQMC approach proposed here is qualitatively different from the meron-cluster methodWiese-95 (); Wiese-99 () and fermion bags methodShailesh-13 (); Shailesh-14 () developed previously, all of which are based on continuous-time QMC (CTQMC)RMP-11 (); Shailesh-14 (); Lei-14 (); Leiwang-14 (). As far as we know, MQMC is the first QMC approach based on auxiliary fields to solve fermion sign problem in a class of spinless (more generally, odd number of species) fermion models. Moreover, MQMC has an important advantage: it is much more efficient than continuous-time QMC in simulating models at low and zero temperatures; the computation-time cost in MQMC scales as while it scales as in continuous-time QMCShailesh-14 () (also see more recent developmenttroyer-14 ()).
As an application of the sign-free MQMC algorithm, we have used it to study the charge density wave (CDW) quantum phase transition of the spinless fermion model with repulsive density interactions on the honeycomb lattice with much larger system size ( sites with up to 24) than previous studies and obtained quantum critical exponents which are in reasonable agreement with renormalization group calculationsunpublished (). We also show that MQMC can solve the fermion-sign problem in a class of models which are beyond the capability of other QMC methods such as the continuous-time QMC.
Majorana quantum Monte Carlo: To explicitly illustrate how MQMC could solve the fermion sign problem in a class of spinless fermion models, we consider the following general Hamiltonian of spinless fermions:
where creates a fermion on site , represents hopping integral and labels density interaction. As we shall show below, the MQMC is fermion-sign-free when the Hamiltonian in Eq. (1) satisfies the following two conditions: (1) only when belong to different sublattices; (2) when belong to different sublattices and when belong to same sublattices. With the first condition, it is clear that the model is invariant under particle-hole transformations: , where has opposite signs for different sublattices and then describes fermions at half-filling. The lattice in question can be any bipartite lattice such as honeycomb and square lattices in 2D as well as cubic and diamond lattices in 3D. For simplify, we hereafter consider the model with only nearest-neighbor (NN) hopping , NN repulsive interaction , and next-nearest-neighbor (NNN) attractive interactions which we call the -- model on the honeycomb lattice (generalizing the MQMC method to models with longer-range hopping/interactions will be straightforward). As shown in Fig. 1, MQMC is fermion sign free in the region where the quantum phase transition between Dirac semimetal and charge density wave (CDW) phases occursLei-14 (). (It is interesting to note that the -- spinless fermion model on the honeycomb lattice feature very interesting phases including quantum anomalous Hall (QAH) phasesRaghu-08 () and pair density wave (PDW) phasesJian-14 ().)
In statistical physics, a key quantity is the partition function. QMC methods are designed to simulate partition functions in a statistical fashion. For the -- model, the partition function after Trotter decomposition is given by
where labels the discrete imaginary time, , and the approximation is good for small or large . HS transformations can be applied to decouple fermion interactions into non-interacting terms interacting with background auxiliary fields. Usual HS decoupling in density channels normally result in minus sign problem in QMC because the Boltzmann weight is a single determinant. However, we observe that the Hamiltonian can be rewritten in terms of Majorana fermions and there are two species of Majorana fermions. In Majorana representation, complex fermions operators are given by:
which enable us to rewrite the Hamiltonian as follows:
where gauge transformations for in only one sublattice were implicitly made so that can be written symmetrically in the two components of Majorana fermions. Now, it is clear that we should perform HS transformations in Majorana hopping channels instead of density channels as done in usual QMC methods. Explicitly, HS transformations for interactions in in MQMC are given by
where and are constants determined through and , respectively. Note that in Eq. (7) the signs of hopping terms are opposite to hopping terms in the HS decompositions of NNN interaction because . The same signs are obtained for decoupling of NN interactions in Eq. (6) because . It is now clear that the free fermion Hamiltonian after the HS transformations is a sum of two parts each of which involves only one component of Majorana fermions. This makes MQMC simulations sign-problem free because the Boltzmann weight can be positive definite, which we shall show below.
Note that auxiliary fields should be introduced independently for each discrete imaginary time . As a result, the partition function is a sum over Boltzmann weight which is a function of auxiliary field configurations in space-time, as given by
Up to an unimportant constant the Boltzmann weight is given by
where represents the transpose of and is a matrix (=the number of lattice sites) is given by
where if are NN sites and 0 otherwise; similarly only if are NNN sites. Now, we can trace out the Majorana fermions since they are free, as shown in Supplemental Material. Because the two components of Majorana fermions are decoupled, tracing out Majorana fermions can be done independently and the Boltzmann weight is a product of two factors:
Note that there is sign ambiguity when taking a square root above, similar to the case of Pfaffian as a square root of determinants.
Fermion sign free: Now we prove that the Boltzmann weight is positive definite by showing that . A key observation is that the Hamiltonian of Majorana fermions can be mapped to a Hamiltonian identical to by the following time-reversal transformation , where is the complex conjugation and is given as below:
Namely, under the time reversal transformation . Because the time-reversal transformation complex conjugates the results of tracing out Majorana fermions, we obtain
which renders the Boltzmann weight for any auxiliary field configuration . Explicitly, it is
where or , which gives rise to the same result. This proves that the MQMC algorithm can solve fermion sign problem in such class of models consisting of spinless fermions. It is the central result in this paper.
Projector MQMC: The MQMC algorithm above simulates finite-temperature partition function in the grand canonical ensemble by computing the trace in Eq. (9). If one is interested in ground state properties, it is of advantage to use the projector algorithm to carry out QMC Sugiyama-86 (); Sorella-89 (); White-89 () since projector QMC is often more efficient than finite-temperature QMC. The expectation value of an operator in the ground state is given by
where is the ground state and is a trial wave function which we assume has a finite overlap with the true ground state. Here, plays the role of usual partition functions and need to be expressed as a sum of Boltzmann weights. In practice, a Slater-determinant wave function describing non-interacting fermions is often chosen as the trial wave function in projector QMC:
where is a matrix ( labels the number of fermions in question). Usually, is an eigenvector of the non-interacting part of the Hamiltonian in question, namely in Eq. (1). In Majarana representation of fermions, and Majorana fermions are decoupled in ; consequently . By introducing similar HS transformations and auxiliary fields as above, the “partition function” is obtained a sum of Boltzmann weight over auxiliary field configurations: . Since and Majorana fermions are decoupled after the HS transformation, we again obtain , where
Similarly, because of the time reversal symmetry . As shown in the Supplemental Material, the Boltzmann weight is given by
where or 2 and is the projection matrix constructed from . Consequently, the projector MQMC is also free from fermion sign problem for a class of spinless fermion models.
Physical observables in MQMC: One important advantage of auxiliary-field QMC algorithms is that physical observables can be obtained conveniently. For instance, time and space dependent Green’s function can be computed directly in DQMC algorithm. We show below that both at finite and zero temperature the computation of physical observables in MQMC is similarly convenient as that in DQMC algorithm.
In QMC, physical observables can be related to single-particle Green’s function: , where the average is done stochastically over auxiliary field configurations. In Majorana representation, it is given by
where we used the results of which is a consequence of the decoupling of the two species of Majorana fermions after the HS transformation. To obtain the Green’s functions, we only need to compute and . Because the two species of Majorana fermions are related by the time reversal symmetry , we obtain . It is straightforward to evaluate the equal-time Majorana Green’s function in finite temperature MQMC:
where the factor 1/2 above comes from the nature of Majorana fermions. Employing Wick’s theorem for each configuration , higher order correlation functions, including density-density and pair-pair correlations, can be obtained from single-particle Green’s functions. For instance, the equal-time density-density correlations are given by .
It is increasingly realized that quantum entanglement could play a key role in understanding quantum many-body systemsCardy-04 (); Kiteav-06 (); Levin-06 (); Haldane-08 (); Eisert-10 (). Quantum entanglement is partially characterized by entanglement entropy, including the von Neumann entropy and Renyi entropy where is the reduced density matrix of subregion . Even though it is still challenging for auxiliary-field QMC algorithms to evaluate von Neumann entropy, it was shown recently that DQMC can provide an efficient way to evaluate Renyi entropy by simulating the reduced density matrix expressed in terms of Green’s functionGrover-13 (); Assaad-14 (). Because MQMC is able to compute Green’s functions efficiently, Renyi entropy can be calculated accurately in MQMC algorithm as long as it is fermion sign free.
Numerical Results: We performed highly-accurate projector MQMC simulations to study the -- model on the honeycomb lattice at zero temperature. For simplicity, we set , , and then vary to find the critical value of , above which the system develops a finite CDW ordering at zero temperature. To measure the CDW order parameter , we calculate CDW structure factor at finite lattice size:
where on sublattice and is the total number of sites. It is obvious that . The simulations are done for lattices up to which is substantially larger than the one in Ref. Lei-14 () indicating that our MQMC algorithm is quite efficient. As shown in Fig. 2(a), we obtain through finite-size scaling of the measured on lattices of . For instance, at . It is clear that the critical value of separating the semimetal and CDW phases is between 1.34 and 1.38. To obtain the critical value of more accurately, we calculate the Binder ratio defined as for various and , where . At the putative critical point, the Binder ratios for different should cross. As shown in Fig. 2(b), the Binder ratios for indeed cross nearly the same point when . Consequently, we conclude that the critical value .
The critical exponents and universality class at the phase transitionsachdevbook (); Herbut-14 () have been analyzed through even larger-scale MQMC simulations by usunpublished () . Because the CPU time-cost scales linearly with , we were able to perform the MQMC simulations on much larger system size ()unpublished () than the one studied by CTQMC ( there)Lei-14 (); consequently the critical exponents obtained by MQMC are reasonably consistent with RG calculations.
Other models sign-free in MQMC: We have shown that MQMC, as a new auxiliary field QMC approach, can solve fermion sign problem in a class of spinless fermion models by utilizing Majorana representation of complex fermions. It will be straightforward to generalize the current MQMC algorithm to solve the fermion sign problem in interacting fermion models with more than one fermion species. Such MQMC fermion-sign free models include the negative- Hubbard model on bipartite lattices whose Hamiltonian is
where and . This model on the honeycomb lattice has a similar semimetal to CDW transition even though the quantum critical exponents can depend on .
More importantly, we can show that the following fermionic model
is sign-free in MQMC when the lattice is bipartite and . It is worth to stress that this class of models are sign-free only in the MQMC method but encounter sign-problem in other QMC methods such as CTQMCShailesh-14 (); Lei-14 (). This shows that the MQMC algorithm discovered by us can solve the fermion-sign of models which go beyond those solvable by CTQMC and other conventional QMC methods.
Acknowledgement: We thank Fakher Assaad, Alexei Kitaev, Ziyang Meng, and Lei Wang for helpful discussions. This work is supported in part by the National Thousand-Young-Talents Program (HY) and the NSFC under grant No. 11474175 (ZXL, YFJ, and HY).
- (1) X.-G. Wen, Quantum Field Theory of Many-body Systems, (Oxford University Press, New York, 2004).
- (2) E. Fradkin, Field Theories of Condensed Matter Physics, Second Edition, (Cambridge University Press, Cambridge, 2013).
- (3) D. J. Scalapino and R. L. Sugar, Phys. Rev. Lett 46, 59 (1981).
- (4) R. Blankenbecler, D. J. scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
- (5) F. Fucito, E. Marinari, G. Parisi, and C. Rebbi, Nucl. Phys. B 180, 369 (1981).
- (6) J. E. Hirsch, D. J. Scalapino, R. L. Sugar, and R. Blankenbecler, Phys. Rev. Lett. 47, 1628 (1981).
- (7) J. E. Hirsch, Phys. Rev. B 31, 4403 (1985).
- (8) E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar, Phys. Rev. B 41, 9301 (1990).
- (9) M. Troyer and U.-J. Wiese, Phys. Rev. Lett 94, 170201 (2005).
- (10) A negative sign in quantum magnets with frustrated superexchange interactions is also a consequence of the fermionic minus sign of the superexchange process.
- (11) J. Hirsch and R. Fye, Phys. Rev. Lett. 56, 2521 (1986).
- (12) S. Hands, I. Montvay, S. Morrison, M. Oevers, L. Scorzato, and J. Skullerud, Eur. Phys. J. C 17,285(2000).
- (13) C. Wu and S.-C. Zhang, Phys. Rev. B, 71, 155115 (2005).
- (14) E. Berg, M. A. Metlitski, and S. Sachdev, Science 338, 1606 (2012).
- (15) M. Hohenadler, Z. Y. Meng, T. C. Lang, S. Wessel, A. Muramatsu, and F. F. Assaad, Phys. Rev. B 85, 115132 (2012).
- (16) D. Zheng, G.-M. Zhang, and C. Wu, Phys. Rev. B 84, 205121 (2011).
- (17) W. Bietenholz, A. Pochinsky, and U. J. Wiese, Phys. Rev. Lett. 75, 24 (1995) .
- (18) S. Chandrasekharan and U. J. Wiese, Phys. Rev. Lett. 83, 16 (1999) .
- (19) S. Chandrasekharan, The European Physical Journal A 49, 90 (2013) .
- (20) E. F. Huffman and S. Chandrasekharan, Phys. Rev. B 89, 111101 (2014) .
- (21) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer and P. Werner, Rev. Mov. Phys, 83, 2011
- (22) L. Wang, P. Corboz and M. Troyer, New J. Phys. 16, 103008 (2014).
- (23) L. Wang and M. Troyer, Phys. Rev. Lett. 113, 110401 (2014).
- (24) A more efficient continuous-time QMC algorithm was intro- duced recently. See M. Iazzi and M. Troyer, Phys. Rev. B 91, 241118 (2015); L. Wang, M. Iazzi, P. Corboz, and M. Troyer, ibid. 91, 235151 (2015).
- (25) Z.-X. Li, Y.-F. Jiang, and H. Yao, arXiv:1411.7383 (2014), New J. Phys. (to be published).
- (26) S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang, Phys. Rev. Lett. 100, 156401 (2008).
- (27) S.-K. Jian, Y.-F. Jiang, and H. Yao, Phys. Rev. Lett. 114, 237001 (2015).
- (28) G. Sugiyama and S. Koogin, Anals of Phys. 168, 1 (1986).
- (29) S. Sorella, S. Baroni, R. Car and M. Parrinello, Europhys. Lett. 8, 663 (1989).
- (30) S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis and R. T. Scalettar, Phys. Rev. B, 40, 506 (1989).
- (31) P. Calabrese and J. Cardy, J. Stat. Mech.: Theory Exp. P06002 (2004).
- (32) A. Kiteav and J. Preskill, Phys. Rev. Lett 96, 110404(2006).
- (33) M. Levin and X.-G. Wen, Phys. Rev. Lett. 96, 110405 (2006).
- (34) H. Li and F. D. M. Haldane, Phys. Rev. Lett 101, 010504 (2008).
- (35) J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mov. Phys, 82, 277 (2010)
- (36) T. Grover, Phys. Rev. Lett 111, 130402 (2013).
- (37) F. F. Assaad, T. C. Lang, and F. P. Toldin, Phys. Rev. B 89, 125121 (2014)
- (38) S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, ed. 2, 2011).
- (39) L. Janssen and I. F. Herbut, Phys. Rev. B 89, 205403 (2014).
Appendix A Supplemental Materials
a.1 Appendix A: Trace involving Majorana fermions
Now we show that the trace of exponentials of bilinear Majorana fermion operators can be expressed as the square root of a determinant (nonetheless, it is formally not a Pfaffian as shown below).
First, we evaluate the trace of a single exponential of bilinear Majorana fermion operators. Suppose where and we need to compute . After diagonalizing , where are eigenvalues of the matrix , it is clear that the trace is given by
which can be reexpressed as the square root of a determinant:
Intuitively, the square root originates from the fact that Majorana fermions carry only half of degrees of freedom of corresponding Hamiltonian in terms of complex fermions.
Then we show that the product of exponentials of bilinear Majorana fermion operators can be grouped into a single exponential of a bilinear Majorana fermion form. Suppose , where , and we would like to evaluate . By observing that , we obtain
where a matrix is defined. Accordingly, we introduce bilinear Majorana fermion operators: . Now, we can show that the trace of the product of exponentials of the Majorana fermions bilinear operator is given by the square root of a determinant:
which proves the result in Eq. (11) of the main text.
a.2 Appendix B: Projector QMC in Majorana representation
We now prove the result in Eq. (19) of the main text. To compute , we compute its square first as follows:
where are the “ghost Majorana fermions” which are independent from but have the same Hamiltonian and ground state wave function as . When combining and into complex fermions , we obtain
where is a projector matrix defined through . Because , we prove Eq. (19) as follows:
a.3 Appendix C: Proof of fermion-sign free in a class of fermionic model with bond interactions
We now prove in details that the fermionic model with fermion species described by Eq. (24) does not encounter fermion-sign problem in our MQMC algorithm. Both the hopping term and the interaction term in Eq. (24) can be rewritten in terms of Majorana fermions:
where gauge transformations for in one sublattice are implicitly made. Then, we can perform a similar Hubbard-Stratonovich (HS) transformation on the bond interactions:
where is a constant satisfying . For , is a real number. After the HS transformation, it is clear that Hamiltonian is a sum of two parts each of which involves only one component of Majorana fermions ( and ):
So Boltzmann weight can be decoupled as the product of identical parts:
where is the Boltzmann weight obtained through tracing out one component of the Majorana fermions, say . Moreover, each part of Hamiltonian is invariant under this anti-unitary time-reversal transformation where is given as below:
As a result, the Boltzmann weight is real and consequently . This proves that the fermionic model with fermion species described by Eq. (24) can be fermion-sign free in MQMC simulations.