A Monte Carlo algorithm

Critical behavior in the cubic dimer model at nonzero monomer density

Abstract

We study critical behavior in the classical cubic dimer model (CDM) in the presence of a finite density of monomers. With attractive interactions between parallel dimers, the monomer-free CDM exhibits an unconventional transition from a Coulomb phase to a dimer crystal. Monomers acts as charges (or monopoles) in the Coulomb phase and, at nonzero density, lead to a standard Landau-type transition. We use large-scale Monte Carlo simulations to study the system in the neighborhood of the critical point, and find results in agreement with detailed predictions of scaling theory. Going beyond previous studies of the transition in the absence of monomers, we explicitly confirm the distinction between conventional and unconventional criticality, and quantitatively demonstrate the crossover between the two. Our results also provide additional evidence for the theoretical claim that the transition in the CDM belongs in the same universality class as the deconfined quantum critical point in the JQ model.

pacs:
64.60.De, 64.60.F-, 75.40.Mg

I Introduction

According to the Landau paradigm, the critical behavior occurring near a second-order phase transition can be understood through the long-wavelength fluctuations of the order parameter.(2) Cases where such a picture does not correctly describe the critical properties, termed non-Landau or simply “unconventional” transitions, present an interesting extension of the established theory of thermal and quantum criticality.

Recent interest in this idea was sparked by the theoretical prediction of such a quantum phase transition in certain two-dimensional antiferromagnets.(3) In these systems, a continuous transition of the Landau paradigm is excluded by the presence of distinct order parameters in the two adjacent phases. It has been argued that the critical behavior is instead governed by “spinons”, quasiparticles carrying spin of , that are deconfined only at the critical point.(3) Evidence for the existence of such “deconfined quantum criticality” has been provided by quantum Monte Carlo studies of certain spin models.(4); (5); (6); (7); (8)

A conventional Landau transition is also excluded for some thermal transitions, when the “disordered” high-temperature phase has some kind of topological order.(9) This can occur in constrained statistical systems, and leads to a second set of unconventional transitions that are, on the surface, only remotely connected to deconfined quantum criticality. Such examples are particularly interesting when the low-temperature phase has broken symmetry, and so the transition separates phases with distinct types of order, topological above the critical temperature and symmetry-breaking below.

It is remarkable that unconventional criticality of the latter variety is exhibited by one of the simplest conceivable constrained classical systems, the close-packed dimer model on the cubic lattice.(9) This model describes the statistics of “dimers”, objects that occupy pairs of neighboring sites of a lattice. In a close-packed dimer model, the allowed configurations are constrained so that every site is covered by exactly one dimer; we also consider here the case when “monomer” defects, where a site’s occupation number deviates from unity, are permitted with a large but finite energy cost .

In the constrained limit, , and with no other interactions, the cubic dimer model (CDM) exhibits a liquid-like phase with strong correlations but no static order.(10); (11) This can be interpreted as the Coulomb phase of an effective lattice gauge theory, in which the monomers act as charges (“monopoles”). A pair of test monomers inserted into the constrained system is subject to an effective interaction, induced by the fluctuations of the dimers. In the Coulomb phase, this interaction obeys a Coulomb law for large separation, implying that the defects can be separated to infinity with finite energy cost, and are hence deconfined.

In the presence of additional interactions between the dimers, a transition can occur into a phase where monomers are instead confined, by an effective potential proportional to their separation. With attractive interactions between neighboring parallel dimers, the phase at low temperature has columnar order [see Fig. 1(a)] that breaks the rotation and translation symmetries of the lattice.(9); (12) While a conventional order parameter can be identified, the critical behavior at this “confinement transition” is modified by the fluctuations of the gauge theory, in ways that cannot be described by a conventional Landau theory.(11); (13); (14); (15)

In this work, we study the critical properties of this transition beyond the constrained limit. Reducing from allows a finite density of monomers in the dimer model, screening the long-range interactions between test monomers and rendering the confinement criterion invalid. The confinement transition at is therefore replaced by a conventional phase transition, described by Landau theory. For sufficiently large , however, the behavior is still strongly influenced by the presence of the nearby confinement critical point. We study this behavior using large-scale Monte Carlo simulations, to test predictions of scaling theory.

In particular, previous work(16); (17) has presented general theoretical arguments that the fugacity is a relevant scaling field under the renormalization group (RG) at the fixed point corresponding to the confinement transition. This result provides quantitative predictions for the behavior at nonzero monomer density, including for example, the shape of the phase boundary at . These applications of scaling theory have previously been tested in certain confinement transitions in spin ice,(16); (17) in which spontaneous symmetry breaking was absent. The present case is more demanding computationally, but also displays considerably richer phenomena, with symmetry breaking coinciding with deconfinement at the transition, and the phase boundary extending to .

This alternative perspective is complementary to those exploited in previous Monte Carlo studies of the CDM, which have been restricted to the fully constrained limit.(9); (12); (14); (15); (18); (19) Our numerical results are entirely consistent with both the qualitative and quantitative predictions of scaling theory. They demonstrate a clear distinction between the critical behavior at and , and thereby provide a particularly stringent test of the interpretation of the transition through confinement.

An additional interesting feature of the CDM transition is that its universality class is expected to be the same as that of the original deconfined critical point in quantum spin models.(3) This allows us to compare our results for this transition with those for the critical point in the JQ model.(4) (Limited comparisons have already been made using results at .(19)) Agreement is again found, within the limits of accuracy that are inherent in such comparisons. Since the evidence for a continuous, non-Landau transition in the CDM is substantial, these results provide quite strong support for the theoretical scenario of deconfined quantum criticality.

Outline

In Section II, we review the cubic dimer model and its phase structure, first in terms of the original dimer degrees of freedom and then through the language of an effective gauge theory and the associated concept of confinement. We then focus, in Section III, on the phase transitions from the dimer crystal, reviewing the critical theories for the transitions and then, in Sections III.2.1III.2.8, presenting a number of predictions for behavior in their neighborhood. These are tested in detail in Section IV, using numerical results from large-scale Monte Carlo simulations of the dimer model. We conclude in Section V with a restatement of our main results and a brief discussion of their significance for the dimer model and for unconventional transitions more broadly. Some details of the Monte Carlo algorithm are given in an Appendix.

Ii Model and phase diagram

ii.1 Cubic dimer model

We consider a classical statistical model in which the elementary degrees of freedom are dimers, defined on the links of a cubic lattice. Labeling each site by a vector of integers, one can define an occupation number , giving the number of dimers on the link joining and , where and is a unit vector. The number of dimers touching site is given by

 n(r)=∑μ[dμ(r)+dμ(r−δμ)]. (1)

In a close-packed hard-core dimer model, one constrains on every site ; instead, we also allow monomers, where a site has either no dimer or multiple dimers. We are interested in the case where the elementary defects, with or , cost a large energy , and so occur at low density. Larger deviations, where , are expected to be negligible near the transition, and are excluded. Defining the on-site energy , we therefore set , , and otherwise.

The transition to an ordered state is driven by an attractive interaction between pairs of nearest-neighbor parallel dimers,(9) whose number is

 N2=∑r∑μ,ν≠μdμ(r)dμ(r+δν). (2)

It is also necessary (see Refs. (18); (19) and the final paragraph of Section II.2) to include additional interactions, and we follow Charrier and Alet (19) by using a term that counts the number of parallel dimers around cubes of the lattice,

 N4=∑r∑μ,ν≠μρ≠μ,νdμ(r)dμ(r+δν)dμ(r+δρ)dμ(r+δν+δρ). (3)

The full configuration energy can then be written as

 E=∑rE(n(r))+Edimer (4) Edimer=v2N2+v4N4.

We are interested in the case of attractive interactions, , and in the following we choose units where .

The partition function is given by

 Z=∑{dμ(r)}e−E/T. (5)

We use a lattice with periodic boundary conditions and an even number of sites in each direction.

ii.2 Phase structure

Considering first the case where , the energy is globally minimized by a columnar dimer crystal, as shown in Fig. 1(a).

The dimers are aligned with a particular cubic axis, and form sheets in which each dimer has four parallel neighbors. There are six such configurations—one example has for and even, and otherwise—with energy .

The broken symmetry in these configurations is characterized by the “magnetization” (per site),

 mμ=2L3∑r(−1)rμdμ(r). (6)

The minimal-energy configurations maximize , having . The order parameter is reduced by fluctuations for but it remains nonzero until a critical temperature , at which point it vanishes and the system enters a phase without order. This phase boundary, between a columnar dimer crystal and a disordered phase, is illustrated in Fig. 2 (which shows results of Monte Carlo simulations; see Section IV).

This description of the phase structure, based on the order parameter, is essentially identical whether is infinite or merely large. (Decreasing from infinity has the quantitative effect of reducing the order–disorder transition temperature , because it favors disorder in the dimer degrees of freedom.) From the point of view of the Landau criterion for phases and transitions, there is no distinction between these two cases.

In fact, the two phases at , illustrated in Fig. 1(b) and (c) and occurring when is, respectively, infinite and finite, are qualitatively distinct. As we will argue, and as our numerical results explicitly demonstrate, there is an important distinction between the transitions in these two cases, based on topological order and deconfinement. To clarify this distinction, we review in Section II.3 the gauge-theoretical description of the dimer model.(10)

The three phases of interest are not qualitatively affected by nonzero , although for large positive the columnar states no longer minimize the energy and an additional ordered phase (“Crystal II”) can appear at low .(19) More importantly, the transition from the columnar dimer crystal into the higher-temperature phases changes from continuous to first order as is decreased below a value .(19) The tricritical point separating these two types of transition apparently strongly influences the effective critical properties near .(9); (19); (18) In our numerical results, we therefore use the value , which allows us to access the critical behavior of the continuous confinement transition, while being small enough to avoid complications from additional ordered phases (see Fig. 2 and Section IV.1).

ii.3 Effective gauge theory

In the limit , the dimer model is constrained to have one dimer per site. This constraint can be rewritten in the form of a Gauss law by defining a “magnetic field”(10); (11)

 Bμ(r)=ηr[dμ(r)−1z], (7)

where is the coordination number, and is on the two sublattices. The lattice divergence, defined by

 divrB≡∑μ[Bμ(r)−Bμ(r−δμ)], (8)

obeys

 divrB=ηr[n(r)−1], (9)

and so monomers act as charges (“magnetic monopoles”), with sign depending on the sublattice.

In these terms, the partition function can be rewritten as

 Z=∑{Bμ(r)}|divrB|≤1e−Edimer/T∏rz(divrB)2, (10)

where is the monomer fugacity. In particular, when , , and so

 Z=∑{Bμ(r)}divrB=0e−E% dimer/T. (11)

The configurations to which the system is constrained at , viz. those of close-packed, hard-core dimers, satisfy a lattice Gauss law, . Such a constraint is preserved by an appropriately chosen RG procedure, and is therefore expected to have important consequences for the long-distance physics, qualitatively distinguishing the cases with and .

More explicitly, the result of such a procedure is to replace by a coarse-grained continuum vector field with a corresponding constraint on its continuum divergence, . In the Coulomb phase at , where continues to fluctuate subject to this constraint, the effective continuum action density is

 Lgauge=12K|B|2=12K|∇×A|2, (12)

where . This is a (Gaussian) fixed point under the RG, at which all other analytic terms are irrelevant, and at which correlations are algebraic.(10); (11)

While higher-order terms are irrelevant, and a mass term for is forbidden by gauge invariance, monomer fugacity is a relevant perturbation at the Coulomb-phase fixed point. For any nonzero , the algebraic correlations are cut off at a length scale set by the monomer separation.

ii.4 Confinement and deconfinement

The magnetization is a local quantity that provides a conventional order parameter for the low-temperature phase. At , the phase transition can alternatively be characterized through the concept of confinement.

Consider first the partition function of Eq. (11), but evaluated in the presence of a set of charges at fixed positions,

 Z[Qr]=∑{Bμ(r)}divrB=Qre−Edimer/T. (13)

The “monopole distribution function” can then be defined, in terms of a test pair of oppositely charged monomers at , as

 Gm(r+,r−)=Z[δr,r+−δr,r−]Z. (14)

The function can be considered as the correlation function corresponding to the monomer fugacity.(17) It is related to the effective interaction, induced between the monomers as a result of the fluctuating dimers.

To understand the behavior of , imagine starting in a defect-free system and taking out a single dimer. This leaves two vacant sites, which are on opposite sublattices and hence carry opposite gauge charge. One can move each monomer by rearranging the dimers; all local rearrangements preserve the signs of the charges. In the ordered phase, attempting to separate the monomers in this way leaves behind a trail of disturbance, costing an energy (at least) proportional to its length, so . In the Coulomb phase, there is no order and so moving the monomers simply scrambles an already-disordered background, leaving no such trail. The scrambling induces only an entropic interaction; the monomers act like point charges in the effective magnetic field, and obeys the Coulomb law, .

For , the large-separation limit of is therefore qualitatively distinct in the phases on either side of . In the ordered phase at , the monomers are confined, costing an unbounded energy to be separated to infinity, and so . In the Coulomb phase at , has a finite limit, and so they are deconfined, with a nonzero limit for .

This distinction applies only for ; for nonzero monomer fugacity, charges are screened in the limit of large separation, and always approaches a finite constant for large separation.

Iii Phase transitions

iii.1 Critical theories

We are interested in the phase transition from the dimer crystal, with order parameter , into phases in which the magnetization vanishes.

At , symmetry-breaking order appears simultaneously with the confinement of monomers, and both of these effects must be incorporated in a critical theory. It has previously been argued(13); (14); (15) that the transition is described by a Higgs theory of a noncompact gauge field and -symmetric matter fields . The action density is

 Lcritical =Lgauge+Lmatter+Lmatter--gauge (15) Lmatter =s|φ|2+u(|φ|2)2 Lmatter--gauge =|(∇−iA)φ|2;

the pure-gauge part of the action is as given in Eq. (12). Note that the matter field is minimally coupled to and so, in the language where is the magnetic field, has electric charge.(20)

The transition occurs when is tuned through its critical value . For , corresponding to the lower-temperature phase, condenses and acquires an effective mass term by the Anderson–Higgs mechanism. This in turn eliminates the algebraic correlations and confines the monomers (magnetic charges) through the Meissner effect. Symmetry arguments show that the magnetization order parameter obeys , where is a Pauli matrix, and so becomes nonzero when condenses. The continuous symmetry of is broken by an eighth-order (in ) term,

 Lcubic=−ucubic∑μm4μ, (16)

which is irrelevant at the critical point, but for selects states where is aligned with one of the 6 cubic directions ().

For , we instead expect the transition, between the ordered dimer crystal and the thermally disordered phase, to be described by a conventional Landau theory. The order parameter is the -component vector , which preferentially aligns with one of the six cubic directions. The Landau action is therefore given by

 LLandau=|∇m|2+sm|m|2+um(|m|2)2+Lcubic, (17)

where is the cubic-anisotropy term given in Eq. (16), again with .(21) This theory can be viewed as the result of adding a nonzero density of magnetic monopoles to Eq. (15), suppressing the algebraic correlations of the gauge field and confining the matter field into the charge-neutral bilinear .

In this case, the anisotropy is only quartic in the critical field , and is known to be relevant. For , the transition would be continuous and in the Heisenberg universality class, but is expected to drive it first order.(22) (In fact, this transition appears to be at most very weakly first-order; see Sections III.2.4 and IV.1.)

iii.2 Scaling

The presence of the continuous confinement transition at and implies that the behavior in its neighborhood is governed by the properties of the corresponding RG fixed point. This includes along the transition line at , and in particular determines the shape of the phase boundary as a function of . In addition, the first-order nature of the transition at is apparently weak enough that, as we argue in Section III.2.4, scaling theory can be applied to this transition. These observations lead to a number of quantitative predictions that we are able to test in our simulations.

Scaling near confinement transition

Previous work(16); (17) has shown that, in addition to the reduced temperature

 t=T−Tc(0)Tc(0), (18)

there is a relevant scaling field at the confinement-transition fixed point given by the monomer fugacity . The reduced free-energy density therefore has a singular part obeying(23)

 fs(t,z,h,L)≈|t|2−αΦ±(z/|t|ϕ,L|t|ν,h/|t|2−α−β)%, (19)

where is a universal function, with the subscript indicating dependence on the sign of . For completeness an applied field coupling to the magnetization has been included; unless stated otherwise we set in the following. The exponents , , , and are related to the RG eigenvalues , , and by

 α =2−dyt (20) β =d−yhyt (21) ν =1yt (22) ϕ =yzyt. (23)

Phase boundary

The phase boundary at is a nonanalyticity of the free energy (for ), and hence of [since ]. According to Eq. (19), this boundary has the form

 |Tc(z)−Tc(0)|∼z1/ϕ. (24)

While this result is asymptotically exact, there are substantial corrections for nonzero that must be incorporated for a correct interpretation of the numerical data. The most important comes from the leading irrelevant scaling variable at the fixed point, which we denote , with RG eigenvalue . Incorporating the unknown, nonuniversal constant into Eq. (19), while setting for simplicity, we can write

 fs(t,z)≈|t|2−αΦ±(z/|t|ϕ,uzωϕ/ν). (25)

The value of at which has a nonanalyticity is now dependent on . Assuming that a Taylor expansion for exists around and dropping corrections of higher order in gives

 |Tc(z)−Tc(0)|∼z1/ϕ(1+Czωϕ/ν), (26)

where is a nonuniversal constant.

A functional-RG study(24) of the Higgs theory, Eq. (15), found a small value of the correction exponent, . This implies that the second term in Eq. (26), as well as omitted higher-order corrections, should be significant, and our results are qualitatively consistent with this. It has in fact been argued that there may be logarithmic corrections to scaling in the JQ model,(6); (7) which would imply that at this fixed point; we are unable to exclude this possibility.

Scaling of monopole distribution function

The monopole distribution function defined in Eq. (14) is the correlation function corresponding to the scaling field , and so for has scaling form

 Gm(R;T,L)≈R−2(d−yz)Γm±(R|t|ν,R/L), (27)

where is the separation(25) and is a universal function. This expression allows one to determine the RG eigenvalue , and hence the exponent , by measuring the monopole correlation function at the critical point, and . Together with Eq. (24), this provides a quantitative prediction for the phase boundary at .

Scaling near order–disorder transition

The order–disorder transition at is described by the Landau theory of Eq. (17), an vector model with cubic anisotropy, . At the Heisenberg-class fixed point, is relevant with RG eigenvalue , causing flow to a discontinuity fixed point, describing a first-order transition.(22); (23)

An argument based on RG theory, however, implies that, for small , this first-order behavior requires inaccessibly large . For small and , the trajectory of the RG flow is initially through the neighborhood of the confinement-transition fixed point. The cubic anistropy is irrelevant at this point, with RG eigenvalue , and is consequently renormalized to a smaller value , where , the “RG time” spent in this part of the trajectory, depends on as . The RG flow subsequently passes through the neighborhood of the Heisenberg fixed point, requiring time to reach the outflow trajectory leading to the discontinuity fixed point. The required length scale to reach the latter fixed point, and hence to see a first-order transition, is therefore . (This additional length scale is a consquence of the dangerous irrelevant scaling variable at the confinement fixed point. The same result can equivalently be derived using only scaling theory, following similar logic to that in Section III.2.7.)

Cubic anistropy is apparently only very weakly relevant at the Heisenberg fixed point, with ,(22) and so the length scale grows extremely rapidly with decreasing . In our numerical results (see Section IV.1), we indeed find no clear signatures of a first-order transition at small , and instead see results consistent with scaling up to the largest accessible system sizes.

On moderate length scales, one expects behavior governed by the Heisenberg fixed point, at which the relevant scaling fields are the reduced temperature , defined with respect to , and the applied field . Since this transition belongs in a different universality class from the confinement transition, the corresponding RG eigenvalues and critical exponents are different. We define () as the correlation-length exponent at the order–disorder transition.

Near the order–disorder fixed point, the fugacity determines only nonuniversal quantities, such as the values of irrelevant scaling fields. As is reduced, one expects a crossover to behavior described by the confinement transition, and hence a breakdown of the order–disorder criticality. On general grounds, this should occur when the characteristic length scale for monomers,(17) , exceeds the system size.

Binder cumulant

The magnetization provides an order parameter for symmetry breaking at both the confinement transition () and the order–disorder transition (). The corresponding Binder cumulant,(26)

 B=1−⟨|m|4⟩3⟨|m|2⟩2, (28)

therefore provides a particularly useful quantity for studying the crossover between the two critical behaviors.

The scaling dimension of vanishes at the fixed points corresponding to both transition classes, giving

 B ≈ΦB(tL1/ν,z/|t|ϕ) (confinement) (29) B ≈ΨB(g(z)[T−Tc(z)]L1/ν′) (order–disorder), (30)

where the function incorporates the dependence on of the nonuniversal constant of proportionality. Near the transition for small nonzero , both Eqs. (29) and (30) apply, and the requirement that they be consistent fixes

 g(z)∼z(ν−ν′)/(ϕν′). (31)

Binder cumulant crossings

According to Eqs. (29) and (30), is independent of at the transition, implying that plots of the Binder cumulant versus temperature for different cross at . Corrections to scaling in fact imply that the crossing points are weakly dependent on . Including leading-order corrections at the order–disorder transition replaces Eq. (30) by

 B≈ΨB(g(z)[T−T%c(z)]L1/ν′)+L−ω′u(z)~ΨB(g(z)[T−Tc(z)]L1/ν′), (32)

where is an unknown function, is the RG eigenvalue of the leading irrelevant scaling operator, and expresses the dependence of its coefficient on . Again requiring consistency with Eq. (29), we find

 u(z)∼z−ω′ν/ϕ. (33)

At this order, the Binder cumulants at and have a crossing at given by(26); (27)

 T×(L1,L2)−Tc(z)∼u(z)g(z)L−ω′2−L−ω′1L1/ν′1−L1/ν′2. (34)

In applying this result to numerical data, it is convenient to fix the ratio ; incorporating the dependence of and then gives

 T×(L,ρL)−Tc(z)∼z−(1+ω′ν−ν/ν′)/ϕL−ω′−1/ν′. (35)

For the confinement transition, one similarly finds

 T×(L,ρL)−Tc(0)∼L−ω−1/ν(z=0). (36)

Slope of Binder cumulant at crossing

The derivative of the Binder cumulant, taken with respect to and evaluated at , can be used to determine the correlation-length exponents and . At , one finds the scaling result

 ∂B∂T∣∣∣T=Tc(0)∼L1/ν. (37)

Similarly, for , the slope is

 ∂B∂T∣∣∣T=Tc(z)≈g(z)Tc(0)[Ψ′B(0)L1/ν′+u(z)~Ψ′B(0)L1/ν′−ω′], (38)

where the correction term from Eq. (32) has been included. In both cases, the slope is proportional to in the limit of large .

The corrections in Eq. (38) have magnitude proportional to and therefore become more significant for smaller . In agreement with the general considerations of Section III.2.4, breakdown of order–disorder scaling occurs for system size below , when the “correction” becomes larger than the zeroth-order term.

Flux stiffness

An additional quantity of interest at the confinement transition is the flux stiffness,(9); (28)

 K−1=13L2∑r⊥,μ⟨[ϕμ(r⊥)]2⟩, (39)

where

 ϕμ(r⊥)=∑rrμ=r⊥Bμ(r) (40)

is the net “magnetic flux” through a closed surface, spanning the periodic boundaries, normal to . The flux stiffness vanishes in the thermodynamic limit for and , while it approaches a constant in the Coulomb phase.

At the confinement phase transition, the quantity has zero scaling dimension,(9) and hence shows a similar crossing point to the Binder cumulant and a scaling form

 LK−1(t,z,L)≈ΦK(tL1/ν,z/|t|ϕ). (41)

The scaling dimension of the flux stiffness is not similarly fixed at the order–disorder transition and so, in contrast to the case of the Binder cumulant, the crossing points of do not converge for .

Iv Numerical results

We have studied the dimer model defined in Section II.1 using a Monte Carlo method based on the directed loop algorithm(29) but adapted to allow for violations of the dimer constraint. An update within the set of constrained dimer configurations can be performed by shifting dimers along a closed loop. To change the number of monomers, one shifts dimers along an open path, creating or annihilating monomers at each end. (This corresponds to adding a string of magnetic flux with a monopole at each end.) Details of the algorithm are given in the Appendix.

Simulations were performed on lattices with periodic boundary conditions and sites, with system sizes up to .

iv.1 Phase structure

The phase structure of the model at has previously been studied extensively using Monte Carlo simulations.(9); (12); (14); (15); (18); (19) For , one finds only the Coulomb phase and the columnar dimer crystal, while for an additional crystalline phase, denoted Crystal II, occurs at low .(19) We are mainly interested in the phase boundary between the Coulomb and columnar phases, which is shifted by nonzero but qualitatively unaffected.

The phase structure at nonzero and is shown in Fig. 2. The transition temperature to the columnar dimer crystal is reduced as increases, while the transition into Crystal II is largely unchanged. At larger , the two phase boundaries merge into a single first-order transition between the disordered phase and Crystal II. This sets an upper limit of on the values of at which the transition to the columnar dimer crystal can be studied.

There is strong evidence that the transition at is continuous for but first order for .(19) As noted in Sections III.1 and III.2.4, we expect that the transition at is first order, but only very weakly so for small .

Energy histograms close to the critical temperature, shown in Fig. 3, exhibit no signs of the double-peak structure that would indicate a first-order transition. The latent heat, if nonzero, should grow with , but only a single peak is resolvable (for our largest available system sizes) up to the maximum accessible .

By contrast, for , a first-order transition is clearly seen even for . In this case, the effects of the tricritical point at and are likely important(19) and our scaling conclusions are less reliable.

iv.2 Confinement transition

We begin by reporting our results for the transition at (and ), which has previously been studied in Ref. (19). The Binder cumulant and flux stiffness at are shown in Fig. 4 for various system sizes. The crossover between the low- and high-temperature values becomes steeper for larger , with a sharply defined crossing point appearing.

The crossing temperatures for sizes in the ratio are shown in Fig. 5, along with fits to Eq. (36) and the corresponding expression for the stiffness.

These fits give values for the critical temperatures of (Binder) and (stiffness).

The close agreement between the two methods of calculating provide evidence that the ordering, measured by the Binder cumulant, and confinement, measured by flux stiffness, occur at a single transition. The remaining discrepancy, , suggests either that our estimates of statistical errors are somewhat optimistic or that small systematic errors remain. (For each , the crossing temperatures of and bound the true critical temperature.(30))

The value of can be determined from the derivative of the Binder cumulant at the critical temperature, as described in Section III.2.7, or equivalently using the flux stiffness. The slope of and at are shown in Fig. 6, as functions of , along with fits to Eq. (37) and a similar expression for the stiffness. We find using the Binder cumulant and using the stiffness.

Our results at are consistent with previous results in the CDM with : Charrier and Alet(19) found , (Binder) and (stiffness).

iv.3 Monopole distribution function

As discussed in Section III.2.3, the scaling behavior of the monopole distribution function at can be used to determine . To calculate , defined in Eq. (14), we require the partition function in the presence of a monopole–antimonopole pair, . The directed-loop Monte Carlo algorithm that we use is well-suited to calculating this quantity, as the intermediate states appearing during the construction of a loop involve a pair of inserted monomers.(29)

Since the partition function can be calculated in Monte Carlo simulations only up to an -dependent factor, we calculate the ratio

 G(L)=Gm(Rmax,L)Gm(Rmin% ,L), (42)

with and both fixed and of order unity. This ratio, which is asymptotically proportional to , is shown in Fig. 7, along with a fit giving .

Using our estimate of from the Binder cumulant gives a crossover exponent (with the error dominated by that on ).

iv.4 Order–disorder transition

Figs. 8 and 9 show the Binder cumulant across the order–disorder transition for several values of .

Consistent with an apparently continuous transition, the crossing points converge for the largest values of ; the convergence is slower (requiring larger ) for smaller , in agreement with Eq. (35).

For each value of , polynomial interpolations are used to find the crossing points for pairs of system sizes in the ratio . These are then fit to Eq. (35) in order to find ; examples are shown in Fig. 10.

Note that the fits use fixed values of and appropriate to the Heisenberg universality class.(31) (The fitted values of are not particularly sensitive to these values; consistent results with larger error bars are found by fitting to the exponent.)

The dependence of the transition temperature on monomer fugacity is shown in Fig. 11.

We find that extremely small values of (on the order of ) are required to give a power-law dependence of , indicating that the corrections to scaling are substantial. A power-law fit to the smallest three values, with , gives , which is in good agreement with the value found using calculations at .

To incorporate leading-order corrections, we use Eq. (26), with the value of the correction exponent found in a functional-RG study(24) of the critical action, Eq. (15). This increases the range of over which a reasonable fit is found; we find a consistent value, , by including . (While the nominal error of the fit is reduced by including data at larger , the quoted value is not necessarily reliable: The small magnitude of implies that significant higher-order corrections should also be expected.)

iv.5 Effective correlation-length exponent

The effective value of the correlation-length critical exponent , determined using the slope of at , is shown as a function of in Fig. 12.

While a clear distinction is seen between the exponents for and , confirming the different universality classes for these two cases, it is difficult to draw further conclusions about the nature of the transition at .

According to the scenario outlined in Sections III.2.4 and III.2.7, one would expect to take the value(31) appropriate to the 3D Heisenberg universality class, at least for intermediate values of : For very small , the influence of the confinement fixed point gives corrections (for even quite large ), described by Eq. (38); for larger values the true first-order nature of the transition should become apparent, causing complete breakdown of scaling. The large uncertainty in , caused by slow convergence of the crossings of the Binder cumulant, makes refinement of difficult.

V Conclusions

Our central result, illustrated in Fig. 11, is the quantitative relationship, derived using the scaling theory and confirmed in large-scale Monte Carlo simulations, between the monopole distribution function at and the phase boundary at . The connection between the two follows from their common dependence on the RG eigenvalue of the monomer fugacity at the fixed point describing the confinement transition. There is no reason to expect such a relationship for a conventional order–disorder transition, and so this result provides clear and direct evidence for the claim of a non-Landau transition characterized by confinement.

As illustrated in Fig. 12, we also directly observe a difference between the correlation-length critical exponents at the confinement () and order–disorder () transitions. This is a direct demonstration of the distinct critical behavior in the generic case, where Landau theory is expected to apply, and in the constrained limit, where it is violated. Related previous examples include the unconventional critical behavior observed in a topologically constrained Heisenberg model(32) and recent studies of the square-lattice dimer model at finite hole density.(33)

Consistently with previous work on the CDM, we see no sign of first-order behavior for and .(9); (12); (19); (14); (15); (18) As discussed in Section III.2.4, we expect a very weak first-order transition for small , with Heisenberg-like critical behavior on moderate length scales. Our results at and are indeed broadly consistent with a continuous transition in the Heisenberg universality class, including (1) the absence of a double-peak structure in the energy histogram (Fig. 3), (2) well-defined crossing points of the Binder cumulant (Fig. 8), and (3) good fits to the finite-size behavior of the Binder crossings, using the exponents of the Heisenberg universality class (Fig. 10). We also observe reasonable finite-size scaling of the temperature-derivative of the Binder cumulant, but we are unable to determine the exponent accurately enough to identify the universality class (Fig. 12).

The scaling arguments of Section III.2.4 indicate that the putative first-order transition at may never be visible with realistic system sizes. Given the substantial corrections to scaling exhibited by this transition, however, such arguments should not be considered definitive. Indeed, at , where corrections due to the tricritical point at are significant, clear first-order behavior is seen at . More work, including studies at intermediate , are needed to clarify this picture.

Universality and deconfined criticality

The zero-temperature ordering transition in certain D quantum antiferromagnets, such as the [] JQ model,(4) is claimed(3) to be described by the same critical theory as given in Eq. (15), and should therefore belong in the same universality class. Universality is one of the most remarkable features of critical behavior, and a clear demonstration in the context of non-Landau criticality would be a striking result. Since the confinement transition in the dimer model necessarily occurs at zero monomer density, it would also provide direct evidence for deconfinement at the quantum critical point.

In the JQ model, the formation of a valence-bond solid (VBS) is a consequence of the condensation of monopoles.(3) The exponent describing the appearance of VBS order is therefore related to the monopole scaling dimension by . Values for and for the critical exponent , displayed in Table 1, show reasonable agreement, albeit with large error bars.

While quantitative comparison of exponents is always challenging, it is made even more so in this case by the presence of large corrections to scaling. Significant corrections are observed in all cases,(6) providing a further, qualitative indication of universality. This indicates the presence of a weakly irrelevant scaling variable at the fixed point, an observation consistent with a recent functional-RG study(24) of the critical action Eq. (15), which found a small correction-to-scaling exponent, .

Although scaling behavior is visible over a range of length scales,(34) the possibility remains that the transition in the JQ model is first order.(35) (This should not be confused with the weak first-order transition expected in the CDM at . Berry phases suppress singly-charged monopoles in the JQ model,(3) so the transition is effectively at zero fugacity, .) Indeed, a drift in the effective exponent values with increasing system sizes was observed by Harada et al.,(8) consistent with a weak first-order transition. Universality between different models and different lattices, as suggested by Table 1, would nonetheless imply that the RG flow is through the neighborhood of a critical fixed point. In this case, one might expect that, as in the CDM,(18); (19) additional perturbations could drive the JQ transition truly continuous, or at least extend the length scale over which scaling can be observed.

Acknowledgements.
The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputing Centre (NSC) and High Performance Computing Center North (HPC2N). *

Appendix A Monte Carlo algorithm

In this Appendix, we give a brief description of the numerical methods used for the calculation.

The Monte Carlo routine is based on the directed loop algorithm,(29) adapted to allow for violations of dimer constraints. The system consists of a three dimensional cubic lattice with sites, or nodes, in each direction and periodic boundary conditions. The dimers can occupy the links of the lattice subject to the dimer constraints described in Section II.1. The Monte Carlo algorithm samples the equilibrium thermodynamic distribution of dimer configurations where the energy of a configuration is given by Eq. (4).

Two kinds of updates are used for sampling the configurations; the first has a high acceptance rate but does not alter the number of monomers, while the second can generate and annihilate monomers.

The first kind of update begins by randomly picking a monomer-free node . The node is accepted with a probability (specified below) that depends only on the local configuration around . If accepted, the system transitions into a configuration in which there are two monomers—one on the node and one on a link connected to . This is accompanied by addition or removal of half a dimer as shown in the example in Fig. 13.

The link monomer has a “direction” attribute, which serves to identify one of the two nodes connected to the link as the node “ahead”, and which initially points away from .

The link monomer can hop to one of the six links connected to the node ahead of it, erasing or creating dimers in the process as shown in Fig. 14.

Any such hopping that does not result in a change in the number of monomers on is allowed. If the link monomer hops back to the same edge, the direction of propagation is flipped. If has a monomer, the link monomer can annihilate with the one on the node, thereby terminating the update process.

The probability of a transition from configuration to is given by

 p(k→q)=¯wq−δkqmin(¯w)∑¯w−min(¯w), (43)

where with the energy of the configuration without considering the two monomers created by the update process. The energy of the half dimer is taken to be half that of a full dimer at the same edge.

The second kind of update is similar to the first, but can start on any node, including ones with monomers. Such an update can terminate with the absorption of the link monomer into a node, creating or annihilating a node monomer. The propagation of the link monomer proceeds similarly to the first update, with probability of transition from to of where . Since the number of monomers can change after the update, the energy is calculated with the energy of every node monomer added by the process.

The relative frequency of the two kinds of update was chosen to minimize the convergence time but otherwise did not affect the outcomes. Thermalization of the system from a monomer-free, fully ordered configuration was achieved by performing the updates until approximately forward-propagation steps of the link monomer had been achieved. Configuration samples were subsequently taken once every forward propagation steps. Averages and error estimates of thermal expectation values and their combinations (e.g., Binder cumulants) were obtained by the blocking method.

The crossing point of Binder cumulants and for systems of size and was obtained by first finding an approximate crossing point from a piecewise linear interpolation of the data. A better estimate was obtained by fitting the data in the temperature interval to quadratics,

 B(L1,T) =B0+M1(T−T×)+A1(T−T×)2 (44) B(L2,T) =B0+M2(T−T×)+A2(T−T×)2

(with fit parameters , , , and ).

The -derivative of the binder cumulant was found using a fluctuation–response relation,

 (∂B∂T)z=−13T2⟨|m|4⟩⟨|m|2⟩2[⟨(Edimer−⟨Edimer⟩)|m|4⟩⟨|m|4⟩−2⟨(Edimer−⟨Edimer⟩)|m|2⟩⟨|m|2⟩], (45)

where is the energy excluding the contribution from monomers.

The Monte Carlo algorithm and its implementation were tested using comparisons with previously published results and, for small system sizes, simple alternative methods. For , results were compared with the exact averages obtained from the enumeration of all possible microscopic states. For system size and parameters and , the energy and susceptibility were compared with an elementary Metropolis algorithm. In addition, the self consistency of the slope of the Binder cumulants estimated using Eq. (45) and from curve fitting through Binder cumulants supports the validity of the sampling algorithm.

References

1. L. D. Landau and E. M. Lifshitz, Statistical Physics, Butterworth–Heinemann, New York (1999).
2. T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, Science 303, 1490 (2004); T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P. A. Fisher, Phys. Rev. B 70, 144407 (2004).
3. A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).
4. J. Lou, A. W. Sandvik, and N. Kawashima, Phys. Rev. B 80, 180414(R) (2009).
5. A. W. Sandvik, Phys. Rev. Lett. 104, 177201 (2010).
6. S. Pujari, K. Damle, and F. Alet, Phys. Rev. Lett. 111, 087203 (2013).
7. K. Harada, T. Suzuki, T. Okubo, H. Matsuo, J. Lou, H. Watanabe, S. Todo, and N. Kawashima, arXiv:1307.0501 (unpublished).
8. F. Alet, G. Misguich, V. Pasquier, R. Moessner, and J. L. Jacobsen, Phys. Rev. Lett. 97, 030403 (2006).
9. D. A. Huse, W. Krauth, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 91, 167004 (2003).
10. C. L. Henley, Ann. Rev. Cond. Matt. Phys. 1, 179 (2010).
11. G. Misguich, V. Pasquier, F. Alet, Phys. Rev. B 78, 100402(R) (2008).
12. S. Powell and J. T. Chalker, Phys. Rev. Lett. 101, 155702 (2008); Phys. Rev. B 80, 134413 (2009).
13. D. Charrier, F. Alet, and P. Pujol, Phys. Rev. Lett. 101, 167205 (2008).
14. G. Chen, J. Gukelberger, S. Trebst, F. Alet, and L. Balents, Phys. Rev. B 80, 045112 (2009).
15. S. Powell, Phys. Rev. Lett. 109, 065701 (2012).
16. S. Powell, Phys. Rev. B 87, 064414 (2013).
17. S. Papanikolaou and J. J. Betouras, Phys. Rev. Lett. 104, 045701 (2010).
18. D. Charrier and F. Alet, Phys. Rev. B 82, 014429 (2010).
19. The dual language is used in Ref. (15), so our electric charges are there referred to as magnetic monopoles.
20. A term is also consistent with symmetry, but is irrelevant at the Heisenberg fixed point.(22)
21. J. M. Carmona, A. Pelissetto, and E. Vicari, Phys. Rev. B 61, 15136 (2000).
22. J. Cardy, Scaling and Renormalization in Statistical Physics (Cambridge University Press, Cambridge, UK, 1996).
23. L. Bartosch, arXiv:1307.3276 (unpublished).
24. This result is asymptotically correct for large compared to the lattice scale. The form of Eq. (15) implies that is then isotropic and sublattice independent.
25. K. Binder, Z. Phys. B 43, 119 (1981).
26. A. M. Ferrenberg and D. P. Landau, Phys. Rev. B 44, 5081 (1991).
27. For , the lattice Gauss law implies that is independent of and so the sum over is redundant.(9)
28. A. W. Sandvik and R. Moessner, Phys. Rev. B 73, 144504 (2006).
29. Reduced size (with periodic boundaries) tends to increase the ordering temperature, and hence the Binder crossing point, because it restricts fluctuations to those commensurate with the system period.(23) The stiffness is sensitive to only those fluctuations that change the total magnetic flux through the system. Threading a unit of flux is easier for a smaller system, and so one needs lower temperature to suppress such fluctuations for smaller .
30. M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi, and E. Vicari, Phys. Rev. B 65, 144520 (2002).
31. O. I. Motrunich and A. Vishwanath, Phys. Rev. B 70, 075104 (2004).
32. S. Papanikolaou, D. Charrier, and E. Fradkin, arXiv:1310.4173 (unpublished).
33. K. Chen, Y. Huang, Y. Deng, A.âB. Kuklov, N.âV. Prokof’ev, and B.âV. Svistunov, Phys. Rev. Lett. 110, 185701 (2013).
34. A.âB. Kuklov, M. Matsumoto, N.âV. Prokof’ev, B.âV. Svistunov, and M. Troyer, Phys. Rev. Lett. 101, 050405 (2008).
72328