Phases of quantum dimers from ensembles of classical stochastic trajectories

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.


I. 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 [1,2] as representations of singlet pairings of quantum spins. The simplest QDMs [3] 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 [4,5].
At zero temperature, the phase diagram of fully packed QDMs is controlled by the ratio ∕ of the energy per plaquette to the flipping frequency (see next section for definitions). The case ∕ = 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" [6]), separating states with different kinds of order at either side of ∕ = 1: for ∕ > 1 this order is known to be "staggered", extending all the way to ∕ = ∞; for ∕ < 1, in contrast, the nature of the order is less clear, the only certainty being that for ∕ = −∞ it is "columnar" [4,5]. The richness of this phase diagram highlights the complexity that * tom.oakes@nottingham.ac.uk 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) [7] 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 [8]. Fully packed classical dimer models (CDMs) also have rich dynamical behaviour [9,10]. And even the simplest constraint of hard-core repulsion as in simple exclusion processes can give rise to interesting nonequilibrium dynamics [11]. 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) [12], which provides a statistical mechanics of trajectories [13][14][15] 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 [16,17]. 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 arXiv:1802.09576v2 [cond-mat.stat-mech] 17 Aug 2018 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.

II. 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 × square lattice with periodic boundary conditions and define ( ) ∈ {0, 1} as the dimer occupation number on the link joining sites and + , where is a lattice vector and ∈ { , }. 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,

A. 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 [3], which can be written schematically as 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 .
Besides the total number of dimers, the Hamiltonian ℍ , in a lattice with periodic boundary conditions has a further conserved quantity, the flux ⃗ Φ. The components of the flux vector are defined by where = (−1) + = ±1 on the two sublattices. The flux ⃗ Φ is conserved by any local dynamics within the space of fully packed dimer configurations [4,5]. Since there are 1 2 2 dimers and each contributes ±1 to one component of ⃗ Φ, possible values of the latter satisfy |Φ | + |Φ | ≤ 1 2 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 [18], and conventional ordered phases ("dimer solids"), which break lattice symmetries. The ordered states, illustrated in Fig. 1, can be further divided according to their value of the flux, which vanishes in the columnar, Fig. 1(a), plaquette and mixed phases Fig. 1(b), and is maximized in the staggered phase Fig. 1(c).
A liquid phase of the dimer model is one that preserves the full symmetry of the lattice, and is characterized by topological order [18] and fractionalized monomer excitations [5]. According to a result of Polyakov [19], 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, ∕ = 1 (and , > 0), can be found exactly, and is a U(1) quantum liquid [20]. (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 and construct the set ℭ of all configurations that can be reached from by successive plaquette flips. By writing ℍ , as a sum of projectors [5], it is straightforward to show that the equal-amplitude superposition ∑ ′ ∈ℭ | ′ ⟩ 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 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 ∕ ), and can be shown to be ground states for all ∕ > 1 [5]. Such states saturate the bound on the flux, with ⃗ Φ = { 1 2 2 , 0} in the configuration shown in Fig. 1(c), and others related by symmetry. These states break rotation and translation symmetries of the lattice.
For ∕ < 1, the energy is instead minimized by a nontrivial superposition of dimer configurations, and so the ground state depends on ∕ . One therefore expects a discontinuity in the first derivative of the ground-state energy at ∕ = 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 Section IV, although the presence of the RK point leads to more subtle behavior on the side with ∕ < 1.) Analytical arguments do not determine the ground state for ∕ < 1, however, and numerical approaches must instead be used. The candidates are those with zero flux, illustrated in Fig. 1(a) and (b).
In the limit → −∞, the ground state maximizes the number of flippable plaquettes f . The maximal value f = 1 2 2 is achieved by the four configurations with dimers arranged in columns; one is shown in Fig. 1(a), and the others are related by symmetry. For large negative ∕ , 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 where the product is over a set of plaquettes tiling the lattice without overlapping, such as those labeled A in Fig. 1.
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 [21] and hence increasing the potential energy. They may therefore compete with the columnar state for ∕ of order unity.
An order parameter for these phases is provided by the magnetization ⃗ [22, 23], which points along one of the square axes, ± or ± , in the columnar phase and along one of the four diagonals, ± ± , 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 used by Banerjee et al. [24], 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 [25].
The phase diagram for ∕ < 1 has been extensively studied using exact diagonalization (ED) [22,26] and quantum Monte Carlo (QMC) [24,27,28], 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 ∕ = −∞ all the way to the RK point ( ∕ = 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 ∕ , was reached by Syljuåsen [27] using projector QMC methods. Ralko et al. [28], combining QMC and ED, agreed with the presence of an intermediate phase, but argued that it showed mixed order. Finally, Banerjee et al. [24], 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.

B. 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 where | ⟩ is the probability vector, with {| ⟩} the configuration basis, and ( ) the probability of configuration at time . The general form of the generator (or master operator) is The positive terms are off-diagonal and encode the possible transitions → ′ and their rates → ′ . The negative terms are diagonal, with the escape rate from configuration , = ∑ ′ ≠ → ′ . The form (9) guarantees probability conservation: the largest eigenvalue of is zero and its left eigenvector is the uniform (or "flat") state: For the specific case of the CDM, a given plaquette, when flippable, flips according to a Poisson process with rate constant [10]. The generator for the dynamics then reads 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 (11) 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, where the normalisation is given by As in the case of the Hamiltonian (2), the classical dynamics generated by conserves the flux for periodic boundary conditions, and thus each flux sector is an irreducible partition of the dynamics. At the RK point, ∕ = 1, the classical generator (11) and the Hamiltonian (2) become identical up to a sign and an overall factor that is determined by the transition rate [3], 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: 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.

C. Trajectory ensembles and dynamical large deviations
The dynamics generated by (11) is realised in terms of stochastic trajectories. For a continuous time Markov chain such as we are considering, a trajectory of overall time extension is a sequence of configurations and jumps between them, where ( = 1, … , ) indicate the times at which jumps between configurations occur. Between jumps the configuration remains unchanged, so that from the time of the last jump, , to the final time, , the configuration in trajectory (15) would be . The dynamics generates an ensemble of trajectories, defined as the set of all possible trajectories (15) 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 ( , ′ ). In particular, the latter is obtained from by summing over all stochastic trajectories (i.e., by contraction), where ′ ( ) 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 [13][14][15]29] defined as the number of configuration changes in a trajectory, which we will denote by the symbol when acting on trajectories. For example, for the trajectory in Eq. (15) we havê ( ) = , as there are a total of jumps 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. Ref. [8].
Associated with the ensemble of trajectories is a corresponding distribution for trajectory order parameters such as the activity, which has probability distribution 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) [12][13][14][15]. At these long times the statistics of are determined by the LD rate function ( ) which can be thought of as an entropy density in the space of trajectories [12][13][14][15]. Equivalent information to that found in ( ) is contained in the moment generating function (MGF) [12][13][14][15], This function generates the moments of , via ⟨ ⟩ = (−1) ( )| =0 . Just as for the probability, the MGF has a LD form at long times, given by the approximate equality in Eq. (18). The function ( ) is the scaled cumulant generating function (SCGF; its derivatives evaluated at = 0 give the cumulants of scaled by time), and plays the role of a free-energy density for trajectories [13][14][15].
One way to interpret Eqs. (17) and (18) is in terms of conditioned and biased trajectory ensembles [30]. The delta function in (17) restricts the sum to trajectories which have total activity . This is analogous to a microcanonical ensemble which restricts configurations to a fixed energy. In turn, in (18) the sum is over all trajectories but the probabilities are exponentially biased (or exponentially tilted). Here is only indirectly controlled by . 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 [12][13][14][15],

D. -ensemble
While both trajectory ensembles encode the same information about the dynamics, the "canonical" ensemble, Eq. (18), is more practical to study. This is sometimes called theensemble [31]. The probabilities of trajectories are exponentially tilted from the natural ones of the dynamics as We denote averages of trajectory observableŝ in this ensemble by ⟨̂ ⟩ , which in terms of the averages of the original dynamical ensemble read The power of the -ensemble is that it allows a full characterisation of the dynamics beyond typical behaviour by tuning away from = 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 the -ensemble is more tractable is that the dynamical partition sum (18) can be written in "transfer matrix" form, where the probability vector |i⟩ represents the distribution from which the initial state is drawn. The operator is a deformation or tilt of the original dynamical generator that reads (for the case of tilting against the activity) [13][14][15] The equivalence between Eqs. (18) and (22) can be proved in the following way. If we write the master operator as = ∑ − ℝ, where denotes the off-diagonal parts of that are responsible for all the possible transitions , and ℝ is the diagonal part of with the escape rates, cf. Eq. (9), the probability of ( ) of a trajectory such as Eq. (15) can be written as (18) we have the exponentially tilted probability Summing over all trajectories to obtain Eq. (18) then gives |i⟩, which is Eq. (22) expressed as a Dyson series for the tilted generator Eq. (23).
In contrast to , defined in Eq. (9), the tilted operator does not define a probability conserving dynamics for ≠ 0. In fact, its largest eigenvalue is ( ), and thus the problem of computing the dynamical partition sum reduces to that of maximising (23). The long-time average activity, obtained from the SCGF as serves as the dynamical order parameter that helps classify the dynamical phase behaviour, with associated susceptibility,

E. Connection to QDM
For the specific case of the dimer model, the tilted operator reads We see from Eq. (2) that this coincides with the general QDM Hamiltonian, if we identify = − and = : Changing is equivalent to changing the ratio ∕ , and so the properties of the -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 weighting. The latter gives ⟨ℚ where is given by Eq. (16) with the replacement → . Using Eq. (20), and applying the same steps that led from Eq. (18) to Eq. (22), one finds The expectation value of ℚ in |gs⟩, the ground state of ℍ − , , is therefore given by the limit of large , Note that this is true for arbitrary initial distribution |i⟩, as long as the overlap ⟨gs|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.

A. Transition path sampling
The main difficulty in sampling -ensemble trajectories is the usual one associated with calculating exponential averages [32]: the trajectories that are easy to generate with the normal dynamics at = 0, Eq. (11), are not the relevant ones for the biased ensemble at ≠ 0, Eq. (20), and the latter are exponentially rare in the original dynamics. Since the tilted operator (23) 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 [33]. In the case of trajectories one such importance sampling scheme is Transition Path Sampling (TPS) [34]. 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, Eq. (11).
The basic idea behind TPS is similar to that of Markov chain Monte Carlo but applied to trajectories [34]. 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 Eq. (20), the key quantity is the change in overall activity Δ between the old trajectory and the new one. As in standard Metropolis, if Δ < 0 the new trajectory is always accepted; if Δ > 0, it is only accepted with probability − Δ . The procedure is repeated until the ensemble of trajectories thus generated converges to the -ensemble.
The non-trivial step is (i). Various methods for generating trajectories have been proposed [34] that guarantee ergodicity in trajectory space and are reasonably efficient. In particular, we use the shifting method [34], described in Fig. 2: 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 [34].
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 Δ 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 -ensemble much more efficiently than brute force sampling (and indeed has been used successfully in this context before [31,35,36]), the exponential cost of sampling long times may render it impractical. An alternative to TPS is the cloning method [37] adapted from quantum diffusion Monte Carlo, which however also suffers from a similar exponential cost [38]. Below we discuss how to parially overcome this limitation for TPS by exploiting umbrella sampling techniques to the trajectory context [38][39][40][41].

B. TPS and trajectory umbrella sampling
The idea of umbrella sampling is the following [38][39][40][41]. We wish to compute exponential averages of the form 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 Eq. (31) as where This simply means that the average of the trajectory observable − ̂ over the original dynamics is the same as the average of the trajectory observable  − ̂ 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 (32) by an empirical average over sp sample trajectories, The sampling error is given by the variance of the average squared of the empirical average, where we have used Eq. (32) 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, where in the last line we have used Eq. (18) for long times. The convexity of the SCGF function implies that (2 ) ≥ 2 ( ) 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 sp . The aim is therefore to find an alternative reference dynamics which makes the convergence of (34) more efficient.

Ideal reference dynamics: generalised Doob transformation
The reweighting factor (33) for a trajectory such as (15) reads, where Δ = − 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 (35). This is known to be given by the generalised Doob transform [30,42,43], which maps the tilted generator (23) to a new stochastic operator̃ whose natural trajectories are those of the -ensemble. For long observation times the transformation is obtained in the following way.
From the components of the left eigenstate of , we construct a diagonal matrix , such that ⟨−| = ⟨ |. We then definẽ where is the identity operator. is stochastic, and for long times is guaranteed to generate the same trajectories as those of the -ensemble. The transition rates iñ are given bỹ while the escape rates coincide, up to a shift, with the original ones,̃ = − ( ).
If the reference dynamics is the one generated bỹ the reweighing factor then reads, We see that this form of  cancels the exponential averaging in the numerator of (35) and the error is no longer exponential in time. This means that an -tilted expectation value like (21) can be computed by simply running the dynamics with̃ [38][39][40][41].

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 first, which amounts to solving the problem exactly. Nevertheless, the form of the ideal transition rates (41) can help guide the definition of convenient approximations for the reference dynamics that are practical [38][39][40][41].
We will consider the transition rates for the reference dynamics that have the form of (41) The aim is to find a vector ⟨ | that approximates the exact ⟨ | and is also tractable numerically. The associated escape rates, in general will not be a uniform shift from the original ones as in (42). The reweighing factor is (up to boundary terms) which together with (35) 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 ⟨ | improves sampling significantly, as we now show. The exact (38) are functions of the whole configuration and in general may have correlations at large distances. A simple approximation is to assume a short-range correlated form for 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 Appendix A. This leads to depending on the configuration only through the total number f ( ) of flippable plaquettes, where the constant parametrises the function. Putting this all together, the sampling proceeds as follows: The reference dynamics we use is given by the transition rates and escape rates Comparison of TPS acceptance rates for = 6 as a function of . The orange curve corresponds to the TPS acceptance when using the original dynamics. The green curve is for the reference dynamics with = from the 2 × 2 approximation of Appendix A. The blue curve corresponds to the optimal value of found from exploring the acceptance rate landscape. The data shown is for trajectories of length = 50 for > 0, and = 5 for ≤ 0 (convergence in time is much faster on the active side ≤ 0). Each point shown corresponds to 5 × 10 6 attempted TPS moves. The optimised values used are = 0.25 for < 0, and = for > 0. Inset: Acceptance rates for = 12 in the region where the active-inactive transition occurs for this system size (see Fig. 4).
In order to sample the original dynamics tilted by − ̂ we have to sample the reference dynamics tilted by  − ̂ , see (32), which we write aŝ In order to account for this tilting we use TPS with the reference dynamics, with an acceptance rate for trajectories given by

Optimization of reference dynamics
The reference dynamics above is parametrised by the constant . While the dynamics is based on the Doob transformation for the open 2 × 2 problem, there is no reason why the value of 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 that maximises the trajectory acceptance rate Γ acc in the TPS simulations.
Finding the optimal value of is a case of exploring the landscape of acceptance rates Γ acc . Figure 3 compares the acceptance rate as a function of when using the original dynamics to that obtained when using the effective reference dynamics, Eqs. (47)(48)(49) with obtained as in Appendix A, and with a value of that optimises even further the acceptance. This latter optimal value of is obtained for each by starting from = and progressively changing until a maximum of Γ acc is reached. The effective reference dynamics defined by this optimal is then the one used for obtaining the corresponding -ensemble. Figure 4 shows the mean activity rate ⟨ ⟩ = ⟨ ⟩∕ evaluated across a range of values, with the results for the smallest system size, = 6, compared with exact diagonalization. The agreement between these results confirms that the method has converged, at least for this , and demonstrates that it is able to resolve detailed features of the dynamics for > 0. For larger system sizes, there is an increasingly sharp drop in the activity at = 0. The activity is related, by Eq. (25), to the first derivative of the thermodynamic potential ( ) (or equivalently, the quantum ground-state energy), and so this indicates a first-order transition at = 0, as expected from the analytical arguments in Section II A. (Note that our simulations are restricted to the zero-flux sector, but this is not expected to change the thermodynamic properties in the limit → ∞.) The activity histogram, shown in the inset, shows the characteristic broadening, compared with a Gaussian distribution, expected for such a transition.

IV. RESULTS
As discussed in Section II A, the ground state for positive ( ∕ > 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 > 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 . For large positive , the system is mostly restricted to the minimally flippable zero-flux configurations, which, as illustrated in Fig. 4(a), have precisely 4 flippable plaquettes for any system size [44]. Second-order perturbation theory gives the ground-state energy as − ( ) = −2 −2 , which leads, using Eq. (25), to a mean activity of ⟨ ⟩∕ = 4 −2 in this limit.
Our main results regarding the phase structure of the QDM are displayed in Fig. 5. Panels (a-d) show the ground-state order-parameter distribution, where ⃗ is the operator defined in Eq. (6) and the expectation value in the QDM ground state |gs⟩ is calculated as described in Section II E. For all negative ( ∕ < 1), the maxima of the distribution ( ⃗ ) occur for ⃗ aligned with the square axes, indicating that the ground state has columnar order. Particularly for small | |, though, the selection of this order is weak, and the distribution has approximate SO(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 ⟨gs| | ⃗ | 2 |gs⟩ and ⟨gs| cos(4 ) |gs⟩ where tan = ∕ . As panels (e-g) of Fig. 5 show, both of these quantities decrease as the RK point ( = 0) is approached, in qualitative agreement with the results of Banerjee et al. [24].
The microscopic model is only symmetric under the discrete rotations of the lattice, and so the approximate SO(2) symmetry of the order parameter ⃗ is emergent. As argued by Fradkin et al. [20], this can be understood qualitatively by considering the renormalization group (RG) flow structure. Properties near the RK point are described by an effective action written in terms of a continuum height field ℎ, where 2 , 4 , and are real parameters, and and 0 denote the space and imaginary-time derivatives, respectively. The coefficient 2 ∼ − ∼ − 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 Eq. (53), it is related to the coarse-grained height by The term therefore breaks SO(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 = −12. Standard scaling arguments in the presence of a dangerously irrelevant perturbation [45,46] therefore imply the existence of an additional length scale, ∝ | | −3 for small negative , with selection between columnar and plaquette order occurring only beyond this large scale [20]. This is qualitatively consistent with the weak columnar ordering observed for small | | 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 rms = ⟨gs| | ⃗ | 2 |gs⟩ 1∕2 , shown in Fig. 5(g). We find in Appendix B that rms ∼ 2 | | 1∕2 for small negative and large , corresponding to a critical exponent = 1 2 . Finite-size corrections, however, cause deviations from this scaling form and saturation at rms ∼ √ ln when = 0. At fixed , crossover between these two forms, ∼ 1∕2 and ∼ 0 (saturation), leads to a reduced effective exponent eff . As shown in Fig. 5(g), we find a reasonable fit to eff = 0.254 for = 16 and eff ≃ 0.1 using exact results for = 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 = 0, which is thermodynamically first-order but shows critical behavior on the negative-side ( ∕ < 1), is a common feature of constrained systems such as dimer models. A classic example is the Kasteleyn transition in 2D classical dimer models [47].

V. 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 [48][49][50], and the one-dimensional Ising model with Glauber dynamics and the transverse field Ising chain [42].
For the QDM, our main result (see Fig. 5) is that the ground state is the columnar phase for all ∕ < 1. Our results in this regime also show an approximate emergent SO(2) symmetry of the order parameter. Both of these findings agree with the observations of Ref. [24], and contradict earlier results [27,28]. 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 = 0, see Fig. 4(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 Fig. 4(b). Furthermore, the two competing phases display different kinds of structural order: while the inactive phase ( > 0) is staggered, the active phase ( < 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 [50][51][52].
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 theensemble Eq. (20)], and to overcome the numerical difficulty of accessing exponentially suppressed trajectories we supplemented TPS with a version of umbrella sampling in trajectory space [38][39][40][41]. 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 [40] 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 [53]. 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 = .
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 twodimer configurations at rate , between single-and doubledimer configurations at rate , and between the no-dimer and single-dimer configurations at rate . The former kind of tran-sition corresponds to a plaquette flip within the 2 × 2 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 2×2 to obtain a SCGF for the number of flips from the largest eigenvalue of the deformed operator 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 2 × 2 regions, as the other transitions correspond to flips in neighbouring regions.
As an approximation to the Doob transform for the full × system, we replace the exact vector ⟨ | by the product of the left eigenvectors ⟨ | of 2×2 for each 2 × 2 region. The components can be expressed in the form = , where is a (dimensionless) "energy" associated to configuration . (Both and depend on ; we suppress this for clarity.) We can characterize each configuration by the number of plaquettes with dimers (i.e., the number in each class in Fig. 6); note that 2 ≡ f ( ), in the notation of Section III B 2. The total number of plaquettes is 0 + 1 + 2 = 2 , and, the fact that the total number of dimers is constrained to be 1 2 2 implies that 2 2 + 1 = 2 . Together these give 0 = 2 = 1 2 ( 2 − 1 ) and allow us to express the dependence of the eigenvector on (as well as on the rates , , and ) as = 2 , using a single parameter obtained from the diagonalisation of 2×2 .
This value of specifies a dynamics that generates the exact -ensemble for the open 2 × 2 problem and that we use as a starting point when optimizing the reference dynamics for the full lattice (see Section III B 3).