Topological sectors, dimer correlations and monomers from the transfer-matrix solution of the dimer model

We solve the classical square-lattice dimer model with periodic boundaries and in the presence of a field $\boldsymbol{t}$ that couples to the (vector) flux, by diagonalizing a modified version of Lieb's transfer matrix. After deriving the torus partition function in the thermodynamic limit, we show how the configuration space divides into 'topological sectors' corresponding to distinct values of the flux. Additionally, we demonstrate in general that expectation values are $\boldsymbol{t}$-independent at leading order, and obtain explicit expressions for dimer occupation numbers, dimer-dimer correlation functions and the monomer distribution function. The last of these is expressed as a Toeplitz determinant, whose asymptotic behavior for large monomer separation is tractable using the Fisher-Hartwig conjecture. Our results reproduce those previously obtained using Pfaffian techniques.


I. INTRODUCTION
The dimer model is a paradigmatic example of a stronglycorrelated system, in which dimers cover the edges of a lattice subject to a close-packing constraint, i.e., each vertex touches exactly one dimer. It was first solved independently by Kasteleyn [1,2] and Temperley and Fisher [3,4] in 1961 using a combinatoric method, in which the partition function is expressed as the Pfaffian of a signed adjacency matrix known as the Kasteleyn matrix. Because the dimer model is exactly solvable, it offers a useful setting for the study of novel phenomena in geometrically frustrated systems [5] or, more specifically, 'Coulomb-phase' physics [6]. In particular, its extensive entropy reflects macroscopic ground-state degeneracy, while the configuration space splits into topological sectors labeled by horizontal and vertical 'flux' components, reflecting topological order [7]. The Pfaffian method can be used to calculate partial partition functions for these sectors, as demonstrated by Boutillier and de Tilière [8].
Moreover, a dimer can be replaced by a pair of monomers, which can be separated by subsequent dimer updates and thus play the role of fractionalized excitations. Fisher and Stephenson's Pfaffian calculation of the monomer distribution function in 1963 [9] implies that, due to the entropy of the background dimer configuration, the monomers interact through an effective Coulomb potential, which is logarithmic in two dimensions. They have also shown that dimer-dimer correlations are long-range with algebraic, rather than exponential, dependence on separation. This is despite the absence of long-range order, and instead a consequence of the close-packing constraint.
Perhaps a more elegant solution of the dimer model is Lieb's transfer-matrix method [10], analogous to the well-known solution of the Ising model by Schultz et al. [11], which maps the problem to free fermions. In this approach, the partition function is expressed in terms of a transfer matrix, which, given a configuration on a row of vertical bonds, generates all dimer configurations compatible with the close-packing constraint on the subsequent row of horizontal and vertical bonds. This can be expressed in terms of spin- 1 2 operators and mapped to fermions through a Jordan-Wigner transformation.
This method has been used in the literature to derive the par-tition function [10] and to determine its vertical-flux decomposition [12,13]. In this work, we show how Lieb's transfer matrix can be modified in order to calculate the full flux-sector decomposition. We also provide a general framework for the calculation of expectation values and explicitly calculate dimer occupation numbers, dimer-dimer correlation functions and the monomer distribution function. For the last of these, we show how the asymptotic dependence for large monomer separation, which was deduced by numerical means in Ref. [9], can be evaluated exactly by applying the Fisher-Hartwig conjecture [14].

Outline
In Sec. II we define the model before showing how it can be formulated in terms of a transfer matrix in Sec. III. We then diagonalize the two-row transfer matrix in Sec. IV, whose spectrum is used to calculate the partition function, including its flux-sector decomposition, in Sec. V, and various expectation values in Sec. VI. We conclude in Sec. VII.

II. MODEL
We consider the close-packed dimer model on an × square lattice with periodic boundary conditions (PBCs), assuming both , even. In the following, we define the flux along with the weights that appear in the partition function.
Denoting by , the dimer occupation number (equal to zero or one) on the bond joining sites and + , with a unit vector in direction ∈ { , }, the flux is given by where = (−1) + = ±1 depending on the sublattice. Due to the close-packing constraint, this is equivalent to the sum of , on links crossing a surface normal to (to see this, one usually defines an effective 'magnetic field' [6,15]). The latter definition highlights that Φ is integer valued, and can only be changed by shifting dimers around a loop encircling the whole system [16]. The flux thus plays the role of a topological invariant. To each configuration, we assign weight e i ⋅ . In the first factor > 0 and are the 'activity' and number of horizontal dimers, respectively. (The total number of dimers + = 1 2 is fixed, so the activity of vertical dimers is set to unity without loss of generality.) Hence, for ≠ 1, the model is anisotropic, with horizontal (vertical) dimers favored for > 1 ( < 1). In the second factor is a field, with components ∈ (− , ], that couples to the flux . An example configuration is shown in Fig. 1.
The partition function is where ℭ 0 denotes the set of all close-packed dimer configurations, and can be thought of as a moment-generating function for Φ . Similarly, expectation values of a function of the dimer occupation numbers , are given by

III. TRANSFER MATRIX
We construct the partition function, Eq. (2), by modifying Lieb's transfer matrix [10] to include the Φ weighting (the Φ weighting can be included without modifying the transfer matrix).
We first define a vector space whose basis vectors |̄ ⟩ correspond to all possible configurations̄ of the dimer occupation numbers on a single row of vertical bonds. As illustrated in Fig. 2, the transfer matrix is defined so that wherē ′ is the configuration on the subsequent row of vertical bonds and ℭ(̄ ,̄ ′ ) is the (possibly empty) set of configura-tions̄ of the intermediate row of horizontal bonds that are compatible with̄ and̄ ′ . The weight function is chosen to give the correct weights for and Φ in the partition function of Eq. (2). On even rows, where = (−1) in Eq. (1), it is given by where while on odd rows is defined in the same way, but with replaced by * . (Here,̄ , denotes the occupation number of the bond between sites = and + 1 in the configuration of the horizontal bonds.) It is convenient to split the action of into two steps:  Fig. 2). The effect on̄ ′ is that an adjacent pair of dimers is removed.
In order to reproduce the weight function , a horizontal dimer on the bond between sites and + 1 in step 2 comes with a factor ( * ) on even (odd) rows. An explicit operator expression for the transfer matrix is obtained by representing occupied and empty vertical bonds by spin up |↑⟩ and down |↓⟩ states, respectively [i.e., eigenstates of , where = ( , , ) are the Pauli matrices]. The above steps are easy to formulate in the spin language. As shown in Fig. 2, step 1 is equivalent to flipping all spins, which is achieved by the operator since ± = 1 2 ( ± ) satisfy + |↓⟩ = |↑⟩ and − |↑⟩ = |↓⟩. In step 2, pairs of neighboring up spins are flipped, so the operator effectively generates a horizontal dimer between sites and + 1, with the correct weight on even rows. Because ( − ) 2 = 0, the operator ( !) −1 ∑ =1 , generates horizontal dimers (PBCs require − +1 = − 1 ), and hence generates an arbitrary number of horizontal dimers. To obtain the correct weights on odd rows, one should instead use the operator * 3 . It is therefore necessary to define two transfer matrices, on even rows and * = † on odd rows. 1 (Note that = because + = − .) We also define the two-row transfer which is manifestly Hermitian. The Φ weighting is included in the transfer-matrix formalism as follows: The operator for the dimer occupation number on a vertical bond is simply since spin up (down) corresponds to an occupied (empty) bond. In terms of this, the vertical flux component on even rows is [see Eq.
construct mutual eigenstates of the two-row transfer matrix and Φ . The partition function, Eq. (2), is then given by (the trace arises due to PBCs in the vertical direction). Similarly, the operator analog of Eq. (3), in the case of the correlation function between observables and ′ in rows 1 ≤ ≤ ′ ≤ , is given by where ( ) = ( ) −1 ( ) and To compute expectation values of dimer observables, it is necessary to find operators that correspond to these quantities. While a suitable operator for the dimer occupation number on vertical bonds has already been defined in Eq. (12), no such operator exactly represents the dimer occupation number on horizontal bonds, since the vector space on which the transfer matrix acts contains only dimer configurations on vertical bonds.
One can nonetheless calculate expectation values involving horizontal dimers using an appropriately constructed operator. From Eqs. (4) and (5), one finds whereas Eqs. (8)- (10) give the operator identity since [ , , ′ , ] = 0. Comparing the right-hand sides, we therefore interpret , as the operator corresponding to the horizontal dimer occupation number̄ , on an even row, but only when appearing in the combination 3 , . Similarly, * , acts as the horizontal dimer occupation number on an odd row in the combination * , † . Setting equal to , ( * , ) on even (odd) rows in Eq. (15) gives the correct combination , ( * , † ) in ( ), allowing one to calculate expectation values involving the horizontal dimer number.

IV. DIAGONALIZATION OF THE TWO-ROW TRANSFER MATRIX
To calculate Eq. (14) it is sufficient to diagonalize the tworow transfer matrix . We do so in this section through a series of transformations.
We map between spins and spinless fermions using the Jordan-Wigner transformation [18][19][20] which identifies spin up and down with filled and empty fermion orbitals, respectively, while preserving the usual (anti)commutation relations In terms of fermions, Eqs. (8) and (12) become while the condition − +1 = − 1 is equivalent to 3 This means that, for example, 2 , does not give the square of the horizontal dimer number; in fact 2 , = 0, whereas̄ 2 , =̄ , .
We now define projectors into the subspaces with even ( = 0) or odd ( = 1) Φ , which satisfy ∑ Π = 1 and (−1) Φ Π = (−1) Π . Then, since (−1) Φ commutes with any quadratic form in fermions, and the fermion operator +1 depends implicitly on through the boundary condition More generally, for any operator containing +1 of Eq. (26), we define an operator that instead only contains +1 of Eq. (32) (and thus depends on ), such that the action of both operators on a state with Φ parity yields the same result, i.e., = ∑ Π . (For operators that do not contain +1 , such as Φ , one has = .) For later reference (see Sec. VI) we note that, after the Jordan-Wigner transformation, the single-row transfer matrix is given by where the operators in the product should be ordered from right to left. We now make a Fourier expansion with and which ensure the correct boundary condition on +1 in Eq. (32) [10]. 4 The fermions obey standard anticommu- 4 As an alternative to the approach in Sec. III, one could instead implement the Φ weighting using = and twisted boundary conditions − +1 = e i − 1 in place of Eq. (6) and − +1 = − 1 , respectively [see Eq.
(1) and text thereafter]. However, a Fourier expansion of the new set of fermions is no longer useful because of the absence of translation symmetry [21,22]. Instead, one would have to perform the gauge transformatioñ = e −i (−1) ∕ back to fermions, before proceeding as in the main text. tation relations, as follows from Eq. (23). Using the result valid for both and ′ in either 0 or 1 , the operator appearing in the exponential of Eq. (31) can be written as Restricting the sum to 0 ≤ ≤ 2 , this becomes where the quadratic form Here, [its Hermitian conjugate means the row vector † = The additional factor of 1 2 for ∈ {0, 2 } prevents double counting of these terms in Eq. (39), and ensures the commutation relation is valid for all 0 ≤ ≤ 2 . 5 Since † ( ) = ( † ), and all quadratic forms in Eq. (39) commute by Eq. (44), the two-row transfer matrix, Eq. (31), is given by which can be reordered as the following product of commuting terms: To proceed, we map to the corresponding one-dimensional quantum Hamiltonian  through Then, by Eq. (30), we have where since the projectors satisfy [Π , ′ ] = 0 and Π Π ′ = Π ′ . After inserting Eq. (48), this implies The Baker-Campbell-Hausdorff formula [23] states that the logarithm in Eq. (52) can be expressed in terms of nested commutators of ( ) and ( † ). Using Eq. (44), these can be expressed in terms of nested commutators of and † , giving 5 For with ⊗ denoting the Kronecker product, Eq. (44) is only true if satisfies the condition = − (or the same for ). However, it is always possible to symmetrize to meet this condition: Using ( † ) = , , one can show where ′ = 1 2 ( − ) is a matrix that satisfies the condition. The matrix ( ) in Eq. (42) has been constructed in this way.
The problem is thus reduced to diagonalization of the 4 × 4 matrix e e † for each .
In order to solve the eigenvalue problem we expand e as a power series and use 2 = 0 to obtain e e † = + + † + † .
After substituting Eq. (42) and writing = ( 1 2 ) , Eq. (54) reduces to a pair of simultaneous equations which, on rearrangement, read The latter is a 2×2 eigenvalue problem, which is easily solved. The result implies and is a unitary matrix whose columns are the eigenvectors of e e † .
By inserting Eq. (58) into Eq. (53), we obtain the freefermion Hamiltonian with dispersion where the and fermions are related by the Bogoliubov transformation for 0 ≤ ≤ ∕2. Both sets of fermions obey standard anticommutation relations. The transformation of Eq. (62) may be expressed as a single transformation valid for all : with Combining Eqs. (34) and (63), the transformation relating the and fermions is This makes it clear that the annihilation operator removes a fermion (or equivalently, removes a vertical dimer) on even sites or adds one on odd sites. According to Eq. (27), it therefore reduces Φ by one. We now construct the spectrum of . As discussed in Sec. III, one can find simultaneous eigenstates of  and Φ . After substituting Eq. (65) into Eq. (27), the latter is given by in terms of fermions, which counts the number of occupied states relative to half filling [the number of available -states is by Eqs. (35) and (36)]. 6 The occupation-number states of the fermions with ∈ form a complete set of mutual eigenstates of  and Φ . From Eq. (50), the complete set of eigenstates of  is given by the union of all eigenstates of  0 that have even Φ eigenvalue and all eigenstates of  1 that have odd Φ eigenvalue. We will denote |Φ ⟩ as the th excited eigenstate of  with vertical flux Φ , and (Φ ) as its eigenenergy. The spectrum of the two-row transfer matrix follows from that of  through Eq. (49): |Φ ⟩ is also an eigenstate of , but with eigenvalue e −2 (Φ ) .
As illustrated in Fig. 3 (top-left panel), the ground-state is half filled and thus denoted by |0⟩ 0 . Formally, it is defined by where ∈ 0 , and has energy which arise when extending the integration bounds, as well as the remaining terms in Eq. (70), can be calculated using the Taylor expansion where F = is the Fermi velocity. The final result is and a similar calculation for the lowest-energy state in the Φ = 1 sector gives Note that the dependence of 0 (0) is the standard result [24] for the ( −1 ) correction to the ground-state energy of fermions with a twist in their boundary conditions (see Footnote 4). Relating this to the (effective) central charge [25], we have and so In particular, for = 0, this gives the expected result of = 1 for a theory containing a single complex fermion. We note, however, that the value of the central charge in the dimer model is controversial, with other arguments suggesting instead = −2 when = 0 [13].

V. PARTITION FUNCTION
In this section, we write down the partition function ( ) using Eq. (14) and eigenvalues of the two-row transfer matrix, before taking the thermodynamic limit.
By Eqs. (30) and (51), one can split ( ) into contributions from each parity sector, giving The projector Π can be expanded using Eqs. (28) and (67) as and hence with̃ The partition function, Eq. (79), can therefore be written as where , = Tr e −  , . Because the trace of an operator is equivalent to the sum of its eigenvalues, one has which reduces to Lieb's partition function for = [see Ref. [10], Eq. (3.14)].
We now take the thermodynamic limit, retaining leadingorder corrections to the free-energy density. To do so for 0,± , we factor out ±e − ( − ∕ ) e i for all terms in the product with < 0 and restrict the product to 0 < ≤ ∕2, which gives where = (2 − 1) by Eq. (35). In the limit , → ∞, we can replace ( ± ∕ ) by its leading-order dependence ( ± ∕ ) [see Eq. (74)], since the next-order terms will eventually be of order ∕ 3 . Hence, Eq. (86) becomes where = e e i , = e −2 and = ∕ . This can be expressed in terms of Jacobi theta functions using the first equality of Eqs. (A4) and (A5): where ( ) is the Dedekind eta function defined in Eq. (A1), and the same for 0,− but with 3 → 4 . An analogous calculation for 1,± yields which is consistent with Eq. (8.41) of Ref. [12] when = 0. When = , 1 (1| ) = 0 and the partition function is in agreement with Ref. [26]. The first factor in Eq. (90) grows exponentially with system volume, and represents the weight of dimer configurations in the bulk, i.e., it specifies the bulk free-energy density [12] As one might expect, bulk does not depend on the choice of boundary conditions, although we note that this is not true in the case of the honeycomb lattice [27]. The remaining terms in ( ), which give leading finite-size corrections to the free-energy density, are boundary dependent and, in the case of PBCs, encode information about topological flux sectors (see subsection below). Previously, these terms have also been evaluated (for = ) with closed [26] and cylindrical [28] boundaries, as well as embeddings on the Möbius strip and Klein bottle [29]. In general, one obtains terms in the free energy proportional to the edge of the system [e.g., 2( + ) for closed boundaries] and of order ∕ . However, with PBCs (i.e., a torus) the edge is zero and we only observe the latter.
Using the modular identities given in Appendix A, one can confirm that the partition function ( ) behaves as expected under 90 • rotations, in spite of the asymmetry between and in the transfer-matrix method. Such a rotation takes where ′ = e − ∕ e i and ′ = e −2 ∕ correspond to and under rotation. Since 2 (i∕ ) = 2 (i ) − i 2 log [30, Sec. 25.12], the remaining (bulk) factor in Eq. (90) is replaced by giving the expected transformation of ( ).

Flux sectors
We now show how the partition function, Eq. (91), divides into topological sectors labeled by the flux. By construction, ( ) is periodic in (with period 2 ), so can be expressed as a Fourier series Comparison of Eqs. (2) and (96) implies where the set ℭ 0 ( ) contains all close-packed dimer configurations with flux . In other words, the Fourier coefficient can be interpreted as the partial partition function, or total weight, of flux sector .
To calculatẽ , we use the second equality of Eqs. (A2)-(A5) to rewrite Eq. (90) as [8] ( ) = e − bulk ∑ ∈ℤ e − ( −2 ) 2 ∕2 ∑ ∈ℤ e i e − 2 ∕2 2 ( ) (98) (the periodicity in is now apparent). The sum over can be written in the same form as the sum over through the Poisson summation formula, giving This result has previously been obtained for the honeycomblattice dimer model using Pfaffian methods [8], while Ref. [12] has used the transfer matrix to calculate the partial partition function of flux sector Φ , equivalent to ∑ Φ ̃ [see their Eqs. (8.19) and (8.36)]. Knowledge of̃ can be used to calculate flux moments. The probability of flux is given by which implies that Φ and Φ are independent variables. This form is known from effective field theories [31,32]. The mean flux vanishes by symmetry, while the mean-square flux is given by and the same for Φ but with → 1∕ .

VI. EXPECTATION VALUES
In this section, we compute various expectation values in the thermodynamic limit, using the spectrum of the two-row transfer matrix.

A. Two-point correlation functions of fermions
For an operator given by a product of fermions, the corresponding time-evolved operator ( ) can also be expressed as a product of ( ) , with the same for each. For example, when = , one has Here, ( ) is defined by extending the definition in Eq. (106) to , even though it does not conserve parity and so does not obey Eq. (105). An expectation value ⟨ ′ ( ′ ) ( )⟩ , can then be expressed in terms of a product of an even number of ( ) operators. Because this is a time-ordered product and , is a free-fermion Hamiltonian, Wick's theorem [33] applies, which allows us to write ⟨ ′ ( ′ ) ( )⟩ , as a sum over products of two-point ( ) correlators in each ( , ) sector. [We similarly extend the definition Eq. (109) to include = , even though Eq. (108) is not valid in this case.] We calculate these two-point correlators in this section.
To do so, we first use Eqs. (16) and (106) for all . Since , , defined in Eq. (82), is a free-fermion Hamiltonian with dispersioñ , and Eq. (109) describes a thermal distribution with effective temperature 1∕ , the two-point cor-relation functions of the fermions are given by where F ( ) = (e +1) −1 is the Fermi-Dirac distribution function.
Hence, denoting = ( , ), the ( ) correlators are These results are exact, with the correct (anti)periodicity in the horizontal direction, and could be used to calculate expectation values for finite system sizes as a function of flux sector. Instead, we take the thermodynamic limit , → ∞, keeping the ratio ∕ and the separation | | finite. In this limit, F ( ) can be replaced by a step function (− Re ) and the discrete values become continuous, giving Some values of these integrals for small | | are shown in Table I, expressed in terms of the quantities which satisfy + = 1 2 . For large | |, the asymptotic behavior is obtained by integrating by parts repeatedly, treating the cases ≫ 1 [where Eq. (74) can be used] and of order unity separately.
These expressions are independent of and , i.e., all four ( , ) sectors make equal contributions in the thermodynamic limit. Hence, Eq. (108) is redundant to this order, and we simply have ⟨ ′ ( ′ ) ( )⟩ = ⟨ ′ ( ′ ) ( )⟩ 0,+ for operators that are products of an even number of fermions. We therefore drop the ( , ) indices from now on. Furthermore, they are independent of , whose leadingorder dependence is ( −1 , −1 ). This implies that expectation values are the same in any fixed flux sector in the thermodynamic limit (but note that that we have taken , → ∞, so this does not apply for ∼ , ). To see this we rewrite Eq.
where and integrating over , one finds that ⟨ ⟩ = ⟨ ⟩ when the latter is independent of .
In subsequent sections we use Eqs. (119)-(122) to calculate various observables in the dimer model in the thermodynamic limit. We expect our results to reproduce those of Ref. [9] in this limit, since the choice of boundary conditions (PBCs versus closed) becomes irrelevant. We also note that asymptotic behavior of correlation functions can be predicted using effective field theories, although the results depend on phenomenological parameters known as the stiffnesses [32].

C. Dimer-dimer correlation functions
Due to the close-packing constraint, the occupation of a given bond by a dimer is influenced by dimers far away. Hence, dimer-dimer correlations are non-trivial even in the absence of interactions. In this section, we show how they can be calculated by extending the above discussion to two-point correlators of , and , .
The results in this section are in agreement with Sec. 7 of Ref. [9].

D. Monomer distribution function
Finally, we characterize the (entropic) interaction between a pair of inserted test monomers by calculating the monomer distribution function where the set ℭ( + , − ) contains all configurations with monomers at sites ± . For simplicity, we consider the case of two monomers on the same row, though the formalism can be extended to the general case. Because − inserts a monomer on site , in the transfermatrix formalism one has which becomes In terms of these, Eq. (152) is a sum of four 2 -point correlators, each of which can be expressed as a sum of products of two-point correlators through Wick's theorem [33]. Then, by Eqs. (155) and (156), the two correlators containing an unequal number of and vanish, while the remaining two are where denotes the symmetric group of order , and (−1) −1 ( , ). Inserting Eq. (157) and using the relation 7 In the case of two monomers on different rows, the operator on each row has an odd number of operators and so does not commute with (−1) Φ . To treat this case, we would not be able to use Eq. (108) and would instead require the analogous expression for anticommuting with (−1) Φ . which can be expressed as a Toeplitz determinant where is an × matrix with elements ( ) , ′ = −2Γ(1− ( − ′ ), 0).
From Table I, the first two non-zero values [cf. Eqs. (11.1) and (11.3) of Ref. [9]], where, up to a factor of , the former is equivalent to the occupation probability of a horizontal bond as calculated in Sec. VI B.
A consistent result was found by Hartwig [39] for the case of monomers separated along a diagonal (i.e., = ) using the Pfaffian method. Note that the algebraic dependence on , stemming mathematically from the discontinuity in , contrasts with the exponential behavior on the triangular lattice [35,36]. As noted by Au-Yang and Perk [40], the decrease with −1∕2 can be understood by relating the dimer model to two uncoupled Ising models at the critical point.

VII. CONCLUSIONS
We have expressed Lieb's transfer matrix for the classical square-lattice dimer model in terms of a free-fermion Hamiltonian, and used its spectrum to rederive some useful results. Although these can equally be derived using Pfaffian techniques, the second quantized approach presented in this paper is perhaps more elegant.
Specifically, our results include the torus partition function which, by including a field , can be interpreted as a moment-generating function of the flux. We have also shown how expectation values can be expressed in terms of the fermionic operators, and evaluated dimer occupation numbers, dimer-dimer correlation functions and the monomer distribution function in the thermodynamic limit, all of which are independent of flux sector for not-too-large flux. Finally, we have derived a new result, namely the asymptotic behavior of the monomer distribution function for large monomer separation along the same row.
The results in this paper are also relevant to the corresponding quantum dimer model at its Rokhsar-Kivelson point [41], while the transfer-matrix method can be extended to other twodimensional lattices. Indeed, the straightforward generalization of Lieb's transfer matrix to the (bipartite) honeycomb and square-octagon lattices, which can both be viewed as a square lattice with certain horizontal bonds removed [i.e., certain terms omitted from the sum in 3 ; see Eq. (9)], has already been demonstrated in Ref. [42].
One advantage of the transfer-matrix method is that dimerdimer interactions can be easily included in the operator formalism, in terms of products of the dimer occupation numbers , and , . For example, on a row of vertical bonds, the operator ∑ , +1, describes interactions between parallel pairs of nearest-neighbor dimers, as studied in Refs. [31,43]. This is a four-fermion interaction, which is non-integrable [31] but could be included perturbatively using standard diagrammatic perturbation theory.
Furthermore, as we will show in a forthcoming publication [34], the well-known height field theory [44,45] of the two-dimensional classical dimer model can be rigorously derived from the fermionic Hamiltonian of Eq. (60). This can be achieved by taking a long-wavelength limit and using the technique of bosonization [46] to express the theory in terms of a single free bosonic field. Interaction operators included perturbatively in this context manifest themselves through renormalization of the 'stiffness' as well as the introduction of (cosine) potential terms consistent with symmetry requirements [34].