Interacting double dimer model on the square lattice

We study phases and transitions of the square-lattice double dimer model, consisting of two coupled replicas of the classical dimer model. As on the cubic lattice, we find a thermal phase transition from the Coulomb phase, a disordered but correlated dimer liquid, to a phase where fluctuations of the two replicas are closely synchronized with one another. Surprisingly, and in contrast to the cubic case, the phase boundary includes the noninteracting point, as we establish using a symmetry-based analysis of an effective height theory, indicating that infinitesimal coupling is sufficient to synchronize the double dimer model. In addition, we observe a novel antisynchronized phase when the coupling between replicas is repulsive, and use Monte Carlo simulations to establish the full phase diagram, including (anti)synchronized columnar and staggered phases, with interactions between parallel dimers in each replica.


I. INTRODUCTION
Symmetry-breaking phase transitions can be described using an 'order parameter,' a local observable that transforms nontrivially under the symmetries, and which is strictly zero in the disordered phase but nonzero in the ordered phase [1]. This precise distinction occurs only in the thermodynamic limit, since for a finite system thermal or quantum fluctuations necessarily restore the symmetry.
Among phase transitions to which this description does not apply, a prominent example is the Berezinskii-Kosterlitz-Thouless (BKT) transition in the two-dimensional (2D) XY model [2][3][4]. Across the BKT transition, no symmetry is broken, as required for a 2D system at nonzero temperature by the Mermin-Wagner-Hohenberg theorem [5,6], and there is hence no local order parameter. One can instead understand the BKT transition as an example of a 'topological' phase transition [7], where the phases are distinguished by their topological properties. An appropriate criterion is the response to a twist applied across the boundaries of the system: the associated energy cost, referred to as the helicity modulus (or phase stiffness), decreases exponentially with system size above the transition, but is nonzero in the thermodynamic limit below it [8,9].
In this work, we investigate a superficially distinct type of topological order in the square-lattice double dimer model, illustrated in Fig. 1, consisting of two coupled replicas of the standard dimer model. The dimer model [10,11] is a paradigmatic example of a strongly correlated statistical system, in which the elementary degrees of freedom are dimers that occupy pairs of adjacent lattice sites subject to the constraint of close packing, i.e., that every site is occupied by exactly one dimer. This constraint, which can be interpreted as a Gauss law for an effective 'magnetic field' [12], has profound consequences for the properties of the system.
In particular, the local Gauss law implies that the effective flux (defined by the lattice analogue of the standard magnetic flux) through any closed surface vanishes, while the flux across FIG. 1. An example configuration of the double dimer model on the square lattice, in which two replicas of the close-packed dimer model (shown in black and white) are defined on the same lattice. According to Eq. (2), parallel pairs of nearest-neighbor dimers within each replica (marked with a star, ⋆) contribute + to the energy, and overlapping dimers contribute + to the energy. Hence, the energy of this configuration is = 15 + 9 .
This allows for the possibility of phases distinguished by the flux variance, which can be either suppressed exponentially or nonzero in the thermodynamic limit. The standard 'single' dimer model on the square lattice is indeed known to exhibit a transition between phases showing these two behaviors, and at which spontaneous breaking of lattice symmetries occurs simultaneously [13,14]. Such transitions, even when accompanied by symmetry breaking, do not fit within the standard Landau paradigm [1], and the transition in the dimer model is in fact known to belong to the BKT universality class, with the flux variance playing the role of the helicity modulus [13,14].
Our model contains interactions, of either sign, both between parallel dimers within each replica, , as studied previously in the single dimer model [13,14], and between dimers that coincide (or 'overlap') in the two replicas, . We show that it exhibits a particularly rich phase structure, including a number of phases with distinct types of order, both topological and symmetry-breaking.
As for the same model on the cubic lattice [15], we find a transition without any symmetry breaking, between a standard 'Coulomb' phase [12] and a 'synchronized' phase, where both replicas remain disordered but their relative fluctuations are suppressed. More precisely, the difference between the fluxes in the two replicas has variance that becomes suppressed exponentially with system size, while each flux separately does not.
In contrast to the three-dimensional (3D) case, however, we find that this transition occurs for infinitesimal coupling = 0 − between replicas, reflecting the critical nature of the noninteracting double dimer model [16]. We also find a novel 'antisynchronized' phase, where the overlap between replicas is minimized, which meets the Coulomb and synchronized phases at the zero-interaction point. Both of these features, along with the other transitions in the phase diagram, can be understood in terms of effective field theories based on 'height models' [17,18], which we derive based on symmetry.
For sufficiently strong repulsive (i.e., > 0) interactions within each replica, we confirm the presence in the single dimer model of a transition into a staggered phase, as noted in previous works [13,19,20], and determine the critical coupling at which it occurs. We also demonstrate the existence of phases in the double dimer model that are simultaneously staggered and (anti)synchronized.
Previous work on the double dimer model on the square lattice has addressed the noninteracting case [16], as well as models that correspond to the limit → +∞ [21] (see Sec. II B 3) and that include nonlocal interactions [22], motivated by a mapping from the quantum dimer model. Other related work has demonstrated the possibility of phase locking transitions in 2D superfluids [23] and the XY model [24].

Outline
In Sec. II we define the interacting double dimer model and present its phase diagram, which is calculated using the methods detailed in the subsequent sections: In Sec. III, we use symmetries to write down height field theories that describe the various phases and transitions, before presenting, in Sec. IV, the Monte Carlo results that underlie our phase diagram and establish the critical properties. We conclude in Sec. V.

II. MODEL
We consider a classical statistical model of dimers on two replicas of an × square lattice with periodic boundaries. To each link of each replica ∈ {1, 2}, we assign a dimer occupation number ( ) which takes values 0 or 1. The closepacking constraint applies separately for each replica and requires that at each site , where the sum is over links connected to .
To each configuration, we assign an energy where and are, respectively, interaction strengths between parallel dimers within each replica and between overlapping dimers in the two replicas (see Fig. 1), and ( ) ∥ counts the number of parallel pairs of nearest-neighbor dimers in replica . The partition function is given by = ∑ − ∕ , where the sum is over all close-packed dimer configurations in both replicas. (We set B = 1 throughout.)

A. Magnetic field and height picture
On the square lattice, it is useful to define a (fictitious) 'magnetic field' [12,25] ( ) on the link joining sites and + , where is a unit vector in direction ∈ { , } (and the lattice spacing is set to 1). Here, = (−1) + = ±1 depending on the sublattice and = 4 is the coordination number. (A similar construction applies to other bipartite lattices such as the honeycomb lattice.) The close-packing constraint for the dimers is then equivalent to the condition that the 'magnetic charge', given by the lattice divergence of ( ) , , is zero on every site. In two dimensions, this divergence-free constraint is resolved by defining a scalar 'height' ( ) on each plaquette, in terms of which where is the two-dimensional Levi-Civita symbol and Δ denotes the lattice derivative [17,18]. (This is the twodimensional analog of = × .) Together, Eqs. (3) and (5) define a one-to-one mapping between dimer configurations and their height representations, which is usually expressed as the following set of rules [14]: One first chooses a plaquette to be the zero of height. Then, moving anticlockwise around sites on sublattice A (B), the height increases (decreases) by 1−1∕ when an occupied bond is crossed. If, instead, an empty bond is crossed, the height decreases (increases) by 1∕ . Example height representations are shown in Fig. 2.
The flux ( ) for each replica can be defined by which, because of the divergence-free constraint, is equivalent to the sum of the magnetic fields on links crossing a surface normal to . The latter definition implies that Φ ( ) takes integer values and furthermore that it can be changed only by shifting dimers around a loop encircling the whole system [26].

B. Phase diagram
Our phase diagram for the square-lattice double dimer model, Eq. (2), is shown in Fig. 3. The fact that the Coulomb, synchronized and antisynchronized phases meet at = = 0 is determined solely from an RG analysis in Sec. III B; the remaining points are obtained numerically using a MC worm algorithm [15,27], as we describe in Sec. IV. In the remainder of this section, we define the phases appearing in Fig. 3.

Independent replicas
For = 0 the two replicas are independent and behave as single dimer models with interactions that favor ( < 0) or disfavor ( > 0) parallel dimers. For = 0, this model exhibits a Coulomb phase [12], where no symmetries are broken and the connected dimer-dimer correlation function ⟨ ′ ⟩ c decreases algebraically with separation. This phase extends to small nonzero , but gives way to ordered phases for sufficiently large | |∕ .
For negative , there is a transition to a phase with columnar order, as illustrated in Figs. 4(a) and (b), breaking translation and rotation symmetries [13,14]. An appropriate order parameter for this phase is the 'magnetization' which takes the values = ± in the four columnar states that maximize ∥ .
Besides the symmetry-breaking order parameter, the two phases are also distinguished by the probability distribution ( ) for the flux . In the thermodynamic limit, the Coulomb phase has ( ) ∝ − 2 | | 2 , where is a function of ∕ (see Appendix A 1). In the ordered phase, by contrast, ( ) is suppressed exponentially with system size for nonzero , since changing the flux requires shifting a row of dimers that spans the whole system, with energy cost proportional to . The mean square flux ⟨| | 2 ⟩ therefore changes its behavior across the transition, being independent of in the Coulomb phase but vanishing in the thermodynamic limit in the columnar phase [14].
We also make use of a third diagnostic of the transition, which is based on confinement of 'monomers', empty sites in the otherwise close-packed configuration. Removing one dimer leaves a pair of monomers on adjacent sites, with unit charges = ±1, which one can separate by locally rearranging the remaining dimers [12]. We define the monomer distribution function m ( + − − ) = m ( + , − )∕ , where m is the sum of Boltzmann weights of all configurations with a pair of monomers fixed at + and − .
In the Coulomb phase, m decreases algebraically with separation (see Appendix A 3), corresponding to a logarithmic effective potential m ( ) ≡ − ln m ( ) ∼ ln| | [26]. In the columnar phase, separating the monomers disturbs the ordered configuration, causing a linear potential m ( ) ∼ | |, and so m ( ) decreases exponentially with | |. The potential m therefore grows without limit in both phases; this is in contrast with the 3D case, where the potential is bounded in the Coulomb phase, and the monomers are said to be 'deconfined' [12]. The different asymptotic behaviors nonetheless allow the phases to be distinguished, and we refer to the 2D Coulomb phase as 'quasideconfined' by analogy with quasi-long-range order in the low-temperature phase of the XY model [4].
For large positive , the system instead reduces the number  [13,19,20]. We treat this phase and the transition in detail in Sec. IV C.
The Coulomb, columnar, and staggered phases of the single dimer model are shown in Fig. 3 on the vertical line ∕ = 0. Note that the Coulomb phase extends to < 0 along the line ∕ = 0.

Coupled replicas
For ≠ 0, the two replicas are coupled, with overlapping dimers favored for < 0 and disfavored for > 0.
The columnar and staggered phases at large | |∕ have order parameters, ( ) and ( ) respectively, in each replica . In the thermodynamic limit, any nonzero coupling fixes the relative values in the two replicas in order to maximize or minimize the overlap, as illustrated for the ground states in Figs. 4(a)-(d). We refer to the resulting phases as columnar/staggered & (anti)synchronized.
For smaller values of | |∕ , phase transitions occur that do not involve symmetry breaking, but can be characterized through the flux distribution and monomer confinement.
The flux distribution in the double dimer model can be described by the 2 × 2 covariance matrix ⟨ ( ) ⋅ ( ′ ) ⟩, but symmetry under replica exchange means we need only consider (2) are the total and relative fluxes corresponding to the fields (±) = (1) ± (2) . Note that the double dimer model can be viewed as a model of directed loops in (−) , formed from the overlap of the two replicas [15], and that (−) is the loop flux in this picture.
For < 0, one can postulate a phase where both replicas remain disordered and their flux variances ( ) remain nonzero, but where their fluctuations are strongly correlated so that the variance of the relative flux ⟨| (−) | 2 ⟩ vanishes in the thermodynamic limit [see Fig. 4(e)]. We have previously identified such a phase, which we call 'synchronized', in the double dimer model on the cubic lattice [15], and we demonstrate in the following that it also occurs on the square lattice. For > 0, we similarly identify an 'antisynchronized' phase, where fluctuations are correlated between the replicas in order to reduce the amount of overlap [see Fig. 4(f)]. The relative flux ⟨| (−) | 2 ⟩ also vanishes in the antisynchronized phase, as we argue in Sec. III B.
The monomer-confinement criterion can also be applied in the double dimer model, where we define m using a pair of monomers of opposite charge in the same replica, say = 1. Each monomer then has nonzero charge for (1) and hence for both (−) and (+) . They are therefore confined, with m ( ) decreasing exponentially with | | and hence ∼ 0 , in the (anti)synchronized phases, where fluctuations of (−) are suppressed.
To provide a signature of the synchronized phase, one can instead insert pairs of monomers in both replicas simultaneously. Two monomers, one in each replica, on the same lattice site form a double charge for (+) , but have zero net charge for (−) . We therefore expect such objects to be quasideconfined even when the relative flux variance is suppressed. Explicitly, we define the double monomer distribution function as d ( + − − ) = d ( + , − )∕ , where d is the sum of Boltzmann weights of all configurations with a pair of monomers fixed at + and − in each replica. Both in the (unsynchronized) Coulomb phase and in the (anti)synchronized phases, d ( ) decreases only algebraically with | |. (We show this directly in Appendix A 5 using an effective field theory.)

Infinite coupling between replicas
The point = 0, ∕ → +∞ corresponds to the dimer loop model [21], which is equivalent to a fully-packed loop model with fugacity = 2. The latter is known to be nonintegrable on the square lattice [28,29] but solvable on the honeycomb lattice, where it is equivalent to a three-coloring model [30].
In the opposite limit ∕ → −∞, the two replicas are perfectly aligned, and so act as a single dimer model with coupling 2 between parallel dimers. The values of ∕ at the columnar and staggered phase boundaries in this limit are therefore exactly half their values at = 0. (For ∕ → +∞, the critical couplings lie in between these two values, because the dimer loop model has higher entropy than the single dimer model.)

III. FIELD THEORIES AND CRITICAL PROPERTIES
Using the height mapping, the long-wavelength properties of the Coulomb phase can be described in terms of a free field theory. In this section, we use symmetry to determine the perturbations to this action that are most relevant under the renormalization group (RG), and hence establish the critical properties at each transition.

A. Single dimer model
To construct a continuum theory we replace the effective magnetic field , and height by coarse-grained fields ( ) and ℎ( ) obeying ( ) = ℎ( ). For the single dimer model, the Coulomb phase has action where is the stiffness, plus irrelevant higher-order terms. In the non-interacting limit ( = 0) the stiffness is from comparison of observables, for example Eqs. (A12) and (A22), with exact results obtained using Pfaffian methods [31,32].
To study the columnar-ordering transition in the single dimer model, we include additional terms in Eq. (8). We require that any action is local, and invariant under both ∕2 rotations and translation of dimers; as discussed in Ref. [14], this imposes constraints on the form of allowed additional terms, which are summarized in Table I. For example, translation of dimers by one lattice constant in the direction maps the height field ℎ( ) → −ℎ( − ) − 1 4 , so the action must be invariant under this change. The critical theory, which includes the most relevant term (in the RG sense) consistent with all requirements, is a sine-Gordon model: Note that if the symmetry of the single dimer model is reduced [33], such as in the case of anisotropic interaction strengths, i.e., ≠ , between parallel dimers [20], the form of the allowed cosine term is modified.
A standard perturbative RG calculation [14] applied to the general sine-Gordon action with an integer, leads to the following conclusions: There is a BKT phase transition at a critical value of the stiffness When < c the cosine term is irrelevant, i.e., it renormalizes to zero in the long distance theory, which is thus a free Coulomb phase. When > c it is relevant and locks the height field to discrete values.
In the case of the columnar-ordering transition where the action, Eq. (10), has = 4, we have In the columnar phase ( > c ) the cosine term locks the height field to values ℎ ∈ 1 8 , 3 8 , 5 8 , 7 8 , which correspond to the average values of the height in the four columnar ground states [see, for example, Fig. 2 [14].

B. Double dimer model
In the double dimer model, each replica has height field ℎ ( ) with identical stiffness , and replicas are coupled by the term (2) , with ∼ . The resulting action for the unsynchronized Coulomb phase may be written where ℎ (±) = ℎ (1) ± ℎ (2) and Note that in the non-interacting limit, i.e., = = 0, one has = 0, = ∞ and  (1) and ℎ (2) .) Here, ∕2 rotation symmetry refers to rotations about a plaquette center, and translation symmetry refers to translations by one lattice constant in the direction. For the SDM ' ' means SDM,col. , which describes the columnar-ordering transition [see Eq. (10)]. For the coupled DDM ' ' means DDM,sync. , which describes the synchronization transition [see Eq. (17)], or DDM,col. , which describes columnar ordering of coupled replicas [see Eq. (19)]. The SDM constraints on [ℎ] are discussed in detail in Ref. [14], and can be used to deduce the coupled DDM constraints on [ℎ (±) ].
We now construct field theories that describe phase transitions in the double dimer model. For independent replicas, rotation and translation symmetry constraints apply separately to both ℎ (1) and ℎ (2) , so each replica has an action given by Eq. (10). Therefore, when = 0, one expects a columnarordering transition with the same critical properties as the single-replica case.
For the double dimer model with non-zero coupling , we require an action local in both replicas, but now invariant under simultaneous ∕2 rotations, and translations, of both replicas. To study the synchronization transition, we focus on the relative height ℎ (−) [ℎ (+) remains non-critical] and include additional terms in Eq. (14). The constraints on allowed terms are easily derived using results for the single dimer model, and are included in Table I. For example, simultaneous translation of dimers by one lattice constant in both replicas maps the height fields ℎ ( ) ( ) → −ℎ ( ) ( − ) − 1 4 , so that the relative height ℎ (−) ( ) → −ℎ (−) ( − ), which must be a symmetry of the action. In this case, the critical theory is a sine-Gordon model with = 1: where, since the cosine term is forbidden by symmetry constraints when = 0, we require (−) ∼ to leading order. The constraints imposed by rotation and translation symmetry are identical (see Table I and ordering occurs when − > −,c . The ordered phase is synchronized (antisynchronized) in regions of the phase diagram with negative (positive) coupling , because the cosine term locks ℎ (−) = 0 ( 1 2 ) in order to minimize the action. (The relative height of any synchronized ground state is clearly ℎ (−) = 0.) We identify the (high-temperature) Coulomb phase in the double dimer model with the low-temperature phase of the XY model [34]. Hence, the synchronization transition is a BKT transition but with an inverted temperature axis.
To locate the phase boundary at fixed ∕ , we measure − as a function of ∕ using MC simulations and, from the crossing with its critical value −,c , identify a critical coupling ( ∕ ) c . However, in the absence of interactions within replicas ( = 0) MC simulations are not necessary, because −,c precisely coincides with the non-interacting limit ( = 0) of Eq. (16). Hence, in this case, the critical coupling ( ∕ ) c = 0, and replicas synchronize for infinitesimal < 0 [using Eq. (15) with ∼ ].
In our phase diagram, ℎ (−) is locked in the vicinity of columnar-ordering transitions when ≠ 0, and columnar ordering of coupled replicas is thus described by a critical theory in ℎ (+) . Adding to Eq. (14) the most relevant term consistent with the constraints on [ℎ (+) ] in Table I, which is a sine-Gordon model with = 2.
The critical stiffness for columnar-ordering of coupled replicas is therefore (In principle, ℎ (+) could lock before ℎ (−) if + > +,c while − < −,c , but we do not observe this.) The ordered phases, for which + > +,c , are columnar & (anti)synchronized. In the columnar & synchronized phase, for example, where ℎ (−) = 0, the cosine term locks the total height to values ℎ (+) = 2ℎ (1) ∈ 1 4 , 3 4 . This is consistent with average values of the height for a single dimer model in the columnar phase (see Sec. III A).

C. Honeycomb lattice
In passing, we consider the double dimer model defined on the honeycomb lattice, which is also bipartite and thus amenable to a height description. As we outline in the following, in the absence of interactions within replicas, i.e., = 0, synchronization on the honeycomb lattice occurs at a critical coupling ( ∕ ) c = 0, as for the square lattice.
The Coulomb phase action for the single dimer model on the honeycomb lattice is given by Eq.

to
= by exact calculations [26]. This is the same as for the square lattice, Eq. (9), and it follows that the double dimer model is again specified by Eqs. (14)-(16)  This finding may be interpreted in the context of a simple geometrically frustrated magnet, the triangular lattice Ising antiferromagnet (TLIAFM), which has Hamiltonian where ⟨ , ⟩ denotes nearest neighbor pairs of sites,  < 0 and = ±1. The TLIAFM has an extensive number of ground states and, as illustrated in Fig. 5, each ground state is in oneto-one correspondence with a close-packed dimer configuration on the honeycomb lattice [26].
In the limit  ∕ → −∞, the double dimer model on the honeycomb lattice is equivalent to a bilayer TLIAFM with Hamiltonian up to additive constants, where the four-spin interaction [35][36][37] derives from the term that counts overlapping dimers in Eq. (2). Hence, in this limit, spins in both replicas are either all aligned ( (1) = (2) ∀ ) or all antialigned ( (1) = − (2) ∀ ) for infinitesimal ∕ < 0, according to our height analysis. In fact, for general  , Eq. (22) is the Hamiltonian of the Ashkin-Teller model on the triangular lattice. The phase diagram of this model, obtained using MC simulations in Fig. 7 of Ref. [38], includes a BKT critical point at ( , ) = (−∞, 0) and is thus consistent with our conclusion.

IV. NUMERICAL RESULTS
In this section we use MC results, obtained using the worm algorithm [15,27], to map out the phase diagram shown in Fig. 3 and study the nature of each transition. There are three types of phase boundaries, which we consider in turn: synchronization, columnar ordering, and staggered ordering.

A. Synchronization transitions
In Sec. III B we identified, when ∕ = 0, a synchronization transition at infinitesimal coupling between replicas, i.e., ( ∕ ) c = 0. The transition is BKT type, where the Coulomb and synchronized phases correspond to the low-and hightemperature phases of the XY model, respectively. In this section, we first provide MC evidence to support this finding, and then describe how the phase boundary, which divides the Coulomb and (anti)synchronized phases, is located in the case ∕ ≠ 0. MC data for the synchronization transition when ∕ = 0 are shown in Fig. 6. According to theoretical arguments, the mean-square flux difference ⟨| (−) | 2 ⟩, shown in the topleft panel, is system-size independent in the Coulomb phase (see Appendix A 2) and decreases exponentially with in the synchronized phase. Hence, the extent of the -independent region in ⟨| (−) | 2 ⟩ provides a rough bound |( ∕ ) c | ≲ 0.2 on the critical coupling. In the thermodynamic limit, however, we expect that ( ∕ ) c scales to zero (see Sec. III B), while ⟨| (−) | 2 ⟩ jumps discontinuously to zero across the transition. The latter is typical of a BKT transition; for example, in the XY model there is a universal jump in the helicity modulus at the critical point [9,34].
As expected for a BKT transition, the synchronization transition [at ( ∕ ) c = 0] is not accompanied by a peak in the heat capacity per site , as shown in the top-right panel of Fig. 6. Instead, near the transition, theory predicts a non-divergent essential singularity, which is unobservable [34,39]. In the XY model, the main feature of the heat capacity per site is a broad peak, which is above the critical temperature and does not diverge with system size. We observe this in the synchronized phase of the double dimer model, i.e., when ∕ < 0, consistent with the correspondence between the phases in the two models.
As discussed in Sec. II B 2, the Coulomb and synchronized phases may be distinguished through the monomer confinement criterion. In the bottom-left panel of Fig. 6, we show the confinement length , defined by which is equivalent to the root-mean-square separation of the test monomers. In the Coulomb phase, where monomers are One also has ∼ 0 in the synchronized phase, where monomers are confined. In our MC data, the region with ∼ at small | ∕ | is thus a signature of a quasideconfined phase. The behavior for large | ∕ | is consistent with a confined phase or quasideconfined monomers with > 4; we have checked that m decays exponentially in this region (not shown), implying the former. Note that for quasideconfined monomers we observe saturation at 2 ∕ 2 ≃ 0.15, which is less than the value expected for fully-deconfined monomers 2 ∕ 2 ≈ 1∕6 [using the result ( 2 + 2)∕6 for the mean-square separation of free monomers hopping on an empty lattice].
In the synchronized phase, where the replicas become strongly correlated, fluctuations in the relative flux (−) are suppressed. However, both replicas remain disordered so fluctuations in the total flux (+) are large in both phases, as shown in the bottom-right panel of Fig. 6. In particular, at For general ∕ , we locate the phase boundary between the (anti)synchronized and Coulomb phases as follows. In the Coulomb phase, the mean-square total and relative flux are given by [14] ⟨| (±) where ± = 1 ± 2 , as derived in Appendix A 2 starting from the continuum theory of Eq. (14). In the MC simulations we measure both ⟨| (±) | 2 ⟩, and solve these equations numerically for the stiffnesses ± using the Newton-Raphson method [40]. As shown in Fig. 7 for ∕ = 0.2, the phase boundary is then located by scanning through ∕ until − crosses its critical value −,c = 2 [see Eq. (18)].
To accurately determine the critical coupling, we use quadratic fits to measure a crossing point ( ∕ ) × for each system size (bottom-left inset of Fig. 7). For a BKT transition, the 16  appropriate finite-size scaling form is [41,42] where and 0 are constants. From our fit for ∕ = 0.2 (top-right inset of Fig. 7), we obtain ( ∕ ) c = −0.293(3). Ten further critical points located in this way are shown in the phase diagram of Fig. 3; this includes transitions into the antisychronized phase at > 0 which, notably, all scale onto the line ∕ = 0.

B. Columnar-ordering transitions
Next, we consider transitions into all columnar-ordered phases. For the case of independent replicas, i.e., when = 0, columnar ordering at < 0 separates the columnar and Coulomb phases. This transition has been studied in detail by Alet et al. in Refs. [13,14] for the single dimer model, where the critical temperature is determined using an order parameter. We first review this approach.
The magnetization , defined in Eq. (7), breaks both translation and rotation symmetry in the columnar phase. Denoting by the number of dimers with orientation , a simpler choice of order parameter is the dimer rotation symmetry breaking a scalar that is sensitive only to rotation symmetry breaking. This is still sufficient to indicate a columnar-ordering transition: In the Coulomb phase, by symmetry one expects ⟨ ⟩ = ⟨ ⟩ so that ⟨ ⟩ is small. In the columnar phase, rotation symmetry is broken and all dimers are either horizontal or vertical. Hence, one expects ⟨ ⟩ = 1 (the total number of dimers is 2 ∕2). Alet et al. observe this behavior in Fig. 9 of Ref. [14].
The critical temperature may be determined accurately using the dimer rotation symmetry breaking Binder cumulant In the vicinity of the critical point, the th moment of the dimer rotation symmetry breaking has scaling form [43] ⟨ ⟩ ∼ ( ∕ ) , where, for a BKT transition, the correlation length diverges as [34] ∼ exp Here, and are unknown constants, is a universal function, and = ( − c )∕ c the reduced temperature. Hence, the Binder cumulant has zero scaling dimension, i.e., where is a new universal function, because Eq. (28) has equal powers of in both numerator and denominator. At the critical temperature = 0, the correlation length diverges and, to leading order, the Binder cumulant has no system size dependence. Hence, depending on the finite-size behavior either side of = 0, MC data for different system sizes may cross at the critical temperature. This is observed for in Fig. 11 of Ref. [14], from which Alet et al. report c = 0.65(1) when = −1, but not for the Binder cumulant of [14,44]. We now generalize this method to locate the phase boundary when ∕ ≠ 0, i.e., for coupled replicas. In this case, columnar-ordering transitions separate the (anti)synchronized phases from the columnar & (anti)synchronized phases. Since translation and rotation symmetry are broken in all columnarordered phases, we again expect a sharp drop in the mean dimer rotation symmetry breaking ⟨ ( ) ⟩ for both replicas, as well as a peak in the corresponding susceptibility in the vicinity of a transition. This is shown in Fig. 8 (left and middle panels) for the transition at ∕ = −0.2, between the columnar & synchronized and synchronized phases.
Columnar ordering in the limits = 0, studied in Refs. [13,14], and ∕ → −∞, equivalent to columnar ordering of a single dimer model with eff = 2 (see Sec. II B 2), is known to be a BKT transition with an inverted temperature axis. We expect the whole phase boundary to share the same critical properties as these points.
We now use the field theory and RG analysis of Sec. III to verify our results. In Fig. 9 (top panel), we measure the monomer distribution function m ( ) at the columnarordering transition for independent replicas ( = −1, = 0 and = c = 0.65), counting only monomers on the same row, i.e., = ( , 0). Each MC simulation can only construct m up to an arbitrary multiplicative constant, so we fix m (1, 0) = 1. The Coulomb phase monomer distribution function has asymptotic form  Fig. 31 of Ref. [14]).
In the case of coupled replicas, one requires the asymptotic form of m ( ) in the (anti)synchronized phases, which is less straightforward. Instead, it is convenient to consider the double monomer distribution function d ( ) (see Sec. II B 2) which, as derived in Appendix A 5, has asymptotic form [45] d ( , 0) ∼ −2 + ∕ .

C. Staggered-ordering transitions
We begin by describing the nature of the staggered phase in the single dimer model. The simplest staggered ground states (which contain no parallel pairs of dimers) have all dimers horizontal, such as Fig. 10(a), or vertical, such as Fig. 10(d). More complicated ground states are obtained by shifting dimers along diagonal loops, or 'staircases', that span the periodic boundaries. For example, Fig. 10(b) is a staggered ground This simple representation of the staggered ground states is specific to two dimensions, and cannot be generalized to the cubic lattice. There is only one ground state, shown in Fig. 10  |Φ | ]. Using this binomial distribution, one may calculate observables deep within the staggered phase. For example, the total number of ground states is which corresponds to a subextensive entropy log ≈ 2 log 2.
We also infer that, since the quantity |Φ |∕ is distributed around 1∕4 with standard deviation ∝ −1∕2 , the flux takes one of the four values ∕ = (± 1 4 , ± 1 4 ) in the thermodynamic limit, thus spontaneously breaking rotation and translation symmetries. At finite there are fluctuations out of these extremal states, but the symmetry-breaking transition remains.
For the double dimer model, the above discussion allows one to write down the partition function exactly, as a function of ∕ , in the limit ∕ → ∞. For example, for both replicas consider only ground states in the first quadrant of Fig. 10, i.e., 0 ≤ Φ , ≤ ∕2. A staircase can be covered by either horizontal or vertical dimers within each replica, and so, including a field that couples to the flux difference (−) , has four possible Boltzmann weights: exp(− ∕ ) when both replicas have the same orientation (i.e., both horizontal or vertical), and exp[± ⋅ (1, −1)] when both replicas have different orientations. Since, in total, there are ∕2 staircases, the contribution of these ground states to the partition function is Because ++ contains all configurations with maximal overlap, we expect that, for = 0 − , the full partition function asymptotically approaches ++ in the thermodynamic limit. By taking suitable derivatives with respect to , one finds that the flux difference (−) is distributed around with variance ≈ −| | ∕ . Hence, as illustrated in Fig. 3, in the staggered phase infinitesimal negative coupling is sufficient to synchronize the two replicas.
We now use MC results to examine transitions into all staggered-ordered phases. To begin, we focus on the case = 0, where staggered ordering at > 0 separates the staggered and Coulomb phases. One expects the same critical properties as for the single-replica case so, for simplicity, we consider a single dimer model with = +1 and vary the temperature.
By analogy with the columnar-ordering transitions (cf. Fig. 8), we use the staggered order parameter [46] to determine the critical temperature. At low temperatures, deep within the staggered phase, one has ⟨ ⟩ = 1 by definition of the ground state manifold, Eq. (36), whereas in the Coulomb phase ⟨ ⟩ is small because the flux distribution ( ) is peaked at = with width ∼ 0 [see Eq. (A11)]. Between these regimes, the sharp drop in ⟨ ⟩ and peak in the corresponding susceptibility shown in Fig. 11 (left and middle panels), are characteristic of a phase transition. In Fig. 11 (right panel), we obtain the critical temperature from the crossing point in the staggered order parameter Binder cumulant [46] Our estimate, c = 0.475 (5), is close to existing results c = 0.449(1) and c = 0.51 of Refs. [13,20], respectively [see also Ref. [19], which reports c = 0.72 (5)]. MC data for the staggered-ordering transition exhibits both first-order and continuous features: ⟨ ⟩ jumps discontinuously at the critical temperature, which is characteristic of a firstorder transition, but the heat capacity per site (not shown) does not diverge strongly with system size (i.e., not ∼ 2 ), which suggests a continuous transition. In the height picture, meanwhile, the transition occurs when the stiffness = 0 [20,46] and hence, from Eq. (A11), all flux sectors become equiprobable at the critical point. Similar behaviors are observed in the quantum spin-1 2 XXZ chain [20], and spin ice subjected to uniaxial pressure [47], whose phase transitions occur at multicritical points of order infinity [48].
To locate the full phase boundary our approach is straightforwardly extended to the case of coupled replicas, using crossing points in the Binder cumulant where = ( (1) , (2) ). Eleven such points are included in our phase diagram, Fig. 3, obtained for system sizes = 64 and = 96. We again infer the critical properties of the whole phase boundary from the limits = 0 and ∕ → −∞.

V. CONCLUSIONS
Our central result is the phase diagram of the classical double dimer model on the square lattice, shown in Fig. 3. As on the cubic lattice, we find a synchronization phase transition at which fluctuations between the two replicas become more strongly correlated, with signatures in the variance of the relative flux and in the monomer distribution function, but no symmetry breaking. The critical properties at this transition are of the BKT type, as expected for such a transition in 2D.
In addition, we find an antisynchronized phase, which was not observed on the cubic lattice, where overlaps between the two replicas are reduced. Our numerical results indicate that the phase boundary with the Coulomb phase runs along the line ∕ = 0 for positive ∕ (except possibly close to = 0, where the finite-size scaling becomes more difficult), as has previously been conjectured [21].
Remarkably, and in contrast with the 3D case, we find that these three phases meet at the noninteracting point = = 0, implying that an infinitesimal coupling between replicas is sufficient to drive the synchronization transition. This conclusion is supported both by our numerical results and by theoretical considerations based on a height field theory.
In forthcoming work [49] we will apply bosonization to the transfer-matrix solution of the dimer model [50]. This provides an alternative perspective on the fact that the synchronization transition is at infinitesimal coupling, because it can be understood as a pairing transition for fermions at zero temperature in 1D. It also allows one to predict the asymptotic form of the phase boundary exactly, based on perturbation theory in terms of the couplings.
While we have focused here on the case of the square lattice, we have also shown that the double dimer model on the honeycomb lattice similarly synchronizes for infinitesimal attractive coupling. Since this model is equivalent to a bilayer triangular lattice Ising antiferromagnet [26], this suggests a possible experimental realization. The synchronization transition would manifest as a highly unusual thermal phase transition between distinct correlated paramagnets, though establishing clear experimental signatures would likely be challenging. A natural extension of the system studied here would involve multiple replicas. With sufficiently strong coupling between adjacent pairs, this could be interpreted as a trajectory either of the classical dimer model imbued with dynamics or of a quantum dimer model in imaginary time. Alternatively, coupling one replica to others and taking the limit → 0 [34] provides a way to introduce a quenched disorder potential on the links of the single dimer model.

ACKNOWLEDGMENTS
The simulations used resources provided by the University of Nottingham High-Performance Computing Service. We are grateful to F. Alet and J. P. Garrahan for helpful discussions.
The probability of flux is obtained by integrating out all other Fourier modes with ≠ , so

DDM d ( ) in the (anti)synchronized phases
Finally, we calculate the double monomer distribution function in the DDM (anti)synchronized phases. In these phases, the cosine term in Eq. (17) is relevant and locks the relative height to values ℎ (−) = 0 ( 1 2 ). Hence, from the continuum version of Eq. (5), the corresponding magnetic field (−) = 0. For the total magnetic field, this implies (+) = 2 (1) , since the cosine term in Eq. (19) is irrelevant. In this case, Eq. (A13) reduces to which is the correct action for the (anti)synchronized phases.
In terms of this, the continuum version of the double monomer distribution function is where is the partition function in the close-packed case and (1) = (2) is the magnetic field in the presence of a pair of test monomers, i.e., ⋅ (1) = ( ). The derivation proceeds as in Appendix A 3 but with → 4 + , and the asymptotic behavior is d ( ) ∼ | | −2 + ∕ . (A27)