Suppression of topological MottHubbard phases by multiple charge orders in the honeycomb extended Hubbard model
Abstract
We investigate the competition between chargedensitywave (CDW) states and a Coulomb interactiondriven topological Mott insulator (TMI) in the honeycomb extended Hubbard model. For the spinful model with onsite () and nextnearestneighbor () Coulomb interactions at half filling, we find two peculiar sixsublattice chargedensitywave insulating states by using variational Monte Carlo simulations as well as the HartreeFock 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 001122type CDW phase, we find a magnetic transition driven by an emergent coupleddimer antiferromagnet on an effective square lattice of singly occupied sites. Possible realizations of the found states are discussed.
I Introduction
Correlation effects in electron systems on a twodimensional honeycomb lattice have been the subject of intensive scrutiny [1], 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 [4], but due to the complex hopping that breaks timereversal symmetry [5; 6]. A related state, the quantum spin Hall (QSH), for which the spin Hall conductance is quantized, emerges in analogy when spinorbit interactions lead to robust spindependent transport [7; 8; 9; 10]. The spin quantum Hall state has also been studied in meanfield theory for chargedordered triangular lattices with spinorbit interactions [11].
Raghu et al. proposed that a topological Mott insulator (TMI) could be stabilized on the extended honeycomb lattice [12]. The QHE would emerge in this scenario from pure Coulomb interactions, with an effective spinorbit interaction being dynamically generated via spontaneous symmetry breaking when the nextnearestneighbor Coulomb interaction is taken into account. Here, in a spinless model, timereversal symmetry is broken spontaneously, while the lattice translational symmetry is preserved; the topological phases are characterized by the Chern number [13]. On the other hand, for QSH in a spinful model, symmetry is broken spontaneously, while timereversal symmetry is preserved; they are characterized by the invariant [14]. Nonlocal Coulomb interactions, which are nonnegligible in graphene [15], may also generate, however, conventional spontaneously symmetrybroken states, such as bond order, charge order, and magnetic order, which then compete with the TMI.
In the spinless Hubbard model with nearestneighbor () and nextnearestneighbor () Coulomb interactions, a Kekulé bond order phase, characterized by a order parameter, has been proposed by meanfield calculations [16]. Exactdiagonalization (ED) [17; 18; 19] and infinite density matrix renormalizationgroup (iDMRG) [20] studies also support the presence of this phase and, further, propose much richer chargedensitywave (CDW) phases. Possible charge instabilities have also been investigated away from half filling [21; 22].
Rich chargedensitywave phases also emerge for the spinful Hubbard model. In addition, when the onsite Coulomb interaction is dominant, collinear antiferromagnetic order appears [23]. 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 [27].
Recent extensive research on the honeycomb extended Hubbard model by ED [17; 18; 19], iDMRG [20], auxiliaryfield quantum Monte Carlo (AFQMC) [28], variational Monte Carlo (VMC) [24], 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 spinorbit 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 onsite and nextnearestneighbor Coulomb interactions, which are relevant for the TMI and CDW states. For simplicity, we drop the nearestneighbor Coulomb interaction . We first use the restricted HartreeFock approximation to get an insight into the plausible CDW states, and then apply the VMC method with a JastrowSlatertype wave function to improve the meanfield state. We find, by the HartreeFock approximation, two types of sixsublattice CDW insulating phases: one is a 220200type CDW phase realized for small and large , and the other is a 001122type CDW phase which appears for larger and . The transition between these two is found to be of first order within the meanfield treatment. When the VMC method is employed, for a fixed large , we find a continuous change from the 220200type to the 001122type CDW phase as is increased. The firstorder transition found in the HartreeFock 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 001122type CDW. Within the meanfield study, the 001122type 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 HartreeFock approximation and the VMC method. In Sec. III, we present the  phase diagram and the properties of the three CDW phases obtained by the HartreeFock 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,
(1)  
Here, denotes the hopping parameter and () denotes the strength of onsite (nextnearestneighbor) Coulomb interaction, as shown in Fig. 1. For simplicity, we only deal with the nextnearestneighbor and neglect the nearestneighbor . Hereafter, we consider repulsive Coulomb interactions () at half filling ().
ii.2 Restricted HartreeFock method
To clarify the plausible ordered phases, we first apply the restricted HartreeFock method. To obtain the meanfield Hamiltonian, we use the HartreeFock decoupling,
(2) 
where are site indices while are spin indices, and determine the order parameters selfconsistently. 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,
(3) 
and magnetic order parameter,
(4) 
per each sublattice . In a chargeuniform 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 nextnearestneighborsite pair . Using this value, the order parameter is defined as
(5)  
where for clockwise (anticlockwise) orientation and is a Pauli spin. This will be a coefficient of dynamically generated spinorbit 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
To investigate the effect of quantum fluctuations beyond the meanfield study, we employ the VMC method [30] by using a JastrowSlatertype wave function [31] given as
(6) 
Here, is an eigenstate of an auxiliary Hamiltonian given by
(7)  
(8)  
(9)  
(10)  
(11)  
(12) 
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 [32] given by
(13) 
This wave function can represent metallic and insulating states with and without charge disproportionation [33]. We optimize the translationalinvariant 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 HartreeFock 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
(14) 
In general, the charge gap in the limit can be estimated as [34; 35]
(15) 
When for , a gap for the particlehole 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 [36]; 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 finitesize charge gap estimated by , which signals the onset of a metallic phase. Hereafter, we first distinguish metallic and insulating phases in a meanfield calculation, and then adopt Eq. (14), as a complementary way, to estimate the size of the charge gap.
Iii Results
iii.1 Meanfield phase diagram
The  phase diagram obtained by the restricted HartreeFock approximation is shown in Fig. 3. In the absence of nextnearestneighbor Coulomb interaction (), a continuous transition occurs at between the nonmagnetic chargeuniform semimetal and antiferromagnetic insulating (AFI) state [37]. These two phases remain for small , consistent with a previous meanfield study [24].
On the other hand, when increases, we find three insulating phases characterized by different charge and magneticorder patterns. As we detail below, two nonmagnetic CDW phases are similar to the ones obtained from meanfield [24], fRG [25; 26], ED [17; 18; 19], and iDMRG [20] 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 HartreeFock approximation.
In the absence of onsite Coulomb interaction (), we find a sixsublattice CDW phase, characterized by three chargerich sites and three chargepoor sites in a unit cell [see Fig. 4(a)]. The charge pattern is like the 220200type, namely, the chargeorder parameters show or . However, there is a charge disproportionation within the three chargerich (chargepoor) sites; one of them is richer (poorer) than the other two. This CDW phase survives for small but nonzero .
At large , we find another sixsublattice CDW phase, characterized by a 001122type charge pattern, namely, , , or [see Fig. 4(b)].
The stability of these CDW phases has been previously tested and they were always found to stabilize in the large region of both spinless [17; 18; 19; 20; 29] and spinful models [24; 25; 26].
Although both CDW phases show sixsublattice orderings, we find a firstorder transition from one to the other. Along a line, the chargeorder 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 meanfield approximation. As we will see later, quantum fluctuations melt the firstorder transition and give a continuous transition or a crossover.
Inside the 001122type sixsublattice 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 firstorder transition between the 220200 and 001122type 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 Mottinsulating state was originally proposed for a large region [12]. To investigate whether the TMI state can be the ground state against the CDW states at the HartreeFock meanfield 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 spinorbit coupling [11]. If such a spontaneous spinorbit coupling is generated dynamically via [12], 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 chargeorder 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 HartreeFock approximation, we find a metallic phase, an antiferromagnetic insulating phase, and both 220200 and 001122type CDW phases. For simplicity, in the present VMC study, we do not consider possible magnetic order in the 001122type CDW phase, namely, the state in Fig. 4(c).
The region of validity of the ordered phases, which are overestimated in the HartreeFock 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 metalinsulator transition point along , estimated to be at , which is in good agreement with the numerically exact value [23]. Furthermore, when the nearestneighbor Coulomb interaction increases, the onsite 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 [28], and our VMC phase diagram agrees with the tendency.
Let us focus on the parameter region along , where the CDW phases are dominant. The chargeorder parameters as a function of are given in Fig. 9. When , or , suggesting that the charge pattern is 220200like. By contrast, when , , , or , suggesting that the charge pattern is 001122like. 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 firstorder transition found in the HartreeFock approximation melts when quantum fluctuations are carefully introduced.
A continuous change of CDW has been reported before in the presence of , which favors twosublattice staggered charge ordering [26]. 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 chargeuniform antiferromagnetic phase is of first order, as in the HartreeFock 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 Discussion
iv.1 Origin of charge order
When is small while is dominant, a chargeuniform state is favorable since prohibits double occupation of electrons. By contrast, for the large and  region, basic chargeorder 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 chargeordering 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 onsite disfavors doubly occupied sites. Consequently, the optimal charge configuration is a 012type 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, 001122type [see Fig. 4(b)] and 021120type [see Fig. 11] charge orderings. Both of them are 18fold degenerate, respectively. In the atomic limit, these two states have exactly the same energy; however, quantum fluctuations lift the degeneracy. The HartreeFock approximation predicts that nonzero hopping favors the 001122type 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 220200type CDW phase.
When , the transition point between the chargeuniform antiferromagnetic collinear state and 001122type CDW state can be estimated by comparing their energies. On a single triangle connected by , the chargeuniform 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 HartreeFock 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 001122type 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 secondorder process. Similarly, for the singly occupied sites connected by the Manhattan distance, , , and are generated by the sixthorder 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 coupleddimer antiferromagnetic model [38; 39; 40; 41] has been previously numerically studied. For , the Néeldimer transition occurs at [39; 40; 41].
In the simple spin system, the nature of this Néeldimer transition is believed to belong to the threedimensional classical Heisenberg universality [40; 41], although its numerical detection is known to be difficult due to the large finitesize 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 spinexchange 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 HartreeFock approximation seems to find these magnetic and nonmagnetic CDW states although the meanfield 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 meanfield magneticnonmagnetic phase boundary significantly inside the CDW phase. Investigating the Néeldimer 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 [3] when the lattice forms a doublelayered triangular structure. This is because the nextnearestneighbor Coulomb interaction on the original honeycomb lattice corresponds to the nearestneighbor one in triangular lattices, and could be dominant when the layers are far enough.
V Conclusions
We have investigated the phase diagram of the Hubbard model on the honeycomb lattice at half filling in the presence of onsite () and nextnearestneighbor () Coulomb interactions. By applying the restricted HartreeFock approximation, we find three sixsublattice CDW insulating phases, namely, (i) a nonmagnetic 220200type CDW phase, (ii) a nonmagnetic 001122type CDW phase, and (iii) a magnetic 001122type CDW phase, as well as the semimetal and antiferromagnetic insulating phases. To investigate the stability of CDW phases beyond the meanfield study, we further apply the VMC method with JastrowSlatertype wave functions. We find that quantum fluctuations destroy the firstorder transition found in the meanfield approximation, and give a continuous change from 220200type CDW to 001122type CDW. On the other hand, when and are large enough, we find a solid 001122type 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 001122type CDW phase are found to behave as in a coupleddimer antiferromagnetic model on the square lattice, where nonmagnetic and magnetic phases are determined by the anisotropy of the spinexchange 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 doublelayered triangular structure.
Acknowledgements.
M.B., R.K., and R.V. acknowledge the support of the German Science Foundation (DFG) through Grant No. SFB/TRR49.References
 Castro Neto et al. [2009] 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. [2010] 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. [2012] L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, Nature (London) 483, 302 (2012).
 Klitzing et al. [1980] K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
 Haldane [1988] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
 Qi et al. [2006] 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 [2006] B. A. Bernevig and S.C. Zhang, Phys. Rev. Lett. 96, 106802 (2006).
 Bernevig et al. [2006] B. A. Bernevig, T. L. Hughes, and S.C. Zhang, Science 314, 1757 (2006).
 König et al. [2007] 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 [2016] Y. Sugita and Y. Motome, J. Phys. Soc. Jpn. 85, 073709 (2016).
 Raghu et al. [2008] S. Raghu, X.L. Qi, C. Honerkamp, and S.C. Zhang, Phys. Rev. Lett. 100, 156401 (2008).
 Thouless et al. [1982] 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. [2011] 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 [2010] C. Weeks and M. Franz, Phys. Rev. B 81, 085105 (2010).
 GarcíaMartínez et al. [2013] N. A. GarcíaMartínez, A. G. Grushin, T. Neupert, B. Valenzuela, and E. V. Castro, Phys. Rev. B 88, 245123 (2013).
 Daghofer and Hohenadler [2014] M. Daghofer and M. Hohenadler, Phys. Rev. B 89, 035103 (2014).
 Capponi and Läuchli [2015] S. Capponi and A. M. Läuchli, Phys. Rev. B 92, 085146 (2015).
 Motruk et al. [2015] J. Motruk, A. G. Grushin, F. de Juan, and F. Pollmann, Phys. Rev. B 92, 085147 (2015).
 Castro et al. [2011] 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. [2013] 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. [2012] S. Sorella, Y. Otsuka, and S. Yunoki, Sci. Rep. 2, 992 (2012).
 Kurita et al. [2016] M. Kurita, Y. Yamaji, and M. Imada, Phys. Rev. B 94, 125131 (2016).
 Volpez et al. [2016] Y. Volpez, D. D. Scherer, and M. M. Scherer, Phys. Rev. B 94, 165107 (2016).
 de la Peña et al. [2017] 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 [2015] M. Golor and S. Wessel, Phys. Rev. B 92, 195154 (2015).
 Scherer et al. [2015] D. D. Scherer, M. M. Scherer, and C. Honerkamp, Phys. Rev. B 92, 155137 (2015).
 [30] 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 [1989] C. Gros, Ann. Phys. 189, 53 (1989).
 Jastrow [1955] 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 [1954] R. P. Feynman, Phys. Rev. 94, 262 (1954).
 Tocchio et al. [2011] L. F. Tocchio, F. Becca, and C. Gros, Phys. Rev. B 83, 195138 (2011).
 Otsuka et al. [2016] Y. Otsuka, S. Yunoki, and S. Sorella, Phys. Rev. X 6, 011029 (2016).
 Sorella and Tosatti [1992] S. Sorella and E. Tosatti, Europhys. Lett. 19, 699 (1992).
 Matsumoto et al. [2001] M. Matsumoto, C. Yasuda, S. Todo, and H. Takayama, Phys. Rev. B 65, 014407 (2001).
 Wenzel et al. [2008] S. Wenzel, L. Bogacz, and W. Janke, Phys. Rev. Lett. 101, 127202 (2008).
 Jiang [2012] F.J. Jiang, Phys. Rev. B 85, 014414 (2012).
 Yasuda and Todo [2013] S. Yasuda and S. Todo, Phys. Rev. E 88, 061301 (2013).
 Tosatti and Anderson [1974] E. Tosatti and P. Anderson, Jpn. J. Appl. Phys. 13, 381 (1974).
 Carpinelli et al. [1996] J. M. Carpinelli, H. H. Weitering, E. W. Plummer, and R. Stumpf, Nature (London) 381, 398 (1996).