Phases of quantum dimers from ensembles of classical stochastic trajectories
Tom Oakes Stephen Powell Claudio Castelnovo Austen Lamacraft Juan P. Garrahan





























ifstarifnextchar[beqslabelbeqsnolabel
ifnextchar[beqlabelbeqnolabel
beqlabel[#1]
    beqnolabel
        beqslabel[#1]
            beqsnolabel
                ifstartom.oakes@nottingham.ac.uk
                School of Physics and Astronomy
                and 
                
                Centre for the Mathematics and Theoretical Physics of Quantum Non-Equilibrium Systems, 
                
                University of Nottingham, Nottingham NG7 2RD, United Kingdom
                
                TCM Group, Cavendish Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom
                School of Physics and Astronomy
                and 
                
                Centre for the Mathematics and Theoretical Physics of Quantum Non-Equilibrium Systems, 
                 
                University of Nottingham, Nottingham NG7 2RD, United Kingdom
                We study the connection between the phase behaviour of quantum dimers and the dynamics of classical stochastic dimers. At the so-called Rokhsar–Kivelson (RK) point a quantum dimer Hamiltonian is equivalent to the Markov generator of the dynamics of classical dimers. A less well understood fact is that away from the RK point the quantum–classical connection persists: in this case the Hamiltonian corresponds to a non-stochastic “tilted” operator that encodes the statistics of time-integrated observables of the classical stochastic problem. This implies a direct relation between the phase behaviour of  quantum dimers and properties of ensembles of stochastic trajectories of classical dimers. We make these ideas concrete by studying fully packed dimers on the square lattice. Using transition path sampling – supplemented by trajectory umbrella sampling – we obtain the large deviation statistics of dynamical activity in the classical problem, and show the correspondence between the phase behaviour of the classical and quantum systems. The transition at the RK point between quantum phases of distinct order corresponds, in the classical case, to a trajectory phase transition between active and inactive dynamical phases. Furthermore, from the structure of stochastic trajectories in the active dynamical phase we infer that the ground state of quantum dimers has columnar order to one side of the RK point. We discuss how these results relate to those from quantum Monte Carlo, and how our approach may generalise to other problems. 
                Phases of quantum dimers from ensembles of classical stochastic trajectories
                    Juan P. Garrahan
                    December 30, 2023
                ============================================================================
                
                
                
                § INTRODUCTION
                
                
                Dimer models are prototypical examples of systems where the degrees of freedom are subject to strong local constraints. Quantum dimer models (QDMs) were initially conceived by Anderson and collaborators <cit.> as representations of singlet pairings of quantum spins. The simplest QDMs <cit.> correspond to systems where  dimers are the basic degrees of freedom, they fully pack a lattice, and have a Hamiltonian with kinetic terms flipping groups of neighbouring dimers – for example pairs of parallel dimers on the square lattice – and potential terms counting the number of such flippable clusters. In spite of their apparent simplicity – and of how old these models are – we lack a full understanding the ground state phase behaviour of such QDMs <cit.>. 
                
                At zero temperature, the phase diagram of fully packed QDMs is controlled by the ratio v/ of 
                the energy per plaquette v to the flipping frequency  (see next section for definitions). The case v/=1, usually called the Rokhsar–Kivelson (RK) point, is of special importance. Here the ground state can be found exactly and is given by an equal superposition of all dimer configurations. The RK point often also delimits different ground state regimes. For example, for the square or honeycomb lattices, the ground state at the RK point is a spin liquid with algebraic decaying correlations (a “Coulomb phase”<cit.>), separating states with different kinds of order at either side of v/=1: for v/>1 this order is known to be “staggered”, extending all the way to v/ = ∞; for v/<1, in contrast, the nature of the order is less clear, the only certainty being that 
                for v/ = -∞ it is “columnar”<cit.>. The richness of this phase diagram highlights the complexity that can emerge from the simple ingredient of a constrained Hilbert space.
                
                Constraints can also play a significant role in classical stochastic many-body systems. For example, kinetically constrained models (KCMs)  <cit.> are simple lattice systems with local constraints in their dynamics that mimic steric restictions, and are one of the paradigms for the slow relaxation characteristic of glassy systems <cit.>. Fully packed classical dimer models (CDMs) also have rich dynamical behaviour <cit.>.  And even the simplest constraint of hard-core repulsion as in simple exclusion processes can give rise to interesting non-equilibrium dynamics <cit.>. In order to fully capture the properties of such complex dynamics it is necessary to study the statistical properties of their trajectories. The natural framework is that of dynamical large deviations (LDs) <cit.>, which provides a  statistical mechanics of trajectories<cit.> which is the dynamical analog of the standard equilibrium ensemble method for static configurations. In this approach, dynamical properties are classified according to time-integrated observables, and the long-time limit plays the role of the thermodynamic limit in the static case. In the long-time regime the statistics of dynamical observables is encoded in LD functions that play the role of dynamical entropies and free-energies. 
                
                In this paper we aim to connect the static quantum phase behaviour of QDMs with the dynamical large deviation properties of the stochastic dynamics of CDMs. This connection starts with the quantum dimers at the RK point, where the Hamiltonian of the quantum system is the same as (minus) the generator of continuous time Markov dynamics of the classical dimers <cit.>. Away from the RK point the Hamiltonian is no longer a stochastic generator for the classical dynamics but is instead a deformation thereof, which corresponds to the  tilted generator that encodes the statistics of the  dynamical activity (number of configuration changes in a trajectory) of the classical system. This connects the statistical properties of classical trajectories to the spectral properties of the quantum problem. In what follows we establish these connections for the case of dimers on the square lattice. 
                
                
                
                
                § MODEL AND DEFINITIONS
                
                
                In both our classical and quantum models, the elementary degrees of freedom are dimers, which occupy the links of a lattice. We consider an L × L square lattice with periodic boundary conditions and define d_μ() ∈{0,1} as the dimer occupation number on the link joining sites  and + _μ, where _μ is a lattice vector and μ∈{x,y}. The allowed configurations are those where the dimers are fully packed, i.e., where every site of the lattice is touched by precisely one dimer,
                ∑_μ[ d_μ() + d_μ( - _μ) ] = 1.
                
                
                 §.§ Quantum dimer model
                
                
                A complete basis for the Hilbert space of the QDM is given by all fully packed dimer configurations. We will use the RK Hamiltonian <cit.>, which can be written schematically as
                [QDMH]
                _̋,v = - ∑(|⟩⟨| + |⟩⟨|) + v∑(|⟩⟨| + |⟩⟨|),
                
                where the sums are over the plaquettes (squares) of the lattice. The first summation is the kinetic energy, where each term flips a pair of dimers around a plaquette with frequency . The second summation is the potential energy which counts the number of flippable plaquettes, each plaquette carrying an energy v. 
                
                Besides the total number of dimers, the Hamiltonian _̋,v in a lattice with periodic boundary conditions has a further conserved quantity, the  flux. The components of the flux vector are defined by
                Φ_μ = ∑_ϵ_ d_μ(),
                
                where ϵ_ = (-1)^r_x + r_y = ± 1 on the two sublattices. 
                The flux  is conserved by any local dynamics within the space of fully packed dimer configurations <cit.>. Since there are 1/2L^2 dimers and each contributes ± 1 to one component of , possible values of the latter satisfy |Φ_x | + |Φ_y |≤1/2L^2.
                
                
                
                
                Possible ground states of the square-lattice QDM can be divided into two kinds: dimer liquids, which are topologically ordered phases that break no symmetries <cit.>, and conventional ordered phases (“dimer solids”), which break lattice symmetries. The ordered states, illustrated in Dimer_phases, can be further divided according to their value of the flux, which
                vanishes in the columnar, Dimer_phases(a), plaquette and mixed phases Dimer_phases(b), and is maximized in the staggered phase Dimer_phases(c).
                
                A liquid phase of the dimer model is one that preserves the full symmetry of the lattice, and is characterized by topological order <cit.> and fractionalized monomer excitations <cit.>. According to a result of Polyakov <cit.>, a quantum U(1) liquid phase, i.e., a phase whose long-wavelength description is a U(1) gauge theory, is not possible in 2D. It can be shown, however, that the ground state at the RK point, v/=1 (and v,>0), can be found exactly, and is a U(1) quantum liquid <cit.>. (This fine-tuned liquid, existing only an isolated point in the phase structure, is consistent with Polyakov's argument.)
                
                To see this, start from an arbitrarily chosen dimer configuration c and construct the set _c of all configurations that can be reached from c by successive plaquette flips. By writing _̋, as a sum of projectors <cit.>,
                _̋, = ∑ (|⟩ - |⟩ )( ⟨| - ⟨|),
                
                it is straightforward to show that the equal-amplitude superposition ∑_c' ∈_c|c'⟩ is an eigenstate of _̋,, with eigenvalue zero, and so, by the Perron–Frobenius theorem, a ground state. Because plaquette flips cannot change the flux, there is such a zero-energy ground state in each flux sector.
                
                We will mostly be concerned with the RK state at zero flux, constructed by choosing for c any dimer configuration in this flux sector. This state clearly preserves all symmetries of the Hamiltonian, and is an example of a quantum dimer liquid.
                
                The staggered states are given by those configurations of the dimer model that have no flippable plaquettes. They are therefore trivially eigenstates of the Hamiltonian with zero energy (for any v/), and can be shown to be ground states for all v/>1<cit.>. Such states saturate the bound on the flux, with = {1/2L^2,0} in the configuration shown in Dimer_phases(c), and others related by symmetry. These states break rotation and translation symmetries of the lattice.
                
                For v/ < 1, the energy is instead minimized by a nontrivial superposition of dimer configurations, and so the ground state depends on v/. One therefore expects a discontinuity in the first derivative of the ground-state energy at v/=1, and hence a first-order quantum phase transition, according to the standard thermodynamic classification. (This will indeed by confirmed by our results, presented in SecResults, although the presence of the RK point leads to more subtle behavior on the side with v/ < 1.) Analytical arguments do not determine the ground state for v/ < 1, however, and numerical approaches must instead be used. The candidates are those with zero flux, illustrated in Dimer_phases(a) and (b).
                
                In the limit v → -∞, the ground state maximizes the number of flippable plaquettes Nf. The maximal value Nf = 1/2L^2 is achieved by the four configurations with dimers arranged in columns; one is shown in Dimer_phases(a), and the others are related by symmetry. For large negative v/t, the ground state is therefore expected to be an ordered state continuously connected to this limit (and so with the same symmetry), referred to as the columnar phase.
                
                In the same way, the plaquette and mixed phases are continuously connected to product states
                [EqProductState]
                |θ⟩ = ∏_p∈(|⟩cosθ + |⟩sinθ),
                
                where the product is over a set  of plaquettes tiling the lattice without overlapping, such as those labeled A in Dimer_phases. The plaquette state has θ = π/4, while the mixed state continuously interpolates between the plaquette and columnar (θ = 0) states. The plaquette and mixed states include resonances that reduce the kinetic energy, at the expense of decreasing the number of flippable plaquettes <cit.> and hence increasing the potential energy. They may therefore compete with the columnar state for v/ of order unity.
                
                An order parameter for these phases is provided by the  magnetization<cit.>,
                [EqMagnetizationDefinition]
                N_μ = ∑_ (-1)^r_μ d_μ(),
                
                which points along one of the square axes, ±_x or ±_y, in the columnar phase and along one of the four diagonals, ±_x ±_y, in the plaquette phase, while interpolating between the two in the mixed phase. It vanishes by symmetry at the RK point and also in the staggered states. Note that the magnetization is distinct from the order parameter M_X used by Banerjee et al. <cit.>, being naturally defined in terms of the dimers, rather than through height fields. The magnetization  is well established as an order parameter in the context of classical dimer models <cit.>.
                
                The phase diagram for v/ < 1 has been extensively studied using exact diagonalization (ED) <cit.> and quantum Monte Carlo (QMC) <cit.>, leading to a variety of contradictory conclusions. Sachdev performed exact diagonalization on lattices up to 6 × 6 and found that the columnar phase extends from v/ = -∞ all the way to the RK point (v/ = 1). Leung et al. extended the calculations to 8 × 8 and concluded that there is an intermediate phase consistent with plaquette order. The same conclusion, though with a different critical value for v/, was reached by Syljuåsen <cit.> using projector QMC methods. Ralko et al. <cit.>, combining QMC and ED, agreed with the presence of an intermediate phase, but argued that it showed mixed order. Finally, Banerjee et al. <cit.>, who used a height representation to access larger system sizes than previous MC studies, concluded that there is no intermediate phase, with the columnar phase extending as far as the RK point.
                
                
                
                 §.§ Classical dimer model
                
                
                The stochastic CDM that we consider is one with continuous-time Markov dynamics 
                within the same set of close-packed dimer configurations. The master equation for the evolution of the probability over configurations can be written in general as
                [ME]
                ∂_|P_⟩ = |P_⟩ ,
                
                where |P_⟩ is the probability vector, 
                [Pt]
                |P_⟩ = ∑_c P_(c) |c⟩ ,
                
                with {|c⟩} the configuration basis, and P_(c) the probability of configuration c at time . The general form of the generator (or master operator) is
                [W]
                 = ∑_c,c'≠ c w_c → c'|c'⟩⟨c| - ∑_c R_c |c⟩⟨c| .
                
                The positive terms are off-diagonal and encode the possible transitions c → c' and their rates w_c → c'. The negative terms are diagonal, with R_c the escape rate from configuration c, R_c = ∑_c'≠ c w_c → c'. The form (<ref>) guarantees probability conservation: the largest eigenvalue of  is zero and its left eigenvector is the  uniform (or “flat”) state: 
                
                    ⟨-|   = 0,
                       ⟨-|   = ∑_c ⟨c|.
                
                
                For the specific case of the CDM, a given plaquette, when flippable, flips according to a Poisson process with rate constant γ<cit.>. 
                The generator for the dynamics then reads
                [MOp]
                 = γ∑(|⟩⟨| + |⟩⟨|) - 
                γ∑(|⟩⟨| + |⟩⟨|) .
                
                The terms in the first summation correspond to transitions due to plaquette flips. All allowed transitions have the same rate γ. The second sum is over the escape rates and ensures conservation of probability. The generator (<ref>) is Hermitian (and hence bistochastic), which means that the uniform state is also the right eigenstate of the zero eigenvalue – and thus the stationary state of the dynamics, 
                
                    | ss⟩   = 0 ,    | ss⟩   =  N^-1|-⟩,
                
                where the normalisation is given by ⟨-| ss|=⟩ 1 ⇒ N = ⟨-|-|$⟩. As in the case of the Hamiltonian (<ref>), the classical dynamics generated byconserves the flux for periodic boundary conditions, and thus each flux sector is an irreducible partition of the dynamics. 
                
                At the RK point,v/=1, the classical generator (<ref>) and the Hamiltonian (<ref>) become identical up to a sign and an overall factor that is determined by the transition rate <cit.>, 
                [RK]
                  = -_̋γ,γ .
                
                It follows that the ground state of the QDM at the RK point coincides (up to normalisation) with the stationary-state probability of the CDM, i.e., an equal superposition of all dimer configurations: 
                [gs]
                |⟩_ RK =  N^-1/2|-⟩ .
                 
                This shows that it is possible to probe the ground state properties of quantum dimers at the RK point from the stationary state dynamics of classical dimers.
                
                
                
                
                 §.§ Trajectory ensembles and dynamical large deviations
                
                
                The dynamics generated by (<ref>) is realised in terms of stochastic trajectories. For a continuous time Markov chain such as we are considering, a trajectoryω_of overall time extensionis a sequence of configurations and jumps between them, 
                [traj]
                ω_ = (c_0, c__1, …, c__K),
                    
                0 < _1 < ⋯ < _K < ,
                
                where_i(i=1,…,K) indicate the times at which jumps between configurations occur. Between jumps the configuration remains unchanged, so that from the time of the last jump,_K, to the final time,, the configuration in trajectory (<ref>) would bec__K.
                
                The dynamics generates an  ensemble of trajectories, defined as the set of all possible trajectories (<ref>) and their probabilities to occur,π(ω_). The ensemble of trajectories contains the information about all possible time correlations, and thus encodes more information about the dynamics than the probability P(c,'). In particular, the latter is obtained fromπby summing over all stochastic trajectories (i.e., by  contraction),
                [cont]
                P(c,') = ∑_ω_π(ω_) δ(c_'(ω_) - c),
                
                where c_'(ω_) is the configuration at time ' in trajectory ω_.
                
                The properties of trajectories can be catalogued by trajectory observables. The simplest of these is the  dynamical activity<cit.> defined as the number of configuration changes in a trajectory, which we will denote by the symbolK̂when acting on trajectories. For example, for the trajectory in traj we haveK̂(ω_) = K, as there are a total ofKjumps in that trajectory.
                
                The activity is a trajectory order parameter – it is extensive in both system size and observation time – and is the natural quantifier of the overall “amount of motion”. It does this quantification of motion in a “structure-agnostic” way; that is, it does not assume any particular configurational property underlying fast or slow relaxation. In the course of a trajectory the activity simply increases by one unit every time the system changes its configuration, irrespective of the nature of those changes. It is particularly well suited for glassy systems where there is no obvious structural order parameter associated to glassiness; see e.g. Chandler2010.
                
                
                Associated with the ensemble of trajectories is a corresponding distribution for trajectory order parameters such as the activity, which has probability distribution
                [PK]
                P_(K) = ∑_ω_π(ω_) δ[K̂(ω_) - K] ∼^-φ(K/).
                
                The approximate equality is the large deviation (LD) form of the probability that is applicable at long times (as long as the correlation times of the dynamics remain finite) <cit.>. At these long times the statistics ofKare determined by the LD  rate functionφ(k)which can be thought of as an entropy density in the space of trajectories <cit.>. 
                
                Equivalent information to that found inP_(K)is contained in the moment generating function (MGF) <cit.>, 
                [Zs]
                Z_(s) 
                = ∑_K P_(K) ^-s K
                = ∑_ω_π(ω_) ^- s K̂(ω_)∼^θ(s) .
                
                This function generates the moments ofK, via⟨K^n|=⟩ (-1)^n ∂_s^n Z_(s)|_s=0. Just as for the probability, the MGF has a LD form at long times, given by the approximate equality in Zs. The functionθ(s)is the scaled cumulant generating function (SCGF; its derivatives evaluated ats = 0give the cumulants ofKscaled by time), and  plays the role of a free-energy density for trajectories <cit.>.
                
                One way to interpret Eqs. (<ref>) and (<ref>) is in terms of  conditioned and  biased trajectory ensembles <cit.>. The delta function in (<ref>) restricts the sum to trajectories which have total activityK. This is analogous to a microcanonical ensemble which restricts configurations to a fixed energy. In turn, in (<ref>) the sum is over all trajectories but the probabilities are exponentially biased (or  exponentially tilted). HereKis only indirectly controlled bys. This is analogous to a canonical ensemble of configurations controlled by an inverse temperature. As in the static case for large volume, the two trajectory ensembles are equivalent for long times. In particular, the rate function and the SCGF are related by a Legendre transform <cit.>,
                [LT]
                θ(s) = - min_k[ φ(k) + s k ] .
                
                
                 §.§ s-ensemble
                
                
                While both trajectory ensembles encode the same information about the dynamics, the “canonical” ensemble, Zs, is more practical to study. This is sometimes called the s-ensemble<cit.>. The probabilities of trajectories are exponentially tilted from the natural ones of the dynamics as
                [sens]
                π_s(ω_) = π(ω_) ^- s K̂(ω_)/Z_(s) . 
                
                We denote averages of trajectory observablesÂin this ensemble by⟨Â|_⟩s, which in terms of the averages of the original dynamical ensemble read
                [avs]
                ⟨Â|_⟩s = ⟨Â^- s K̂|⟩/Z_(s) . 
                
                The power of thes-ensemble is that it allows a full characterisation of the 
                dynamics beyond typical behaviour by tuningsaway froms=0. In particular, at long times, the analytic structure of the SCGF determines the phase structure of the dynamics. This is analogous to what occurs with the free-energy in static problems. 
                
                One reason that thes-ensemble is more tractable is that the dynamical partition sum (<ref>) can be written in “transfer matrix” form,
                [Zss]
                Z_(s) = ⟨-|^𝕎_s|i⟩,
                
                where the probability vector|i⟩represents the distribution from which the initial state is drawn. The operator𝕎_sis a deformation or tilt of the original dynamical generator that reads (for the case of tilting against the activity) <cit.>[Ws]
                _s = ^-s∑_c,c'≠ c w_c → c'|c'⟩⟨c| - ∑_c R_c |c⟩⟨c|.
                
                
                The equivalence between ZsZss can be proved in the following way. If we write the master operator as = ∑_μ_μ - , where_μdenotes the off-diagonal parts ofthat are responsible for all the possible transitionsμ, andis the diagonal part ofwith the escape rates, cf. W, the probability ofπ(ω_τ)of a trajectory such as traj can be written as 
                
                    π(ω_τ)
                    =
                    ⟨ c_τ_K |
                    e^- (t-t_K)  _μ_K e^- (t_K-t_K-1) ⋯_μ_1 e^- t_1 | c_0 ⟩ 
                    ,
                
                In Zs we have the exponentially tilted probability
                
                    π(ω_) ^- s K̂(ω_)   =
                    ⟨ c_τ_K |
                    e^- (t-t_K)  
                    e^-s_μ_K
                       × 
                    e^- (t_K-t_K-1) ⋯
                    e^-s_μ_1 e^- t_1 | c_0 ⟩ 
                    .
                
                Summing over all trajectories to obtain Zs then gives
                
                    Z_(s)    = ∑_K=0^∞ ∑_μ_1⋯∑_μ_K ∫_0^ d_1 
                    ⋯∫__K-1^ d_K
                    
                       ⟨ - |
                    e^- (-_K)  
                    e^-s_μ_K e^- (_K-_K-1) ⋯
                    e^-s_μ_1 e^- _1 |  i⟩ 
                    ,
                
                which is Zss expressed as a Dyson series for the tilted generator Ws.
                
                In contrast to, defined in W, the tilted operator _s does not define a probability conserving dynamics fors ≠ 0. In fact, its largest eigenvalue isθ(s), and thus the problem of computing the dynamical partition sum reduces to that of maximising (<ref>). The long-time average activity, obtained from the SCGF as 
                [Ks] 
                lim_→∞⟨K̂|_⟩s/ = - θ'(s) .
                
                serves as the dynamical order parameter that helps classify the dynamical phase behaviour, with associated susceptibility,
                [chi] 
                χ_s = lim_→∞⟨K̂^2|_⟩s - ⟨K̂|_⟩s^2/ = θ”(s) .
                
                
                 §.§ Connection to QDM
                
                
                For the specific case of the dimer model, the tilted operator reads 
                [MOps]
                _s = ^-sγ∑(|⟩⟨| + |⟩⟨|) 
                - γ∑(|⟩⟨| + |⟩⟨|) .
                
                We see from QDMH that this coincides with the general QDM Hamiltonian, if we identify = ^-sγandv = γ:
                [WsH]
                _s = - _̋^-sγ,γ.
                
                Changingsis equivalent to changing the ratiov/t, and so the properties of thes-ensemble of classical trajectories of the CDM are directly related to those of the QDM.
                
                To be specific, consider the ground-state expectation value of a quantum observable represented by an operator  that is diagonal in the basis of dimer configurations. To find this, we evaluate  in the configuration at the midpoint of each trajectory and average over trajectories with s weighting. The latter gives ⟨(/2) ⟩_s = ∑_c P_s(c,/2) ⟨c||c⟩, where P_s is given by cont with the replacement π→π_s. Using sens, and applying the same steps that led from Zs to Zss, one finds
                ⟨(2) ⟩_s = ⟨-|^/2𝕎_s^/2𝕎_s|i⟩⟨-|^𝕎_s|i⟩.
                
                The expectation value of  in |⟩, the ground state of _̋^-sγ,γ, is therefore given by the limit of large ,
                ⟨||⟩ = lim_→∞⟨(2) ⟩_s.
                
                Note that this is true for arbitrary initial distribution |i⟩, as long as the overlap ⟨|i|⟩ is nonzero.
                
                In what follows we exploit the relationship between the quantum and classical models, and connect the ground state phase diagram of the QDM to the properties of CDM trajectories explored numerically via path-sampling methods.
                
                
                
                § TRAJECTORY SAMPLING OF CLASSICAL DIMERS
                
                
                 §.§ Transition path sampling
                 
                
                The main difficulty in sampling s-ensemble trajectories is the usual one associated with calculating exponential averages <cit.>: the trajectories that are easy to generate with the normal dynamics at s=0, MOp, are not the relevant ones for the biased ensemble at s ≠ 0, sens, and the latter are exponentially rare in the original dynamics. Since the tilted operator (<ref>) is not a dynamical generator, there is no simple way to generate the relevant trajectories directly. 
                
                As in a static context (think for example of sampling the equilibrium of a spin system at finite temperature), this is resolved by importance sampling <cit.>. In the case of trajectories one such importance sampling scheme is Transition Path Sampling (TPS) <cit.>. TPS is a set of numerical techniques developed for generating rare trajectories that are infrequent enough that their observation through brute force simulations is unfeasible. Using original dynamics to generate trial trajectories, TPS performs a biased random walk through trajectory space towards the region of rare trajectories that exhibit the desired behavior. TPS is particularly appropriate for sampling dynamics with detailed balance, as in the case of the CDM, MOp.
                
                The basic idea behind TPS is similar to that of Markov chain Monte Carlo but applied to trajectories <cit.>. In its simplest form the procedure is as follows: (i) Given a trial trajectory, a new trajectory is proposed (as we describe below in detail). (ii) The proposed trajectory is accepted or rejected according to a Metropolis criterion. Since in our case we want to sample sens, the key quantity is the change in overall activity Δ K between the old trajectory and the new one. As in standard Metropolis, if Δ K < 0 the new trajectory is always accepted; if Δ K > 0, it is only accepted with probability ^-s Δ K. The procedure is repeated until the ensemble of trajectories thus generated converges to the s-ensemble.
                
                The non-trivial step is (i). Various methods for generating trajectories have been proposed <cit.> that guarantee ergodicity in trajectory space and are reasonably efficient. In particular, we use the  shifting method <cit.>, described in TPS cartoon: Given an initial trajectory, a new trajectory is proposed by choosing an arbitrary cut time, _ cut (chosen uniformly between 0 and ), keeping either the portion of the trajectory  after the cut, > _ cut, and shifting  back to time 0; or keeping the portion of the trajectory  before the cut, < _ cut, and shifting   forward to . These two options are chosen with equal probability.
                
                The remaining part of the old trajectory is discarded and has to be replaced by a new partial trajectory. In the case of a backward shift, the missing part is that from τ_ cut onwards; the new portion is obtained by shooting a new trajectory with the original dynamics from the configuration at _ cut (after the shift, i.e., the final configuration of the original trajectory) up to the final time .  In the case of a forward shift, the missing part is between time 0 and _ cut; to fill it one shoots a new trajectory with the original dynamics  forwards starting from the configuration at _ cut (after the shift, i.e., the initial configuration of the original trajectory) for a length _ cut and then inverts time. This is a valid procedure in the case of detailed balance dynamics <cit.>.
                
                For step (ii) we need the difference in activities between the trajectories. The current and proposed trajectory share the portion that has been shifted, and so the difference in activity comes only from the newly generated part. Since Δ K is a time-extensive quantity, acceptance will be suppressed exponentially. The fundamental limitation of this version of TPS for sampling long trajectories is that ergodicity in the dynamics implies that proposed trajectories diverge exponentially fast from their seed. This should be contrasted with sampling a spin model, for example, where new configurations can be proposed by flipping just a single spin, thus preventing the energy difference from growing with system size.
                
                The above means that, while simple TPS can sample the s-ensemble much more efficiently than brute force sampling (and indeed has been used successfully in this context before <cit.>), the exponential cost of sampling long times may render it impractical. An alternative to TPS is the  cloning method<cit.> adapted from quantum diffusion Monte Carlo, which however also suffers from a similar exponential cost <cit.>. Below we discuss how to parially overcome this limitation for TPS by exploiting umbrella sampling techniques to the trajectory context <cit.>. 
                
                
                
                
                
                 §.§ TPS and trajectory umbrella sampling
                
                
                The idea of umbrella sampling is the following <cit.>. We wish to compute exponential averages of the form
                [expav]
                ⟨^- s K̂|=⟩∑_ω_π(ω_) ^- s K̂(ω_),
                
                where π(ω_) is the probability at which trajectories are generated by the original dynamics. Consider now an alternative (and proper stochastic)  reference dynamics where the same trajectories are generated with probability π_ ref(ω_). We can rewrite expav as
                [expav2]
                ⟨^- s K̂|=⟩∑_ω_π_ ref(ω_) 
                π(ω_)/πref(ω_)^- s K̂(ω_) 
                = 
                ⟨^- s K̂|⟩ref,
                
                where 
                [R]
                (ω_) = π(ω_)/πref(ω_),
                
                This simply means that the average of the trajectory observable ^- s K̂ over the original dynamics is the same as the average of the trajectory observable ^- s K̂ over the reference dynamics, ⟨⋯|⟩ref. The “umbrella” compensates for the change of probability. 
                
                This can be useful in the following way. Given a reference dynamics, we would estimate (<ref>) by an empirical average over N_ sp sample trajectories,
                [sp]
                ⟨^- s K̂|=⟩⟨^- s K̂|_⟩ ref≈1/N_ sp∑_α=1^N_ sp(ω^α) ^- s K̂(ω^α) .
                
                The sampling error is given by the variance of the average squared of the empirical average,
                
                    ε^2 
                       = 
                     var_ ref(  
                    1/N_ sp∑_α=1^N_ sp(ω^α) ^- s K̂(ω^α))
                    /1/N_ sp∑_α=1^N_ sp(ω^α) ^- s K̂(ω^α)^2_ ref
                       =
                    1/N_ sp(
                    ⟨^2 ^-2 s K̂|_⟩ ref/⟨^- s K̂|_⟩ ref^2
                    
                    -1
                    )
                    
                       =
                    1/N_ sp(
                    ⟨^-2 s K̂|⟩/⟨^- s K̂|^⟩2
                    
                    -1
                    ) ,
                
                where we have used expav2 between the second and third lines to recast ε^2 in terms of the original averages. 
                
                Consider the case where the reference dynamics is just the original one, =1. The sampling error reads,
                [e22]
                ε^2 
                =
                1/N_ sp(
                ⟨^-2 s K̂|⟩/⟨^- s K̂|^⟩2
                
                -1
                ) 
                ≈^ [θ(2s)-2θ(s)]/N_ sp , 
                
                where in the last line we have used Zs for long times. The convexity of the SCGF function implies that θ(2s) ≥ 2θ(s) always, and the error diverges exponentially with time. This is why sampling with the original dynamics is inefficient, and accurate estimation requires exponentially many trajectories N_ sp. The aim is therefore to find an alternative reference dynamics which makes the convergence of (<ref>) more efficient. 
                
                
                
                
                  §.§.§ Ideal reference dynamics: generalised Doob transformation
                
                
                The reweighting factor (<ref>) for a trajectory such as (<ref>) reads,
                [RR]
                (ω_) = 
                ^-_1 Δ R_c_0w_c_0 → c__1/w_c_0 → c__1^ ref^-(_2-_1) Δ R_c__1⋯ ,
                
                where Δ R_c = R_c - R_c^ ref. It contains exponential factors for all the time periods between jumps and ratios of the transition probabilities for each jump in the trajectory. 
                The ideal choice for a reference dynamics would be one that cancels the exponential growth of the numerator in the first term of (<ref>). This is known to be given by the  generalised Doob transform<cit.>, which maps the tilted generator (<ref>) to a new stochastic operator  whose natural trajectories are those of the s-ensemble. For long observation times the transformation is obtained in the following way. 
                
                From the components l_c of the left eigenstate of _s, 
                [ls]
                ⟨l_s|_s = ⟨l_s|θ(s) ,
                
                we construct a diagonal matrix L_s, such that ⟨-| L_s = ⟨l_s|.  
                We then define 
                [Doob]
                 = L_s (_s - θ_s 𝕀) L_s^-1,
                
                where 𝕀 is the identity operator. _s is stochastic, 
                [Doob2]
                ⟨-| = 0 ,
                
                and for long times is guaranteed to generate the same trajectories as those of the s-ensemble. 
                
                The transition rates in  are given by 
                [WD]
                w̃_c→ c' = ( l_c'/l_c) ^-sw_c→ c' ,
                
                while the escape rates coincide, up to a shift, with the original ones,
                [RD]
                R̃_c = R_c - θ(s) .
                
                If the reference dynamics is the one generated by  the reweighing factor then reads,
                [RRD]
                 = ^τθ(s)^s K̂l_c_0/l_c_ .
                
                We see that this form of  cancels the exponential averaging in the numerator of (<ref>) and the error is no longer exponential in time. 
                This means that an s-tilted expectation value like (<ref>) can be computed by simply running the dynamics with <cit.>.
                
                
                
                
                  §.§.§ Effective reference dynamics
                
                
                While the ideal reference dynamics is provided by the Doob transformed generator , this is not a useful solution in practice, as one needs to diagonalise _s first, which amounts to solving the problem exactly. 
                Nevertheless, the form of the ideal transition rates (<ref>) can help guide the definition of convenient approximations for the reference dynamics that are practical <cit.>.
                
                We will consider the transition rates for the reference dynamics that have the form of (<ref>)[Wref]
                w_c→ c'^ ref = ( ℓ_c'/ℓ_c) 
                ^-sw_c → c' .
                
                The aim is to find a vector ⟨ℓ| that approximates the exact ⟨l_s| and is also tractable numerically. The associated escape rates, 
                [Rref]
                R_c^ ref = ^-s∑_c' ≠ c( ℓ_c'/ℓ_c) w_c → c' , 
                
                in general will not be a uniform shift from the original ones as in (<ref>).
                The reweighing factor is (up to boundary terms)
                [gref]
                 = e^s K̂ e^ - ∫_0^ d' Δ R_c(') .
                
                which together with (<ref>) indicates that, in contrast to the Doob-transformed dynamics, in general sampling will be exponentially difficult. Despite this, a judicious choice of ⟨ℓ| that reasonably approximates ⟨l_s| improves sampling significantly, as we now show.
                
                The exact l_c(<ref>) are functions of the whole configuration c and in general may have correlations at large distances. A simple approximation is to assume a short-range correlated form for ℓ_c and write them as a products of local factors. These local factors could in turn be of the exact form for small-enough local regions. We pursue this approach by considering 2×2 neighbourhoods with open boundary conditions, as described in app_open_2by2. 
                This leads to ℓ_c depending on the configuration only through the total number N_ f(c) of flippable plaquettes,
                [ellD]
                ℓ_c = e^D N_ f(c), 
                
                where the constant D parametrises the function.
                
                Putting this all together, the sampling proceeds as follows: The reference dynamics we use is given by the transition rates 
                [WrefD]
                w_c→ c'^ ref = e^D [N_ f(c')-N_ f(c)] 
                e^-sw_c → c' , 
                
                and escape rates 
                [RrefD]
                R_c^ ref = e^-s∑_c' ≠ c e^D [N_ f(c')-N_ f(c')]w_c → c' . 
                
                In order to sample the original dynamics tilted by e^-s K̂ we have to sample the reference dynamics tilted by e^-s K̂, see (<ref>), which we write as
                [g]
                ĝ = e^-D [N_ f(c_)-N_ f(c_0)]- ∫_0^ d' ( R_c(') - R_c(')^ ref).
                
                In order to account for this tilting we use TPS with the reference dynamics, with an acceptance rate for trajectories given by 
                [min]
                Γ_ acc(ω→ω') = min(1,ĝ(ω')/ĝ(ω)).
                
                
                  §.§.§ Optimization of reference dynamics
                
                
                The reference dynamics above is parametrised by the constant D. While the dynamics is based on the Doob transformation for the open 2 × 2 problem, there is no reason why the value of D that corresponds to the exact solution for the small system will provide the optimal dynamics for the large system. The reference dynamics can be optimised by choosing the value of D that maximises the trajectory acceptance rate Γ_ acc in the TPS simulations. 
                
                
                Finding the optimal value of D is a case of exploring the landscape of acceptance rates Γacc. Figure <ref> compares the acceptance rate as a function of s when using the original dynamics to that obtained when using the effective reference dynamics, Eqs. (<ref>–<ref>) with D_s obtained as in app_open_2by2, and with a value of D that optimises even further the acceptance. This latter optimal value of D is obtained for each s by starting from D=D_s and progressively changing D until a maximum of Γ_ acc is reached. The effective reference dynamics defined by this optimal D is then the one used for obtaining the corresponding s-ensemble.
                
                
                
                § RESULTS
                
                
                Figure <ref> shows the mean activity rate ⟨ k⟩ = ⟨ K ⟩/ evaluated across a range of s values, with the results for the smallest system size, L=6, compared with exact diagonalization. The agreement between these results confirms that the method has converged, at least for this L, and demonstrates that it is able to resolve detailed features of the dynamics for s>0. For larger system sizes, there is an increasingly sharp drop in the activity at s=0. The activity is related, by Ks, to the first derivative of the thermodynamic potential θ(s) (or equivalently, the quantum ground-state energy), and so this indicates a first-order transition at s=0, as expected from the analytical arguments in SecQDM. (Note that our simulations are restricted to the zero-flux sector, but this is not expected to change the thermodynamic properties in the limit L→∞.) The activity histogram, shown in the inset, shows the characteristic broadening, compared with a Gaussian distribution, expected for such a transition.
                
                As discussed in SecQDM, the ground state for positive s (v/>1) is fully staggered and so has maximal flux. Our simulations are, however, restricted to the zero-flux sector, and so the long-time dynamics is dictated by the ground state within this sector. For s>0, the activity decreases in a sequence of rounded steps, which are strongly system-size dependent. This is apparently a consequence of commensuration effects, as different low-flippability configurations are favored depending on the precise value of s. For large positive s, the system is mostly restricted to the minimally flippable zero-flux configurations, which, as illustrated in FigActivity(a), have precisely 4 flippable plaquettes for any system size <cit.>. Second-order perturbation theory gives the ground-state energy as -θ(s) = -2^-2s, which leads, using Ks, to a mean activity of ⟨ K ⟩ /  = 4^-2s in this limit.
                
                Our main results regarding the phase structure of the QDM are displayed in FigMagnetization. Panels (a–d) show the ground-state order-parameter distribution,
                
                p(') = ⟨|δ( - ')|⟩,
                
                where  is the operator defined in EqMagnetizationDefinition and the expectation value in the QDM ground state |⟩ is calculated as described in SecConnectionToQDM. For all negative s (v/ < 1), the maxima of the distribution p() occur for  aligned with the square axes, indicating that the ground state has columnar order. Particularly for small | s|, though, the selection of this order is weak, and the distribution has approximate (2) symmetry under continuous rotations, with the largest probabilities occurring on a ring of fixed ||.
                
                
                To characterize quantitatively the degree of selection of columnar order, we consider the quantities ⟨|||^2|⟩ and ⟨|cos (4ϕ)|⟩ where tanϕ = N_y/N_x. As panels (e–g) of FigMagnetization show, both of these quantities decrease as the RK point (s=0) is approached, in qualitative agreement with the results of Banerjee et al. <cit.>.
                
                The microscopic model is only symmetric under the discrete rotations of the lattice, and so the approximate (2) symmetry of the order parameter  is emergent. As argued by Fradkin et al. <cit.>, this can be understood qualitatively by considering the renormalization group (RG) flow structure. Properties near the RK point are described by an effective action
                [EqAction]
                 = 1/2(∂_0 h)^2 + 1/2ρ_2( h)^2 + 1/2ρ_4(^2 h)^2 + λcos(2π h) + …,
                
                written in terms of a continuum height field h, where ρ_2, ρ_4, and λ are real parameters, and  and ∂_0 denote the space and imaginary-time derivatives, respectively. The coefficient ρ_2 ∼ - v ∼-s is tuned through zero at the RK point, which, in spite of the first-order nature of the phase transition, corresponds to a critical fixed point of the height field theory.
                
                While the magnetization  does not appear explicitly in EqAction, it is related to the coarse-grained height by
                [EqMagnetizationAndHeight]
                N_x +  N_y ∼exp[- π/2(h + 1/2)].
                
                The λ term therefore breaks (2) down to the discrete subgroup of lattice symmetries, and determines the ultimate direction of the RG flow, towards columnar order for positive λ. At the RG fixed point corresponding to the RK point, however, λ is strongly irrelevant, with RG eigenvalue y_λ = -12. Standard scaling arguments in the presence of a dangerously irrelevant perturbation <cit.> therefore imply the existence of an additional length scale, ∝| s|^-3 for small negative s, with selection between columnar and plaquette order occurring only beyond this large scale <cit.>. This is qualitatively consistent with the weak columnar ordering observed for small | s | at the system sizes accessible in our MC simulations.
                
                The same effective action can be used to calculate the critical behavior of the the root-mean-square magnetization Nrms = ⟨|||^2|⟩^1/2, shown in FigMagnetization(g). We find in AppRMSN that Nrms∼ L^2| s|^1/2 for small negative s and large L, corresponding to a critical exponent β = 1/2. Finite-size corrections, however, cause deviations from this scaling form and saturation at Nrms∼ L√(ln L) when s = 0. At fixed L, crossover between these two forms, ∼ s^1/2 and ∼ s^0 (saturation), leads to a reduced effective exponent βeff. As shown in FigMagnetization(g), we find a reasonable fit to βeff = 0.254 for L = 16 and βeff≃ 0.1 using exact results for L = 6, consistent with such a scenario. Further results at much larger system sizes would likely be needed to confirm the expected scaling behavior, and we leave this to future work.
                
                
                Note that the unusual nature of the phase transition at s=0, which is thermodynamically first-order but shows critical behavior on the negative-s side (v/ < 1), is a common feature of constrained systems such as dimer models. A classic example is the Kasteleyn transition in 2D classical dimer models <cit.>.
                
                
                
                § CONCLUSIONS
                
                
                Our results here provide an example of the connection between the 
                statistical properties of long-time trajectories of a classical system and the properties of the low-lying spectrum of a related quantum system. Here we have focused on the classical fully packed dimer model on the square lattice and, correspondingly, the quantum dimer model. The connection works both ways as we have illustrated: from the known existence of a quantum phase transition in the QDM at the RK point we infer the existence of a transition – which we confirm numerically – between active and inactive dynamical phases in the CDM. Conversely, from the statistics of atypical trajectories of the CDM we learn about the ground state properties of the QDM away from the RK point. Other examples of this classical–quantum connection include classical exclusion processes and  XXZ chains <cit.>, and the one-dimensional Ising model with Glauber dynamics and the transverse field Ising chain <cit.>.
                
                For the QDM, our main result (see FigMagnetization) is that the ground state is the columnar phase for all v/ < 1. Our results in this regime also show an approximate emergent (2) symmetry of the order parameter. Both of these findings agree with the observations of Banerjee, and contradict earlier results <cit.>. We note that the method we use here is closer in spirit to that of the earlier work, based upon projector MC. For the CDM, the main result is the non-trivial structure of fluctuations in the dynamics away from typical behaviour. The first-order transition at s=0, see FigActivity(a), implies a coexistence in the equilibrium dynamics of space–time regions of high and low activity, and therefore a broad distribution of the dynamical order parameter, see FigActivity(b). Furthermore, the two competing phases display different kinds of structural order: while the inactive phase (s>0) is staggered, the active phase (s<0) is columnar (with both plaquette and mixed order being metastable due to their stability over shorter length scales). This change in the nature of configurations in order to optimise large dynamical fluctuations is reminiscent of what occurs in other systems, such as simple exclusion processes where – even in a state where the typical activity and current are featureless – rare inactive trajectories are associated with phase separated states and atypical large currents to hyperuniform (super-homogeneous) states <cit.>. 
                
                A consequence of the large fluctuations in the dynamics is that sampling rare trajectories is difficult. This is more so in a system like the CDM with periodic boundary conditions due to the constrained nature of configuration space and the conservation of the flux. To sample trajectory space we used transition path sampling [a Monte Carlo meta-dynamics in the space of trajectories guaranteed to converge to the s-ensemble sens], and to overcome the numerical difficulty of accessing exponentially suppressed trajectories we supplemented TPS with a version of  umbrella sampling in trajectory space <cit.>. TPS is well suited to our problem as the CDM dynamics obeys detailed balance (and is in fact bi-stochastic). Our umbrella sampling could be improved by obtaining the reference dynamics in an adaptive manner, as is done in <cit.> for cloning dynamics. Other interesting avenues to pursue include considering open boundary conditions (where we expect exploration of dynamics to be easier due to the absence of flux conservation), and to study in a similar manner as here dimer coverings in other geometries including higher dimensions.
                
                The method we have presented can easily be generalized to other geometries and other systems. All that is required is that the system of interest has an RK point, i.e., that for certain values of the parameters the Hamiltonian is equivalent to a stochastic generator <cit.>. If that is the case, properties of the ground state away from the RK point can be recovered from the rare fluctuations of the stochastic system, just as we have done here for the QDM away from v=t.
                
                
                
                 The simulations used resources provided by the University of Nottingham High-Performance Computing Service. This work was supported by EPSRC Grant Nos. EP/M019691/1 (SP), EP/P034616/1 (CC & AL), EP/K028960/1 (CC), and EP/M014266/1 (JPG).
                
                
                
                
                
                § L=2 CDM WITH OPEN BOUNDARIES
                
                
                We can obtain an approximation to the Doob-transformed dynamics (see SecTPSUmbrella) of a large system by focusing on the properties of a local region of size . This corresponds to a dimer model with four sites and  open boundary conditions: while a dimer is connected to each of the sites, these dimers may be directed outwards and so not contained within the  region. We can thus think of each site as occupied by a either a dimer or a monomer. As shown in FigOpenTwoByTwo, there are seven configurations in this open L=2 problem: two configurations with two dimers within the region, four configurations with one dimer and two monomers, and a single configuration with four monomers. The dynamics of the larger CDM, SecCDM, induces a dynamics between these seven states of the local  region.
                
                
                The dynamical generator in the reduced system has the form
                
                    ^  = [     -4b       b       b       b       b       0       0;       b  -(a+b)       0       0       0       a       0;       b       0  -(a+b)       0       0       0       a;       b       0       0  -(a+b)       0       0       a;       b       0       0       0  -(a+b)       a       0;       0       a       0       0       a -(2a+c)       c;       0       0       a       a       0       c -(2a+c);         ],
                
                where the components correspond to configurations with, respectively, four monomers (first row), a single dimer and two monomers (rows 2–5), and two dimers (rows 6–7). The above generator has three kinds of transitions: between the two-dimer configurations at rate c, between single- and double-dimer configurations at rate a, and between the no-dimer and single-dimer configurations at rate b. The former kind of transition corresponds to a plaquette flip within the  region, while the latter two are when the flip occurs at its boundary. The values of the rates depend on the size of the system on which the smaller region is embedded and can be obtained numerically from simulations. 
                
                As explained in the main text, we can deform ^ to obtain a SCGF for the number of flips from the largest eigenvalue of the deformed operator
                
                    ^_s  = [     -4b       b       b       b       b       0       0;       b  -(a+b)       0       0       0       a       0;       b       0  -(a+b)       0       0       0       a;       b       0       0  -(a+b)       0       0       a;       b       0       0       0  -(a+b)       a       0;       0       a       0       0       a -(2a+c)  e^-s c;       0       0       a       a       0  e^-s c -(2a+c);         ].
                
                We count only the transitions between the two-dimer configurations in the region, to avoid over-counting when we reconstruct the large system by overlaying  regions, as the other transitions correspond to flips in neighbouring regions.
                
                
                As an approximation to the Doob transform for the full L× L system, we replace the exact vector ⟨l_s| by the product of the left eigenvectors ⟨ℓ| of ^_s for each  region. The components ℓ_c can be expressed in the form ℓ_c = e^ε_c, where ε_c is a (dimensionless) “energy” associated to configuration c. (Both ℓ_c and ε_c depend on s; we suppress this for clarity.)
                
                We can characterize each configuration c by the number N_n of plaquettes with n dimers (i.e., the number in each class in FigOpenTwoByTwo); note that N_2 ≡ Nf(c), in the notation of SecEffectiveReferenceDynamics. The total number of plaquettes is N_0 + N_1 + N_2 = L^2, and, the fact that the total number of dimers is constrained to be 1/2L^2 implies that 2N_2 + N_1 = L^2. Together these give N_0 = N_2 = 1/2(L^2 - N_1) and allow us to express the dependence of the eigenvector on s (as well as on the rates a, b, and c) as ε_s = D_s N_2, using a single parameter D_s obtained from the diagonalisation of ^_s.
                
                This value of D_s specifies a dynamics that generates the exact s-ensemble for the open  problem and that we use as a starting point when optimizing the reference dynamics for the full lattice (see SecOptimizationOfReferenceDynamics).
                
                
                
                § ROOT-MEAN-SQUARE MAGNETIZATION NEAR RK POINT
                
                
                As argued in SecResults, close to the RK point and below the length scale for columnar ordering, one can set λ = 0 and drop higher-order terms in EqAction, leaving a quadratic action . Expressing the magnetization  in terms of the height h using EqMagnetizationAndHeight, we then have
                [EqRMSN]
                Nrms^2 ∼ L^2 ∫^2  ^-π^2/8 hrms^2(),
                
                where hrms^2() = ⟨|[h() - h(0)]^2 |⟩ is the mean-square height difference for positions separated by .
                
                By writing  in terms of the Fourier transform of h, one can express this ground-state expectation value as an integral over wavevector  and frequency ω. Integrating over ω and the angle between  and  gives
                [EqhrmsIntegral]
                hrms^2() = 1/2π∫_0^Λ k 1-J_0(k ||)/√(ρ_2 + ρ_4 k^2),
                
                where k = ||, Λ is an ultraviolet cutoff of order the inverse of the lattice spacing, and J_0 is a Bessel function.
                
                For ||≫Λ^-1, the integral in EqhrmsIntegral can be evaluated analytically, giving
                
                Nrms^2 ∼ L^2 [Ψ(L√(ρ_2/ρ_4)) - Ψ(Λ^-1√(ρ_2/ρ_4))],
                
                in terms of the function
                Ψ(x) = ∫_1^x/2 u   u ^2I_0(u)K_0(u),
                
                where I_0 and K_0 are modified Bessel functions of the first and second kind, respectively. The behavior for smaller || (i.e., of order the lattice spacing) is not well described by the continuum action, and, according to EqRMSN, will make a contribution of order L^2Λ^2 to Nrms^2.
                
                The results quoted in SecResults for the root-mean-square magnetization Nrms follow from the asymptotic behavior of the function Ψ. For large x, Ψ(x) ∼ x^2, and so in the large-L limit with fixed nonzero ρ_2 ∼| s|, Nrms∼ L^2| s|^1/2. (The magnetization is therefore extensive, as expected in an ordered phase.) For small x, Ψ(x) ∼|ln x|, and so Nrms∼ L√(ln L) at s = 0. In both cases, the contributions from lattice-scale effects are of lower order.
                
                99Anderson1973 P. W. Anderson, Resonating valence bonds: A new kind of insulator?Mater. Res. Bull.81531973.
                
                Fazekas1974 P. Fazekas and P. W. Anderson, On the ground state properties of the anisotropic triangular antiferromagnetPhilosophical Mag.304231974.
                
                RK1988 D. S. Rokhsar and S. A. Kivelson, Superconductivity and the quantum hard-core dimer gas6123761988.
                
                Moessner2011
                R. Moessner and K. S. Raman,  Quantum dimer models, in C. Lacroix, P. Mendels, and F. Mila (eds),  Introduction to Frustrated Magnetism, Springer Series in Solid-State Sciences, Vol. 164 (Springer, New York, 2011). [10.1146/10.1007/978-3-642-10589-0_17]. 
                
                Chalker J. T. Chalker,  Spin liquids and frustrated magnetism, in C. Chamon, M. O. Goerbig, R. Moessner, and L. F. Cugliandolo (eds),  Topological Aspects of Condensed Matter Physics, Lecture Notes of the Les Houches Summer School, Vol. 103, August 2014 (Oxford University Press, Oxford, 2017) [10.1093/acprof:oso/9780198785781.001.0001].
                
                HenleyReview C. L. Henley, The “Coulomb phase” in frustrated systemsAnnu. Rev. Cond. Matt. Phys.11792010 [10.1146/annurev-conmatphys-070909-104138].
                
                Ritort2003
                F. Ritort and P. Sollich,
                Glassy dynamics of kinetically constrained modelsAdv. Phys.522192003.
                
                Chandler2010
                D. Chandler and J. P. Garrahan,
                Dynamics on the Way to Forming Glass: Bubbles in Space-TimeAnnu. Rev. Phys. Chem.611912010.
                
                Henley1997 C. L. Henley, Relaxation time for a dimer covering with height representationJ. Stat. Phys.894831997 [10.1007/BF02765532].
                
                Oakes2016 T. Oakes, J. P. Garrahan, and S. Powell, Emergence of cooperative dynamics in fully packed classical dimers930321292016.
                
                Derrida2007
                B. Derrida,
                Non-equilibrium steady states: fluctuations and large deviations of the density and of the currentJ. Stat. Mech.P070232007.
                
                Touchette2009
                H. Touchette,
                The large deviation approach to statistical mechanicsPhys. Rep.47812009. 
                
                DeformationPaper1 J. P. Garrahan, R. L. Jack, V. Lecomte, E Pitard, K. van Duijvendijk, and F. van Wijland, Dynamical first-order phase transition in kinetically constrained models of glasses981957022007.
                
                DeformationPaper2 V. Lecomte, C. Appert-Rolland, and F. van Wijland, Thermodynamic formalism for systems with Markov dynamicsJ. Stat. Phys.127512007.
                
                DeformationPaper3 J. P. Garrahan, R. L. Jack, V. Lecomte, E Pitard, K. van Duijvendijk, and F. van Wijland, First-order dynamical phase transition in models of glasses: an approach based on ensembles of historiesJ. Phys. A420750072009.
                
                Henley2004 C. L. Henley, From classical to quantum dynamics at Rokhsar–Kivelson pointsJ. Phys.: Condens. Matter16S8912004 [10.1088/0953-8984/16/11/045].
                
                Laeuchli A, M, Läuchli, S. Capponi, and F. F. Assaad, Dynamical dimer correlations at bipartite and non-bipartite Rokhsar–Kivelson pointsJ. Stat. Mech.P010102008 [10.1088/1742-5468/2008/01/P01010].
                
                Balents2010 L. Balents, Spin liquids in frustrated magnetsNature4641992010.
                
                Polyakov1977 A. Polyakov, Quark confinement and topology of gauge theoriesNucl. Phys. B1204291977.
                
                Fradkin2004 E. Fradkin, D. A. Huse, R. Moessner, V. Oganesyan, and S. L. Sondhi, Bipartite Rokhsar–Kivelson points and Cantor deconfinement6922442004.
                
                FootnoteProductState The product state in EqProductState has ⟨θ| Nf|θ⟩ = 1/4L^2(1 + sin^4θ + cos^4θ).
                
                Sachdev1989 S. Sachdev, Spin-Peierls ground states of the quantum dimer model: A finite-size study4052041989.
                
                FootnoteOrderParameter The order parameter defined by Sachdev <cit.> is related to ours by Ψ_col^1 = N_x +  N_y.
                
                Banerjee D. Banerjee, M. Bögli, C. P. Hofmann, F.-J. Jiang, P. Widmer, and U.-J. Wiese, Interfaces, strings, and a soft mode in the square lattice quantum dimer model902451432014.
                
                Alet2006 F. Alet, Y. Ikhlef, J. L. Jacobsen, G. Misguich, and V. Pasquier, Classical dimers with aligning interactions on the square lattice740411242006.
                
                Leung1996 P. W. Leung, K. C. Chiu, and K. J. Runge, Columnar dimer and plaquette resonating-valence-bond orders in the quantum dimer model54129381996.
                
                Syljuasen2006 O. F. Syljuå sen, Plaquette phase of the square-lattice quantum dimer model: Quantum Monte Carlo calculations732451052006.
                
                Ralko2008 A. Ralko, D. Poilblanc, and R. Moessner, Generic Mixed Columnar-Plaquette Phases in Rokhsar-Kivelson Models1000372012008.
                
                Baiesi2009
                M. Baiesi and C. Maes and B. Wynants
                Fluctuations and response of nonequilibrium statesPhys. Rev. Lett.1030106022009.
                
                Chetrite2015
                R. Chetrite and H. Touchette, 
                Nonequilibrium Markov processes conditioned on large deviationsAnn. Henri Poincaré1620052015.
                
                Hedges2009 L. O. Hedges, R. L. Jack, J. P. Garrahan, and D. Chandler, Dynamic order-disorder in atomistic models of structural glass formersScience32313092009.
                
                Gingrich2015
                T. R. Gingrich and P. L. Geissler,
                Preserving correlations between trajectories for efficient path samplingJ. Chem. Phys.1422341042015.
                
                
                Chandler1987
                D. Chandler, 
                 Introduction to Modern Statistical Mechanics (Oxford University Press, 1987).
                
                TPSreview P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Transition path sampling: Throwing ropes over rough mountain passes, in the darkAnnu. Rev. Phys. Chem.532912002.
                
                Elmatad2010
                Y. S. Elmatad,  R. L. Jack,  D. Chandler, and J. P. Garrahan,
                Finite-temperature critical point of a glass transitionProc. Natl. Acad. Sci. USA107127932010.
                
                Speck2012
                T. Speck, A. Malins, and C. P. Royall,
                First-Order Phase Transition in a Model Glass Former: Coupling of Local Structure and DynamicsPhys. Rev. Lett.1091957032012.
                
                Giardina2011
                C. Giardina, J. Kurchan, V. Lecomte, and J. Tailleur,
                Simulating rare events in dynamical processesJ. Stat. Phys.1457872011.
                
                Ray2017
                U. Ray, G.K. Chan, and D.T. Limmer,
                
                Importance sampling large deviations in nonequilibrium steady states: Part 1 arXiv:1708.00459 (unpublished).
                
                Ray2017b
                U. Ray, G.K. Chan, and D.T. Limmer,
                 Exact fluctuations of nonequilibrium steady states from approximate auxiliary dynamics, arXiv:1708.09482 (unpublished).
                
                Nemoto2016
                T. Nemoto, F. Bouchet, R. L. Jack, and V. Lecomte,
                Population dynamics method with a multi-canonical feedback controlPhys. Rev. E930621232016.
                
                Klymko2017
                K. Klymko, P. L. Geissler, J.P. Garrahan, and S. Whitelam,
                 Rare behavior of growth processes via umbrella sampling of trajectories, arXiv:1707.00767 (unpublished).
                
                Jack2010
                R. L. Jack and P. Sollich,
                Large Deviations and Ensembles of Trajectories in Stochastic ModelsProg. Th. Phys. Supp.1843042010.
                
                Garrahan2016
                J.P. Garrahan,
                Classical stochastic dynamics and continuous matrix product states: gauge transformations, conditioned and driven processes and equivalence of trajectory ensemblesJ. Stat. Mech.0732082016.
                
                OakesThesis T. Oakes, unpublished.
                
                Senthil2004 T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P. A. Fisher, Quantum criticality beyond the Landau-Ginzburg-Wilson paradigm701444072004.
                
                Sreejith2014 G. J. Sreejith and S. Powell, Critical behavior in the cubic dimer model at nonzero monomer density890144042014.
                
                Kasteleyn1963 P. W. Kasteleyn, Dimer Statistics and Phase TransitionsJ. Math. Phys.42871963.
                
                Appert2008
                C. Appert-Rolland, B. Derrida, V. Lecomte and F. van Wijland,
                Universal cumulants of the current in diffusive systems on a ringPhys. Rev. E780211222008.
                
                Lecomte2012
                V. Lecomte, J. P. Garrahan, and F. van Wijland,
                Inactive dynamical phase of a symmetric exclusion process on a ringJ. Phys. A451750012012.
                
                Karevski2017
                D. Karevski and G.M. Schutz,
                Conformal invariance in driven diffusive systems at high currentsPhys. Rev. Lett.1180306012017.
                
                Jack2015
                R. L. Jack, I. R. Thompson, and P. Sollich, 
                Hyperuniformity and Phase Separation in Biased Ensembles of Trajectories for Diffusive SystemsPhys. Rev. Lett.1140606012015.
                
                Carollo2017
                F. Carollo, J.P. Garrahan, I. Lesanovsky, and C. Pérez-Espigares, 
                Fluctuating hydrodynamics, current fluctuations and hyperuniformity in boundary-driven open quantum chainsPhys. Rev. E960521182017.
                
                Castelnovo2005
                C. Castelnovo, C. Chamon, C. Mudry, and P. Pujol, From quantum mechanics to classical statistical physics: Generalized Rokhsar–Kivelson Hamiltonians and the “Stochastic Matrix Form” decompositionAnn. Phys.3183162005.
