Synchronization transition in the double dimer model on the cubic lattice

TL;DR

This paper investigates a classical cubic-lattice double dimer model, revealing a continuous synchronization transition with strong correlations between replicas, characterized by a novel universality class, and explores the phase diagram with additional interactions.

Contribution

It introduces the concept of a synchronization transition in coupled dimer models, characterizes its universality class, and provides an exact Bethe lattice solution that captures the qualitative phase structure.

Findings

Identifies a continuous synchronization transition without symmetry breaking.

Classifies the transition within the 3D inverted-XY universality class.

Provides an exact Bethe lattice solution that reproduces the phase diagram qualitatively.

Abstract

We study the classical cubic-lattice double dimer model, consisting of two coupled replicas of the close-packed dimer model, using a combination of theoretical arguments and Monte Carlo simulations. Our results establish the presence of a 'synchronization transition' at a critical value of the coupling, where both replicas remain disordered but their fluctuations become strongly correlated. We show that this unconventional transition, which has neither external nor spontaneous symmetry breaking, is continuous and belongs to the 3D inverted-XY universality class. By adding aligning interactions for dimers within each replica, we map out the full phase diagram including the interplay between columnar ordering and synchronization. We also solve the coupled double dimer model exactly on the Bethe lattice and show that it correctly reproduces the qualitative phase structure, but with mean-field critical behavior.

Figures

Tables

Table 1. Table 1: Behavior in different phases of the double dimer model of: the variance of the flux difference, 𝚽 ( − ) = 𝚽 ( 1 ) − 𝚽 ( 2 ) superscript 𝚽 superscript 𝚽 1 superscript 𝚽 2 \bm{\Phi}^{(-)}=\bm{\Phi}^{(1)}-\bm{\Phi}^{(2)} ; the variance of the total flux, 𝚽 ( + ) = 𝚽 ( 1 ) + 𝚽 ( 2 ) superscript 𝚽 superscript 𝚽 1 superscript 𝚽 2 \bm{\Phi}^{(+)}=\bm{\Phi}^{(1)}+\bm{\Phi}^{(2)} ; and the confinement length, ξ 𝜉 \xi . ‘Small’ means Var ⁡ 𝚽 ( ± ) Var superscript 𝚽 plus-or-minus \operatorname{Var}\bm{\Phi}^{(\pm)} decreases exponentially with linear system size L 𝐿 L and ξ 𝜉 \xi is independent of L 𝐿 L , while ‘Large’ means Var ⁡ 𝚽 ( ± ) ∼ L similar-to Var superscript 𝚽 plus-or-minus 𝐿 \operatorname{Var}\bm{\Phi}^{(\pm)}\sim L and ξ ∼ L similar-to 𝜉 𝐿 \xi\sim L . In the columnar & (anti)synchronized phases, these observables behave as in the columnar phase.

	Coulomb	Columnar	Synchronized
$\mathrm{Var}{𝚽}^{(-)}$	Large	Small	Small
$\mathrm{Var}{𝚽}^{(+)}$	Large	Small	Large
$\xi$	Large	Small	Small

Equations

$$E=J\left[N_{\parallel }^{(1)}+N_{\parallel }^{(2)}\right]+K\sum _l{d}_l^{(1)}{d}_l^{(2)}\text{,}$$

$$M_{\mu }^{(\alpha )}=\frac{2}{L^3}\sum _{\mathbf{r}}(-1{)}^{r_{\mu }}{d}_{\mathbf{r},\mu }^{(\alpha )}\text{,}$$

$$n_{\mathbf{r}}^{(\alpha )}=\sum _{\mu }\left[d_{\mathbf{r},\mu }^{(\alpha )}+d_{\mathbf{r}-{\mathbf{\delta }}_{\mu },\mu }^{(\alpha )}\right]$$

$$Z=\sum _{\begin{array}{c}{c}^{(1)}\in {\mathfrak{C}}_0\\ c^{(2)}\in {\mathfrak{C}}_0\end{array}}{\mathrm{e}}^{-E/T}$$

$$B_{\mathbf{r},\mu }^{(\alpha )}={ϵ}_{\mathbf{r}}\left[d_{\mathbf{r},\mu }^{(\alpha )}-\frac{1}{q}\right]\text{,}$$

$$Q_{\mathbf{r}}^{(\alpha )}=\sum _{\mu }\left[B_{\mathbf{r},\mu }^{(\alpha )}-B_{\mathbf{r}-{\mathbf{\delta }}_{\mu },\mu }^{(\alpha )}\right]\text{,}$$

$${\Phi }_{\mu }^{(\alpha )}$$

$$=\frac{1}{L}\sum _{\mathbf{r}}{ϵ}_{\mathbf{r}}{d}_{\mathbf{r},\mu }^{(\alpha )}\text{,}$$

$$G_{\text{m}}({\mathbf{r}}_{+}-{\mathbf{r}}_{-})=\frac{1}{Z}\sum _{\begin{array}{c}{c}^{(1)}\in \mathfrak{C}({\mathbf{r}}_{+},{\mathbf{r}}_{-})\\ c^{(2)}\in {\mathfrak{C}}_0\end{array}}{\mathrm{e}}^{-E/T}\text{.}$$

$${\xi }^2=\frac{{\sum }_{\mathbf{R}}|\mathbf{R}{|}^2{G}_{\text{m}}(\mathbf{R})}{{\sum }_{\mathbf{R}}{G}_{\text{m}}(\mathbf{R})}\text{,}$$

$${\mathcal{L}}_{\text{SDM}}=\frac{1}{2}\kappa |\mathbf{B}{|}^2=\frac{1}{2}\kappa |\nabla \times \mathbf{A}{|}^2\text{,}$$

$${\mathcal{L}}_{\text{DDM}}=\frac{1}{2}{\kappa }_{+}|{\mathbf{B}}^{(+)}{|}^2+\frac{1}{2}{\kappa }_{-}|{\mathbf{B}}^{(-)}{|}^2\text{,}$$

$${\mathcal{L}}_{\text{SDM,crit.}}={\mathcal{L}}_{\text{SDM}}+|(\nabla -\mathrm{i}\mathbf{A})\mathbf{z}{|}^2+s|\mathbf{z}{|}^2+u|\mathbf{z}{|}^4\text{,}$$

$${\mathcal{L}}_{\text{DDM,crit.}}={\mathcal{L}}_{\text{DDM}}+|(\nabla -\mathrm{i}{\mathbf{A}}^{(-)})\phi {|}^2+s_{-}|\phi {|}^2+u_{-}|\phi {|}^4\text{,}$$

$$w\bm{\big{(}}(i,j)^{(\alpha)}\bm{\big{)}}P\bm{\big{(}}(i,j)^{(\alpha)}\rightarrow(j,k)^{(\alpha)}\bm{\big{)}}=\\ w\bm{\big{(}}(j,k)^{(\alpha)}\bm{\big{)}}P\bm{\big{(}}(j,k)^{(\alpha)}\rightarrow(i,j)^{(\alpha)}\bm{\big{)}}\,{\text{.}}$$

$$w\bm{\big{(}}(j,k)^{(\alpha)}\bm{\big{)}}\propto\exp\left[-\left(JN_{(j,k)}^{(\alpha)}+Kd_{(j,k)}^{(\bar{\alpha})}\right)/T\right],$$

$$P\bm{\big{(}}(i,j)^{(\alpha)}\rightarrow(j,k)^{(\alpha)}\bm{\big{)}}=\frac{w\bm{\big{(}}(j,k)^{(\alpha)}\bm{\big{)}}-\textrm{min}(\bm{w})\delta_{i,k}}{\sum_{k}w\bm{\big{(}}(j,k)^{(\alpha)}\bm{\big{)}}-\textrm{min}(\bm{w})},$$

$$w\bm{\big{(}}(j,k)^{(1)},(m,n)^{(2)}\bm{\big{)}}\\ =w\bm{\big{(}}(j,k)^{(1)}\bm{\big{)}}w\bm{\big{(}}(m,n)^{(2)}\bm{\big{)}}\exp(K\delta_{j,m}\delta_{k,n}/T),$$

$$\mathrm{Var}{\mathbf{\Phi }}^{(-)}\sim f_{\Phi }(L^{1/\nu }t)$$

$${\xi }^2/L^2\sim f_{\xi }(L^{1/\nu }t)\text{,}$$

$$\mathrm{Var}{\mathbf{\Phi }}^{(-)}\sim f_{\Phi }(L^{1/\nu }t)+u{L}^{-|y_u|}{\tilde{f}}_{\Phi }(L^{1/\nu }t)\text{,}$$

$$T_{\times }(L_1,L_2)-T_{\text{c}}\sim \frac{L_2^{-|y_u|}-L_1^{-|y_u|}}{L_1^{1/\nu }-L_2^{1/\nu }}\text{.}$$

$$T_{\times }(L_1,\rho L_1)-T_{\text{c}}\sim L_1^{-|y_u|-1/\nu }\text{,}$$

$${\frac{\partial }{\partial T}\mathrm{Var}{\mathbf{\Phi }}^{(-)}|}_{T=T_{\text{c}}}\sim L^{1/\nu },$$

$$c=c_0+A{L}^{\alpha /\nu },$$

$$G(L)=\frac{G_{\text{m}}({\mathbf{R}}_{\text{max}};L)}{G_{\text{m}}({\mathbf{R}}_{\text{min}};L)},$$

$$G(L)\sim L^{-2\left(d-\frac{\varphi }{\nu }\right)},$$

$${\sigma }_{12}=\{\begin{array}{rl}-1& \text{ObD, antiparallel}\\ 0& \text{ObD, perpendicular}\\ -0.2& \text{no ObD.}\end{array}$$

$$Z=\sum _{c\in {\mathfrak{C}}_0}1\text{.}$$

$$⟨d_0⟩=\frac{1}{Z}\sum _{c\in {\mathfrak{C}}_0}{d}_0\text{.}$$

$$Z$$

Full text

I Introduction

Confinement phase transitions bring together a number of important concepts in modern condensed matter physics. These include fractionalized excitations Balents (2010); Henley (2010), i.e., emergent objects that cannot be constructed from finite combinations of the elementary constituents; non-Landau phase transitions Senthil et al. (2004a, b), where a description in terms of a symmetry-breaking order parameter is not sufficient to describe the critical behavior; and topological order Castelnovo et al. (2012), loosely defined as structure that is detectable only through global measurements.

Such transitions are known to exist in classical dimer models Kasteleyn (1961); Temperley and Fisher (1961), in which the degrees of freedom are dimers, objects that occupy the bonds of a lattice subject to a close-packing constraint: Each site is touched by exactly one dimer. Introducing a pair of empty sites – defects in the constraint referred to as monomers – one can distinguish phases by the effective interactions induced between the monomers by the background dimer configurations. A confinement transition separates a phase where these interactions are confining, i.e., the potential grows without bound with increasing distance, from one where they are not, and so the monomers can be separated to infinity with finite free-energy cost Henley (2010).

Known examples of confinement transitions in classical dimer models can be divided into two types. The first, where confinement occurs simultaneously with spontaneous symmetry breaking, includes the well-studied case of the columnar-ordering transition in the dimer model on the cubic lattice Alet et al. (2006a) and the corresponding square-lattice transition Alet et al. (2005, 2006b). The second type is where symmetry is broken externally, and includes the Kasteleyn transition in the honeycomb-lattice dimer model Kasteleyn (1963) and its generalizations to three dimensions (3D) Bhattacharjee et al. (1983); Jaubert et al. (2008), as well as the ‘1GS’ variant of the cubic dimer model studied by Chen et al. Chen et al. (2009).

Here, we demonstrate the existence of a ‘pure’ confinement transition, with neither spontaneous nor external symmetry breaking, in the double dimer model Raghavan et al. (1997); Damle et al. (2012); Kenyon (2014), comprising two coupled replicas of the close-packed dimer model. In this paper we focus on the cubic lattice, and confirm using Monte Carlo (MC) simulations the presence of a synchronization transition at a finite ratio of temperature $T$ to the coupling between replicas. Subsequent work will address the case of the square lattice.

To understand the transition, consider overlaying any pair of dimer configurations and deleting all dimers that coincide. As illustrated in Fig. 1, the result is a gas of directed loops, referred to as the ‘transition graph’ and corresponding to the set of dimer rearrangements that take one configuration to the other Kasteleyn (1963); Krauth and Moessner (2003). A coupling that favors overlapping dimers then amounts to an energy cost per unit loop length. This makes possible a loop proliferation (or ‘condensation’) transition, between a phase at low $T$ with only sparse short loops and one at high $T$ with a finite density of boundary-spanning loops, as a result of competition between energy and entropy Banks et al. (1977); Dasgupta and Halperin (1981). (See Sec. II for a precise definition.) In terms of the dimers, the proliferation of loops amounts to a phase transition between replicas being mostly synchronized and mostly independent, which we will refer to as a ‘synchronization transition’.

An equivalent characterization is provided by the concept of confinement: Consider first removing a single dimer in one replica, leaving an adjacent pair of monomers, and then rearranging dimers to separate the monomers to arbitrary distance $R$. In terms of the overlap between the two replicas, this necessarily produces an open string joining the pair. The condition that loops proliferate is equivalent to the entropy of a long string overcoming its energy cost ($\propto R$) Savit (1980) and hence that such a monomer pair is deconfined Henley (2010). Note that a double pair of monomers, one in each replica, is instead deconfined in both phases, including in the limit of perfect synchronization between the layers.

As is apparent from the nature of the two phases, the synchronization transition on the cubic lattice is an unusual example of a classical transition in 3D with no local order parameter. We argue theoretically that the transition, if continuous, belongs in the inverted-XY universality class, and demonstrate using MC results that this is indeed the case. It can therefore be seen as an unusual example of a scalar Higgs transition, similar to the 1GS model Chen et al. (2009) and the helical-field transition in spin ice Powell (2012, 2013), but without any external symmetry breaking.

Previous studies of the double dimer model have focused on the square lattice, and particularly on the cases where the coupling between replicas is either infinite Raghavan et al. (1997) or vanishing Kenyon (2014). Damle et al. Damle et al. (2012) have showed that certain matrix elements in a quantum dimer model can be related to the partition function of an interacting double dimer model.

A related transition between low- and high-overlap phases has been observed in paired replicas of plaquette spin models Garrahan (2014); Turner et al. (2015). While there are many similarities between these transitions, it should be noted that the plaquette models have no analogue of the confinement criterion, and the boundary is a first-order line terminating at a critical endpoint. Recent work on the bilayer Kitaev model Tomishige et al. (2018) has also observed a transition from a spin liquid to a phase where spins on the two layers pair into singlets.

Outline

In Sec. II we define the double dimer model, including couplings between and within replicas, and present theoretical arguments for its phase structure and critical properties. We then describe, in Sec. III, the MC method that we use, which extends the standard worm algorithm. Our numerical results, including a detailed study of the critical properties of the synchronization transition, are presented Sec. IV. Finally, in Sec. V, we show that the double dimer model, including coupling between replicas, can be solved exactly on the Bethe lattice. We conclude in Sec. VI with a brief discussion of analogous synchronization transitions in other systems and of possible experimental realizations.

II Model

We consider the double dimer model, consisting of two replicas $\alpha\in\{1,2\}$ of the dimer model, both defined on an $L\times L\times L$ cubic lattice with periodic boundaries. Denoting by $d_{l}^{(\alpha)}$ the dimer occupation number (equal to zero or one) of replica $\alpha$ on bond $l$, the configuration energy can be written as

[TABLE]

where $N_{\parallel}^{(\alpha)}$ denotes the number of parallel pairs of nearest-neighbour dimers Alet et al. (2005, 2006b, 2006a) in replica $\alpha$. An example configuration is shown in Fig. 1.

When $J<0$, the first term favors parallel alignment of dimers within each replica. It is minimized by configurations with columnar order, as illustrated in Figs. 2(a)–(c), which break translation and rotation symmetries. A suitable order parameter is the ‘magnetization’ $\bm{M}^{(\alpha)}$, a vector for each replica $\alpha$ with components

[TABLE]

where $\bm{r},\mu$ denotes the bond between sites $\bm{r}$ and $\bm{r}+\bm{\delta}_{\mu}$ (with $\bm{\delta}_{\mu}$ a unit vector in direction $\mu\in\{x,y,z\}$). On the cubic lattice, there are six such configurations for each replica with $\bm{M}^{(\alpha)}=\pm\bm{\delta}_{\mu}$.

The second term in Eq. (1) couples the two replicas, by counting the number of bonds that are occupied in both. We mainly focus on the case $K<0$, where the effect is to favor overlapping dimers; a fully synchronized example, minimizing this term, is shown in Fig. 2(d).

The dimer configurations are subject to the close-packing constraint $n_{\bm{r}}^{(\alpha)}=1$, where

[TABLE]

is the number of dimers at site $\bm{r}$. This constraint is applied separately to each replica $\alpha$, and so the partition function

[TABLE]

sums, for each $\alpha$, over the set $\mathfrak{C}_{0}$ of all configurations that obey $n_{\bm{r}}^{(\alpha)}=1$ for all $\bm{r}$. (We set $k_{\text{B}}=1$.)

II.1 Loop picture

On a bipartite lattice, such as the cubic lattice, configurations may be represented by a ‘magnetic field’, defined on the bonds of the lattice Huse et al. (2003); Henley (2010). Specifically, bond $\bm{r},\mu$ is assigned the field

[TABLE]

where $\epsilon_{\bm{r}}=(-1)^{r_{x}+r_{y}+r_{z}}=\pm 1$ depending on the sublattice and $q$ is the coordination number.

The close-packing constraint for the dimers is equivalent to the condition that the (magnetic) ‘charge’, given by the lattice divergence of $B_{\bm{r},\mu}^{(\alpha)}$,

[TABLE]

is zero on every site. The normalization of $B_{\bm{r},\mu}^{(\alpha)}$ is chosen so that removing a dimer (and so breaking the close-packing constraint) leaves a pair of monomers on opposite sublattices with $Q_{\bm{r}}=\pm 1$.

When the two replicas are overlaid, the result can be interpreted as a set of directed loops. To see this, consider the relative magnetic field $B_{\bm{r},\mu}^{(-)}=B_{\bm{r},\mu}^{(1)}-B_{\bm{r},\mu}^{(2)}$, which takes values on each bond of $\pm 1$ or [math]. The former is interpreted as a loop element directed along $\pm\bm{\delta}_{\mu}$, while the latter means that the dimers overlap and is interpreted as the absence of a loop element. Since the relative field is clearly also divergenceless, these elements indeed form a set of closed loops. Note that swapping the two replicas switches the direction of each loop.

Given the magnetic field $B_{\bm{r},\mu}^{(\alpha)}$, it is possible to define a vector flux $\bm{\Phi}^{(\alpha)}$ by

[TABLE]

which is equivalent to the sum of the magnetic fields on links crossing a surface normal to $\bm{\delta}_{\mu}$. It is therefore invariant under local dimer rearrangements; in fact, with periodic boundary conditions (PBCs), changes in flux are only possible by shifting dimers around a loop encircling the whole system. For example, adding a loop in $B_{\bm{r},\mu}^{(-)}$ that spans the system once in direction $\mu$ increases the relative flux $\Phi_{\mu}^{(-)}=\Phi_{\mu}^{(1)}-\Phi_{\mu}^{(2)}$ by one.

II.2 Phase diagram

The phase diagram of the double dimer model, Eq. (1), obtained using the MC method detailed in Sec. III, is shown in Fig. 3. In this section we define the phases shown and explain how the phase structure can be understood theoretically. In Sec. IV we describe how the phase boundaries, as well as the critical properties at each, are determined in the simulations.

II.2.1 Independent replicas

We first review the phase structure for $K=0$, where the two replicas act as independent (single) dimer models. For $K=J=0$, the cubic dimer model exhibits a Coulomb phase Huse et al. (2003); Henley (2010), in which the dimers are disordered and their correlations take a dipolar form. A single phase transition at $(J/T)_{\text{c}}=-0.597$ separates this from a low-temperature phase with nonzero columnar order parameter $\langle\bm{M}\rangle\neq\bm{0}$ Alet et al. (2006a). The transition is apparently continuous with critical exponents compatible with the tricritical universality class. Theoretical arguments Powell and Chalker (2008); Charrier et al. (2008); Chen et al. (2009), however, suggest that the critical properties should be described by the so-called noncompact $\mathbb{CP}^{1}$ theory (see Sec. II.3), and additional interactions Charrier and Alet (2010); Papanikolaou and Betouras (2010) indeed modify the exponents to values compatible with this universality class Sreejith and Powell (2014).

Besides the order parameter $\langle\bm{M}\rangle$, the Coulomb and columnar phases can be distinguished either through the flux or through the effective interactions between monomers. Consider first the latter, which involves introducing a test pair of monomers with opposite charge into the background of close-packed dimers. In the case of the double dimer model, the pair can be inserted into either replica; we choose replica $1$, and define the monomer distribution function as

[TABLE]

The set $\mathfrak{C}(\bm{r}_{+},\bm{r}_{-})$ contains all configurations with charge distribution $Q_{\bm{r}}=\delta_{\bm{r},\bm{r}_{+}}-\delta_{\bm{r},\bm{r}_{-}}$, i.e., with monomers of opposite charge at sites $\bm{r}_{\pm}$ (e.g., holes on opposite sublattices). The effective interaction $U_{\text{m}}$ is defined by $G_{\text{m}}(\bm{R})=\mathrm{e}^{-U_{\text{m}}(\bm{R})/T}$, and can be interpreted as the free-energy cost associated with the imposed charge distribution. (Translation symmetry, assuming PBCs, implies $G_{\text{m}}$ and $U_{\text{m}}$ depend only on the displacement $\bm{R}=\bm{r}_{+}-\bm{r}_{-}$.)

For $K=0$, the replicas are independent, and so $G_{\text{m}}$ reduces to the monomer distribution function for the single dimer model. In the Coulomb phase for small $\lvert J\rvert/T$, $U_{\text{m}}$ is a Coulomb potential at large separation, $U_{\text{m}}(\bm{r}_{+}-\bm{r}_{-})\sim U_{\text{m}}(\infty)-\kappa/(4\pi\lvert\bm{r}_{+}-\bm{r}_{-}\rvert)$, where $\kappa$ (the ‘flux stiffness’) and $U_{\text{m}}(\infty)$ are finite (positive) constants Henley (2010). In the low-temperature phase, separating the monomers necessarily disrupts the columnar order along a string joining them [consider Figs. 2(a)–(c)], and so has a free-energy cost proportional to distance Henley (2010). This qualitative distinction, between a confining interaction (preventing separation to infinite distance) at low $T$ and deconfinement at high $T$, provides an alternative characterization of the phase transition.

In practice, it is convenient to use the confinement length $\xi$, defined by

[TABLE]

which represents the root-mean-square separation of the test monomers. In the Coulomb phase, $G_{\text{m}}(\bm{R})\rightarrow\mathrm{e}^{-U_{\text{m}}(\infty)/T}>0$ for large separation, and so $\xi\sim L$. In the columnar phase, by contrast, $U_{\text{m}}(\bm{R})$ grows without limit as $\lvert\bm{R}\rvert\rightarrow\infty$, $G_{\text{m}}(\bm{R})\rightarrow 0$, and so $\xi$ is an $L$-independent constant.

A related criterion for the phases can be expressed in terms of the flux $\bm{\Phi}^{(\alpha)}$. The mean flux vanishes by symmetry in both phases, while the variance, $\operatorname{Var}\bm{\Phi}^{(\alpha)}=\langle\lvert\bm{\Phi}^{(\alpha)}\rvert^{2}\rangle$, scales differently with system size in the two: In the Coulomb phase, flux fluctuations are large, $\operatorname{Var}\bm{\Phi}^{(\alpha)}\approx L/\kappa$ Huse et al. (2003). In the columnar phase, the variance is exponentially small in $L$, because shifting dimers along a loop spanning the system disturbs the columnar order and hence costs energy $E\sim JL$ 111To see the connection to the monomer-confinement criterion, imagine changing the flux by removing a dimer to create a monomer pair, winding one monomer around the periodic boundaries, and then recombining the pair Hermele et al. (2004).. Because the two replicas are independent, the variances of the total and relative flux $\bm{\Phi}^{(\pm)}=\bm{\Phi}^{(1)}\pm\bm{\Phi}^{(2)}$ are identical, and equally serve to distinguish the two phases.

II.2.2 Coupled replicas

Consider now $J=0$ and nonzero coupling $K<0$ between replicas. In the ensemble defined by Eq. (9), $B_{\bm{r},\mu}^{(2)}$ is divergenceless while $B_{\bm{r},\mu}^{(1)}$ has nonzero divergence at $\bm{r}_{+}$ and $\bm{r}_{-}$. This implies the presence of an open string in $B_{\bm{r},\mu}^{(-)}$ that runs between these two sites, and along which the two replicas differ. In the limit $K/T\rightarrow-\infty$, the string will take the shortest possible path, resulting in an energy proportional to separation and hence a confining effective interaction $U_{\text{m}}$. Comparing this limit with the case where $K=0$ (and $J=0$), it follows that there must be a confinement transition, a qualitative change in the large-separation form of $U_{\text{m}}$, between the two. In our results, shown in Fig. 3, we indeed find such a transition at a critical coupling $(K/T)_{\text{c}}=-1.400$.

At temperatures above this point, where the entropy of the open string overcomes its energy cost, closed loops of $B_{\bm{r},\mu}^{(-)}$ are also free-energetically favorable. As a result, loops spanning the system boundaries, which cost an energy $E\sim KL$ and are hence suppressed exponentially in $L$ at low temperatures, ‘proliferate’ at the transition. Because these loops change the relative flux $\bm{\Phi}^{(-)}$, the high-temperature phase has $\operatorname{Var}\bm{\Phi}^{(-)}\approx L/\kappa_{-}$, as at $K=0$ but with a modified flux stiffness $\kappa_{-}$. By contrast, the variance of the total flux, $\operatorname{Var}\bm{\Phi}^{(+)}$, is large ($\approx L/\kappa_{+}$) in both phases, because identical loops in both replicas costs zero energy [consider Fig. 2(d)]. The behavior of the flux variances in the different phases is summarized in Tab. 1.

While confinement and the flux variance thereby provide precise definitions of the phases, we also expect loop proliferation to reduce the overlap $\sum_{l}d_{l}^{(1)}d_{l}^{(2)}$ between the replicas. We therefore refer to the low-$T$, high-overlap phase as ‘synchronized’ and the high-$T$, low-overlap phase as ‘unsynchronized’, although the overlap is nonzero in both phases and so does not provide an order parameter in the strict sense.

It should be noted that the energy $\propto K$ associated with each element of a directed loop (or open string) is not the only contribution to its free-energy cost. In regions devoid of loops, overlapping dimers can be rearranged without changing the loop configuration. (For example, in the configuration of Fig. 1, flipping the parallel pair of overlapping dimers around the bottom-left plaquette in both replicas does not create a new directed loop). To leading order, this results in an entropy that scales with the number of overlapping dimers. Since the introduction of a directed loop reduces this number, and so the entropy, at finite $K/T$ we expect an additional free-energy cost per unit length of loops, which can be thought of as renormalizing $K/T$ towards more negative values 222The significance of this effect can be estimated by comparing with a simple approximation that neglects it and treats the loops as simply costing energy ${\frac{1}{2}\lvert K\rvert}$ per unit length. (A loop of length $\ell$ reduces the number of overlapping dimers by ${\frac{1}{2}\ell}$.) This model has a proliferation transition at ${T_{\text{c}}\simeq 0.33\lvert K\rvert}$ Dasgupta and Halperin (1981). In fact, our MC simulations give ${T_{\text{c}}=0.714\lvert K\rvert}$ (see Fig. 7) – the additional free-energy cost of loops due to the entropy of overlapping dimers, neglected in our approximation, means a higher temperature is needed for them to proliferate..

The arguments for the phase structure can be straightforwardly extended to include both $J/T$ and $K/T$. At large negative $J/T$ and $K=0$, both replicas are columnar ordered, but the relative orientations of the two magnetizations $\bm{M}^{(1)}$ and $\bm{M}^{(2)}$ are arbitrary. Infinitesimal negative $K/T$ is sufficient to split this degeneracy extensively and therefore to synchronize the two replicas, giving the columnar & synchronized phase, illustrated in Fig. 2(a), with $\langle\bm{M}^{(1)}\rangle=\langle\bm{M}^{(2)}\rangle$.

For positive $K/T$, any pair of columnar configurations with distinct magnetizations has zero overlap and hence minimal energy. There is an accidental degeneracy, between antiparallel ($\bm{M}^{(1)}=-\bm{M}^{(2)}$) and perpendicular ($\bm{M}^{(1)}\cdot\bm{M}^{(2)}=0$) magnetizations in the two replicas, illustrated in Figs. 2(b) and (c) respectively, which can be resolved by order by disorder Henley (1989); Chalker (2017). The elementary fluctuations, which involve flipping a single pair of parallel dimers around a plaquette, cost energy $+6\lvert J\rvert$ in both cases, but additionally may cost energy $+2K$ in the case of perpendicular magnetization. The free energy is therefore lower in the antiparallel arrangement, suggesting that this is selected by order by disorder. Our MC results (see Sec. IV.4) are indeed consistent with a phase where $\langle\bm{M}^{(1)}\rangle=-\langle\bm{M}^{(2)}\rangle$, which we refer to as columnar & antisynchronized.

Comparison with the single dimer model further allows some quantitative details of the phase boundaries to be inferred: The critical point separating the Coulomb and columnar phases for $K=0$ is clearly $(J/T)_{\text{c}}=-0.597$ as for the single dimer model. Similarly, when $K/T\rightarrow-\infty$, the two replicas are perfectly synchronized, and behave as a single dimer model with effective interaction $J_{\text{eff}}=2J$. The critical temperature for columnar ordering is therefore given by $\frac{1}{2}(J/T)_{\text{c}}=-0.299$ in this limit. As shown in Fig. 3, the critical value of $J/T$ closely approximates this limiting value already for $K/T=-2$.

II.3 Field theories and critical properties

A continuum description for the Coulomb phase in the single dimer model is given by replacing the effective magnetic field $B_{\bm{r},\mu}$ by a continuum vector field $\bm{B}$ Huse et al. (2003); Henley (2010). The latter is subject to the constraint $\bm{\nabla}\cdot\bm{B}=0$, inherited from the close-packing constraint on the dimers, and hence can be expressed as $\bm{B}=\bm{\nabla}\times\bm{A}$ in terms of the vector potential $\bm{A}$. The continuum (Euclidean) action density is then given by

[TABLE]

where $\kappa$ is the flux stiffness introduced in Sec. II.2.1, plus irrelevant higher-order terms.

In the double dimer model, one can similarly introduce a continuum magnetic field $\bm{B}^{(\alpha)}$ for each replica, with the same stiffness $\kappa$ for each. The coupling $K$ leads to a term $\lambda\bm{B}^{(1)}\cdot\bm{B}^{(2)}$, with $\lambda\sim K$, and so an effective action for the unsynchronized Coulomb phase can be written as

[TABLE]

where $\bm{B}^{(\pm)}=\bm{B}^{(1)}\pm\bm{B}^{(2)}$ and $\kappa_{\pm}=\frac{1}{2}(\kappa\pm\lambda)$. The synchronization transition, at which fluctuations of $\bm{B}^{(-)}$ are suppressed, occurs when $K<0$ and hence $\kappa_{-}>\kappa_{+}$.

Confinement transitions from the Coulomb phase, such as the synchronization and columnar-ordering transitions, can be described by introducing ‘matter’ fields to enforce the restriction to discrete values Chen et al. (2009); Powell (2011). Condensation of these fields then leads, by the Higgs mechanism, to suppression of magnetic-field fluctuations. The structure of the critical theory is determined by considering representations of the projective symmetry group (PSG) Wen (2002) under which the matter fields transform.

In the case of the columnar-ordering transition in the single dimer model, the critical theory is Powell and Chalker (2008); Charrier et al. (2008); Chen et al. (2009)

[TABLE]

where $s$ and $u$ are real parameters and $\bm{z}$ is a two-component complex vector (which is said to be ‘minimally coupled’ to $\bm{A}$). The PSG analysis shows that the field $\bm{z}$ transforms as a spinor under real-space rotations and allows one to express the magnetization as $\bm{M}\sim\bm{z}^{\dagger}\bm{\sigma}\bm{z}$. In this description, the ordering transition occurs when $s$ is reduced below its critical value and $\bm{z}$ condenses, giving a nonzero magnetization and also suppressing fluctuations of the magnetic field via the Higgs mechanism.

In the double dimer model, the matter field should couple identically to both replicas. We therefore expect the critical properties at the columnar-ordering transition in the double dimer model to be the same as in the single-replica case. While some theoretical aspects of this transition remain unresolved Sreejith et al. (2019), its properties have been well characterized numerically Alet et al. (2006a); Charrier and Alet (2010); Papanikolaou and Betouras (2010).

To describe the synchronization transition, one must similarly include a matter field $\varphi$ whose role is to restrict $B_{\bm{r},\mu}^{(-)}$ to $\pm 1$ or [math]. Because these values are integers, the PSG is trivial in this case 333In the notation of Ref. Chen et al. (2009), the background flux is zero, and so the static gauge configuration $\bar{\alpha}$ vanishes., and so the result is a scalar Higgs theory,

[TABLE]

where $\bm{B}^{(-)}=\bm{\nabla}\times\bm{A}^{(-)}$. We have not included a field coupling to $\bm{B}^{(+)}$, which would remain noncritical across the synchronization transition.

In 3D, the scalar Higgs theory has a continuous transition in the XY universality class but with an inverted temperature axis Dasgupta and Halperin (1981). (A more direct route to the same critical theory starts from the loop picture and uses the standard mapping from integer loops to the XY model Banks et al. (1977).) We therefore expect the synchronization transition to belong to the inverted-XY universality class.

III Worm algorithm

We obtain numerical results using the MC worm algorithm Sandvik and Moessner (2006), in which non-local loop updates are performed. With PBCs, loops can span the boundaries, allowing changes in flux.

III.1 Single loops

We begin by reviewing the standard implementation of the worm algorithm, using single loops. As illustrated in Fig. 4, the process is broken down into a series of local steps:

Choose a lattice site $i=i_{0}$ and a replica $\alpha$ at random. 2. 2.

In the current configuration of replica $\alpha$, site $i$ is connected by a dimer to a neighboring site $j$. Delete this dimer, denoted by $(i,j)^{(\alpha)}$. 3. 3.

Select a neighbour of $j$, called $k$, using a local detailed balance rule (described below), and insert a new dimer $(j,k)^{(\alpha)}$. 4. 4.

If $k=i_{0}$, close the loop. Otherwise return to step 2, using $i=k$.

Since all loops are performed without rejection, the worm algorithm is highly efficient.

The transition probability P\bm{\big{(}}(i,j)^{(\alpha)}\rightarrow(j,k)^{(\alpha)}\bm{\big{)}}, with which the site $k$ is selected in step 3, is determined as follows. The requirement for global detailed balance translates into the local detailed balance condition Sandvik and Moessner (2006)

[TABLE]

Here, w\bm{\big{(}}(j,k)^{(\alpha)}\bm{\big{)}} is the equilibrium probability of the configuration obtained on insertion of the dimer $(j,k)^{(\alpha)}$ in step 3. For the double dimer model, Eq. (1) implies

[TABLE]

where $N_{(j,k)}^{(\alpha)}\in\{0,1,2,3,4\}$ is the number of nearest-neighbour dimers parallel to $(j,k)^{(\alpha)}$ in the same replica, whilst $d_{(j,k)}^{(\bar{\alpha})}$ is the dimer occupation number of the other replica $\bar{\alpha}$ on bond $(j,k)$. A solution for the transition probabilities, chosen to reduce backtracks (where $k=i$), is then

[TABLE]

where $\bm{w}$ is a $6$-component vector containing elements w\bm{\big{(}}(j,k)^{(\alpha)}\bm{\big{)}} for all $k$ Sreejith and Powell (2014).

Step 2 produces configurations containing two test monomers: a stationary monomer at site $i_{0}$, and a moving monomer at site $j$ (see Fig. 4). Hence, single-loop updates may be used to construct the monomer distribution function $G_{\text{m}}(\bm{R})$, by tallying the monomer separation $\bm{R}$ after each step 2. Since the local detailed balance rule correctly samples only configurations produced by step 3, it is necessary to tally an amount 1/\sum_{k}w\bm{\big{(}}(j,k)^{(\alpha)}\bm{\big{)}}, rather than unity Alet et al. (2006b).

III.2 Double loops

As discussed in Sec. II.2.2, single loops are suppressed in the synchronized phase, whereas simultaneous loops in both replicas are not. Therefore, to avoid problems with ergodicity, it is necessary to perform double-loop updates.

Double loops are constructed as follows. In step 1, loops begin in both replicas from the same randomly chosen site $i=l=i_{0}$. Both loops perform step 2 as for a single-loop update, deleting dimers $(i,j)^{(1)}$ and $(l,m)^{(2)}$. Step 3 now corresponds to $36$ choices, with $6$ in each replica. The insertion of new dimers $(j,k)^{(1)}$ and $(m,n)^{(2)}$ is associated with a configuration probability

[TABLE]

where the factor $\exp(K\delta_{j,m}\delta_{k,n}/T)$ prevents double-counting of dimer overlap when the bonds $(j,k)$ and $(m,n)$ are identical. Eq. (17) is then used to obtain the transition probabilities, with $\bm{w}$ a $36$-component vector containing elements w\bm{\big{(}}(j,k)^{(1)},(m,n)^{(2)}\bm{\big{)}}. In step 4, the process terminates when both loops close simultaneously, i.e. $k=n=i_{0}$. Otherwise we return to step 2, using $i=k$ and $l=n$.

Double-loop updates are performed without rejection, but their efficiency is poor at higher temperatures. This is because the probability of simultaneous closure is small, and so updates are unnecessarily long. To reduce this problem, we use double loops only for large $\lvert K\rvert/T$. We also define a spring potential V\bm{\big{(}}(j,k)^{(1)},(m,n)^{(2)}\bm{\big{)}}=\frac{1}{2}k_{\text{s}}\lvert\bm{r}_{k}-\bm{r}_{n}\rvert^{2}, where $k_{\text{s}}$ is a spring constant and $\bm{r}_{k}$ denotes the position vector of site $k$. This is imposed by multiplying the equilibrium probabilities of Eq. (18) by $\exp(V)$, and favours the selection of sites $k$ and $n$ with smaller separation in step 3. The potential only modifies equilibrium probabilities of configurations during the double-loop construction, and so does not affect detailed balance for the close-packed dimer configurations.

III.3 Simulation parameters

We performed simulations on systems up to a maximum linear size $L=96$ with PBCs. Following equilibration, data points are typically obtained by averaging over $10^{6}$ MC sweeps, where a sweep is defined such that all lattice bonds are visited once on average. Statistical errors are estimated using a jackknife resampling method. A spring constant $k_{\text{s}}=2$ is used for double loops.

We have checked the MC data by comparing with exact results for the case $L=2$ and with the limits discussed in Sec. II.2.1.

IV Numerical results

We have used the worm algorithm to find the phase diagram shown in Fig. 3, and to determine the critical properties of each transition. In this section, we present our results for each of the phase boundaries in turn, and also for the nature of the ordered phases at large negative $J/T$.

IV.1 Synchronized $\longleftrightarrow$ Coulomb

We first focus on the synchronization transition, between the synchronized and Coulomb phases. In particular, we consider the case $J=0$, $K=-1$, and vary the temperature. In the vicinity of the critical point, double-loop updates are not required.

Our data for the flux difference variance $\operatorname{Var}\bm{\Phi}^{(-)}$ and normalized confinement length $\xi/L$ are shown in Fig. 5. Both quantities are small (large) at low (high) temperatures, indicating a phase transition between synchronized and Coulomb phases (see Sec. II.2.2). In particular, the high-temperature limit $\xi^{2}/L^{2}\simeq 0.25$ is observed, which closely matches the mean-square separation of $(L^{2}+2)/4$ for free monomers hopping on an empty lattice Chen et al. (2009). This is evidence for deconfined monomers in the Coulomb phase.

In contrast to $\operatorname{Var}\bm{\Phi}^{(-)}$, the variance of the total flux, $\operatorname{Var}\bm{\Phi}^{(+)}$, is large in both phases, as shown in Fig. 6. This confirms that the dimers in each replica remain disordered, even though relative fluctuations between the two replicas are suppressed. In fact, as $\lvert K\rvert/T$ increases and the replicas become more synchronized, $\operatorname{Var}\bm{\Phi}^{(+)}$ becomes larger. In the limit of perfect synchronization, $K/T=-\infty$, $\bm{\Phi}^{(1)}=\bm{\Phi}^{(2)}$ and so $\operatorname{Var}\bm{\Phi}^{(+)}=4\operatorname{Var}\bm{\Phi}^{(1)}$, double the value at $K=0$, where $\bm{\Phi}^{(1,2)}$ are independent and their variances add.

In order to classify the phase transition, we use finite-size scaling arguments Cardy (1996). At the transition of interest, both $\operatorname{Var}\bm{\Phi}^{(-)}$ and $\xi^{2}/L^{2}$ have zero scaling dimension Alet et al. (2006a); Powell (2013), and so for a continuous transition at critical temperature $T_{\text{c}}$, obey the scaling forms

[TABLE]

and

[TABLE]

where $t=(T-T_{\text{c}})/T_{\text{c}}$ is the reduced temperature, $\nu$ is the correlation-length exponent, and $f_{\Phi}$ and $f_{\xi}$ are universal functions. At the critical temperature $t=0$, Eqs. (19) and (20) become independent of system size, predicting a distinct crossing point in MC data at $T=T_{\text{c}}$. This is observed (see Fig. 5, insets), indicating that the transition is continuous.

In reality, we observe a weak dependence on system size at the critical point, which may be explained by corrections to scaling. Including the leading-order correction, Eq. (19) becomes

[TABLE]

where $u$ is a constant, $-\lvert y_{u}\rvert$ is the renormalization group (RG) eigenvalue of the leading irrelevant scaling operator, and $\tilde{f}_{\Phi}$ is a universal function. For two system sizes $L_{1}$ and $L_{2}$, this implies a crossing temperature $T_{\times}$ scaling as Binder (1981); Ferrenberg and Landau (1991)

[TABLE]

Fixing the ratio $\rho=L_{2}/L_{1}$ gives

[TABLE]

with an identical result applying to the $\xi^{2}/L^{2}$ crossing point. We determine $T_{\text{c}}$ by fitting to these expressions with $\rho=2$, using quadratic fits in the vicinity of the crossing point to measure $T_{\times}$ for each $L_{1}$ (insets of Fig. 5). From our results, shown in Fig. 7, we obtain consistent critical temperatures $T_{\text{c}}=0.71447(4)$ (flux) and $T_{\text{c}}=0.714444(2)$ (confinement length).

In order to determine the correlation length exponent $\nu$, we evaluate the temperature derivative of Eqs. (19) and (20) at the critical point. For $\operatorname{Var}\bm{\Phi}^{(-)}$ this gives

[TABLE]

and one finds the same result for $\xi^{2}/L^{2}$. The system size dependence of the slope at $T_{\text{c}}$ is extracted from quadratic fits. The results are shown in Fig. 8, and a fit to Eq. (24) yields consistent estimates $\nu=0.671(8)$ (flux) and $\nu=0.677(3)$ (confinement length). These values are compatible with the 3D XY universality class, for which $\nu_{\text{3DXY}}=0.6717(1)$ Campostrini et al. (2006).

Now equipped with estimates for $T_{\text{c}}$ and $\nu$, we replot the data of Fig. 5 against $L^{1/\nu}t$ in Fig. 9. Near the critical temperature, a good data collapse is obtained for all but the smallest system size. The curves, which represent the universal functions $f_{\Phi}$ and $f_{\xi}$, are consistent (up to normalization) with those in Fig. 6 of Ref. Chen et al. (2009).

As shown in Fig. 10 (left panel), a single peak in the heat capacity per site $c$ is observed at the transition temperature, indicating a single phase transition between the synchronized and Coulomb phases. To measure the specific heat exponent $\alpha$, we consider its scaling at the critical point,

[TABLE]

where $c_{0}$ represents the regular part, and $A$ is a constant. A fit to this form in Fig. 10 (bottom right panel) yields $\alpha/\nu=0.13(11)$, and using $\nu=0.677(3)$ (flux) gives a rough estimate $\alpha=0.09(7)$. In the 3D XY universality class, the corresponding value is $\alpha_{\text{3DXY}}=-0.0151(3)$ Campostrini et al. (2006). Our results satisfy hyperscaling $\alpha=2-d\nu$, where $d=3$ is the spatial dimension.

We next measure the crossover exponent $\phi$, which can be found by considering the monomer distribution function $G_{\text{m}}$ Powell (2013). Each MC simulation can only construct $G_{\text{m}}$ up to an arbitrary multiplicative constant, so we define the ratio

[TABLE]

where $\lvert\bm{R}_{\text{max}}\rvert\sim L$, $\lvert\bm{R}_{\text{min}}\rvert=1$, and the system size dependence of $G_{\text{m}}$ has been shown explicitly. At the critical point, this has scaling form Sreejith and Powell (2014)

[TABLE]

for sufficiently large systems. A fit to this form in Fig. 11 yields $\phi/\nu=2.4820(6)$, and using $\nu=0.677(3)$ (confinement length) gives $\phi=1.680(8)$. This value is compatible with the 3D XY universality class, for which $\phi_{\textrm{3DXY}}=d\nu_{\textrm{3DXY}}-\beta_{\textrm{3DXY}}=1.6665(3)$, using the exponents reported in Ref. Campostrini et al. (2006).

Finally, we consider the same phase boundary, between the synchronized and Coulomb phases, at points where $J\neq 0$. The critical point $(K/T)_{\text{c}}$ for each $J/T$, plotted in Fig. 3, has been obtained from the crossing point of $\operatorname{Var}\bm{\Phi}^{(-)}$ for system sizes $L=16$ and $L=24$. We expect that the universality class is the same for each point along the boundary, and have confirmed this for the points $J/T=-0.16$ and $J/T=-0.24$.

IV.2 Columnar & (Anti)synchronized $\longleftrightarrow$ Coulomb

As discussed in Sec. II.2.1, independent replicas ($K=0$) exhibit a continuous transition between the columnar and Coulomb phases Alet et al. (2006a). We now consider columnar ordering of coupled replicas ($K\neq 0$), i.e., the transition between the columnar & (anti)synchronized and Coulomb phases. Our results indicate that columnar ordering is driven first-order when replicas are coupled. (Certain other additional interactions have previously been shown to have this effect Charrier and Alet (2010); Papanikolaou and Betouras (2010).)

According to Eq. (20), a continuous (confinement) transition is characterized by a crossing point in $\xi^{2}/L^{2}$, at the critical temperature. We plot this quantity in Fig. 12, in the vicinity of a transition between the columnar & synchronized and Coulomb phases. A distinct crossing point is not observed [cf. Fig. 5 (insets)], and hence the transition is not continuous. Similar behaviour is obtained for $\operatorname{Var}\bm{\Phi}^{(\pm)}$.

Instead, the transition must be first-order. One thus expects a bimodal energy histogram in the vicinity of the critical point, which can be seen in Fig. 13 (red). The same behavior is also obtained for transitions between the columnar & antisynchronized and Coulomb phases (blue). In contrast, a single peak is observed for columnar ordering of independent replicas (green), as expected for a continuous transition.

Six points along this first-order phase boundary are included in the phase diagram of Fig. 3. These have been located by identifying peaks in the heat capacity per site, using system size $L=32$.

IV.3 Columnar & Synchronized $\longleftrightarrow$ Synchronized

We next consider the transition between the columnar & synchronized and synchronized phases. In the limit $K/T\rightarrow-\infty$, this phase boundary corresponds to columnar ordering of a single dimer model with $J_{\textrm{eff}}=2J$ (see Sec. II.2.2). This is known to be an (apparently) continuous transition in the tricritical universality class, and we expect the whole phase boundary to share the same critical properties as this point.

Since the flux difference variance $\operatorname{Var}\bm{\Phi}^{(-)}$ and confinement length $\xi$ are small in both (synchronized) phases, we locate the phase boundary using crossing points in the total flux variance $\operatorname{Var}\bm{\Phi}^{(+)}$, for system sizes $L=16$ and $L=24$. We have analyzed the points $K/T=-2.0$ and $K/T=-1.2$ in greater detail (not shown), and verified the expected critical properties.

IV.4 Columnar & (Anti)synchronized phases

Finally, we consider the different possible columnar-ordered phases at negative $J/T$ and both signs of $K/T$. To classify these, it is convenient to use the covariance of the replica magnetizations $\sigma_{12}=\langle\bm{M}^{(1)}\cdot\bm{M}^{(2)}\rangle$, which, deep within the columnar-ordered region, indicates the relative orientations of the magnetizations.

In the columnar phase at $K=0$, the two replicas are independent, and so $\sigma_{12}=\langle\bm{M}^{(1)}\rangle\cdot\langle\bm{M}^{(2)}\rangle=0$, since the mean magnetization vanishes by symmetry. Deep within the columnar & synchronized phase for $K<0$, the 6 ground states with $\bm{M}^{(1)}=\bm{M}^{(2)}=\pm\bm{\delta}_{\mu}$ [see Fig. 2(a)] dominate, giving $\sigma_{12}=1$.

For positive $K$, there are two sets of configurations that minimize the energy: 6 where the magnetizations are antiparallel, $\bm{M}^{(1)}=-\bm{M}^{(2)}=\pm\bm{\delta}_{\mu}$ [see Fig. 2(b)], and $6\times 4=24$ where they are perpendicular, $\bm{M}^{(1)}\cdot\bm{M}^{(2)}=0$ [see Fig. 2(c)]. Because the degeneracy between the two sets is accidental (i.e., not required by symmetry), it is liable to be resolved by order by disorder (ObD). There are, a priori, three possibilities: ObD favoring antiparallel magnetizations; ObD favoring perpendicular magnetizations; and no ObD, leaving all orientations equally likely. For large negative $J/T$, where columnar order is well established and so $\langle\lvert\bm{M}^{(\alpha)}\rvert\rangle\simeq 1$, these give limiting values of

[TABLE]

MC results for $\sigma_{12}$ are shown in Fig. 14. The expected behaviour is obtained in the columnar phase ($\sigma_{12}=0$ when $K/T=0$), and the columnar & synchronized phase ($\sigma_{12}\rightarrow 1$ for $K/T<0$). In the columnar & antisynchronized phase, the data appear to converge towards $\sigma_{12}=-1$ as $J/T$ becomes more negative. This indicates that ObD selects an arrangement with antiparallel magnetizations, in agreement with consideration of the elementary fluctuations, as in Sec. II.2.2.

V Bethe lattice

To gain further insight into the synchronization transition, we consider the double dimer model on the Bethe lattice, which, we will show, can be solved exactly. This provides an approximation to the model on the cubic lattice that is in the spirit of mean-field theory. In particular, we expect it to reproduce the qualitative behavior correctly, with a critical temperature that approximates the true value, but to fail to predict the critical properties.

We first consider a ‘Cayley tree’, illustrated in Fig. 15, a graph where each vertex has $q$ neighbors, except for those at the boundaries, and where there are no closed loops. To avoid contributions from the boundaries, which in the thermodynamic limit constitute a finite fraction of the vertices, we define the Bethe lattice as the part of the Cayley tree that is far away from all boundaries. A dimer model on the Bethe lattice has dimers occupying the bonds of the lattice (i.e., the edges of the graph).

Statistical mechanics problems with nearest-neighbour interactions are often exactly solvable on the Bethe lattice, since the absence of circuits allows one to formulate a recurrence relation for the partition function. This method has been used for the Ising model Baxter (1982), whilst a similar calculation has been performed for spin ice on the Husimi tree Jaubert et al. (2008); Jaubert (2009); Powell (2013). Here we apply it to the synchronization transition on the Bethe lattice.

V.1 Noninteracting dimers

To illustrate the method, we begin with a simpler calculation. Consider a single close-packed dimer model, with no interactions, on the Cayley tree. In this case the partition function is simply

[TABLE]

A quantity of interest is the mean dimer occupation number for the central bond, or root, of the Cayley tree, given by

[TABLE]

To begin, the partition function is written as

[TABLE]

In the first line, the quantity $Z_{n}(d_{0})$ is the ‘partial partition function’ of the left, or equivalently right, branch of the Cayley tree, when the root dimer occupation number is fixed to $d_{0}$. The index $n$ enumerates the branch depth. The same logic may be applied to Eq. (30), and results in

[TABLE]

A branch with root $d_{0}$ and depth $n$ consists of $(q-1)$ ‘subbranches’, rooted at $d_{1}$ and with depth $n-1$ (see Fig. 15). This observation allows the construction of recurrence relations which connect the partial partition functions of branches of depth $n$ and $n-1$. By allowing for all consistent configurations of the subbranches, while applying (at the roots) the constraint that each site should be covered by exactly one dimer, one finds

[TABLE]

It is convenient to introduce the variable

[TABLE]

for which Eqs. (34) and (35) imply

[TABLE]

Next we consider only sites on the Bethe lattice, deep within the Cayley tree, by taking the thermodynamic limit $n\rightarrow\infty$. Here, the solution is a fixed point satisfying $x_{n}=x_{n-1}=x$, so that Eq. (33) may be re-written

[TABLE]

whilst Eq. (37) becomes

[TABLE]

This has solution

[TABLE]

and substitution into Eq. (38) yields

[TABLE]

This is given, as expected, by the ratio of the number of dimers to the number of bonds.

In reality, the recurrence relation of Eq. (37) does not converge towards its fixed point in the thermodynamic limit, but instead oscillates indefinitely. To perform a more rigorous treatment, one can permit monomers with a small nonzero fugacity $z$, modifying Eq. (37) to

[TABLE]

This recurrence relation does converge in the thermodynamic limit, and Eq. (41) is easily retrieved by subsequently taking the limit $z\rightarrow 0$.

Using the same approach, one can also calculate the response to monomer insertion, which, as discussed in Sec. II.2, allows one to distinguish confined and deconfined phases. The distribution function $G_{\text{m}}$, defined in Eq. (9), involves a pair of monomers and cannot easily be calculated using the recurrence relation. Instead, we consider the corresponding quantity for a single monomer,

[TABLE]

where $Z_{\text{m}}=\sum_{c\in\mathfrak{C}(\bm{r}_{+})}\mathrm{e}^{-E/T}$ is the partition function with a monomer inserted at $\bm{r}_{+}$. This takes the form of an expectation value (specifically, of a monomer insertion operator Sreejith et al. (2019)) and we therefore refer to it as the ‘monomer expectation value’. While it vanishes due to the requirement of charge neutrality when periodic boundary conditions are applied, it can be nonzero with open boundary conditions, including on the Bethe lattice.

Suppose $\bm{r}_{+}$ is taken as the left side of the root $d_{0}$. Then the left (right) branch of the Cayley tree ‘sees’ an occupied (unoccupied) root, and the system has partition function $Z_{\text{m}}=Z_{n}(1)Z_{n}(0)$. The partition function without monomers $Z$ is again given by Eq. (32). From the definition of Eq. (36), and its solution in Eq. (40), one obtains

[TABLE]

The result is nonzero, indicating that an isolated monomer can occur with finite free-energy cost $\Delta F_{\text{m}}=-T\ln\Psi_{\text{m}}$. Monomers are therefore deconfined, as expected in the noninteracting dimer model.

V.2 Synchronization transition

Now consider the double dimer model of Eq. (1) on the Cayley tree. Since parallel pairs of dimers cannot be defined, we set $J=0$, leaving configuration energies

[TABLE]

and a partition function given by Eq. (4). The quantity of interest is the mean energy per site deep within the interior of the tree, which we take as its value on the root bond,

[TABLE]

assuming translational symmetry (at least on average). We therefore require the correlation function

[TABLE]

on the same bond. This may be calculated in analogy with Sec. V.1, although the algebra is more involved.

The partition function is written as

[TABLE]

where the reduced coupling $k=K/T$ has been introduced for convenience. In the first line, the quantity $Z_{n}\left(d_{0}^{(1)},d_{0}^{(2)}\right)$ is the ‘partial partition function’ of the left, or equivalently right, branch of the Cayley tree, when the root dimer occupation numbers are fixed to $d_{0}^{(1)}$ and $d_{0}^{(2)}$. Similarly, Eq. (47) may be written

[TABLE]

In order to construct recurrence relations, one must again allow for all possible configurations of the subbranches, while applying (at the roots) the constraint that each site should be covered by exactly one dimer in each replica. The results are

[TABLE]

Next we define the variables

[TABLE]

and take the thermodynamic limit. The solutions are again fixed points, and Eq. (50) may be rewritten

[TABLE]

whilst Eqs. (51)–(54) translate into a system of coupled, nonlinear equations given by

[TABLE]

The solutions to this system depend on the value of the reduced coupling $k$. For $k\leq k_{\text{c}}=-\log(q-1)$, there is a single solution

[TABLE]

whereas, for $k>k_{\text{c}}$, we find additionally

[TABLE]

which, as we have confirmed by a linear stability analysis, is the only stable solution. (Note that the critical value $k_{\text{c}}$ is negative – as expected, the transition occurs for attractive coupling $K<0$.)

Substitution of this result into Eq. (56), and then into Eq. (46), yields the final result for the mean energy per site of

[TABLE]

In line with $x_{n}$ in Sec. V.1, the recurrence relations for $u_{n}$, $v_{n}$ and $w_{n}$ derivable from Eqs. (51)–(54) may oscillate in the thermodynamic limit. Again, convergence is achieved by allowing monomers with fugacity $z$ and taking the limit $z\rightarrow 0$.

In Fig. 16 (top panel) we show the temperature dependence of the mean energy per site for a Bethe lattice with the same coordination number as the cubic lattice, $q=6$, and with $K=-1$. There is a second-order phase transition at $T_{\text{c}}=1/\log(5)\simeq 0.62$, characterized by a kink in the mean energy per site. Note that on the cubic lattice, our corresponding result (with $J=0$) is $T_{\text{c}}\simeq 0.71$. The low-temperature phase is always perfectly synchronized, since there is an energy $K$ for every dimer in a given replica. The high-temperature phase, which we identify with the Coulomb phase, is unsynchronized. In particular, when $k=0$, the mean energy per bond is $K/q^{2}$. This is sensible, because in this limit the replicas are independent, and from Eq. (41) the probability of double bond occupation is $1/q^{2}$.

To confirm our identification of the high-temperature solution with the (unsynchronized) Coulomb phase, we return to the monomer expectation value $\Psi_{\text{m}}$ defined in Eq. (43). In this case, we consider the partition function with a monomer in a single replica (again on the left side of the root $d_{0}$), which is

[TABLE]

while the partition function without monomers $Z$ is given by Eq. (49). From the definitions of Eq. (55), and their solution in Eqs. (60) and (61), one obtains

[TABLE]

The result is shown in Fig. 16 (bottom panel), using the same parameters as for the mean energy per site. In the low-temperature synchronized phase ($k<k_{\text{c}}<0$), $\Psi_{\text{m}}=0$ and the free-energy cost for an isolated monomer, $\Delta F_{\text{m}}=-T\ln\Psi_{\text{m}}$, is infinite, while in the high-temperature unsynchronized phase, $\Psi_{\text{m}}>0$ and $\Delta F_{\text{m}}$ is finite. This qualitative distinction, equivalent to the criterion based on $G_{\text{m}}$ introduced in Sec. II.2, implies that the synchronization transition on the Bethe lattice is a bona fide confinement transition.

While the model on the Bethe lattice with $q=6$ gives a reasonable approximation to the critical temperature on the cubic lattice, it does not reproduce the correct critical behavior. This is directly evident for the heat capacity, $\frac{\partial}{\partial T}\langle E\rangle$, which, according to Eq. (62), has a discontinuity at $T=T_{\text{c}}$, as expected for a mean-field theory. For the monomer expectation value, the duality mapping to the XY model Powell (2013) gives $\Psi_{\text{m}}\sim t^{\beta}$ for $t>0$, where $\beta$ is the magnetization exponent. From Eq. (64), we find the mean-field value $\beta=\frac{1}{2}$.

VI Conclusions

We have studied the phases of the double dimer model on the cubic lattice using a combination of theoretical arguments and MC simulations. Our results demonstrate the presence of a synchronization transition at a critical coupling between the replicas, which has no local order parameter but can be characterized through the confinement of monomers.

In a subsequent publication, we will address the same model on the square lattice; extensions to dimer models on other bipartite lattices are straightforward. By adapting the synchronization criterion introduced here, analogous transitions can also be expected in other systems consisting of two coupled replicas of a fractionalized phase. These include the Coulomb phase in ice models Henley (2010); Castelnovo et al. (2012), where a pair of monopoles in one replica would similarly become confined upon synchronization.

Experimental realizations of such transitions could be possible in various frustrated systems. In 2D, these include magnetic materials with a bilayer structure as well as nanomagnet arrays Nisoli et al. (2013), which have been used to simulate ice models with a variety of geometries, constructed in a double-layer configuration. A 3D synchronization transition could be possible between magnetic moments of two types, for example, in pyrochlore oxides with magnetic ions on both the A and B sites of the crystal structure Gardner et al. (2010). In these cases, one expects a thermodynamic phase transition (see, e.g., Fig. 10), but with no magnetic ordering. We leave the detailed study of possible experimental signatures to future work.

A natural extension of the double dimer model that we have treated here is to consider the case of multiple replicas $\alpha\in\{1,2,\dotsc,n\}$. With strong coupling between ‘adjacent’ replicas $\alpha$ and $\alpha+1$, this can be interpreted as the Suzuki–Trotter decomposition of the partition function for a quantum dimer model Moessner and Raman (2011). The double-loop algorithm introduced in Sec. III.2 could be extended to the case of multiple replicas, giving a method that is similar (at least in spirit) to the membrane algorithm Henry and Roscilde (2014) previously applied to quantum ice.

Acknowledgements.

The simulations used resources provided by the University of Nottingham High-Performance Computing Service. We are grateful to J. P. Garrahan for helpful discussions.