Fermion Bag Approach to Lattice Hamiltonian Field Theories

Using a model in the Gross-Neveu Ising universality class, we show how the fermion bag idea can be applied to develop algorithms to Hamiltonian lattice field theories. We argue that fermion world lines suggest an alternative method to the traditional techniques for calculating ratios of determinants in a stable manner. We show the power behind these ideas by extracting the physics of the model on large lattices.


Introduction
Traditionally lattice field theories have been studied by the lattice community using the Lagrangian path integral approach. However, the Hamiltonian formalism does offer some compelling reasons for its application. The first reason has to do with additional symmetries that can be found in Hamiltonian theories and may be relevant to phase transitions. As an example theory, the spectrum of Dirac fermions can be obtained through the following Hamiltonian, when placed on a honeycomb lattice. Here x, y denotes nearest neighbor sites, and σ is a flavor index. Now specifically taking σ = {↑, ↓}, and H = H 0 + H int with H int as the Hubbard interaction, the model has an S U (2) symmetry that is broken at a critical coupling U = U c when the fermions transition from a massless Dirac fermion phase to an antiferromagnetic phase. Its critical behavior belongs to the Gross-Neveu Heisenberg universality class, corresponding to the Heisenberg SU (2) symmetry broken in the effective continuous Gross-Neveu theory. [1] However, it is not clear how to reach this same universality class using a four-fermion Lagrangian lattice theory. When staggered fermions are combined with a four-fermion term, as in [2], the same symmetry is not broken at the critical coupling and the universality class is different from that of the Hubbard model, as evidenced by different critical exponents in the calculation. e-mail: emilie.huffman@duke.edu. Work done in collaboration with Shailesh Chandrasekharan Another related reason to use the Hamiltonian formalism is that sometimes a Hamiltonian is signproblem-free, but the proposed Lagrangian path integral analog is not. For the following Hamiltonian model, where x is a spatial lattice site,d labels the directions such that x, d labels a unique nearest neighbor bond, and η xy is defined as η x,1 = 1, η x,2 = (−1) x 1 , the bilinear terms in the Hamiltonian describe one flavor of four-component Dirac fermions at low energy. We can get one-flavor of massless Dirac fermions using reduced staggered fermions in the Lagrangian approach, but the interaction term will cause a sign problem.
On the other hand, the Hamiltonian version does not have a sign problem, as seen in [3,4]. This model exhibits a second-order phase transition from massless fermions to a massive charge-densitywave (CDW) ordered phase, with a broken Z 2 symmetry. This puts it in the Gross-Neveu Ising universality class.
The model in (3) will be our focus in the demonstration of a new continuous-time Hamiltonian Quantum Monte Carlo algorithm, which we call the Hamiltonian fermion bag algorithm. Since the solution was found, multiple calculations have been done to determine the critical exponents for the model (3), but there have been discrepancies between them, as summarized in Table 1. (Some of these calculations are actually for the model (3) on a honeycomb lattice with η x,d = 1, which has the same critical behavior.) While the continuous time algorithms seem to have produced η-values of around .3, the discrete time MQMC algorithm is significantly different with an η value of .45. One might be tempted to assume the difference was due to discrete versus continuous time, but the discrete MQMC algorithm also used larger lattices than did the three continuous time results. In our application of the Critical Exponents from Quantum Monte Carlo Past Work ν η Critical Coupling (V/t) c Largest lattice (sites) CT-INT (GS) [5,6] .80 (3) .302(7) 1.304 225 MQMC (GS) [7] .77 (3) .45(2) 1.296 576 CT-INT (Finite-T) [8] .74 (4) .275(25) -484 fermion bag algorithm, we manage to go to even larger lattices (our largest lattice size used has 2304 lattice sites) than the other calculations and find an η value that is larger (.54 (6)) than that of the other lattice calculations. We conclude that larger lattices are necessary to get the critical behavior of model (3), and demonstrate that we have a new algorithm capable of doing calculations at larger lattices than before. Our calculations for lattices with 2304 lattice sites are on par with the largest lattices ever used for these spatially two-dimensional Hamiltonian models (2592 sites have been reached for the Hubbard model [1,9]), and additionally we show stable behavior for the algorithm at lattices of 4096 sites.

Defining the Fermion Bags
A review of the application of fermion bags to Lattice Field Theories can be found in [10]. A key idea is to expand out the exponentiated fermionic terms in the path integral sum, and for each term to divide the lattice space into smaller groups of sites, known as fermion bags, that correspond to the configuration. For example, an exponentiated fermion-interaction term can be expanded as whereψ x and ψ x are Grassmann fields, and n xy represents a specific configuration of n xy = ±1 values. In this way there would be no more fermion interaction terms in the exponent.
The fermion bags could then be defined for each term in the expansion. For example, suppose we had a term in the expansion for some model that was of the form [10] For this particular term there would be four groups of sites that could be treated independently: (x 1 ), (x 2 ), (x 3 , x 4 ), and all the remaining sites. Grassmann integration could be performed over each of the four groups independently. The fermion bag idea has two main potential advantages: 1. sometimes sign problems in these alternative expansions are solved, and 2. these groupings of fermions into bags can be helpful in terms of allowing for fast algorithmic updates [11]. We now turn our attention back to the Hamiltonian formalism, and develop a Hamiltonian fermion bag algorithm based on similar principles. We assume a Hamiltonian that may be written in as a sum of terms of the following form where x is a spatial lattice site, andd labels the directions such that x, d labels a unique nearest neighbor bond. The operators c a x † and c a x are fermionic creation and annihilation operators at the site x with a flavor a = 1, 2.., N f . The couplings of the model are defined through the real constants δ x,d > 0 and α x,d . To get model (3) from this model, up to a constant, we can set ω [12]. Terms of the form (6) contain two key properties that are essential to this algorithm. First, they are local, since each only pertains to two site bond. Second, they are exponentiated fermionic bilinears, which will be important for the application of the Blankenbecler-Scalapino-Sugar (BSS) formula, given later on.
Assuming that H = s,d H x,d we use an expansion similar to the Stochastic Series Expansion (SSE). Here the partition function is expanded as  However, once the temperature is low enough (and thus β is large enough), we should expect to have so many bonds that all the sites end up entangled into a single bag, leaving us with no ability to take advantage of small bags. We can circumvent this problem by dividing the low temperature space into several high temperature spaces. Since the space-time density of bonds is a physical quantity related to the energy density of the system [12], we expect a fixed density of bonds for a particular set of model parameters. At high temperatures we will have fewer bonds and many small fermion bags, and as mentioned before when the temperature is lowered fermion bags will begin to merge to form a single large fermion bag. This suggests that at some optimal temperature the fermion bags may efficiently break up the system into smaller regions that do not depend on the system size. Even at low temperatures, we may be able to divide the imaginary time axis into many time-slices and update a single time-slice efficiently. This is illustrated in the right image of Figure 1, where the imaginary time extent is divided into four-time slices and in the shaded time-slice there are eight fermion bags, instead of the four shown in the middle image.
A nice property is that the sizes of these fermion bags appear independent of the lattice size. In Figure 2 we took model 3 near its critical point, and taking β = 4.0 we divided the imaginary time direction into 16 time-slices and studied the fermion bag size as a function of the lattice size. For equilibrated configurations of L = 48, 64 and 100, the average maximum fermion bag size within a time slice was about 30-independent of L (see Figure 2). Further tests suggests that the optimal temperature is roughly 0.25. Since bond insertions in different fermion bags commute with each other, we can efficiently update fermion bags in space-time blocks (shown as a box in the shaded time slice in the right image of Figure 1) involving 30 to 60 spatial sites within each time slice.  Figure 2. The maximum fermion bag size in a timeslice of .25 as a function of lattice size. We see that this size remains constant even for large lattices so long as the timeslice remains the same size.

Mathematical Details
To study the critical behavior of a model using the expansion (7), we need a mathematical expression for Here 0 ≤ t 0 ≤ β is a time where the operator C n is introduced. In the partition function sector C 0 = I (the identity operator) and in the observable sector C 1 =Ô, whereÔ is operator for the order parameter we would like to measure. The operator O must be written in terms of exponentiated bilinear operators. Because every factor in (8) is an exponentiated bilinear factor, we can apply the BSS formula to (8) and get where B x,d and O n are N × N matrices, where N is the number of sites, and now the matrices have a commutation property analogous to that of the H x,d operators in that B x,d , B x ,d = 0 if the corresponding bonds share no sites in common. For model (3) in particular, B x i ,d i is the identity matrix except in a 2 × 2 block labeled by the rows and columns of the sites that touch the bond x i , d i . Within this block, B x i ,d i takes the form Also, the matrix O 0 is the identity, whereas O 1 comes from the order parameter operator C = (−1) L/2 n (0,0) − 1/2 n (L/2,0) − 1/2 and is given by Now we have what we need to make fast updates according to the fermion bags. As mentioned earlier, we divide the configuration space into smaller timeslices-for model (3) we used timeslices of width 0.25, with t 0 chosen to be at the beginning of the first time slice. We then update bonds within each time-slice sequentially. During the update of a time-slice we define two N × N matrices: the background matrix M B (which is a product of all of the B x,d matrices outside the selected time-slice and O n ), and the time-slice matrix M T , which is the product of all the B x,d matrices within the timeslice being updated. The right image of Figure 1 shows what contributes to M B and M T . When the configuration of bonds within the time-slice is changed then only M T changes to M T . The ratio R is given by where we have defined two new N × N matrices: Since the bond matrices B x,d in different fermion bags commute, it is easy to verify that ∆ is non-zero only within a block which contains spatial sites connected to fermion bags that change. If we randomly choose a spatial block containing about 30 − 60 sites and focus on updating the bonds only within that block, during such a block-update the size of the matrix ∆ cannot be greater than the sum of the sites in the fermion bags that touch the sites within the block. We refer to this set of sites, which can be larger than the block size, as a super-bag and denote its size as s. Since ∆ is non-zero only in an s × s block, it is easy to show that the computation of R (the ratio of the weight of the current configuration with that of the background configuration that existed at the time when the block update began) using (11), reduces to the computation of the determinant of an s × s matrix.
Since G B and M T are fixed matrices during the entire block-update, they can be computed and stored. All proposals to update the current configuration within the block reduce to taking the determinant of an s × s matrix, independent of system size.
Beyond this main savings, there is one other aspect where the fermion bags allow us to save time, and that is for the updates to G B . First note that given two matrices, M 1 and M 2 , we have the identity  d 1 to B x f ,d f matrices that belong to a specific fermion bag. Thus the non-identity block of each G f is a matrix with f rows and f columns corresponding to the fermion bag sites. We can then combine the G f matrices into a G T matrix, which has distinct blocks according to the fermion bags. Thus, while G B is an N × N matrix, we can use the idea of fermion bags along with the identity (12) to reduce the number of O N 3 operations. While this does not reduce the scaling of the algorithm, it does significantly reduce the prefactor. These are the main aspects of the algorithmic updates that are helped by the fermion bags. The full details of the algorithm are given in [13].

Results
The results at the time of Lattice 2017 are summarized by the plots in Figure 3. The left image demonstrates the stability of the algorithm by showing the high-temperature equilibration of 100×100 lattices (10,000 sites) at β = 4.0. Now low temperature equilibrations of 64×64 lattices and 100×100 lattices can be found in [13]. For β = L sweeps on 48×48, 64×64, and 100×100 lattices on a single core are 4 hours, 30 hours, and 30 days respectively.
The middle figure and left figure give our calculations that were done at the couplings V = 1.304t and V = 1.296t. These were the critical couplings found by [6] and [7]. At the critical coupling, the power law behavior of C = a L 1+η , is expected. Our calculations include lattices up to size 48×48 (2304 sites) and fit the data to (13). For V = 1.304t, the power law fit was poor and the η value was only .149(9). However, when we  Data for the couplings V = 1.304t and V = 1.296t, which were the critical couplings found in [6,7]. dropped the larger lattices and only kept up to lattice size L = 20, we could achieve a good fit with η = .31(2), consistent with the previous result of η = .318(8) that went up to L = 20 for specifically model (3). [5,6,8] For V = 1.296t, the power law fit was again poor with an η value of .35 (1). Again dropping the larger lattices and this time keeping lattices up to L = 24 led to a good fit with η = .41(4), consistent with the previous result of η = .43(2) that went up to L = 22 for specifically model (3). [7] Both of these results suggest that larger lattices make a substantial difference in the calculation and therefore are important in determining the critical behavior. In [13] we combine our data from several different couplings (using lattices up to L = 48) and use a combined fit to calculate exponents of η = .54(6), ν = .88 (2), and V c = 1.279(3)t.