Suppression of topological Mott-Hubbard phases by multiple charge orders in the honeycomb extended Hubbard model
We investigate the competition between charge-density-wave (CDW) states and a Coulomb interaction-driven topological Mott insulator (TMI) in the honeycomb extended Hubbard model. For the spinful model with on-site () and next-nearest-neighbor () Coulomb interactions at half filling, we find two peculiar six-sublattice charge-density-wave insulating states by using variational Monte Carlo simulations as well as the Hartree-Fock approximation. We observe that conventional ordered states always win with respect to the TMI. The ground state is given in the large- region by a CDW characterized by a 220200 (001122) charge configuration for smaller (larger) , where 0, 1, and 2 denote essentially empty, singly occupied, and doubly occupied sites. Within the 001122-type CDW phase, we find a magnetic transition driven by an emergent coupled-dimer antiferromagnet on an effective square lattice of singly occupied sites. Possible realizations of the found states are discussed.
Correlation effects in electron systems on a two-dimensional honeycomb lattice have been the subject of intensive scrutiny , both in metallic and insulating materials [2; 3]. A special focus has been the possible realization of the quantum Hall effect (QHE), i.e., the quantization of the Hall conductivity in two dimensions, not as the result of an external magnetic field , but due to the complex hopping that breaks time-reversal symmetry [5; 6]. A related state, the quantum spin Hall (QSH), for which the spin Hall conductance is quantized, emerges in analogy when spin-orbit interactions lead to robust spin-dependent transport [7; 8; 9; 10]. The spin quantum Hall state has also been studied in mean-field theory for charged-ordered triangular lattices with spin-orbit interactions .
Raghu et al. proposed that a topological Mott insulator (TMI) could be stabilized on the extended honeycomb lattice . The QHE would emerge in this scenario from pure Coulomb interactions, with an effective spin-orbit interaction being dynamically generated via spontaneous symmetry breaking when the next-nearest-neighbor Coulomb interaction is taken into account. Here, in a spinless model, time-reversal symmetry is broken spontaneously, while the lattice translational symmetry is preserved; the topological phases are characterized by the Chern number . On the other hand, for QSH in a spinful model, symmetry is broken spontaneously, while time-reversal symmetry is preserved; they are characterized by the invariant . Nonlocal Coulomb interactions, which are non-negligible in graphene , may also generate, however, conventional spontaneously symmetry-broken states, such as bond order, charge order, and magnetic order, which then compete with the TMI.
In the spinless Hubbard model with nearest-neighbor () and next-nearest-neighbor () Coulomb interactions, a Kekulé bond order phase, characterized by a order parameter, has been proposed by mean-field calculations . Exact-diagonalization (ED) [17; 18; 19] and infinite density matrix renormalization-group (iDMRG)  studies also support the presence of this phase and, further, propose much richer charge-density-wave (CDW) phases. Possible charge instabilities have also been investigated away from half filling [21; 22].
Rich charge-density-wave phases also emerge for the spinful Hubbard model. In addition, when the on-site Coulomb interaction is dominant, collinear antiferromagnetic order appears . Inclusion of and/or gives rise to the competition of magnetic and several CDW states [24; 25; 26]. In this case, charge and spin degrees of freedom are intertwined, especially, off half filling, and the formation of superstructures often enhances the geometrical frustration of spins that allows noncollinear magnetic order to coexist with the CDW .
Recent extensive research on the honeycomb extended Hubbard model by ED [17; 18; 19], iDMRG , auxiliary-field quantum Monte Carlo (AFQMC) , variational Monte Carlo (VMC) , and functional renormalization group (fRG) [29; 25; 26] suggests that TMI is less likely than conventional ordered states, often CDW, at half filling. This fact leads us to ask the following questions:
Which are the plausible CDW patterns for the honeycomb extended Hubbard model at half filling?
Can CDW coexist with magnetic order when spin degrees of freedom are present?
Is TMI always suppressed by CDW states? Could they coexist via the spontaneous formation of both effective spin-orbit interaction and charge ordering?
To answer these questions, we revisit the ground states of the extended Hubbard model on the honeycomb lattice. We focus on the spinful model at half filling, and consider the on-site and next-nearest-neighbor Coulomb interactions, which are relevant for the TMI and CDW states. For simplicity, we drop the nearest-neighbor Coulomb interaction . We first use the restricted Hartree-Fock approximation to get an insight into the plausible CDW states, and then apply the VMC method with a Jastrow-Slater-type wave function to improve the mean-field state. We find, by the Hartree-Fock approximation, two types of six-sublattice CDW insulating phases: one is a 220200-type CDW phase realized for small and large , and the other is a 001122-type CDW phase which appears for larger and . The transition between these two is found to be of first order within the mean-field treatment. When the VMC method is employed, for a fixed large , we find a continuous change from the 220200-type to the 001122-type CDW phase as is increased. The first-order transition found in the Hartree-Fock approximation seems to melt when quantum fluctuations are carefully taken into account. On the other hand, when is sufficiently large, we again find the stable 001122-type CDW. Within the mean-field study, the 001122-type CDW phase shows a magnetic transition at sufficiently large . The magnetic moment appears for two singly occupied sites out of six sites in a unit cell, and these spins align antiferromagnetically via superexchange interaction (see Fig. 4 for the representation of the various phases). Furthermore, for the parameter region that we have studied, the energy of the TMI state is found to always be higher than these CDW states. The TMI does not coexist with CDW, and thus CDW is harmful for stabilizing the TMI phase at half filling, consistent with previous studies.
This paper is organized as follows: In Sec. II, we present the honeycomb extended Hubbard model and introduce the Hartree-Fock approximation and the VMC method. In Sec. III, we present the - phase diagram and the properties of the three CDW phases obtained by the Hartree-Fock approximation. We then show how quantum fluctuations modify the phases by means of the VMC method. In Sec. IV, we comment on the origin of the charge and magnetic orders in CDW phases. Finally, in Sec. V, we draw our conclusions.
Ii Model and methods
ii.1 Extended Hubbard model
We consider the extended Hubbard model on the honeycomb lattice at half filling,
Here, denotes the hopping parameter and () denotes the strength of on-site (next-nearest-neighbor) Coulomb interaction, as shown in Fig. 1. For simplicity, we only deal with the next-nearest-neighbor and neglect the nearest-neighbor . Hereafter, we consider repulsive Coulomb interactions () at half filling ().
ii.2 Restricted Hartree-Fock method
To clarify the plausible ordered phases, we first apply the restricted Hartree-Fock method. To obtain the mean-field Hamiltonian, we use the Hartree-Fock decoupling,
where are site indices while are spin indices, and determine the order parameters self-consistently. In this paper, we consider six independent sites, called A, B, C, D, E, and F, in a unit cell, which give a system size of , with being a linear size (see Fig. 2).
To characterize each phase, we calculate the number of electrons,
and magnetic order parameter,
per each sublattice . In a charge-uniform nonmagnetic phase, and for all . When for a different sublattice pair of , the phase shows charge disproportionation. Similarly, when , the phase is magnetic.
In the TMI phase, the expectation value gives an imaginary number for the next-nearest-neighbor-site pair . Using this value, the order parameter is defined as
where for clockwise (anticlockwise) orientation and is a Pauli spin. This will be a coefficient of dynamically generated spin-orbit interaction, and having a nonzero is a necessary condition of the TMI phase.
To determine whether the phase is metallic or insulating, we calculate the density of states and estimate the size of the charge gap.
ii.3 Variational Monte Carlo method
Here, is an eigenstate of an auxiliary Hamiltonian given by
where , , , , , and are variational parameters. The honeycomb lattice contains two sites ( and ) in a unit cell, and we choose for (). On the other hand, is the charge Jastrow factor  given by
This wave function can represent metallic and insulating states with and without charge disproportionation . We optimize the translational-invariant Jastrow factor and the variational parameters in the auxiliary Hamiltonian. Each ordered phase is characterized by the charge order and magnetic order parameters, as in the Hartree-Fock approximation. Hereafter, we focus on the lattice systems , with and .
To estimate the size of the charge gap, we calculate the charge structure factor defined as
When for , a gap for the particle-hole excitation vanishes and the state is metallic. On the other hand, when for , converges to a nonzero value and the state is insulating. In practical calculations, we choose .
Note that the metallic phase in the honeycomb Hubbard model shows singular behavior due to the presence of a Dirac cone ; therefore, the aforementioned criterion for is not optimal to distinguish metallic and insulating states. However, in practical VMC calculations, in most cases, we observe an abrupt decrease of a finite-size charge gap estimated by , which signals the onset of a metallic phase. Hereafter, we first distinguish metallic and insulating phases in a mean-field calculation, and then adopt Eq. (14), as a complementary way, to estimate the size of the charge gap.
iii.1 Mean-field phase diagram
The - phase diagram obtained by the restricted Hartree-Fock approximation is shown in Fig. 3. In the absence of next-nearest-neighbor Coulomb interaction (), a continuous transition occurs at between the nonmagnetic charge-uniform semimetal and antiferromagnetic insulating (AFI) state . These two phases remain for small , consistent with a previous mean-field study .
On the other hand, when increases, we find three insulating phases characterized by different charge- and magnetic-order patterns. As we detail below, two nonmagnetic CDW phases are similar to the ones obtained from mean-field , fRG [25; 26], ED [17; 18; 19], and iDMRG  calculations. In addition to the two phases, we find a CDW phase with antiferromagnetic order in the large- and - region.
iii.2 CDW states
We now discuss the charge and magnetic properties of the three CDW states for sufficiently large obtained by the Hartree-Fock approximation.
In the absence of on-site Coulomb interaction (), we find a six-sublattice CDW phase, characterized by three charge-rich sites and three charge-poor sites in a unit cell [see Fig. 4(a)]. The charge pattern is like the 220200-type, namely, the charge-order parameters show or . However, there is a charge disproportionation within the three charge-rich (charge-poor) sites; one of them is richer (poorer) than the other two. This CDW phase survives for small but nonzero .
At large , we find another six-sublattice CDW phase, characterized by a 001122-type charge pattern, namely, , , or [see Fig. 4(b)].
Although both CDW phases show six-sublattice orderings, we find a first-order transition from one to the other. Along a line, the charge-order parameters for six sites in a unit cell show a clear jump at very small [see Fig. 5(a)]. However, this might be an artifact of the mean-field approximation. As we will see later, quantum fluctuations melt the first-order transition and give a continuous transition or a crossover.
Inside the 001122-type six-sublattice CDW phase, we find a magnetic transition when is further increased [see Fig. 4(c)]. Magnetic order appears for the sites with and the nearest spins align antiferromagnetically. As shown in Fig. 5(b), this magnetic transition is continuous and is distinguished from the first-order transition between the 220200- and 001122-type CDW phases.
We also calculate the density of states to estimate the size of the charge gap. The number of bands are, at most, six since the unit cell contains six sites. A shown in Fig. 6, all three CDW phases have a gap at half filling, suggesting their insulating nature.
iii.3 Absence of a topological Mott insulator
A topological Mott-insulating state was originally proposed for a large- region . To investigate whether the TMI state can be the ground state against the CDW states at the Hartree-Fock mean-field level, we compare the energy of each state, as shown in Fig. 7. The energy of the TMI state is found to always be higher than the CDW state.
In general, the quantum Hall effect is allowed within the CDW states for a suitable lattice structure in the presence of the spin-orbit coupling . If such a spontaneous spin-orbit coupling is generated dynamically via , the coexistence of CDW and TMI is not prohibited. However, within the parameter region we have studied, we do not find a coexistence of CDW and TMI. When the charge-order parameter is not , we numerically find that the expectation value does not give an imaginary value, and the order parameter for the TMI is zero.
These results suggest that the TMI does not occur in the honeycomb extended Hubbard model at half filling.
iii.4 Variational Monte Carlo results
The phase diagram obtained by the VMC method is given in Fig. 8. As in the Hartree-Fock approximation, we find a metallic phase, an antiferromagnetic insulating phase, and both 220200- and 001122-type CDW phases. For simplicity, in the present VMC study, we do not consider possible magnetic order in the 001122-type CDW phase, namely, the state in Fig. 4(c).
The region of validity of the ordered phases, which are overestimated in the Hartree-Fock approximation, get shrunk and the metallic phase region gets enlarged when taking into account the effect of quantum fluctuations through the Jastrow correlation factor. Indeed, we find a single metal-insulator transition point along , estimated to be at , which is in good agreement with the numerically exact value . Furthermore, when the nearest-neighbor Coulomb interaction increases, the on-site is effectively screened. As a result, when is fixed, a metallic phase is expected to be more stable for larger . Such behavior is already observed in recent AFQMC calculations , and our VMC phase diagram agrees with the tendency.
Let us focus on the parameter region along , where the CDW phases are dominant. The charge-order parameters as a function of are given in Fig. 9. When , or , suggesting that the charge pattern is 220200-like. By contrast, when , , , or , suggesting that the charge pattern is 001122-like. The 001122 pattern is more stable for larger . Surprisingly, there is no clear jump in between, and these two CDW phases are found to be continuously connected. The first-order transition found in the Hartree-Fock approximation melts when quantum fluctuations are carefully introduced.
A continuous change of CDW has been reported before in the presence of , which favors two-sublattice staggered charge ordering . Since and favor CDW states with different ordering wave vectors, it is natural to expect incommensurate CDW states in the presence of both and . However, here we find a continuous change of CDW phases without . The size of the unit cell is always kept to six sites in our case.
When , charge distribution becomes uniform and the collinear antiferromagnetic phase appears. The transition between the CDW phase and the charge-uniform antiferromagnetic phase is of first order, as in the Hartree-Fock approximation.
Finally, we discuss the charge gap obtained by the VMC method. As shown in Fig. 10, when and are relatively larger than , for the CDW phases converges to a nonzero value, suggesting the two phases are insulating. However, at smaller- values, approaches a tiny value, and we cannot firmly conclude that the phases are insulating. Larger system sizes are needed to investigate the size of the charge gap in this region.
iv.1 Origin of charge order
When is small while is dominant, a charge-uniform state is favorable since prohibits double occupation of electrons. By contrast, for the large- and - region, basic charge-order patterns of the CDW phases can be understood by considering the atomic limit ().
Let us first consider the case of small but nonzero . When while is turned on, the original honeycomb lattice decouples into doubled triangular lattices. One has only to minimize the energy for each sublattice or , separately, by choosing the plausible charge-ordering patterns. It is enough to find such configurations by minimizing the energy of a local triangle connected by . At half filling, the number of electrons per triangle should be three. Besides, positive favors empty sites to reduce the energy loss in , while on-site disfavors doubly occupied sites. Consequently, the optimal charge configuration is a 012-type charge pattern. Note that the number of degenerate configurations is finite, but still large. The 012 pattern is sixfold degenerate for each sublattice or , and, therefore, the total number of degeneracy on the honeycomb lattice is .
These configurations can be divided into two types that cannot be transformed to each other by mirroring or rotation; namely, 001122-type [see Fig. 4(b)] and 021120-type [see Fig. 11] charge orderings. Both of them are 18-fold degenerate, respectively. In the atomic limit, these two states have exactly the same energy; however, quantum fluctuations lift the degeneracy. The Hartree-Fock approximation predicts that nonzero hopping favors the 001122-type pattern.
In contrast to the case, in the absence of , the ground state in the atomic limit shows macroscopic degeneracy. Since double occupation is no longer prohibited by , the number of singly occupied sites is flexible. Meanwhile, the total number of electrons should be identical to the number of sites at half filling. Therefore, any two different local charge patterns can be turned into without any additional energy cost. As in the case, nonzero hopping lifts the degeneracy and seems to select the 220200-type CDW phase.
When , the transition point between the charge-uniform antiferromagnetic collinear state and 001122-type CDW state can be estimated by comparing their energies. On a single triangle connected by , the charge-uniform state gives per site when the spins are assumed to be fully polarized, while the CDW state gives per site. As a result, the CDW phase is stable when . The phase boundaries obtained by the Hartree-Fock approximation and the VMC method are in good agreement with this prediction.
iv.2 Origin of magnetic order
Let us consider the origin of magnetic order coexisting with CDW for large and . In the 001122-type CDW phase, spin degrees of freedom survive for singly occupied sites [see Fig. 4(b)], and superexchange interaction is generated via virtual hopping processes. For the singly occupied sites next to each other, is generated by the second-order process. Similarly, for the singly occupied sites connected by the Manhattan distance, , , and are generated by the sixth-order process. These processes realize the square Heisenberg model with , as shown in Fig. 12.
When , this is basically the simple antiferromagnetic square Heisenberg model, and conventional antiferromagnetic Néel order is expected. On the other hand, when , the system nearly decouples into isolated dimers, and staggered dimer order is expected. Such a coupled-dimer antiferromagnetic model [38; 39; 40; 41] has been previously numerically studied. For , the Néel-dimer transition occurs at [39; 40; 41].
In the simple spin system, the nature of this Néel-dimer transition is believed to belong to the three-dimensional classical Heisenberg universality [40; 41], although its numerical detection is known to be difficult due to the large finite-size effect [39; 40; 41]. One may expect the same scenario in the present system; however, spin and charge degrees of freedom are coupled in the Hubbard model. Furthermore, the anisotropic spin-exchange interactions are not static, but generated by spontaneous charge order. These factors possibly modify the nature of the transition and make it much harder to identify the magnetic transition numerically.
The Hartree-Fock approximation seems to find these magnetic and nonmagnetic CDW states although the mean-field method, in principle, cannot represent dimer singlet pairs. Besides, the effective ratio could be a nonmonotonic function of and . Therefore, it is not trivial how dimer and magnetic phases appear in the - phase diagram. Quantum fluctuations may modify the mean-field magnetic-nonmagnetic phase boundary significantly inside the CDW phase. Investigating the Néel-dimer transition point and its nature is beyond the scope of this paper, and it is left as a subject for future study.
iv.3 Possible realizations
In general, Coulomb interactions become smaller as the distance is increased, and hence it is not simple to realize . However, these CDW phases may be found in silicon adatoms [42; 43] or optical lattices  when the lattice forms a double-layered triangular structure. This is because the next-nearest-neighbor Coulomb interaction on the original honeycomb lattice corresponds to the nearest-neighbor one in triangular lattices, and could be dominant when the layers are far enough.
We have investigated the phase diagram of the Hubbard model on the honeycomb lattice at half filling in the presence of on-site () and next-nearest-neighbor () Coulomb interactions. By applying the restricted Hartree-Fock approximation, we find three six-sublattice CDW insulating phases, namely, (i) a nonmagnetic 220200-type CDW phase, (ii) a nonmagnetic 001122-type CDW phase, and (iii) a magnetic 001122-type CDW phase, as well as the semimetal and antiferromagnetic insulating phases. To investigate the stability of CDW phases beyond the mean-field study, we further apply the VMC method with Jastrow-Slater-type wave functions. We find that quantum fluctuations destroy the first-order transition found in the mean-field approximation, and give a continuous change from 220200-type CDW to 001122-type CDW. On the other hand, when and are large enough, we find a solid 001122-type CDW phase.
In contrast to previous studies on the spinful model, where nonmagnetic CDW phases were proposed, we find a magnetic transition within the CDW phase. In the magnetic CDW phase, only the singly occupied sites contribute to magnetism, and they show collinear antiferromagnetic order. To understand the origin of magnetic order, we considered the possible superexchange interactions through virtual hopping processes. The spins in the 001122-type CDW phase are found to behave as in a coupled-dimer antiferromagnetic model on the square lattice, where nonmagnetic and magnetic phases are determined by the anisotropy of the spin-exchange interaction.
For the parameter region we have studied, the TMI phase is found to be less favorable than the CDW phases, and there is no coexisting region of TMI and CDW. Therefore, CDW is harmful for the TMI phase in the honeycomb extended Hubbard model at half filling. These CDW phases may be found in silicon adatoms or optical lattices where sizable on the honeycomb lattice could be realized by favoring a double-layered triangular structure.
Acknowledgements.M.B., R.K., and R.V. acknowledge the support of the German Science Foundation (DFG) through Grant No. SFB/TRR49.
- Castro Neto et al.  A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- McChesney et al.  J. L. McChesney, A. Bostwick, T. Ohta, T. Seyller, K. Horn, J. González, and E. Rotenberg, Phys. Rev. Lett. 104, 136803 (2010).
- Tarruell et al.  L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, Nature (London) 483, 302 (2012).
- Klitzing et al.  K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
- Haldane  F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
- Qi et al.  X.-L. Qi, Y.-S. Wu, and S.-C. Zhang, Phys. Rev. B 74, 085308 (2006).
- Kane and Mele [2005a] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005a).
- Bernevig and Zhang  B. A. Bernevig and S.-C. Zhang, Phys. Rev. Lett. 96, 106802 (2006).
- Bernevig et al.  B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006).
- König et al.  M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007).
- Sugita and Motome  Y. Sugita and Y. Motome, J. Phys. Soc. Jpn. 85, 073709 (2016).
- Raghu et al.  S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang, Phys. Rev. Lett. 100, 156401 (2008).
- Thouless et al.  D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
- Kane and Mele [2005b] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005b).
- Wehling et al.  T. O. Wehling, E. Şaşıoğlu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blügel, Phys. Rev. Lett. 106, 236805 (2011).
- Weeks and Franz  C. Weeks and M. Franz, Phys. Rev. B 81, 085105 (2010).
- García-Martínez et al.  N. A. García-Martínez, A. G. Grushin, T. Neupert, B. Valenzuela, and E. V. Castro, Phys. Rev. B 88, 245123 (2013).
- Daghofer and Hohenadler  M. Daghofer and M. Hohenadler, Phys. Rev. B 89, 035103 (2014).
- Capponi and Läuchli  S. Capponi and A. M. Läuchli, Phys. Rev. B 92, 085146 (2015).
- Motruk et al.  J. Motruk, A. G. Grushin, F. de Juan, and F. Pollmann, Phys. Rev. B 92, 085147 (2015).
- Castro et al.  E. V. Castro, A. G. Grushin, B. Valenzuela, M. A. H. Vozmediano, A. Cortijo, and F. de Juan, Phys. Rev. Lett. 107, 106402 (2011).
- Grushin et al.  A. G. Grushin, E. V. Castro, A. Cortijo, F. de Juan, M. A. H. Vozmediano, and B. Valenzuela, Phys. Rev. B 87, 085136 (2013).
- Sorella et al.  S. Sorella, Y. Otsuka, and S. Yunoki, Sci. Rep. 2, 992 (2012).
- Kurita et al.  M. Kurita, Y. Yamaji, and M. Imada, Phys. Rev. B 94, 125131 (2016).
- Volpez et al.  Y. Volpez, D. D. Scherer, and M. M. Scherer, Phys. Rev. B 94, 165107 (2016).
- de la Peña et al.  D. S. de la Peña, J. Lichtenstein, and C. Honerkamp, Phys. Rev. B 95, 085143 (2017).
- Kaneko et al. [2016a] R. Kaneko, L. F. Tocchio, R. Valentí, and C. Gros, Phys. Rev. B 94, 195111 (2016a).
- Golor and Wessel  M. Golor and S. Wessel, Phys. Rev. B 92, 195154 (2015).
- Scherer et al.  D. D. Scherer, M. M. Scherer, and C. Honerkamp, Phys. Rev. B 92, 155137 (2015).
-  The variational Monte Carlo code is based on a code first developed by Robert Rüger. See R. Rüger, “hVMC”, https://github.com/robertrueger/hVMC (2011) (unpublished).
- Gros  C. Gros, Ann. Phys. 189, 53 (1989).
- Jastrow  R. Jastrow, Phys. Rev. 98, 1479 (1955).
- Kaneko et al. [2016b] R. Kaneko, L. F. Tocchio, R. Valentí, F. Becca, and C. Gros, Phys. Rev. B 93, 125127 (2016b).
- Feynman  R. P. Feynman, Phys. Rev. 94, 262 (1954).
- Tocchio et al.  L. F. Tocchio, F. Becca, and C. Gros, Phys. Rev. B 83, 195138 (2011).
- Otsuka et al.  Y. Otsuka, S. Yunoki, and S. Sorella, Phys. Rev. X 6, 011029 (2016).
- Sorella and Tosatti  S. Sorella and E. Tosatti, Europhys. Lett. 19, 699 (1992).
- Matsumoto et al.  M. Matsumoto, C. Yasuda, S. Todo, and H. Takayama, Phys. Rev. B 65, 014407 (2001).
- Wenzel et al.  S. Wenzel, L. Bogacz, and W. Janke, Phys. Rev. Lett. 101, 127202 (2008).
- Jiang  F.-J. Jiang, Phys. Rev. B 85, 014414 (2012).
- Yasuda and Todo  S. Yasuda and S. Todo, Phys. Rev. E 88, 061301 (2013).
- Tosatti and Anderson  E. Tosatti and P. Anderson, Jpn. J. Appl. Phys. 13, 381 (1974).
- Carpinelli et al.  J. M. Carpinelli, H. H. Weitering, E. W. Plummer, and R. Stumpf, Nature (London) 381, 398 (1996).