# Ising nematic fluid phase of hard-core dimers on the square lattice

###### Abstract

We present a model of classical hard-core dimers on the square lattice that contains an Ising nematic phase in its phase diagram. We consider a model with an attractive interaction for parallel dimers on a given plaquette of the square lattice and an attractive interaction for neighboring parallel dimers on the same row (viz column) of the lattice. By extensive Monte carlo simulations we find that with a finite density of holes the phase diagram has, with rising temperatures, a columnar crystalline phase, an Ising nematic liquid phase and a disordered fluid phase, separated by Ising continuous phase transitions. We present strong evidence for the Ising universality class of both transitions. The Ising nematic phase can be interpreted as either an intermediate classical thermodynamic phase (possibly of a strongly correlated antiferromagnet) or as a phase of a 2D quantum dimer model using the Rokhsar-Kivelson construction of exactly solvable quantum Hamiltonians.

###### pacs:

## I Introduction

Hard-core dimer models have been useful in identifying and handling basic aspects of the configurational frustration in quantum spin systems with non-magnetic ground states such as RVB wavefunctions and valence bond singlet crystals. Hard-core dimer models have been used extensively in the context of quantum condensed matter model constructions, such as the quantum dimer models by Rokhsar and Kivelson,Rokhsar and Kivelson (1988) to characterize topological fluid quantum phases,Moessner and Sondhi (2001a) strongly correlated superconducting phases,Fradkin and Kivelson (1990) as well as crystalline phases and stripe phasesPapanikolaou et al. (2007) (see, e.g Ref.[Fradkin, 2013]).

Such modeling led to useful insights and basic intuitive understanding of aspects of quantum frustration.Moessner et al. (2000); Moessner and Sondhi (2001b) Similar models can be used to understand the high-temperature behavior of such strongly correlated quantum antiferromagnets, assuming that the valence bond spin gap has an onset at much higher energies.

Ising nematic liquid phases and more generally, electronic liquid crystalline phases,Kivelson et al. (1998) have been observed in several physical systems, notably in the two-dimensional electron gas (2DEG) in large magnetic fields, in the bilayer ruthenate SrRuO in magnetic fields, in the iron-based superconductors and most interestingly, and in the cuprate superconductors.Fradkin et al. (2010); Fradkin (2012); Kivelson et al. (2003) In particular, nematic order has been reported in YBaCuO in much of the pseudogap regime,Daou et al. (2010) and in very underdoped YBaCuO over a large temperature range,Hinkov et al. (2008) as well as in BiSrCaCuO,Lawler et al. (2010) and in the iron-based superconductors LaFeAsO and Ba(FeCo)As.de la Cruz et al. (2008); Chu et al. (2010); Chuang et al. (2010); Fang et al. (2008); Xu et al. (2008)

In this paper, we study the emergence of an Ising nematic dimer phase in a classical dimer model, using extensive numerical classical Monte Carlo simulations and analytical arguments. We consider dimer models with two types of interactions: a) an attractive interaction with strength for parallel dimers on the same plaquette of the square lattice (which favors columnar ordered states), and b) an attractive interaction of strength for dimers on next-nearest-neighbor bonds on the same row (viz column) of the square lattice (which favors all tilted columnar states). We investigated the phase diagram of this model, shown schematically in Fig. 2 A, B and C, at full dimer coverage and at a finite density of holes as a function of temperature and for varying dimer interactions and .

In the fully packed regime, 2D dimer models can describe ordered phases, columnar phases with varying degree of translation symmetry breaking,Fradkin et al. (2004); Vishwanath et al. (2004); Alet et al. (2006); Castelnovo et al. (2007); Papanikolaou et al. (2007) as well as quantum critical systemsRokhsar and Kivelson (1988); Ardonne et al. (2004) and topological phases.Moessner and Sondhi (2001a); Fendley et al. (2002) The nematic phase of this system is characterized by spatially disordered dimer configurations that globally break rotational invariance, as shown in Fig. 1, i.e. long range orientational order that breaks spontaneously the point group symmetry of the square lattice but leaves translation symmetry intact. We observe that such a phase cannot be stable on the square lattice against columnar energy perturbations in the absence of holes (the undoped system). Only in the existence of mobile holes, we find that a nematic phase exists for a finite temperature range (see Fig. 2B). Moreover, we find that in the presence of generic perturbations such as a small attractive columnar interaction, this liquid phase only occurs in an intermediate temperature range, comprised between a high temperature isotropic liquid and a low temperature columnar crystal, separated by corresponding (classical) 2D Ising phase transitions.

These phases can be interpreted as occurring in a classical dimer model at various coverings. In the context of quantum antiferromagnets, a bond occupied by a dimer means that it is occupied by a valence bond singlet of two neighboring spinsAnderson (1987); Kivelson et al. (1987) or as a label of an unsatisfied bond in frustrated quantum magnets.Moessner and Sondhi (2001b); Bergman et al. (2007) In doped antiferromagnets with a spin gap, in addition to valence bonds one has to include the charge degrees of freedom in the form of holes. A simple model of such systems is the doped quantum dimer model.Rokhsar and Kivelson (1988); Fradkin and Kivelson (1990) Dimer models of the type we discuss here can also serve to represent, qualitatively, forms of orbital order of the type used in models of the iron superconductors,Lv and Phillips (2011) and in the analysis of the STM data in BiSrCaCuO,Lawler et al. (2010) as well as for the analysis on valley order in certain 2DEGs in magnetic fields.Abanin et al. (2010); Kumar et al. (2013) Although the dimer models do not represent many aspects of the microscopic physics of these systems faithfully, its symmetries are the same as the point groups of the more microscopic models and, hence, should be useful to describe many aspects of these systems, particularly their phase transitions.

Here we show that a classical dimer model with a finite density of holes can lead, with suitable interactions, to a nematic phase in dimer models arising from the melting of a columnar state. However, it has been unknown how to stabilize such a nematic liquid in classical or/and quantum dimer models by using isotropic Hamiltonians with local interactions. By using the RK-construction,Rokhsar and Kivelson (1988); Ardonne et al. (2004); Castelnovo et al. (2004); Papanikolaou et al. (2007) it is straightforward to find an exactly solvable quantum dimer models with a stable quantum nematic phase whose ground-state wave function has local weights of the same form as the Gibbs weights of the classical dimer model we study here, given the configurations are assumed to be orthonormal. Namely,

(1) |

where

(2) |

is the norm of the wave function and also the partition function of a classical system with Gibbs weight in the same dimension.

We fully characterize the behavior of the corresponding classical model through analytic arguments and extensive simulations. First, we show that the key ingredient for stabilizing nematic dimer liquids is a strong same-row or same-column dimer-aligning, spatially isotropic, interaction. In a fully packed lattice this interaction leads to crystalline phases.Alet et al. (2005b); Charrier et al. (2008); Charrier and Alet (2010); Alet et al. (2006); Papanikolaou et al. (2007) Here we will show that at finite hole density it stabilizes an Ising nematic fluid phase. Further, we demonstrate numerically that all relevant transitions at finite hole density are of Ising nature and that the nematic phase, even though it breaks explicitly rotational symmetry, by having most dimers point to the same direction, also has disordered liquid-like correlations along the orientation direction. However, in the close-packed case, and when only the dimer-aligning interaction is present, it is unclear whether the transition is in the Ising or the KT universality class. Theoretical arguments predict the KT transition but the simulations give, as we will discuss, inconsistent results possibly due to finite-size effects.

This work is organized as follows. In Sec. II we present the generalized dimer model and discuss its symmetry properties. In Sec. III we present and analyze the results of our Monte Carlo simulations which lead to the phase diagrams of Fig. 2. In Sec. IV we present results form a mean-field theory based on a formulation of the generalized dimer model in a Grassmann variables formulation. Our conclusions are presented in Sec. V.

## Ii The Model

The classical dimer model is on a square lattice of linear dimension (total number of sites ) covered by hard-core dimers and holes with the hole density being The configuration space consists of all the possible dimer coverings satisfying the hard-core constraint with a given number of holes. The energy of a dimer configuration is

(3) | ||||

where is the dimer occupation number of the link of the square lattice with endpoints at the two nearest neighbor sites and (with ). The first term of the total energy of Eq. (3) is proportional to the number of parallel dimers on the same row (or column) on next-nearest-neighbor sites, while the second term is proportional to the number of plaquettes with two parallel dimers, is the inverse temperature of the classical model and we assume that the interaction of parallel dimers on a plaquette is attractive, . Below, we will refer to the first term of Eq.(3) the dimer-aligning interaction (with strength ) and to the second term the plaquette interaction (with strength ).

The model with configurational energy of Eq. (3) is explicitly invariant under translations (by one lattice spacing) in both directions of the square lattice and under point group symmetries. However, in the special case , the configurational energy is also invariant under the independent rigid shifts of the dimer configurations on the links of each separate rows and columns. This discrete “sliding” symmetry formally relates dimer configurations with different tilts.Fradkin et al. (2004) As it stands, this symmetry formally leads to a sub-extensive (on the square lattice) entropy even at which hinders the equilibration of this system for large lattice sizes. A discrete sliding symmetry of this type also arises in the 2D classical Ising model with only four spin interactions.Xu and Moore (2004) However this symmetry can be broken by fixed boundary conditions to select a particular tilted state. This symmetry is also explicitly broken by the plaquette interactions even for arbitrarily small values of . For this reason, it is important to focus on the case .

The partition function for a dimer model with no aligning interaction was studied in detail both in the fully packed caseAlet et al. (2005b, 2006); Papanikolaou et al. (2007) and with a finite density of holes.Papanikolaou et al. (2007) In the fully packed regime it was found that, for , the system is characterized by a line of fixed points parametrized by . At temperatures below a Kosterlitz-Thouless (KT) transition at , the fully-packed dimer model enters a four-fold degenerate columnar phase which has a spontaneously broken translation and rotational symmetry.Alet et al. (2005b); Papanikolaou et al. (2007) For any density of holes, the line of fixed points is replaced by a disordered phase, but the columnar phase survives, unless the density of holes is large enough.Papanikolaou et al. (2007)

The phase diagram in the presence of the aligning interaction as well as the plaquette interaction at various hole densities is complex. Apart from the already described disordered liquid and columnar solid phases, it is expected that a new phase emerges, which only breaks the rotational symmetry but not the translational one. The existence of this Ising nematic phase has been proven for the case on dilute systems by Heilmann and Lieb Heilmann and Lieb (1979). Here, we will explore the emergence of this phase primarily through numerical simulations, supplemented by theoretical arguments.

Intuition for this model can come from the fact that the close-packed square lattice dimer model has an effective field theoretic description in terms of a continuum Gaussian (free boson) field. The mapping of the square lattice classical dimer model to a height model,Ardonne et al. (2004) proceeds by assigning a height variable on each plaquette. On the even sublattice, going clockwise, the height changes by if a dimer is on a link, and by when there isn’t one. On the odd sublattice, the changes reverse. Moreover, the dimer hard-core constraint implies that there are only local dimers configurations on each site, guaranteed in the dual height model by the compactification . In terms of the rescaled field the actual action is

(4) |

with at the non-interacting dimer point and with a charge perturbation, , biasing the coarse-grained height field to take integer values.

The different observables (and hence perturbations) of this theory consist of the composite operators with units of charge and units of vorticity.Kadanoff and Brown (1979); Kadanoff-1979b Their scaling dimensions are

(5) |

As shown in Ref.[Papanikolaou et al., 2007], the columnar order parameter corresponds to the elementary charge operator , which is relevant with , while the nematic order parameter corresponds to which is irrelevant at the free dimer point. The operator , which naturally appears due to discreteness of the height values is a strongly irrelevant operator, but it is the one that drives the KT transition, when the plaquette attractive interaction is added. Moreover, the hole operator, corresponding to a hole density, is represented by the fundamental vortex operator , a typically strongly relevant operator which drives the system to a high temperature liquid. Finally, by using OPE based on the density operator definitions Fradkin et al. (2004), it is easy to show that the leading effect of the dimer-aligning interaction on the same row or the same column at the non-interacting dimer point is just a renormalization of the free field stiffness. This conclusion is verified in the simulations, since in the close-packed limit (and infinitesimal) the columnar crystalline phase appears. On the dilute lattice, such a picture breaks down and the emergence of the phase diagram is not understood analytically.

## Iii Results from Monte Carlo Simulations

We simulate the partition function of this dimer model with configurational energy of Eq.(3) by means of a Monte Carlo worm-algorithm with a local heat-bath detailed balance condition.Alet et al. (2006) We have set the energy scale to be . The study is made in the canonical ensemble, similar to previous analogous studies.Papanikolaou et al. (2007); Alet et al. (2006); Charrier et al. (2008); Charrier and Alet (2010) The observables, apart from the energy related specific heat, are related to the possible symmetry breaking, the columnar and nematic order. The local columnar order parameter Alet et al. (2005b); Papanikolaou et al. (2007) is defined with respect to the dimer occupation number at each site ,

(6) |

with the global order parameter being

(7) |

and is normalized such that the four columnar states correspond to

(8) |

A columnar state breaks the rotation and the translation symmetries of the lattice.

The rotation symmetry breaking (nematic) order parameter is defined to be

(9) |

where and are, respectively, the number dimers on horizontal and vertical bonds in configuration . One also considers the corresponding susceptibilities . In particular, for a second order transition

(10) |

while for a first order transition

(11) |

Finally, the Binder cumulant

(12) |

is a scale-invariant quantity in the case of a continuous transition, and should thus exhibit a crossing point at criticality as a function of the system size. In addition, the finite size-scaling (FSS) of its derivative with respect to the temperature,

(13) |

provides a direct way to determine the critical exponent of the correlation length for nematic fluctuations.

### iii.1 Case , : attractive dimer-aligning interactions (but no plaquette interactions) on a close-packed dimer lattice

Firstly, we investigate the simplest case, being the case where only the dimer-aligning attractive interaction is present, and , on a close-packed dimer lattice. In this limit the system has a discrete “sliding” symmetry which, for periodic boundary conditions leads to a sub-extensive ground state entropy. This feature indicates that the ground state may be any of the possible tilted phases, such as those discussed in Ref.[Fradkin et al., 2004]. This feature poses serious numerical difficulties. Nevertheless, as we will see in Section III.2, these difficulties are suppressed once , which favors columnar order, no matter how small it is. Since in the case of and the transition from the four-old degenerate columnar order state is in the KT universality class,Alet et al. (2005a); Papanikolaou et al. (2007) here we might expect the same scaling. However the existence of many tilted ground states (each with a four-fold global degeneracy) makes the analysis more difficult and convergence problematic.

Close-packed dimer configurations on the square lattice present critical correlations at as aforementioned. While from the field theory (and renormalization group arguments) we expect, as described, a renormalization of the KT-transition temperature toward a solid phase,which should resemble a columnar solid under appropriate boundary conditions on an infinite lattice, our simulations present a complex picture shown in Figs. 3, 4, 5, and 6. From this data, it is reasonable to conclude that this is an Ising transition, given that the anomalous exponent of correlation function of the rotational order parameter (having a parent Ising symmetry) appears to be , consistent both with KT and with Ising universality, and that the columnar order is absent. On the other hand, the specific heat does seem to diverge is consistent with either a cusp or with a weak a logarithmic divergence, both of which have an exponent . this may be consistent with an Ising transition, which has its prototypical logarithmic divergence. In a KT transition, the specific heat displays a peak clearly above the critical temperature (for example detected through the peak of the rotational susceptibility) and does not diverge; In our case, the specific heat is clearly non-divergent but the peak appears unconventionally close () to the critical temperature for the system sizes considered (as noticed in a comparison of Figs. 5 and 6). In addition, Binder cumulants (not shown) have non-Ising behavior which display a crossing but not at the universal expectations for Ising systems. It may well be that this discrepancy implies that the Ising signatures mask the true nature of the transition. The simulations were performed on double-periodic boundary conditions for reasons of numerical efficiency.

However, in this
situation the sub-extensive entropy discussed above is not suppressed. Hence our simulations cannot distinguish a columnar state from any
state with broken translational invariance but at finite tilt. Since the entropy is
sub-extensive the associated discrete “sliding” symmetry must be spontaneously broken. However,
performing fixed boundary condition simulations appears numerically intractable.
For all these reasons we believe that our simulations with periodic boundary conditions have not converged at low temperatures for .

### iii.2 Case , : Competition between attractive dimer-aligning and plaquette interactions on a close-packed dimer lattice

The next case study is for and present in the system for a close-packed lattice, but with ; since we aim to focus on the non-perturbative effects of the dimer-aligning interaction. It is clear that in the opposite limit the dimer-aligning perturbation would only shift the location of the KT transition with no change in the universality class. However, in the limit we probe , the system’s behavior is rather complex, exactly due to the non-trivial crossover that we find in the case where . While the maximum system size we studied is many times larger than the sizes studied in Ref. [Alet et al., 2005a], it appears still that the system has not converged to a well defined transition. The expectation is that the system should display a KT transition since the sub-extensive entropic degeneracy is lifted for any finite value of . In Figs. 7, 8, 9, 10, and 11 we present strong evidence that this transition has not converged. The rotational order parameter displays a clear transition with its susceptibility displaying a large peak with scaling similar to Ising transitions (), and its Binder cumulant shows a crossing point close to the one expected for Ising transitions (1.13), even though there is an apparent ongoing drift of the crossing point toward higher values. The correlation length exponent , derived from the slope of the cumulant at the crossing, appears to be close to the Ising value , but a drift is observed with increasing system size. But, while these observations appear to be consistent with an Ising transition (similar observations with the case), the specific heat displays a clear transition, with an apparent non-Ising critical exponent , while there is some drift as well toward lower values. Further, the columnar order parameter apparently displays a a very weakly-converging transition, where for it is not even detectable (cf. Fig. 7). Its susceptibility displays a very weak convergence, with a very wide hump for moderate system sizes, and apparently , but again the shape of the susceptibility does not allow us to conclude that this exponent (derived from just the peak’s scaling) can be accurate – given that in such estimates there is always the implicit assumption of a scaling collapse, i.e. the shapes of the curves shall be self-similar in the scaling regime. Further, its Binder cumulant does not display a crossing, another evidence of a strong crossover. We should note that the KT transition is the only known transition where the Binder cumulants of the order parameter do not display a crossing (for XY models), so the observed behavior could signify the onset of KT transition scaling. This collection of results suggests that our data is still inside a strong crossover region even at the largest sizes we studied.

### iii.3 Case , : Only dimer-aligning interaction on a dilute dimer lattice and the emergence of an Ising phase transition towards a nematic fluid

Then, we proceed by considering the effect of finite hole density . From previous work of some of usPapanikolaou et al. (2007) (shown in Fig. 2A), there is concrete intuition on the possible expectations when a plaquette interaction is included. However, it has not been known how a dimer-aligning interaction affects the phase structure. As shown in Fig. 2C, we find a new phase, the dimer Ising Nematic phase, which is characterized by orientational order of the dimers (as seen in Fig. 1) but no translational symmetry breaking, i.e. the system remains a fluid. It is characteristic that the stability of this phase over possible effects from plaquette interactions is guaranteed by the presence of holes, since they provide a source of extensive entropy to the liquid, a feature absent at close packing due to the hard-core dimer frustration.

0 | 0.90(2) | - | 0.25(4) | - | - |

3.125 | 0.84(1) | 0.1(1) | 0.25(4) | 1.0(1) | 1.16(5) |

4.6875 | 0.81(1) | 0.1(1) | 0.29(6) | 1.0(1) | 1.16(4) |

6.25 | 0.78(1) | 0.1(1) | 0.25(4) | 0.95(7) | 1.15(4) |

9.375 | 0.72(1) | 0.1(1) | 0.26(4) | 1.1(1) | 1.13(2) |

12.5 | 0.67(1) | 0.07(7) | 0.22(5) | 1.1(2) | 1.13(2) |

25 | 0.49(1) | 0.03(3) | 0.24(5) | 1.1(2) | 1.14(3) |

- | 0 | 1/4 | 1 | 1.13 |

At finite small dilution but no plaquette interaction , we find a clear transition line that separates the isotropic hole-dimer liquid phase from the Ising nematic phase. This transition is clearly in the Ising universality class, as verified by the values of the anomalous exponents as well as the Binder cumulant crossing values, being in consistency with the Ising universality. The evidence is presented in Figs. 12, 13, 14, 15 and Table 1. While in the figures we show a typical case at dilution and the consistency shown by the order parameters, their susceptibilities, the specific heat and the Binder cumulants, at the table we show our estimations for particular critical exponents and universal cumulant values for a number of hole densities we studied, displaying a strong consistency. For , a clear phase transition is detectable around and the rotational symmetry is spontaneously broken at low temperature as shown by the rotational order parameter, but the translational symmetry is not, as the vanishing, with the system size, columnar order parameter shows. The evidence for the Ising universality class comes from the specific heat which is consistently showing scaling with , the rotational susceptibility which gives and the Binder cumulant ratio which not only displays a crossing point at the expected Ising location, but also its derivative at the crossing scales with system size in a way that shows that .

Finally, in Figs. 16, 17 we show how the nematic correlations evolve with system size and hole density at a low temperature . As can be seen from Fig. 2C, at this temperature and all but the close-packed case, the system is an Ising nematic liquid which is stable to perturbations (such as the plaquette interaction). Through the scaling with the system size, it is clear that at long distances dimers have liquid correlations, albeit double than at infinite temperature, due to the fact that deep in the nematic phase they are primarily oriented in the same direction. The correlations have clearly converged for . As a function of hole density, we observe the dramatic effect of having non-decaying correlations at , signifying the absence of liquid correlations, but they immediately become liquid-like as the system gets dilute. Starting from the perfect close-packed crystal (at zero dilution, cf. Fig. 2C), assuming fixed boundary conditions, then to first order, if some dimers (introducing holes) are removed from the lattice, then there are ways to remove them, therefore contributing an extensive entropy without though destroying the columnar order. Further, to second order, at the expense of energy , each pair of holes has a sub-extensive contribution by moving freely on the line it was created, again winning over the crystalline confining tendency. Since a pair can be formed on each line, it leads to an extensive entropy contribution as well, that can dominate over the energy cost. Therefore, the nematic liquid is also deconfining along the oriented direction.

We should note that the Ising transition line, should end at a tricritical point and coexistence at a higher from the ones we investigated. The investigation of this discontinuous transition is beyond the scope of this work, but should be in the tricritical Ising universality class.

### iii.4 Case , : Competition between dimer-aligning and plaquette interactions in a dimer model at finite hole density displaying two Ising phase transitions

Finally, we investigate the case of a finite hole density, where the system contains both the dimer-aligning and plaquette interactions. In this case, we observe a more complex phase diagram, where two Ising transition lines emerge from the close-packed limit (where the transition is in the KT universality class possibly, as we argued) (shown in Fig. 2B), separating the Ising nematic dimer liquid from the disordered isotropic dimer-hole liquid at high temperatures and the columnar solid phase at low temperatures. The evidence for our conclusion that both transitions are in the Ising universality class are displayed in Figs. 18, 19, 20, 21, 22, 23 for a particular hole density and interaction strength , chosen for clarity. Our study is extended over several choices of and (not shown).

The columnar order parameter displays a transition at , which is remarkably weaker than the transition displayed by the rotational order parameter at . The weakness of the low temperature transition can be further observed in the behavior of the specific heat, which displays a strong peak (with a very weak divergence, as expected for Ising transitions) at but a very weak peak at .

The columnar order susceptibility displays a strong peak at with a peak which shows an anomalous exponent , lower than the Ising-expected , but expected given the weak convergence we observe throughout for this transition with the system size. However, the nematic order susceptibility shows a clear Ising behavior at the corresponding nematic-disordered transition, with a clear Ising anomalous exponent for over two decades of scaling, confirming the quick convergence of that transition.

Finally, regarding the behavior of the Binder cumulants, as possibly expected, the rotational order’s cumulant shows a textbook behavior at the high-temperature nematic-disordered transition, where there is a strong crossing at , which given our definitions, is the universal value expected for the Ising universality class. Further, the cumulant’s derivative scales in a consistent way with the Ising expectations. However, the columnar order’s cumulant displays a complex behavior with two crossings at both transitions, since the columnar order parameter displays fluctuations sensitive at both transitions. The lower transition’s Binder ration displays a crossing around that appears to drift toward lower value and the Binder ratio derivative displays Ising-like scaling (with ). We believe that the Binder ratio has not yet converged to the Ising value for finite-size and finite-anisotropy reasons. Namely, the columnar order parameter is defined in a unit cell while the simulation box has always been square . In studies of the Binder ratio for the square-lattice Ising model on rectangular boxes Selke (2006), it was shown that the critical Binder cumulant value had a strong dependence on the system’s aspect ratio. We expect a strong dependence also in our case, and simulations in thin strips are expected to show a much faster convergence for this transition.

## Iv Mean-Field Theory

The understanding of this Ising nematic phase transition can be elucidated also through a mean-field study, similar to the approach followed in Ref.[Papanikolaou et al., 2007]. Using a Grassmann representation for the dilute dimer model and performing two Hubbard-Stratonovich transformations, with the exact same field definitions, we find the following effective thermodynamic potential,

(14) |

where are order parameters coupled to neighboring same-axis aligned dimers and hole density respectively, while and where and the hole chemical potential. Further, is non-zero and unity when are neighboring sites on the same lattice line, i.e. is non-zero and unity when are nearest neighboring sites. The only difference with the study of Ref.[Papanikolaou et al., 2007] is that the attractive neighboring interaction is in a different direction (aligning instead of plaquette).

We use the ansatz and and we look for solutions of the extremal equation

(15) |

The equations can be either solved numerically or analytically, with the latter being in an expansion in (small dimer density) and then, (high temperature), giving finally the leading term for the critical hole density at high temperatures to be

(16) |

leading to the dashed line shown in Fig. 24.

## V Conclusions

In this paper we examined the phase diagram of a 2D dimer model on the square lattice at zero and finite hole densities. The important difference between this work and previous studies is the consideration of the interplay of attractive interactions on a plaquette (which was considered previously) with attractive (aligning) dimer interactions on the same row (and column) for dimers on next-nearest neighboring bonds of the lattice. By means of large-scale Monte Carlo simulations, combined with a scaling analysis and a mean-field theory we find that the phase diagram of this simple model, shown in Fig. 2, in addition to a line of fixed points at zero hole density, generally contains three phases: a dimer-hole liquid, a columnar phase (valence bond solid) and a novel Ising nematic phase. We presented detailed results for the scaling behavior of the columnar and nematic (orientational) order parameters, specific heat and Binder cumulants at the various transitions. In the general case our results are consistent with two Ising transitions as the interaction and hole density are varied.

These results are of interest to both the understanding of classical dimer models and to quantum dimer models which arise in the context of frustrated quantum magnets. In the latter context the theory we presented here describes the behavior of RVB-type wave functions at finite hole density. All phases described here are non-magnetic as the dimers represent valence bond singlets.

Experimentally relevant models of nematic phases should have a columnar phase dominating at low temperatures, consistent with strong experimental evidence in cuprate and iron pnictides materials. In this work, we considered a framework that includes geometric frustration originating in the hard-core dimer exclusion effect on the square lattice. We find that exactly because of the geometric frustration, a nematic phase at low temperatures is always unstable to columnar fluctuations, due to entropic reasons, implying that negligible amounts of the plaquette interaction yields a columnar phase at low temperatures. It would be very interesting to investigate similar phase diagrams on different lattices such as the kagome, where a spin-liquid has been identified as the ground-state of the Heisenberg antiferromagnetYan et al. (2011) pointing toward useful dimer descriptions for the ground-state wavefunction.Schwandt et al. (2010); Rousochatzakis et al. (2013)

###### Acknowledgements.

We would like to thank Cristian Batista, Steven Kivelson, Jörg Schmalian and Bernard Nienhuis for inspiring discussions. This work is supported in part by the National Science Foundation through the grant DMR-1064319 at the University of Illinois (EF), and through a DTRA grant No. 1-10-0021 (SP) at Yale University. We note that we used the ALPS libraries alps (). This work benefited greatly from the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center as well as the computing facilities of MPI-PKS Dresden## References

- Rokhsar and Kivelson (1988) D. S. Rokhsar and S. A. Kivelson, Phys. Rev. Lett. 61, 2376 (1988).
- Moessner and Sondhi (2001a) R. Moessner and S. L. Sondhi, Phys. Rev. Lett. 86, 1881 (2001a).
- Fradkin and Kivelson (1990) E. Fradkin and S. Kivelson, Mod. Phys. Lett. B 4, 225 (1990).
- Papanikolaou et al. (2007) S. Papanikolaou, E. Luijten, and E. Fradkin, Phys. Rev. B 76, 134514 (2007).
- Fradkin (2013) E. Fradkin, Field Theories of Condensed Matter Systems (Cambridge University Press, Cambridge, UK, 2013), 2nd ed.
- Moessner et al. (2000) R. Moessner, S. L. Sondhi, and P. Chandra, Phys. Rev. Lett. 84, 4457 (2000).
- Moessner and Sondhi (2001b) R. Moessner and S. L. Sondhi, Phys. Rev. B 63, 224401 (2001b).
- Alet et al. (2005a) F. Alet, J. L. Jacobsen, G. Misguich, V. Pasquier, F. Mila, and M. Troyer, Phys. Rev. Lett. 94, 235702 (2005a).
- Fradkin et al. (2004) E. Fradkin, D. A. Huse, R. Moessner, V. Oganesyan, and S. L. Sondhi, Phys. Rev. B 69, 224415 (2004).
- Vishwanath et al. (2004) A. Vishwanath, L. Balents, and T. Senthil, Phys. Rev. B 69, 224416 (2004).
- Kivelson et al. (1998) S. Kivelson, E. Fradkin, and V. Emery, Nature 393, 550 (1998).
- Fradkin et al. (2010) E. Fradkin, S. A. Kivelson, M. J. Lawler, J. P. Eisenstein, and A. P. Mackenzie, Annu. Rev. Condens. Matter Phys. 1, 7.1 (2010).
- Fradkin (2012) E. Fradkin, in Modern Theories of Many-Particle Systems in Condensed Matter Physics, edited by D. Cabra, A. Honecker, and P. Pujol, Proceedings of the Les Houches Summer School, 2009 (Springer, Berlin, Germany, 2012), vol. 843 of Lecture Notes in Physics, pp. 53–116.
- Kivelson et al. (2003) S. A. Kivelson, I. Bindloss, E. Fradkin, V. Oganesyan, J. Tranquada, A. Kapitulnik, and C. Howald, Rev. Mod. Phys. 75, 1201 (2003).
- Daou et al. (2010) R. Daou, J. Chang, D. LeBoeuf, O. Cyr-Choinière, F. Laliberté, N. Doiron-Leyraud, B. J. Ramshaw, R. Liang, D. A. Bonn, W. N. Hardy, et al., Nature 463, 519 (2010).
- Hinkov et al. (2008) V. Hinkov, D. Haug, B. Fauqué, P. Bourges, Y. Sidis, A. Ivanov, C. Bernhard, C. T. Lin, and B. Keimer, Science 319, 597 (2008).
- Lawler et al. (2010) M. J. Lawler, K. Fujita, J. W. Lee, A. R. Schmidt, Y. Kohsaka, C. K. Kim, H. Eisaki, S. Uchida, J. C. Davis, J. P. Sethna, et al., Nature 466, 347 (2010).
- de la Cruz et al. (2008) C. de la Cruz, Q. Huang, J. W. Lynn, J. Li, W. Ratcliff II, J. L. Zaretsky, H. A. Mook, G. F. Chen, J. L. Luo, N. L. Wang, et al., Nature 453, 899 (2008).
- Chu et al. (2010) J.-H. Chu, J. G. Analytis, K. D. Greve, P. L. McMahon, Z. Islam, Y. Yamamoto, and I. R. Fisher, Science 329, 824 (2010).
- Chuang et al. (2010) T.-M. Chuang, M. Allan, J. Lee, Y. Xie, N. Ni, S. Bud’ko, G. S. Boebinger, P. C. Canfield, and J. C. Davis, Science 327, 181 (2010).
- Fang et al. (2008) C. Fang, H. Yao, W.-F. Tsai, J. P. Hu, and S. A. Kivelson, Phys. Rev. B 77, 224509 (2008).
- Xu et al. (2008) C. Xu, M. Müller, and S. Sachdev, Phys. Rev. B 78, 020501(R) (2008).
- Alet et al. (2006) F. Alet, Y. Ikhlef, J. L. Jacobsen, G. Misguich, and V. Pasquier, Phys. Rev. E 74 (2006).
- Castelnovo et al. (2007) C. Castelnovo, C. Chamon, C. Mudry, and P. Pujol, Annals of Physics 322, 903 (2007).
- Ardonne et al. (2004) E. Ardonne, P. Fendley, and E. Fradkin, Annals of Physics 310, 493 (2004).
- Fendley et al. (2002) P. Fendley, R. Moessner, and S. L. Sondhi, Phys. Rev. B 66, 214513 (2002).
- Anderson (1987) P. W. Anderson, Science 235, 1196 (1987).
- Kivelson et al. (1987) S. A. Kivelson, D. Rokhsar, and J. P. Sethna, Phys. Rev. B 35, 865 (1987).
- Bergman et al. (2007) D. L. Bergman, R. Shindou, G. A. Fiete, and L. Balents, Phys. Rev. B 75, 094403 (2007).
- Lv and Phillips (2011) W. Lv and P. Phillips, Phys. Rev. B 84, 174512 (2011).
- Abanin et al. (2010) D. A. Abanin, S. A. Parameswaran, S. A. Kivelson, and S. L. Sondhi, Phys. Rev. B 82, 035428 (2010).
- Kumar et al. (2013) A. Kumar, S. A. Parameswaran, and S. L. Sondhi, Phys. Rev. B 88, 045133 (2013).
- Castelnovo et al. (2004) C. Castelnovo, C. Chamon, C. Mudry, and P. Pujol, Annals of Physics 318, 316 (2004).
- Alet et al. (2005b) F. Alet, J. L. Jacobsen, G. Misguich, V. Pasquier, F. Mila, and M. Troyer, Phys. Rev. Lett. 94 (2005b).
- Charrier et al. (2008) D. Charrier, F. Alet, and P. Pujol, Phys. Rev. Lett. 101, 167205 (2008).
- Charrier and Alet (2010) D. Charrier and F. Alet, Phys. Rev. B 82, 014429 (2010).
- Xu and Moore (2004) C. Xu and J. E. Moore, Phys. Rev. Lett. 93, 047003 (2004).
- Heilmann and Lieb (1979) O. J. Heilmann and E. H. Lieb, J. Stat. Phys. 20, 679 (1979).
- Kadanoff and Brown (1979) L. P. Kadanoff and A. C. Brown, Annals of Physics 121, 318 (1979).
- Selke (2006) W. Selke, J. Stat. Mech. 2007, P04008 (2007).
- Yan et al. (2011) S. Yan, D. A. Huse, and S. R. White, Science 332, 1173 (2011).
- Schwandt et al. (2010) D. Schwandt, M. Mambrini, and D. Poilblanc, Physical Review B 81, 214413 (2010).
- Rousochatzakis et al. (2013) I. Rousochatzakis, Y. Wan, F. Mila, and O. Tchernyshyov, arXiv preprint arXiv:1308.0738 (2013).
- (44) A. F. Albuquerque etal., J. Magn. Magn. Mater. 310, 1187 (2007); M. Troyer, B. Ammon, and E. Heeb, Lect. Notes Comput. Sci. 1505, 191 (1998); see http ://alps.comp-phys.org.