Hamiltonian models of lattice fermions solvable by the meron-cluster algorithm

We introduce a half-filled Hamiltonian of spin-half lattice fermions that can be studied with the efficient meron-cluster algorithm in any dimension. As with the usual bipartite half-filled Hubbard models, the na\"ive $U(2)$ symmetry is enhanced to $SO(4)$. On the other hand our model has a novel spin-charge flip ${\mathbb Z}^C_2$ symmetry which is an important ingredient of free massless fermions. In this work we focus on one spatial dimension, and show that our model can be viewed as a lattice-regularized two-flavor chiral-mass Gross-Neveu model. Our model remains solvable in the presence of the Hubbard coupling $U$, which maps to a combination of Gross-Neveu and Thirring couplings in one dimension. Using the meron-cluster algorithm we find that the ground state of our model is a valence bond solid when $U=0$. From our field theory analysis, we argue that the valence bond solid forms inevitably because of an interesting frustration between spin and charge sectors in the renormalization group flow enforced by the ${\mathbb Z}^C_2$ symmetry. This state spontaneously breaks translation symmetry by one lattice unit, which can be identified with a $\mathbb{Z}_2^\chi$ chiral symmetry in the continuum. We show that increasing $U$ induces a quantum phase transition to a critical phase described by the $SU(2)_1$ Wess-Zumino-Witten theory. The quantum critical point between these two phases is known to exhibit a novel symmetry enhancement between spin and dimer. Here we verify the scaling relations of these correlation functions near the critical point numerically. Our study opens up the exciting possibility of numerical access to similar novel phase transitions in higher dimensions in fermionic lattice models using the meron-cluster algorithm.


I. INTRODUCTION
The connections between lattice models of quantum many body physics and continuum field theories remains a forefront topic of research in both high energy and condensed matter physics. In high energy physics one often starts with a continuum field theory, and then introduces a lattice to enable nonperturbative numerical computations. Even though the lattice is viewed as a calculational artifice, it has lead to many profound physical insights. In condensed matter physics on the other hand, one starts with the lattice, which is derived from the physical crystal structure and the connection to continuum field theories emerges in the long distance physics close to a critical point. Despite these contrasting motivations, there are recurring common themes across the two approaches leading to a wonderful synergistic exchange of ideas.
Since our work here may be of interest to both communities, we explain our motivations from both viewpoints in turn. In particle physics, there is a great interest in mechanisms through which fermions acquire masses. While fermion masses are usually thought to arise from fermion bilinear operators, it is well known that strong four-fermion couplings can also generate masses through spontaneous symmetry breaking [1,2]. The role of four-fermion couplings is still poorly understood since such couplings are irrelevant perturbatively at the free fermion fixed point in two or more spatial dimensions and have to be of finite strength when they generate a mass, requiring a non-perturbative analysis. From a Lorentz-symmetry point of view the interactions constructed with scalar or pseudo-scalar fermion bilinears are referred to as Gross-Neveu (GN) couplings [3], while those constructed with vector or pseudo-vector fermion bilinears are called Thirring couplings [4]. These relativistic four-fermion field theories can be studied through the following Lagrangian where we use the Euclidean signature and we will keep this convention through the paper. The case with two flavors of Dirac fermions that is relevant to our paper is discussed in Appendix B. The fixed point structure of renormalization group (RG) flows in the space of these couplings is rich and reveals that the GN and Thirring couplings can have very different physics at strong couplings depending on the number of fermion flavors especially in 2 + 1 dimensions [5][6][7].
In condensed matter physics one usually works with Hamiltonians instead of Lagrangians, which are motivated more naturally by the physics of the underlying material. For many important systems, such as one dimensional metals and two and three dimensional semi-metals, the long distance physics is described by relativistic four-fermion field theories. These Hamiltonians usually have two parts: an electron hopping term that in certain cases can result in massless Dirac fermions and an electron-electron interaction term that is modeled by a four-fermion interaction. These kind of models are typified by the iconic Hubbard model, where c iα destroys an electron on lattice site i with spin α and n iα = c † iα c iα . While the choice of the matrix elements t ij depends on the material to be modeled, by now many physically motivated cases give Dirac fermions whose masslessness is protected by lattice symmetries or topology [8]. The topic of interest is how symmetries are spontaneously broken at strong interactions that covert the semi-metal to an insulator and the nature of the quantum critical point [9]. These models, although motivated from the physics of electrons in crystals, are natural Hamiltonian discretizations of the fourfermion mass generation problem in Eq. (1). The complexities of the RG flows of the continuum theory manifest themselves here in rich phase diagrams that are currently poorly understood, making this an intensely studied area of research.
Generally speaking the only available unbiased method to study such strongly coupled four-fermion lattice models is the Monte Carlo (MC) method [10,11]. Given the great interest in this topic, there has been a large body of work from both communities that has led to important progress and also a number of unresolved issues. Most work found in the high energy literature has been focused on understanding phase diagrams and fermion mass generation at strong four-fermion couplings [12][13][14][15][16][17]. The value of the critical number of fermion flavors below which fermion bilinear condensates in the Thirring model has been an interesting quantity to compute [18][19][20]. Lattice calculations with a local Lagrangian typically suffer from the fermion doubling problem. Due to this problem simple lattice four-fermion couplings get mapped into a linear combination of many continuum four-fermion couplings, which makes it difficult to understand which continuum model is being explored without careful fine tuning. For example recently it was recognized that the lattice Thirring model and the lattice GN models in 2 + 1 dimensions seem to flow to the same continuum theory at the critical point [21].
The field has become more exciting recently due to new interest from condensed matter physics, which was originally inspired in part by the physics of graphene. Again the focus has been on mass terms generated by four-fermion interactions on the honeycomb lattice Hubbard model. Here the motivation is to study the semi-metal to antiferromagnetic transition and its universality class [22]. While the original work focused on SU (2), generalizations to SU (N ) have been carried out [23]. Interacting spinless lattice Hamiltonian models have also been demonstrated to be free of sign problems, which then allows one to study the simplest fermion mass generation mechanism at a quantum critical point, the chiral Ising transition [24,25]. Perhaps the most striking recent result is that strong four-fermion couplings may also create fermion masses without symmetry breaking [26][27][28][29]. In 2 + 1 dimensions there is growing evidence that this mechanism of symmetric fermion mass generation may be a feature of continuum quantum field theory since it appears to be connected to an exotic quantum critical point [30][31][32][33]. In 3 + 1 dimensions it has been suggested that this mechanism may be helpful to construct lattice chiral gauge theories [34][35][36]. Another interesting scenario that has been proposed is the existence of second order quantum phase transitions between different massive fermion phases [37,38]. While such transitions cannot easily be understood within the Landau-Ginzburg paradigm, there is speculation that they may arise through a deconfined quantum critical point with emergent gauge fields [39]. It is also believed that some of them could even be driven by topological terms in the low energy bosonic theory [30,40]. A four-fermion realization of this transition has also been proposed [41,42], where the existence of symmetry enhancement has also been studied at multi-critical points [43].
Despite the power of MC methods, two bottlenecks quickly arise. The first one is the sign problem that limits the type of fermion models that can be studied. Furthermore, even when sign problems are solved, most traditional MC studies of quantum critical points in fermionic systems can only be performed on rather small lattice sizes, due to the poor scal-ing of the computational time with system size. Given these hurdles it is clearly of great interest to find new four-fermion models that are sign problem free and can be accessed with efficient MC algorithms. Recently, alternate types of fermion MC methods have become available, which are able to reach somewhat larger lattice sizes. For example, the recently proposed fermion-bag approach [44][45][46] can study system sizes involving up to 10, 000 sites [47]. The meron-cluster algorithm, a precursor to the fermion-bag approach, is even more efficient while being applicable only to a more restricted class of models [48][49][50]. The main motivation behind our current work is to design fermionic models that can be studied on large lattices using the meron-cluster algorithm.
Interestingly the new four-fermion model we introduce in this work is not only amenable to efficient simulations, it also has some novel physical features. In brief we demonstrate below a novel mechanism by which the valence bond solid (VBS) state which breaks translations symmetry is realized in a one-dimensional fermionic system at arbitrarily weak coupling. It is well known that the Hubbard model Eq. (2) at half filling on bipartite lattices has an explicit SU (2) s spin symmetry and a hidden SU (2) c charge symmetry [51]. The hopping term of the Hubbard model has an additional a spincharge flip symmetry Z C 2 under which the two SU (2) symmetries are interchanged. Indeed, as is well known when U is repulsive (attractive) the spin (charge) sector is favored resulting in anti-ferromagnetic (superconducting) correlations. It is interesting to ask what is the fate of the Hubbard model when the interactions added preserve Z C 2 ? Here we show by field theoretic arguments and explicit MC simulations that in one spatial dimension, when the interactions preserve Z C 2 , the system releases the frustration between spin and charge sectors by forming a VBS. This is a novel mechanism for the formation of a VBS in the one dimensional Hubbard model. From a field theory point of view, we argue that our model can be considered as a Hamiltonian lattice regularization of the twoflavor chiral-mass GN model with a spontaneously broken Z χ 2 chiral symmetry. Using the meron-cluster representation we observe that our model is also related to the θ = π phase of the CP 3 model where the charge conjugation symmetry is spontaneously broken [52]. When U = 0 a combination of GN and Thirring couplings is introduced. These new couplings induce a quantum phase transition to a phase where Z χ 2 symmetry is restored, which turns out to be the SU (2) 1 Wess-Zumino-Witten (WZW) model perturbed by a marginally irrelevant coupling. The universality class of the transition in our model has been studied earlier by Affleck et al. [53], who has argued that at the quantum critical point the marginal coupling vanishes, enhancing the symmetry of the theory. This transition has also been well studied numerically within the context of quantum spin-half chain [54,55].
Our paper is organized as follows. In Section II we introduce our lattice model and explain its symmetries. In Section III A we argue that our lattice model is naturally mapped into a continuum model with a variety of four-fermion couplings. We discuss the continuum symmetries and relate them to the lattice symmetries. In Section III B we apply nonabelian bosonization to the continuum four-fermion theory and rewrite the low energy physics in terms of bosonic excitations. This helps us uncover the phase diagram of our lattice model. In Section III C, we explain how the symmetries of the lattice model are enhanced at the critical point and derive expressions for the spin and dimer correlation functions near the critical point and in the conformal phase. In Section IV we verify the theoretical analysis against MC results obtained using the meron-cluster algorithm and determine the critical point. We also confirm our estimate for the critical point using an exact diagonalization method. In Section V we summarize our results and provide an outlook for the future.

II. OUR LATTICE MODEL
In this section we introduce our lattice model and discuss its symmetries. We focus on models with two flavors (or spins) of lattice fermions which can be annihilated and created at each lattice site j by the usual fermionic Fock operators c jα and c † jα , where α =↑, ↓. The choice of our model is constrained by the ability to use the meron-cluster algorithm as discussed in [49]. While a variety of models fall in this class, here we focus on a particularly simple one whose Hamiltonian is where The symbol i, j refers to a bond between nearest neighbor sites i and j. Our model above remains solvable by the meroncluster algorithm in the presence of some other carefully chosen interactions, for example, the Hubbard interaction, In this work we will study the Hamiltonian This model and its variants can be defined on a bipartite lattice in any dimension and remain solvable using the meron-cluster approach. The constraints of this solvability typically make them strongly interacting and end up within massive phases. However, as we will show in this work, they can still be useful to study phase transitions to massless phases. Here we will study this interesting quantum phase transition of Eq. (6) in one spatial dimension as a function of U . Our model has a variety of interesting lattice symmetries that become manifest when written in terms of Majorana operators. We define two such operators γ 1 j and γ 2 j for spin-up fermions on each lattice site through the relations for even j, and for odd j. Note that n j↑ = 1 2 (−iγ 1 j γ 2 j + 1) in both cases. Similarly we define two more Majorana operators γ 3 j and γ 4 j using the spin-down fermions. In terms of the four Majorana operators the two parts of the Hamiltonian take the form Let us now argue that The SO(4) symmetry becomes obvious when we realize that the following six operators satisfy the so(4) algebra and commute with H. In fact the four operators γ µ j transform as an SO(4) vector. We will argue in the next section that this symmetry is the vector subgroup of the full chiral symmetry group in the continuum. In addition, the H J has a discrete symmetry generated by C ↑ := i j γ 1 j γ 3 j γ 4 j , which we denote as Z C 2 for convenience. It is easy to verify that C ↑ flips the sign of γ 2 j but not the other three Majorana operators. Hence along with C ↑ the four Majorana operators actually transform under the O(4) group. Note that in the Dirac fermion language, C ↑ c j↑ C ↑ = (−1) j c † j↑ is the familiar particle-hole transformation on the spin-up fermions. Therefore it can also be viewed as the spin-charge flip transformation that is more familiar in condensed matter physics. On a regular lattice H is also invariant under translations by one lattice site T † a c jα T a = c j+1α . As we will argue in the next section, there is a remnant discrete Z χ 2 subgroup of the continuum chiral symmetry group buried within T a . When U = 0, the above lattice symmetries actually lead to an interesting degeneracy in the energy spectrum when the lattice size is a multiple of four, as will be shown in Appendix A. In fact all energy levels are evenly degenerate, and in particular the ground state is doubly degenerate.
In order to compute quantities in our lattice model with the meron-cluster algorithm we use the continuous time formulation of the partition function, Notice that here we use Tr to denote the trace over the full Hilbert space. Later we will use tr to denote the trace over local Hilbert spaces or local fields. The integrals over Euclidean time are always assumed to be time ordered such that β > t k > .... > t 2 > t 1 > 0. The choice of H b leads to a simple formula for the trace in the fermionic Hilbert space. In particular it does not contain any determinants of large matrices, as is the case in the traditional auxiliary field methods. Instead, it can be shown that where { 1 , 2 , ...} is a set of loops and W ( i ) is the weight associated with the loop i [49]. The loops can be identified by introducing two parallel bonds for each H b at the appropriate imaginary time, as illustrated in Fig. 1. As can be seen from the figure, these bonds naturally divide the lattice into disconnected loop clusters. When the trace over the fermionic Hilbert space is performed each cluster gets a weight W ( ) = 2(1 ± e −U t /2 ), where t is the linear temporal size of the loop. The sign associated with each loop is given by (−1) nt+n b /2+1 , where n t is the number of temporal winding of the loops and n b is the number of bonds in the loop. The fermionic nature of the problem is hidden in this sign. Note that when U = 0 clusters with a negative sign (merons) are naturally forbidden. On the other hand when U is very large all clusters are allowed and from the cluster representation our model becomes identical to the Heisenberg spin-half chain [56,57].
FIG. 1. Illustration of a configuration of bonds that naturally divides space-time into loops. The fermionic trace is given as a product of weights associated with each loop as described in the text.

III. CONTINUUM FIELD THEORY ANALYSIS
In this section, we review the results for the non-abelian bosonization of the Hubbard model pioneered by Affleck et al. [53,58], and the field theoretic picture for the phase transition from critical to VBS in one dimension. While this scenario is now standard, a new aspect of our model is the spin-charge flip symmetry Z C 2 which is present in H J at a lattice level even beyond the quadratic level. This extra symmetry appears in an interesting way in the continuum description, the RG flow and the resulting phase diagram.

A. Connection to the Gross-Neveu-Thirring Model
In this section we will identify the 1 + 1 dimensional continuum QFT which is described by our lattice model in one spatial dimension, given in Eq. (6). Since the model is strongly interacting, in principle such an identification can be questionable. Still we can try to perform a tree level analysis by expanding the lattice Hamiltonian in terms of modes near the Fermi points of the free Hamiltonian and then include interactions perturbatively. Such an analysis will teach us how the lattice symmetries are embedded within the continuum ones [59]. In particular we show below that the H J term can be identified with the strongly coupled lattice-regularized two-flavor chiral-mass GN model, while the H U term introduces a combination of the usual GN coupling and a Thirring coupling. More details of these couplings and the mapping are explained in Appendix B.
In order to perform the tree level analysis described above, we introduce a small parameter ε and deform H J → H ε J as follows, Clearly all lattice symmetries discussed in the previous section are maintained order by order in ε. Expanding in powers of ε, we get where the leading term describes free fermions, and is a four-fermion interaction. Higher order terms H are irrelevant perturbatively near the free fermion fixed point and therefore their exact forms will not affect our discussion. By focusing on H we can identify the continuum model and also map the lattice symmetries to the continuum. Our lattice model with ε = 1 will be regarded as merely another lattice discretization of the same continuum model.
In order to analyze the low energy physics near the free fermion fixed point, we expand our lattice fields c jα in terms of smooth fields ψ α,L (aj) and ψ α,R (aj) near the two Fermi momenta ∓k F , where k F = π/2a and a is the lattice spacing. Without loss of generality we choose for both α. Inserting Eq. (18) into H J , we can expand in powers of lattice spacing a to obtain the leading low energy effective Hamiltonian in the continuum as This is the Hamiltonian for two-flavor free massless Dirac fermions. It is easy to verify that the above free Hamiltonian is invariant under O(4) L × O(4) R which we refer to as the full continuum chiral symmetry. The SO(4) transformation in each chiral sector can be decomposed into spin SU (2) s and charge SU (2) c transformations, which becomes explicit if we introduce a 2 × 2 matrix of operators, and define the spin and charge transformations as , where S L(R) and Q L(R) are SU (2) s and SU (2) c matrices in each chiral sector. When written in terms of Ψ L(R) , the Hamiltonian is clearly invariant under the above transformation. However, since (−S L(R) , −Q L(R) ) should be identified with (S L(R) , Q L(R) ), the symmetry group of the continuum Hamiltonian H (2)cont J is actually (SU (2) s × SU (2) c ) Z 2 ∼ = SO(4) in each chiral sector. These transformations are generated on the Hilbert space by the spin current operators and the charge current operators In addition the Hamiltonian is also invariant under two inde- , under which J i sL(R) and J i cL(R) exchange with each other. Including them the free Hamiltonian is indeed invariant under the O(4) L × O(4) R symmetry as stated above.
We know from Eq.
a belongs to the vector subgroup of the chiral symmetry group, and therefore T a is effectively a Z 2 subgroup of the chiral symmetry group, hence denoted by Z χ 2 .
Using a similar analysis as above we can identify the following continuum operators that corresponds to H where In order to understand the symmetries of this term it is useful to express it in terms of spin and charge currents using the In this form it is easy to see that this term breaks the chiral O(4) symmetry of the free theory down to the diagonal O(4) L=R subgroup. However, the Z χ 2 group generated by T a survives. Therefore H Coming to the Hubbard interaction H U in Eq. (5), the corresponding continuum Hamiltonian was already derived by Affleck et al. [60] and is given by up to irrelevant pieces. Note that H U preserves all the symmetries of H J except the spin-charge flip symmetry which means the O(4) × Z χ 2 symmetry is reduced to SO(4) × Z χ 2 . The last two terms in Eq. (27) only renormalize the speed of light and do not introduce any interactions, which is clearer in the language of conformal field theory (CFT). Throwing away these chiral terms, and normalizing the kinetic term to have a coefficient of one, we obtain the final continuum Hamiltonian of our lattice model as where λ s = 2(ε − U 4J ) and λ c = 2(ε + U 4J ). Assuming the H (6,8) terms do not modify the above analysis except perhaps to change the couplings λ s and λ c , we argue that Eq. (6) is a lattice regularization of Eq. (28). In Appendix B we construct the Lagrangian of Eq. (28) and identify it with a linear combination of GN-Thirring couplings.
In summary the above discussion shows that the O(4) symmetry of our lattice model in Eq. (6) can then be identified with the diagonal O(4) subgroup of the continuum Hamiltonian Eq. (28). The translation by one site symmetry T a can be identified as the remnant chiral symmetry Z χ 2 discussed above. While it is known that lattice translation is part of the continuum chiral symmetry, we have shown explicitly how it is embedded along with the flavor symmetries. When U = 0, Our model is on the λs = λc line due to the spin-charge flip symmetry, and is in a massive phase described by a lattice regularized chiral-mass GN model with a spontaneously broken Z χ 2 symmetry. When U increases to Uc, λs = 0 and our model is at the critical point of a second order phase transition, whose low energy physics is described by the SU (2)1 WZW model. When U > Uc, λs < 0 and a marginally irrelevant coupling would modify the low energy WZW theory mildly.

B. Bosonization and Phase Diagram
In order to understand the physics of our lattice model we begin by discussing the RG flows of the continuum model Eq. (28) near the free fermion fixed point. The one-loop βfunctions in the continuum are given by [60] which shows that the spin and charge currents are completely decoupled in the low energy theory, i.e., terms involving both spin currents and charge currents and respecting the SO(4) symmetry are irrelevant. The RG flow diagram based on this β-function is shown in Fig. 2. The couplings are relevant when λ s(c) > 0, and irrelevant when λ s(c) < 0.
Our lattice model falls somewhere on the (λ s , λ c ) plane. For example when U = 0 the model must be on the λ c = λ s line due to the spin-charge flip symmetry, which is identified with the chiral-mass GN model with the well known physics of asymptotic freedom. Here we expect M (x) = 0 due to the spontaneous breaking of the Z χ 2 chiral symmetry [3], which is related to the dimer order parameter and verified by our MC results presented later. When U > 0 our model moves away from the spin-charge symmetric axis towards the conformal phase in the second quadrant where λ s < 0 and λ c > 0. This implies the existence of a quantum phase transition at some critical value U = U c , which we find to be second order.
The most convenient way to understand the physics of the conformal phase is through the language of non-abelian bosonization [61], where Eq. (28) is mapped to the following non-linear sigma model, g and h are SU (2) group-valued fields, describing the spin and charged sectors respectively, and is the action for the Wess-Zumino-Witten (WZW) theory in the Euclidean signature. While the first term in S WZW [g] involves an integration over 2-sphere space-time, the second term involves an integration over a 3-disk whose boundary is the 2-sphere. The currents J i sL(R) and J i cL(R) are given by where ∂ = ∂ x + i∂ t and∂ is the complex conjugate of ∂. Note that although in the Euclidean space those currents are usually called holomorphic and anti-holomorphic and denoted by J andJ , here we keep the convention from the Minkowski space and still call the currents as left and right handed for convenience. The bosonized action Eq. (30) enjoys almost the same symmetries as Eq. (28), except for the fact that (g, h) and (−g, −h) should be identified. This constraint comes from the fact that we have written the bosonized theory with an SU (2) s × SU (2) c symmetry rather than SO(4). This means that observables in the fermionic theory will always be mapped into terms that do not violate this identification. For example tr g itself is not an observable. Besides, the chiral symmetries of S WZW in each sector are SO(4) instead of SU (2) L × SU (2) R . Also note that under spin-charge flip symmetry g ↔ h.
Armed with these bosonization results we see that when U > 0 the charge excitations have higher energies than the spin excitations. Thus, the long distance physics of the conformal phase in the second quadrant must be described by This long distance physics can also be seen directly in the lattice model, where the charge sector becomes energetically less favorable when U > 0 and the low energy physics is described by a Heisenberg spin-half chain. In a famous theorem, Lieb, Shultz and Mattis showed that spin-half chains can only be in one of two possible ground states [62]: a dimerized phase where the translation symmetry is broken or a critical phase. In fermionic viewpoint, the dimerized phase is the chirally broken GN phase while the critical phase is the CFT described by Eq. (36) with λ s ≤ 0. This implies the universality class of the transition at U = U c in our model can also be understood from simpler spin-half chain models. One such example is with both J 1 , J 2 > 0. When J 2 = 0, this model is the usual Heisenberg spin chain and is well known to be described the WZW conformal phase. Also, as Majumdar and Ghosh pointed out long ago, when J 2 /J 1 = 0.5 the ground state is doubly degenerate due to dimerization, suggesting that the theory at those couplings is in a different phase [63,64]. Previous studies have shown that the phase transition between the two phases occurs at J 2 /J 1 = 0.241167(5) [55,65]. Although the phase transition has also been studied without a sign problem in one-dimensional spin-half chain using the J-Q idea [66,67], our model provides a route to the onedimensional WZW critical phase to VBS quantum phase transition starting from a lattice fermionic Hubbard model Hilbert space instead of a spin-half chain Hilbert space. A very narrow region of VBS was also reported in an extended Hubbard model [68].

C. Symmetry Enhancement and Correlation Functions
Since the charge sector decouples when U > 0, the long distance physics of the lattice model seems to only have the diagonal SU (2) spin symmetry and the remnant Z χ 2 chiral symmetry, under which g → iσ 3 giσ 3 . These symmetries can be observed in the continuum low energy theory described by Eq. (36). However, if we set λ s = 0 the continuum theory is invariant under the enhanced chiral transformation g → S L gS † R where S L and S R are two independent SU (2) matrices implementing the left and right spin symmetries of the fermion model, except that S L = S R = −1 needs to be quotiented out. Based on the RG flow diagram shown in Fig. 2 we see that by tuning to the critical point U = U c we can indeed set λ s = 0, and thus we must see the enhanced chiral symmetry there. In our work we look for this symmetry enhancement, although observing it could be non-trivial since there will always be higher dimensional operators that break the symmetry on the lattice. Also, when U = U c but close to it, since λ s is a marginal coupling, the logarithmic corrections due to it would be visible in the lattice MC data at long distances. Here we derive the form of these corrections and use them in our analysis.
It was explained in [58,69] that the symmetry enhancement is visible in correlation functions of spin and dimer operators on the lattice defined through the relations The leading continuum terms of these two operators are the primary fields ϕ S = tr h tr(giσ 3 ) and ϕ D = tr h tr g of the WZW model. In the fermionic language these operators are given by ϕ S ∝ ψ †α L σ 3 αβ ψ β R + h.c. and ϕ D (x) ∝ M (x) defined in Eq. (25). We see that ϕ D can be used as an order parameter on the lattice for Z χ 2 . Since ϕ S and ϕ D transform into each other under the SO(4) chiral transformations, their correlation functions will be related to each other. Using the well known techniques in CFT [70][71][72], their conformal dimensions have been computed to be h S = h D = 1 2 , and the two point correlation functions of the spin and dimer are for large values of r at the critical point U = U c . In order to be able to fit our MC data away from the critical point we have to include the logarithmic corrections to Eq. (40), due to the presence of the marginal coupling λ s . We begin with the β-function of the marginal coupling λ s This is the same as Eq. (29) except that we have replaced the energy scale µ with a length scale r, which reverses the sign of the beta function. Integrating to first order, we obtain where λ 0 is the coupling at some short distance length scale r 0 . Note that the solution λ s (r) above is meaningful for very large values of r only when λ 0 < 0 (i.e., when U > U c and the theory is in the conformal phase), and therefore the coupling λ s (r) is marginally irrelevant. When λ 0 > 0 the coupling λ s (r) diverges at r = r 0 e 2π/λ0 , which is the remnant of the Landau pole of perturbation theory and signals new physics at long distances due to the formation of a fermion mass. This means in the broken phase the above form of λ s (r) will only be valid close to the critical point and small values of r.
With the caveats discussed above, we can use λ s (r) in Eq. (42) in both phases to derive the corrections for G i (r). Noting that it depends on r both explicitly and through λ s (r), we can write G i (r) =: G i (λ s (r), r) to be more precise. We note that it satisfies the following RG equation: where h i is the conformal dimension of ϕ i at the critical point and γ i (λ s (r)) is the anomalous dimension of ϕ i induced by the marginal operator. Integrating this equation we get ln G i (λ s (r), r) G i (λ 0 , r 0 ) = r r0 −2(h i + γ i (λ s (r )))d ln r , (44) from which we can derive the corrections due the presence of the marginal operator once we know the γ i (λ s (r)). First let us focus on the conformal phase and near the critical point, i.e. λ 0 0. In 2d CFT the conformal dimension h i of the operator ϕ i can be calculated using the finite size energy E i = 2πh i /L of the state |ϕ i through the state-operator correspondence [73]. Then the anomalous dimension γ i of the operator is related to the change in this energy due to the presence of the marginal operator. In our case, to leading order in λ s , this leads to the expression [53], where In the above formula s L = s R = 1/2 for both spin and dimer states, while s tot depends on the state |ϕ i . The values of b i can now be calculated for both spin and dimer fields and are tabulated in Table I below.
Lattice Field WZW field stot sL sR hi bi S i ϕS = tr h tr gσ i 1 1/2 1/2 1/2 1/4 D ϕD = tr h tr g 0 1/2 1/2 1/2 −3/4 Similar expressions have been derived earlier in the context of spin chains [55,74] and verified numerically [55]. Note again that for a fixed value of λ 0 but large values of r, the above expressions make sense only when λ 0 ≤ 0, i.e. when we are in the conformal phase as already mentioned above. On the other hand when r is in a fixed range but λ 0 → 0 the above expressions give us the leading corrections to conformal behavior on both sides of the critical point. We could do a similar analysis in the massive phase. However, in the massive phase the anomalous dimension Eq. (45) obtained in the conformal phase cannot be completely correct. Besides, like the β-function, it was derived in perturbation theory and so there could be some non-perturbative corrections to it due to the presence of the mass scale M . In particular if we assume the correlation function to the leading order is of the form G i (r) ∼ e −αiM r , this would contribute additional non-perturbative terms to the anomalous dimension γ i (λ s ). For example if we assume λ s (r) = − 2π ln M r (49) as a possible definition for the RG invariant mass scale M , then it is easy to see that Since we are only deriving corrections to G(r) due to the marginal coupling perturbatively, in this work we ignore such non-perturbative effects. Then, Eqs. (47) and (48) can still be used in the massive phase near the critical point except that we will find λ 0 > 0 and r is constrained to be such that λ 0 /2π ln(r/r 0 ) < 1.

A. Monte Carlo Results
We have used the meron-cluster algorithm to compute the equal time spin and dimer correlation functions defined through the expressions where O i (r), i = S, D are the spin and dimer operators defined in Eqs. (38) and (39). The symbol r is used for the spatial lattice site j at some fixed time slice. In our work we choose β = L where L is the size of our spatial lattice in lattice units. Based on our estimate for the speed of light this leads to a small effective temperature. In particular we have evidence that physical temperatures and physical length scales are related by T phys ≈ 0.25(La) −1 in our lattice model. We computeG i (r) at various lattice sizes and U . In Fig. 3 we illustrate some of our results at U = 0, 0.5, 1 and 1.5 on a lattice size of L = 128. For clarity we focus on the region of 8 ≤ r ≤ 40. The figure clearly shows a non-zero expectation value for the dimer order parameter D j (t) at U = 0 while there is no such expectation value for the spin operator. This is consistent with our expectations that for small values of U the lattice model breaks the Z χ 2 chiral symmetry as explained in Section III B. Another important point to note is that the spin and dimer correlation functions also behave very differently at least for small values of U but slowly begin to become similar by U ≈ 1.5.
Assuming we are in the vicinity of the critical point around U 1.5, we want to fit our MC data to Eqs. (47) and (48). However, since we work on a finite lattice with periodic boundary conditions, the correlation functions receives finite size corrections. Fortunately, in the conformal phase the finite size corrections can be obtained using the map from an infinite plane to a cylinder, which results in replacing r by where L is the spatial size [73,74]. Furthermore, the spin correlation functions clearly shows oscillating behavior due to the higher order operators in the observable with ferromagnetic behavior. Taking these corrections into account, we make the following ansatz for the lattice correlation functions at U U c on a finite lattice, where we have introduced a single fit parameterλ 0 = λ 0 /(1 + λ0 2π lnr 0 ), which is an RG invariant quantity and numerically equals the coupling measured atr 0 = 1 in the lattice unit. The critical point is obtained whenλ 0 = 0.
We have performed a combined fit for both G S (r) and G D (r) to the form Eqs. (53) and (54) at each fixed value of U ≥ 1.5, which is tabulated in the combined fit column of Table II. Each of these is a four parameter fit involving A 1 , A 2 , B,λ 0 and uses all data from L = 64, 80, 96, 128 and 12 ≤ r ≤ 40. As can be seen from these results, our data fit well to the form Eqs. (53) and (54) for all couplings in the range 1.5 ≤ U ≤ 2.0. Since at the critical point we expect λ 0 ≈ 0, applying conservative systematic errors related to our fitting procedures we estimate U c = 1.75 (5). In order to visualize that the correlation functions indeed obeys the power law near the critical point, we plot them in log scale in Fig. 4. We only plot at even r in order to avoid the distracting oscillation.
The quality of the combined fit column in Table II becomes poor for U > 2.0, and changing the fitting range of r does not seem to improve the fits much. In fact the fit seems to have a very bad reduced chi-squared χ 2 ν = 9.68 when U → ∞. This seems a bit surprising since the forms described by Eqs. (53) and (54) must work well in the entire conformal phase. One reason is that our data is very precise at U = ∞, where H is the same as the Heisenberg spin chain, and therefore the results become sensitive to higher order terms in the lattice model which are ignored in the theoretical analysis. In this limit the meron clusters have the same weight as the nonmerons, and therefore we do not need to check merons and the algorithm is much more efficient. However, we believe there is more to the story here. Fortunately, there are precise asymptotic results for the spin and dimer correlation functions in the Heisenberg spin chain [74][75][76]: and indeed they are consistent with Eqs. (53) and (54) in thẽ r → ∞ limit. More precisely in the limit U → ∞, we must observe the following constraints among the fit parameters However results from the combined fit give us 2πA 1 (−λ 0 ) 1/2 ≈ 1.06, 4π 2 A 2 ≈ 1.06 and 8π 3 B(−λ 0 ) −3/2 ≈ 1.18. While the first two constraints are not far from expectations, the last one seems quite a bit off. We want to check how the fits change if we allow forλ 0 to be different in the spin and dimer correlation functions. Therefore we performed separate fits of our data, which are tabulated in the spin fit and dimer fit columns in Table II, and plotted in Fig. 5. Note that in the vicinity of the critical point both values tend to become small as expected. Focusing on the U → ∞ limit, the separate fits now suggest 2πA 1 (−λ s 0 ) 1/2 ≈ 1.01, 4π 2 A 2 ≈ 1.09 and 8π 3 B(−λ d 0 ) −3/2 ≈ 0.98. Although the deviation in A 2 still seems to be large, this is understandable because A 2 is a higher order correction. Now they are in much better agreement with the above constraints, suggesting that for some reason that we do not yet understand the values ofλ 0 for spin  and dimer begin to drift apart as we go away from the critical point.

B. Exact Diagonalization Results
In order to confirm that our estimate for the critical point is reasonable, we also use an alternate idea outlined in [54], based on the low energy physics of the WZW model Eq. (36) that emerges at the critical point. Since the WZW model is invariant under the enhanced chiral transformations with two independent SU (2) ("left" and "right") symmetries, the energy eigenstates can be labeled with spin quantum numbers (s L , s R ). It is known that the ground state has (s L , s R ) = (0, 0) with momentum 0, while the first excited state has (s L , s R ) = ( 1 2 , 1 2 ) with momentum π [53]. However, since the lattice model is only invariant under the diagonal SU (2) subgroup, the energy eigenstates on the lattice will only be la-beled by the total spin s tot . The four-fold degeneracy requires fine tuning to the critical point where the singlet (s tot = 0) and the triplet (s tot = 1) state together form an (s L , s R ) = ( 1 2 , 1 2 ) state. This suggests an alternative method to determine the critical point: we can compute the lowest five energy eigenvalues as a function of U and L using an exact diagonalization method and locate the coupling where the first excited state becomes four-fold degenerate. When L is finite, the critical coupling where the energies of these two total spin sectors cross can be referred to as a pseudo-critical point U c (L). This point turns into the true critical point as L becomes large.
In order to implement the above idea we have computed the first five eigen-energies by diagonalizing the Hamiltonian on small lattices with the coordinate descent method [77,78]. The behavior of the lowest five states as a function of U at L = 14 and L = 16 is plotted in Fig. 6. We observe that in the broken phase (small U ) the ground state and the first excited state turn out to be spin singlets with s tot = 0. The next three degenerate excited states form a triplet with s tot = 1. In contrast when U > U c , the triplet states have lower energy as compared to the singlet state. We thus can determine U c (L) as the coupling where the triplet and singlet states cross, which are tabulated in Table III and plotted in Fig. 7 as a function of 1/L 2 . Based on Fig. 7 we can estimate U c (L) ≈ 1.746 (1), which is consistent with our estimate in the previous section and with the estimate using finite size scaling in the conference proceeding published earlier [59]. Looking more closely at Fig. 6 we observe some peculiarities. For example at U = 0 we note that the ground state is degenerate at L = 16 but not at L = 14. We have explained the reason for this degeneracy in Appendix A. In fact we show that there is an interesting difference in the energy spectrum when the lattice size is L = 4n versus when it is L = 4n − 2, where n = 1, 2, · · · . Moreover, this difference permeates even away from U = 0 and is observed in the dramatically different values of U c (L) when L is small. The difference however decreases rapidly as L increases, and both of them approach to the true critical point as expected. Another peculiarity comes from the momentum quantum number k for the first five energy eigenstates. We note that k flips between π and 0 when comparing L = 4n with 4n + 2. This is as expected because similar phenomena are also observed in the Heisenberg spin chain, with k = 0 for L = 4n and k = π for L = 4n + 2 [79]. The extra π phase in our model compared to the Heisenberg spin chain is due to the fermionic nature of our model, since there is an intrinsic extra minus sign when the system is translated by one site when L = 2n.
While it is not easy to extend our meron-cluster algorithm for the more general Hamiltonian with the parameter ε introduced in Eq. (14), we can extend the above exact diagonalization method to determine the critical point for arbitrary ε. For example, when ε = 0.1 and J = 1, U c (L) is tabulated in Table IV. When ε is small the perturbative analysis should be a good guide and we obtained at leading order U c = 4Jε = 0.4 in Section III A, which is in good agreement with Table IV. Using the exact diagonalization method we have also confirmed that our model at U = 0 is indeed in a massive phase with spontaneously broken Z χ 2 symmetry. In such phase the energy gap to the first excited state is expected to become exponentially small as L becomes large, while the gap to the second excited state remains non-zero at L → ∞. We plot these gaps in Fig. 8, where the solid lines are exponential fits. These features are qualitatively visible, along with the peculiar degeneracy in our model.

V. SUMMARY AND CONCLUSIONS
In this work we constructed a strongly correlated lattice fermion Hamiltonian that was solvable by the meron-cluster  Table III. algorithm. We were able to add the Hubbard coupling U to our model without losing this property. This combined lattice model had SO(4) × Z χ 2 symmetry for all values of U , and an extra spin-charge flip symmetry Z χ 2 at U = 0. While the lattice model can be formulated in any dimension, here we studied it in one spatial dimension. In order to study the phase diagram of the model as a function of U we used non-abelian bosonization to relate our model to two decoupled SU (2) 1 WZW models with marginal couplings that depend on U . In particular we discovered that the model undergoes a quantum phase transition between a massive phase where the Z χ 2 symmetry is spontaneously broken and a conformal phase without such symmetry breaking. The massive phase is the well known asymptotically free chiral-mass GN model with massive fermions and spontaneous symmetry breaking, while the conformal phase is the single SU (2) 1 WZW model. Using the meron-cluster method and exact diagonalization we provided further evidence of this scenario.
We view the current study as just a first step in demonstrating that meron-cluster algorithms can be useful to study phase diagrams of lattice fermion systems where fermions become massive, while the low energy bosonic physics is still interesting either due to frustration or the presence of topological terms. Extending our model to 2 + 1 dimensions would be very interesting. For example in [42] it was argued that in an extended Hubbard model very similar to ours, there is a direct second order phase transition between a VBS phase and an antiferromagnetic phase, and an enhanced SO(5) symmetry is expected to emerge at the transition. The critical exponents at this exotic transition were computed in the earlier studies, but unfortunately they do not match those obtained from a loop gas formulation of the same phase transition [80]. On the other hand the critical exponents of the fermionic realization seem consistent with the bounds obtained from conformal bootstrap [81]. As far as we know, this controversy remains unresolved, partially because the fermionic model was studied on rather small lattices with less than 500 lattice sites, due to difficulties associated with auxiliary field fermion MC algorithm. We are currently exploring if our model and its extensions allow us to study this transition more efficiently using the meron-cluster approach.

ACKNOWLEDGMENTS
We would like to thank Emilie Huffman and Hersh Singh for helpful discussions and pointing us to the literature on in-teresting fermionic quantum critical behavior. SC would like to thank Uwe-Jens Wiese for discussions. HL would like to thank Zhe Wang for sharing his computer codes on the coordinate descend method, and Xin Zhang for helpful discussions. The material presented here is supported by the U.S. Department of Energy, Office of Science, Nuclear Physics program under Award Numbers DE-FG02-05ER41368. RKK was supported in part by NSF award DMR-1611161. This work used computational resources provided by the Extreme Science and Engineering Discovery Environment (XSEDE) [82], which is supported by National Science Foundation grant number ACI-1548562.
This implies that the pair of states with k = ±k 0 but all other quantum numbers being the same will be degenerate as long as k 0 = 0, π/a. This is easily understandable since a state with a fixed lattice momentum will have a partner with a negative momentum and both will have the same energy as long as parity is a symmetry of the theory.
Interestingly, we will now show that there are additional pairs of degenerate states due to the C ↑ operator, which has the following properties: = e −ia k k(c π/a−k↑ c † π/a−k↑ +c † k↓ c k↓ ) = e −i k (π−ak)c k↑ c † k↑ +akc † k↓ c k↓ = e −i k (ak−π)(c † k↑ c k↑ −1)+akc † k↓ c k↓ and Using these relations it is easy to show that C ↑ |E, k, l s , s 3 , l q , q 3 , α ∝ |E, k + (s 3 − q 3 + L/2 + 1) π a , l q , q 3 , l s , s 3 , α . (A8) This relation shows that additional pairs of state can have the same energy. For example when l s = l q or s 3 = q 3 , the two states |E, k, l s , s 3 , l q , q 3 , α and |E, k + (s 3 − q 3 + L/2 + 1) π a , l q , q 3 , l s , s 3 , α are different but have the same energy irrespective of the value of k. Thus, the only situation where an energy eigenstate could remain non-degenerate is when k = 0 or π/a, l s = l q = j and s 3 = q 3 = m. In this case Eq. (A8) mixes the states |E, k, j, m, j, m, α and |E, (k + L/2 + 1)π/a, j, m, j, m, α . Indeed when L = 4n − 2 these two states are identical and C ↑ cannot be used to pair degenerate states. However, when L = 4n, C ↑ mixes |E, 0, j, m, j, m, α with |E, π/a, j, m, j, m, α . This means the whole spectrum has an even degeneracy when L = 4n. The RG flow diagram is shown in Fig. 9.
Using the relations given by Eqs. (B3), (B4) and (B7) we see that these couplings are related to g 1 , g 2 and g 3 via the linear transformation Notice that the pure Thirring coupling g 3 is obtained when λ s = λ c −λ c = 0. This direction is special because the beta function vanishes for all the couplings and we obtain a line of fixed points. In Table V we summarize the symmetries that are preserved by the various couplings. In our lattice model we introduced two independent couplings H ε J and H U . In order to explore the full three dimensional space discussed above we will need one more independent coupling. Using the same symmetry argument above, we need a lattice interaction which breaks SU (2) c down to U (1) c while preserving particle hole symmetry. The most straightforward way to do this is to include the interaction in the lattice Hamiltonian, where n i = n i↑ + n j↓ is the total fermion number operator at the site i. At the tree level in the continuum limit we can show that which maps to the Lagrangian Therefore we see that H U + H V gives the Thirring coupling, while H U − H V gives the usual GN model. It's interesting to observe that when both U and V are positive, H U favors the spin sector while H V favors the charge sector, and the frustration between them gives the conformal Thirring model. Furthermore, H V can also be included in the meron-cluster algorithm for a range couplings.