Rejection-free quantum Monte Carlo in continuous time from transition path sampling

Continuous-time quantum Monte Carlo refers to a class of algorithms designed to sample the thermal distribution of a quantum Hamiltonian through exact expansions of the Boltzmann exponential in terms of stochastic trajectories which are periodic in imaginary time. Here, we show that for (sign-problem-free) quantum many-body systems with discrete degrees of freedom -- such as spins on a lattice -- this sampling can be done in a rejection-free manner using transition path sampling (TPS). The key idea is to converge the trajectory ensemble through updates where one individual degree of freedom is modified across all time while the remaining unaltered ones provide a time-dependent background. The ensuing single-body dynamics provides a way to generate trajectory updates exactly, allowing one to obtain the target ensemble efficiently via rejection-free TPS. We demonstrate our method on the transverse field Ising model in one and two dimensions, and on the quantum triangular plaquette (or Newman-Moore) model. We show that despite large autocorrelation times, our method is able to efficiently recover the respective quantum phase transition of each model. We also discuss the connection to rare event sampling in continuous-time Markov dynamics.


I. INTRODUCTION
Statistical mechanics has facilitated the study of manybody systems through Monte Carlo sampling for several decades.Originally developed for classical Hamiltonians, classical Monte Carlo methods allow the sampling of the Boltzmann distribution by proposing updates to a configuration together with a criterion, such as Metropolis [1,2] or Glauber [3], to accept or reject the changes.At the same time, Monte Carlo sampling can be generalised to study the Boltzmann distribution of quantum manybody Hamiltonians [4].However, a complication here is that the weights of the configurations require the evolution to be computed in imaginary time.
A standard approach is to split the imaginary-time evolution of the quantum partition sum into discrete time steps, giving rise to "trajectories" of configurations in discretised imaginary-time, where the weightings of the trajectories can be estimated for small time steps through a Trotter-Suzuki decomposition [5,6].This procedure can be made exact by considering a continuous-time expansion instead, in terms of a perturbative expansion of the Boltzmann exponential in the interaction picture and to all orders, yielding a sum over classical trajectories in continuous time [7][8][9][10].In either of these approaches, there is added complexity in relation to the classical case, as imaginary-time trajectories have to be sampled by making changes to the temporal degrees of freedom in addition to the spatial ones [11].
In this paper we focus on the continuous-time expansion of the quantum Boltzmann distribution through trajectory ensembles.For systems with discrete local degrees of freedom (and in the absence of a sign problem [12,13]) we show how the trajectories that define the quantum partition sum can be sampled through the realisation of the so-called Doob dynamics (see e.g.Refs.[14][15][16]), which optimally samples rare events.Specifically, we show how this can be achieved efficiently through exact trajectory updates which are local in space but extensive in time.In this way we can generate ensembles of trajectories using a version of transition path sampling (TPS) which is rejection-free (in the sense it does not require any acceptance criteria) [17,18].
We illustrate the method by studying the transverse field Ising model (TFIM) in one and two spatial dimensions.In this case we sample many-body trajectories by updating the trajectories of individual spins, while keeping the other spins as an effective time-dependent background.This approach is similar to the one in Ref. [17], which studied the TFIM on the Bethe lattice, and to the one in Ref. [18], which studied transition events of the classical 2D Ising model.It also generalises the approach adopted in Ref. [19] which studied the ground state properties of the 1D TFIM via time-and space-local TPS updates.Here, we demonstrate that this approach is capable of predicting the well known continuous phase transitions for the 1D and 2D TFIMs, despite the fact that the local updates suffer from large autocorrelation times close to criticality.It is important to note that there exist update schemes which can overcome the difficulties of criticality, for example, cluster updates [20,21] (a generalisation of the Swendsen-Wang algorithm [22] to continuous-time ensembles), and the Worm algorithm [23].Reference [24] demonstrates the use of our approach where these alternative methods are difficult to formulate.
We also show how to use our approach with more complex spin interactions by means of customised local updates.As a concrete example, we study the quantum Newman-Moore model [19,[24][25][26][27][28][29], also known as the quantum triangular plaquette model (QTPM).This model builds on the classical TPM [30][31][32], a system of spins with three-body interactions studied in the context of glasses, making it quantum by adding a trans-arXiv:2305.08935v2[cond-mat.stat-mech]16 Jul 2024 verse magnetic field.For the QTPM we devise local updates by simultaneously changing the trajectories of three neighbouring spins.Our method can effectively recover the first-order quantum phase transition of the QTPM [19,[24][25][26]28].We further show how thermal annealing can be implemented to improve sampling close to the transition point.
The paper is structured as follows.Section II formalises the connection between the quantum Boltzmann distribution and rare-trajectory sampling, explaining how observables can be estimated from trajectories.Section III then explains how the optimal Doob dynamics can be used to sample trajectories from the quantum Boltzmann distribution.Section IV demonstrates how this approach can be used to incorporate a singlespin update scheme for spin models with a simple transverse field, in analogy to Ref. [17], by implementing the method for the 1D and 2D TFIM.Section V generalises the method to Hamiltonians with more complex interaction terms, using the QTPM as a concrete example.We give our conclusions in Sec.VI, where we discuss the possibility of more involved schemes, including updates which are collective in space.Appendix B demonstrates how to use our approach for rare event sampling in classical continuous-time Markov dynamics, where we determine the dynamical large deviations (LDs) of the TFIM and QTPM as an example.

A. The quantum partition function and continuous-time stochastic dynamics
We consider a system with a Hamiltonian Ĥ and a discrete set of configurations that defines a basis {|x⟩} of its Hilbert space, H [33].In this basis, the Hamiltonian can be decomposed into a diagonal and a off-diagonal part, Ĥ = Ĥc − K,

K =
x,y̸ =x Kx→y ≡ x,y̸ =x We will refer to operators which act on the whole Hilbert space by symbols with a hat (e.g.Ĥc ), and their individual matrix elements by symbols without a hat.The statistical properties of Ĥ at some finite temperature are characterised by the partition function, Equation ( 3) can be expressed as a sum over stochastic paths by considering its Dyson series expansion [34], where the integrals over times t m are performed with the limits β ≥ t M ≥ • • • ≥ t 0 .Note that the product and integrals in Eq. ( 4) are omitted for the case of M = 0.Each path in Eq. ( 4) can be interpreted as a classical stochastic trajectory in continuous time with time extent β, where the operator Kxm−1→xm makes the instantaneous transition x m−1 → x m at the time t m .We denote a stochastic trajectory of time extent β with M jumps as Note that the last jump occurs at t M ≤ β.We then represent Eq. ( 4) as a sum over these jumps, where x(t) denotes the configuration of the trajectory x at time t, the delta function δ[x(0), x(β)] allows only for trajectories which are periodic in time, and K(x) is the trajectory observable (which will henceforth be denoted by a calligraphic font) which returns the number of transitions, M , which occur in the trajectory x.The sum over {x} indicates a sum over all paths for any number of jumps, M ∈ [0, ∞).We can then write the probability of the trajectory x as The probability Eq. ( 7) can be related to that for trajectories generated by a continuous-time Markov dynamics in the following way.We define the Markov generator Ŵ = K − R, with the transitions given by K and their associated escape rates, R, R = x,y̸ =x Since we consider quantum systems without a sign problem [12,13], we can assume that the probabilities Eq. ( 7) are positive real numbers.With these definitions, the relation between the Hamiltonian and the stochastic generator is It follows that the trajectories that define the partition function Eq. ( 6) are related to those generated by the stochastic dynamics of Ŵ with the following properties: (i) The trajectory must start and finish in the same configuration, i.e., it should be a stochastic bridge.
A trajectory which does not meet this criterion has zero probability of occurring.
(ii) The probability of a trajectory that satisfies (i) is exponentially biased (with respect to the probability of occurring under Ŵ ) by the time integral of Ĥc − R along the trajectory.That is, its probability is multiplied by the factor e Condition (i) is needed due to the imposition of periodic paths in Eq. ( 4).The bias in (ii) accounts for the difference in the diagonal parts of Ŵ and Ĥ.Many methods have been developed to efficiently sample biased dynamics as in (ii).One popular approach is transition path sampling (TPS), see e.g.Refs.[19,35,36].Here, we will use TPS focusing on techniques for sampling trajectory updates which are local in space, cf.Refs.[17][18][19], and which also account for (i).

B. Calculating observables
The thermal expectation value of an observable Ô is given by, This can be rewritten as an imaginary-time average by noticing that we can arbitrarily move Ô through the exponential due to the properties of the trace, In order to consider how to compute Eq. ( 11) from trajectories we will deal separately with the cases when Ô is diagonal and off-diagonal.
For the case of a diagonal operator, Ôdiag = x O(x) |x⟩ ⟨x|, Eq. ( 11) can be directly written as a sum over all trajectories with the probability given by Eq. ( 7).This leads to where is the time-integrated trajectory observable.That is, we average the value of the observable over all times and over all trajectories.
For an off-diagonal operator, Ôoff−diag = y̸ =x O x→y |y⟩ ⟨x|, we can expand both exponentials in Eq. (11).This expansion in the integrand of Eq. (11) gives a sum over all trajectories which have the transition x → y at the time t, summed over y ̸ = x.Compared to Eq. (6), the jumps at time t appear with factors O x→y rather than K x→y .Thus, if we want to express ⟨ Ôoff−diag ⟩ as a sum over trajectories with the probability given by Eq. ( 7), we need to account for the change K x→y → O x→y at time t.Furthermore, since we are time averaging, we can replace the integration over time by K x→y (x) which counts the number of jumps x → y which occur in trajectory x.We can then write

III. OPTIMAL SAMPLING
We now explain how trajectories can be sampled efficiently and exactly from the partition function Eq. ( 6) by means of a rejection-free form of TPS.A general method for converging to an ensemble of trajectories such as Eq. ( 7) is TPS, a form of Monte Carlo sampling in trajectory space [35].A typical problem with standard TPS is that acceptance of trajectory updates can become exponentially small in system size (and/or the length of the trajectory), thus slowing down convergence.
In the language of stochastic dynamics, the procedure we now describe is sometimes referred to as obtaining the Doob dynamics [14][15][16], a proper (normalised) stochastic dynamics derived from the original generator Ŵ , which generates trajectories with a conditioned/biased probability P β (x), such as Eq.(7).While this approach is exact in theory, it often relies on the computation of expressions which are, in practice, analytically (and often numerically) intractable.However, we will show later in the paper that we can in fact implement this general idea via optimal local updates.

A. Edge configurations
The first obstacle is to sample the initial/final configuration of trajectories, x 0 .While Eq. ( 6) has no explicit probability distribution for the initial configuration x 0 , its probability will be decided by its possible transitions and its diagonal component in Eq. (1).That is, we choose some initial configuration, x 0 , with probability Once the trajectory edges have been selected, we can then sample the remainder of the trajectory from the subset of trajectories which have the required boundary conditions.That is, we want to sample from the reduced dynamics where x i and x f are the initial and final configurations.Note that, although we are only interested in trajectories which have x i = x f , it will be useful to solve the more general problem for later for the TFIM and QTPM.

B. Continuous-time dynamics
The continuous-time Monte Carlo (CTMC) method [37][38][39] is a dynamical protocol which can be employed to sample the dynamics of Eq. (6).While this method is straightforward for a normalised time-homogeneous dynamics, it will be useful to formulate the general case of time-dependent dynamics for what follows.
Consider time-dependent instantaneous transition rates, Kx→y (t; β, x f ) ≥ 0, for transitions from configuration x to configuration y at time t ≤ β.We use the tilde to distinguish these from the off-diagonal entries of the Hamiltonian, cf.Eq. ( 2), and also allow for a dependence on the final configuration, x f , and the overall trajectory time, β.The associated escape rate from configuration x at time t is The escape rate in turn determines the distribution for the waiting time τ , P x , for a jump out of x after t, After a waiting time has been drawn from Eq. ( 17), the transition into y at time t + τ is chosen with probability C. Optimal transition rates Given a dynamics defined by transition rates Kx→y (t; β, x f ), the CTMC method generates trajectories with probability density for trajectories starting from x i (where x f acts as a parameter in the definition of the rates).However, we are interested in sampling trajectories from the (conditioned and tilted) distribution Eq.(7) where the final configuration is x f , that is, In the expression above we allow (for later convenience) the final configuration to be fixed to a different value from the initial one.As written, Eq. ( 20) is not the distribution of trajectories generated with transition rates K xi−1→xi , due to Eq. ( 20) satisfying conditions (i) and (ii) above (the former generalised to some fixed initial and final configurations), and the need for the explicit normalisation Z β (x i , x f ).Furthermore, while Eq. ( 19) appears simpler than Eq. ( 20) due to the apparent absence of the conditioning of the boundary conditions in time, this is accounted for by the fact that the transition rates in Eq. ( 19) are time-dependent.
Our aim is to find a dynamics Kx→y (t; β, x f ) such that Pβ (x) = P β (x); that is, the optimal dynamics (or Doob dynamics) for sampling P β (x) [15].These optimal transition rates can be computed from the probability that the transition x → y occurs at the time t under the original dynamics, conditioned by the probability that the trajectory is in configuration x at time t, where in the second line we have inserted the definition of Kx→y , Eq. ( 2), and in the third line we used Eq. ( 15).From this, we are able to calculate the escape rate from the configuration x at time t, where we have used Eq. ( 2) to account for all possible transitions out of x and the derivative of Eq. ( 15) to obtain the last line.Based on the rates from Eqs. ( 21) and ( 22), this time-dependent dynamics gives Pβ (x) = P β (x) if the distribution of initial configurations is set to be x i [40].Appendix A explains how the optimal transition rates can be used with the CTMC algorithm from the previous section to sample the time-dependent dynamics.

D. Example: single two-level system
As an illustration of the ideas above, we consider a simple problem which will be important later in this paper: a single two-level system with the states σ = {−1, +1}, which we will refer to as a spin.A generic Hamiltonian for such a system reads where X, Ŷ and Ẑ are the usual Pauli operators.By calculating the matrix exponential of Eq. ( 23), is can be shown that the reduced dynamics, Eq. ( 15), takes the form ) where θ t = t 1 + g 2 .Using Eq. ( 24), we are able to define a time-dependent stochastic dynamics which allows us to exactly sample trajectories with probabilities given by Eq. ( 7).The initial/final configuration σ 0 is chosen with probability Z β (σ 0 , σ 0 )/ σ Z β (σ, σ), and the timedependent transition rates for the dynamics are By simply running this dynamics we sample the trajectories we need for computing thermal averages in this problem.The results are shown in Figs.1(a, b) for the average magnetisations ⟨ X⟩ and ⟨ Ẑ⟩ for various inverse temperatures β.The data points show the trajectoryaveraged values using Eqs.( 12) and ( 13).The numerical results coincide with the exact averages,

IV. TRANSVERSE FIELD ISING MODEL IN DIMENSIONS ONE AND TWO
We now show how the trajectory sampling approach can be used to sample from the Boltzmann distribution of the transverse field Ising model (TFIM) with N spins Ẑi Ẑj (27) where Xi and Ẑi are the Pauli operators acting on site i, and nearest neighbours are denoted by ⟨i, j⟩.The method we use in this section can be equally adapted to accommodate for any other potential (the second summation) with the same single-body kinetic term, see App.B for an example with a global diagonal operator.In the next section we demonstrate how to deal with other kinetic terms.
In what follows we set h = 1 and consider periodic boundary conditions in space.Sampling directly from the full partition function Eq. ( 4) is difficult.We instead use a single-spin update scheme, whereby we update the entire trajectory for a single spin, keeping the trajectories of all other spins fixed by sampling directly from a reduced partition function.

A. Redrawing trajectories for an individual spin
We now explain how our formulation can be used to perform single spin updates of the many-body trajectories.These updates adjust the entire trajectory of a single spin within the many-body trajectory, modelling the other spins as an effective time-dependent environment (or heat bath [17]).We stress that this update is equivalent to that presented in Ref. [17], which considered Eq. ( 27) on the Bethe lattice, and similarly that of Ref. [18] for the classical Ising model.Here, we instead use this approach on the TFIM in 1D and 2D to demonstrate how the update can be used to investigate the ground state properties and phase transitions of quantum lattice models.

Factorisation of the partition function
We consider the many-body trajectory ω as a collection of the N individual spin trajectories, ω = {σ 1 , • • • , σ N }, with σ j denoting the time series of transitions for spin j each with a total of M j transitions.The partition function becomes In the above, all K σ→σ ′ = 1, cf.Eq. ( 27), but we keep them in the expression to keep track of the spin transitions.
At this point, we notice that we are able to factorise Eq. ( 29) into a factor which depends on a specific σ l , and a factor which does not, where {i; ⟨i, l⟩} denotes sites i which are nearest neighbours of l, and {j ̸ = l; ⟨j, p⟩} sites j different from l which are nearest neighbours of p.
We can now express Eq. ( 30) as a sum over each spin l and for each l over the partial trajectories that is, the individual trajectories of all the spins other than l.The first factor in Eq. ( 30) depends on the trajectory of spin l and its neighbouring spins through the interactions i;⟨i,l⟩ Jσ i (t)σ l (t), which we rewrite in the form of an effective time-dependent longitudinal field The second factor in Eq. ( 30) only depends on the trajectories of the spins p ̸ = l, and we denote it N (σ (l) ).
The partition function can now be written as Now suppose we have some many-body trajectory ω, and we have randomly chosen to update site l.The partition function for spin l is Z l β [g l ], with the time-dependent longitudinal field g l (t) which is piecewise constant, meaning it remains constant except for a discrete number of times, M , when it instantaneously transitions to some other value.M is the total number of times any of the neighbouring spins transition.The partition function can now be written in the following way: The spin configurations σ m = σ l (τ m ) are the configuration of spin l at the times τ m where the field g l (t) instantaneously changes, and we set σ M +1 = σ 0 .Equation ( 34) is a sum over all possible spin configurations at each of these times.The weight for each sequence of spin configurations {σ m } is given by the product in Eq. (34).The factors Z ∆τm (σ m , σ m+1 ; g l (τ m )) are exactly the partition function of the single spin problem, Eq. ( 24), with time extent ∆τ m = τ m+1 − τ m and longitudinal field g l (τ m ), which is taken to be the value after the mth transition of neighbouring spins, and can be written as In the above equation, σ l m = σ l (τ m : τ m+1 ) denotes the partial trajectory of spin σ l between times τ m and τ m+1 .
It is simple to check that Eq. (34) with Eq. (35) gives We now explain the strategy to update a trajectory ω, making use of the above equations and Fig. 2, which demonstrates the update for the 1D TFIM.We first randomly select a spin l ∈ {1, . . ., N } to update, keeping the trajectories for all other spins fixed.The trajectories of spins neighbouring l define the time-dependent longitudinal field g l (t).The process of selecting spin l and determining the time-dependent field takes us from Eq. ( 29) to Eq. ( 33), and allows us to sample a trajectory for spin l through the partition function We then exploit the fact that the longitudinal field is constant in between the transition times τ m for spins neighbouring l.In particular, we write Z l β [g l ] as a sum over all edge configurations, σ m , which gives Eq. ( 34).The edge configurations can be sampled using transfer matrices, which is explained in Sec.IV A 2, see Fig. 2(b).Finally, we sample stochastic bridges between configurations σ m and σ m+1 as explained in Sec.IV A 3, see Fig. 2(c).

Sampling the edge configurations
Given some time-dependent longitudinal field g l (t), we wish to sample the states of the spin σ m = σ l (τ m ) at the times τ m which the longitudinal field changes value.The probability of observing the sequence of configurations Each of the partition functions Z ∆τm (σ m , σ m+1 ; g l (τ m )) has four possible values.Notice that if we were to think of each σ m as a spin on a periodic 1D lattice with M +1 lattice sites, then Z ∆τm (σ m , σ m+1 ; g l (τ m )) can be thought of as a transfer matrix between lattice sites m and m + 1.
To efficiently sample the sequence of configurations, we start by sampling σ 0 .The probability of observing some σ 0 is given by with σ M +1 = σ 0 .Notice that calculating P (σ 0 ) is equivalent to performing M +1 matrix multiplications (for 2×2 matrices), and can then be calculated numerically.The configurations σ m for m = 1, . . ., M can be calculated recursively using The numerator is calculated as a multiplication over M + 1−m matrices and m real numbers, and the denominator is known from the previous iteration of the calculation.
It is simple to show that sampling from each of these M + 1 distributions yields a sequence of configurations sampled with probability given by Eq. ( 36).

Sampling the full trajectory
After we have sampled the configurations σ m , all that is left to do is to sample from each of the partition functions Z ∆τm (σ m , σ m+1 ; g l (τ m )).That is, we must sample stochastic bridges between configurations σ m and σ m+1 as was done previously for the two-level system.This indicates that the trajectory between times τ m and τ m+1 is sampled from a dynamics initialised in configuration σ m , with transition rates (39) Each partial trajectory between times τ m and τ m+1 is sampled with probability Concatenating the partial trajectories gives the fully sampled trajectory for the spin σ l , which is sampled with probability Notice that the trajectories are defined to be piecewise continuous, and so we can drop the intermediate δ[σ m , σ(τ m )] terms in Eq. ( 40),

Detailed balance
For the trajectory update to be rejection free, it must obey detailed balance with respect to Eq. ( 29).It can be verified that this is true by considering the transition probabilities.That is, given some trajectory ω, the probability to generate a new trajectory ω is where the factor of N comes from the fact that site l is chosen randomly.Notice that since P (σ l | σ (l) ) ∝ P β (ω), it follows that and detailed balance is obeyed.

B. Monte Carlo method
The method of resampling a trajectory for a single spin given in the previous section can now be used to implement a Monte Carlo algorithm to sample the many-body trajectories using the following steps: 1. Create some initial seed trajectory, ω = (σ 1 , • • • , σ N ) with inverse temperature β.

Choose a random lattice site
3. Generate a new trajectory, σ l , for the spin.

Repeat from step 2 until convergence.
There are various ways to generate an initial seed trajectory.The key requirement is periodicity in the time interval [0, β].The simplest way to achieve this is to generate the trajectory for each spin independently with a non-interacting dynamics, but allowing for the possibility of a local longitudinal field (see Sec. III D).The choice of the initial trajectory can be guided by the phase one is trying to target.For example, for the TFIM with J/h > 1 we expect to see ferromagnetic behaviour, and so we could choose a large local longitudinal field to force an initial trajectory with a higher likelihood of magnetic ordering.A second approach, which we discuss in Sec.V B is thermal annealing.
For |J/h| > 1, there is spontaneous breaking of the symmetry at the level of the ground state.Indeed, depending on the initial trajectory seed, the single-spin update approach will break the Z 2 symmetry and will only be ergodic over one of the ground states.In practice, one should randomly perform the global update σ i (t) → −σ i (t) for all spins simultaneously to ensure both states are explored.In order to investigate the effects of the local updates, we choose to not do this here.

C. Continuous phase transition in the TFIM
We now demonstrate the effectiveness of our approach for the 1D and 2D TFIM.For fixed h = 1, both models are known to undergo a continuous phase transition at the critical points J c = 1 [41] and J c ≈ 0.32747 [20,21,23], respectively.We investigate their ground state by using the previously described trajectory sampling algorithm with β = 128 and a variety of system sizes, N .Figures 3(a, b) show the average transverse magnetization, in 1D and 2D, respectively, while Figs.3(c, d) show the two-spin correlator, For 1D, we compare our results for systems of size N = 16 to 1024 to those from infinite matrix product states (iMPS) [42,43].These methods allow us to determine the properties of the 1D chain to high accuracy.In particular, we also use the infinite time-evolving block decimation (iTEBD) method [42,43].For 2D, we show results for systems N = 4 × 4 to 32 × 32.These demonstrate agreement with the 2D density matrix renormalization group (DMRG) [44][45][46] for N = 4 × 4 (dashed line).The characterisation of a continuous phase transition is facilitated through the use of an order parameter which displays a discontinuous derivative at the critical point J c in the large system size limit.The commonly used order parameter for Ising models is the longitudinal magneti-sation, Figures 3(e, f) demonstrate a sharp increase for M |z| at the critical point, J/J c = 1.

D. Trajectory autocorrelations
We now investigate the autocorrelation properties of the sampling dynamics using the single-spin update for the 1D and 2D TFIM.Under TPS we generate a Markov chain of trajectories, {ω 0 → ω 1 → • • • → ω Ntraj }.In order to test the convergence to the target trajectory ensemble, cf.Eq. ( 7), we define the autocorrelation between two trajectories in the Markov chain separated by µ TPS updates, where σ j ν (τ ) indicates the value of the j-spin of the ν-th trajectory at time τ .The trajectory ensemble average can then be estimated from the Markov chain, where ⟨ Ẑ⟩ is the equilibrium expectation value of Ẑi (for any i).With the above definition, this autocorrelator is normalized to be C 0 = 1 and C ∞ = 0 for N traj → ∞.We show C µ as a function of J/J c in Fig. 4 for µ/N ∈ [1,7], for (a) the 1D TFIM and (b) the 2D TFIM.The same data is shown in panels (c) and (d) as a function of µ/N .Far from the critical point the trajectory autocorrelation function decays approximately exponentially.Near the critical point, this decay indicates that trajectories remain correlated even after a considerable number of TPS iterations.This is a manifestation of the expected slow down of Monte Carlo sampling near criticality at the level of trajectory sampling.This is a short-coming of the method when compared to methods which use non-local updates, such as rejection-free cluster updates [20,21] and Worm algorithms [7,8,23].However, it is important to note that it is not always easy to find cluster updates for interaction terms which are more complex than the Ising interaction in Eq. ( 27).On the contrary, the approach described here can be generalized to investigate arbitrary diagonal terms with a transverse field, see Ref. [24] for triangular plaquette interactions and App.B for an example with a global diagonal operator.

V. QUANTUM TRIANGULAR PLAQUETTE MODEL
We now show how the local update scheme can be generalised to models with more complex kinetic terms in their Hamiltonians.As a concrete example we consider the quantum triangular plaquette model (QTPM) [19,[24][25][26][27][28][29], a generalisation of the classical TPM studied in the context of glassy systems [30][31][32].This is a model defined on a triangular lattice with interactions between a subset of the triangular plaquettes and a magnetic field transverse to them.We write the Hamiltonian of the QTPM as where △ indicates the upward pointing triangular plaquettes in the triangular lattice, see Fig. 5.We have chosen a representation of the QTPM where the interactions are off-diagonal in the classical basis and the magnetic field is longitudinal.This corresponds to the dual model of the usual QTPM.
Most numerical studies [19,24,25] indicate that the QTPM has a first-order quantum phase transition at |J c /h c | = 1 between two distinct phases.In our recent paper, Ref. [24], we used the method presented in Sec.IV with single-spin updates on the Hamiltonian dual to Eq. ( 50).The aim of this section is to demonstrate how the update scheme can be adjusted to account for Hamiltonians like Eq. ( 50) using local updates which redraw the trajectories for multiple spins, rather than a single spin.

A. Plaquette updates
The obvious Monte Carlo update is to randomly select a plaquette labelled by △ encompassing three sites {i, j, k} ∈ △.The trajectory of this plaquette is η △ = {σ i , σ j , σ k }.The corresponding transitions are those where the three spins in the plaquette flip simultaneously, in accordance with the kinetic term in Eq. ( 50).Given the eight different possible states of these three spins, it might appear that one needs to consider this eight-level system in the simulations.However, as we will show below, what actually matters is the change in sign in the plaquette magnetisation, which reduces again the analysis to that of a local twolevel system, see Fig. 5.
The time-dependent effective dynamics for the plaquette whose trajectory we choose to update is given by the six neighbouring plaquettes which contain any of the spins i, j or k.Generalising the approach used for the TFIM above, we do not modify the trajectories of these neighbouring plaquettes when updating plaquette △.The neighbouring plaquettes will have a total of M flips at times 0 ≤ τ 1 ≤ τ 2 ≤ • • • ≤ τ M ≤ β.These transi-tions will force a single spin in η △ to flip at each of these times.

Factorisation of the partition function
We start by writing the partition function as a sum over all possible imaginary-time trajectories [47], where a plaquette trajectory is given by and m denotes the corresponding plaquette flip.As in the case of the TFIM, we can write Eq. ( 52) in a factorised form.Here, however, there will be three contributions.The first will contain all degrees of freedom and transitions within the plaquette △.The second will contain all the transitions associated with the six other plaquettes that contain any of the spins {i, j, k} ∈ △, and it will determine the effective dynamics for the plaquette △.The third factor contains all other contributions which do not affect the dynamics of plaquette △.The partition function reads where △ ′ labels the plaquettes which share a site with △, and △ ′′ labels those which do not.We will now single out and propose a trajectory update for just one of the N plaquettes.The last factor does not depend on any of the spins within plaquette △: as before, we denote σ (△) the trajectory of all the spins except from those {i, j, k} ∈ △, and write this factor as N △ σ (△) .The second factor does depend on the spins within plaquette △, but corresponds to plaquette flips which are not for plaquette △.By this we mean, each lattice site {i, j, k} ∈ △ will belong to two other plaquettes.When either of these plaquettes flip, so will the spin σ {i,j,k} .However, the trajectory of each of these plaquettes is fixed, and thus we call this factor M △ ({η △ ′ }), where {η △ ′ } denotes the trajectories of all neighbouring plaquettes.The first factor corresponds to plaquette △ and depends on {η △ ′ }.We can now rewrite Z TPM β as a sum over all plaquettes △ and over all trajectories of plaquettes {η (△) } which are not △, The partition function Z △ β {η △ ′ } is the sum over all possible trajectories for the plaquette △, subject to the transitions which occur in the neighbouring plaquettes at fixed times.
As was done for the TFIM, we can now factorise the partition function into M +1 components, where M is the total number of times any of the neighbouring plaquettes flip.We can write where ∆τ m = τ m+1 − τ m .The plaquette configuration η ′ m is related to η m by a single spin flip, which is predetermined by the neighbouring plaquettes.This is true for all m ̸ = 0: periodic boundary conditions in imaginary time gives the condition η ′ 0 = η 0 = η M +1 .

Determining the trajectory sector and sampling the edge configurations
An important point to notice is that the trajectory space of η △ , for some chosen △, is composed of four disconnected sectors, which we refer to as the trajectory sectors.An initial plaquette configuration η △ 0 = {σ i 0 , σ j 0 , σ k 0 } will belong to one of the four local sectors shown in Fig. 5.At the times τ m when there is a transition in a neighbouring plaquette △ ′ , the sector of the plaquette △ will transition to one of the other three sectors, as shown in Fig. 5; however, the sector it moves to is entirely determined by which plaquette △ ′ transitioned.Thus, the initial sector for η 0 will predetermine which local sector the plaquette will occupy at any time.
This observation implies that the local plaquette effectively behaves as a two-level (rather than an eight-level) system, and the flipping of neighbouring plaquettes at times t m changes the energies of the two levels according to the rules shown in Fig. 5.The configurations η m in Eq. ( 56) can be sampled using the transfer matrix approach as was done for the TFIM, see Sec IV A. While the transfer matrices for the QTPM are 8 × 8 matrices, the previous observation allows them to be treated as four separate 2 × 2 matrices.This allows us to make use of the results in Sec.II, which is computationally easier than solving the original eight-level system.To determine the trajectory sector the plaquette will lie in, we must first calculate Eq. ( 56) for all four sectors, and use this as a weighting to randomly select one of the four sectors.Then, we can sample the edge configurations η m after each time τ m , as was done in Sec.IV.This gives a sequence of configurations {η 0 , • • • , η M } which are sampled with probability

Mapping to a spin and sampling stochastic bridges
The partition function Z △ ∆τm (η ′ m , η m+1 ) describes the ensemble of all possible trajectories for the plaquette △ between the times τ m < t < τ m+1 .Within these times, the plaquette can only hop between two distinct configurations for some chosen η ′ m .These configurations will have opposite magnetisation (c.f.Fig. 5).As such, it is convenient to treat the plaquette as a spin κ △ , with and an effective longitudinal field which is constant between times τ m < t < τ m+1 .We can now write where η △ m = η △ (τ m : τ m+1 ) and κ △ m = κ △ (τ m : τ m+1 ) are partial trajectories between times τ m and τ m+1 for the plaquette and spin respectively.Now that the partition function is expressed as the dynamics of a single spin, it can again be sampled us-ing the results from Sec. III D with a longitudinal field g(τ m ), and then mapped back onto the plaquette.The probability of generating each stochastic bridge is The probability of generating the full plaquette update is the product of sampling the edge configurations and then sampling stochastic bridges between them,

Detailed balance
We can now verify that the single plaquette update obeys detailed balance and is thus rejection free.Given some trajectory ω, the probability to generate some new trajectory ω is where the factor N −1 comes from the fact the plaquette △ is chosen at random.Since P (η △ | σ (△) ) ∝ P β (ω), it follows that and thus detailed balance is obeyed.

B. Thermal annealing
The expected first-order phase transition at J c can slow down the convergence for J ≈ J c and large inverse temperatures, β.In particular, if the initial seed trajectory is chosen to be in one of the two phases, there is the possibility that the update procedure cannot explore the entire trajectory space, thus remaining stuck in the incorrect phase.This is a consequence of the large barriers which need to be overcome to move between phases.Furthermore, making assumptions a priori on which phase the trajectory should belong to for some value of J could bias the results.
A common technique used in Monte Carlo sampling to overcome such metastability is thermal annealing.In this approach, we start from a small inverse temperature, β = 0.1, and gradually increase it to the target inverse temperature, β = 128.Our annealing schedule is to make N updates to the trajectory and then increase β by ∆β = 0.1 for β < 32, and ∆β = 1 for β ≥ 32.When the inverse temperature is increased, we have to modify the trajectory to account for this.In practice, we just stretch the trajectory time by a factor of (β + δβ)/β.After reaching the target inverse temperature, we do N × 10 2 updates, and then restart the process.We repeat this procedure 10 3 times.
We find this approach to work well.Indeed, when sufficiently far from J = 1, the generated trajectories (at the target dynamics) have behaviour corresponding to their correct phase with high accuracy.However, when close to J = 1, the trajectories can have properties which correspond to either phase (in practice, we find that the process picks just one of the phases for each run of the annealing process).While we cannot be certain this approach guarantees the correct amount of mixing between the phases, it provides a less biased way to propose initial trajectory seeds, and still demonstrates the first-order behaviour of the transition point.
C. First-order quantum phase transition of the QTPM We now demonstrate how the plaquette update scheme with temperature annealing can be used to investigate the quantum phase transition of the QTPM.As described in Ref. [24], the finite size study of the QTPM has to be done with care, since, depending on the size and aspect ratio of the system used, Eq. ( 50) can have many or no symmetries.We focus on square lattices N = L × L, with L = {4, 5, 8, 11, 16, 23, 32} chosen such that there are no such symmetries [24].We also compare our results against those from 2D DMRG for N = 4 × 4.
Figure 6(a) shows the average transverse three-spin correlator, The inset shows the results close to the transition point.tion, As we increase the system size, the crossover between the two phases becomes increasingly step-like, suggesting a first-order singularity in the large-size limit.The behaviour of the magnetic susceptibility, is consistent with this interpretation, with its peak getting higher and narrower around J = 1 with the system size, see Fig. 6(c).

VI. CONCLUSIONS
In this paper we have leveraged the connection between the continuous-time expansion of the quantum Boltzmann distribution and the rare-event sampling in stochastic dynamics for systems which have no sign problem.In particular, we have focused on the so-called stoquastic Hamiltonians, where computing the partition sum is equivalent to sampling imaginary-time trajectories of a continuous-time Markov chain.Each trajectory in the ensemble is conditioned to return to the initial configuration, and its probability is exponentially biased (or tilted) due to the difference between the Hamiltonian and the associated stochastic generator.Such trajectories can be accessed with a method like TPS, i.e.Monte Carlo for trajectory ensembles.Specifically, we showed that in systems with finite-state local degrees of freedom (such as spins on a lattice) one can use an approach similar to that of Refs.[17,18] to devise a rejection-free TPS scheme, by means of an exact local generation of trajectory updates which is especially simple in spin-1/2 mod-els.In fact, as we showed above, this gives the optimal, or Doob, dynamics for sampling the rare trajectories.
We illustrated the effectiveness of this approach by studying the quantum phase transitions of two classes of models.The first included the TFIM in 1D and 2D, where the transition is well known to be continuous.We showed that the rejection-free TPS method correctly characterises their quantum phase transition, even in the near-critical regime where trajectories take many TPS iterations to decorrelate.The second class of models we considered are quantum plaquette models.In particular, we studied the QTPM, and showed how to generalise the rejection-free method to local multi-spin updates.While less understood than the TFIM, the QTPM has a quantum phase transition which is firstorder.Again, the rejection-free method efficiently recovered the quantum phase transition.As a by-product we also computed to high accuracy the statistics of dynamical observables (both in the large deviation regime and for finite imaginary times) in the trajectory ensembles that resolve the quantum partition sums of both the TFIM and the QTPM, which is shown in Appendix B.
The method described here was used recently in Refs.[24,48], where the ground state phase transitions of various spin models were characterised for large system sizes.We foresee that our approach will be useful for numerical investigations of glassy models, where nonlocal updates could be difficult to formulate.While the method we used here is based on local updates, the key property is not locality but simplicity, which allows the sampling of exact trajectory moves by solving Eqs. ( 21) and ( 22) that define the Doob dynamics.It would be interesting to find other (perhaps non-local) moves that also solve Eqs. ( 21) and (22).While finding exact updates might prove difficult, one might be able to formulate an approximately optimal dynamics, with the error accounted for in an acceptance test, such as Metropolis in a TPS scheme.The difficulty here is finding a dynamics which provides improvements in sampling while defined in terms of transition rates which remain computationally cheap (perhaps building on the use of tensor networks for optimal sampling [49]).A second proposition could be to implement more advanced sampling methods to overcome local minima, such as parallel tempering.While this could present some technical challenges, it might be advantageous for the sampling of glassy models.

ACKNOWLEDGMENTS
We acknowledge financial support from EPSRC Grant no.EP/V031201/1.LC was supported by an EPSRC Doctoral prize from the University of Nottingham.Calculations were performed using the Sulis Tier 2 HPC platform hosted by the Scientific Computing Research Technology Platform at the University of Warwick.Sulis is funded by EPSRC Grant EP/T022108/1 and the HPC Midlands+ consortium.We acknowledge access to the University of Nottingham Augusta HPC service.
The MGF has a form similar to the partition function, Eq. (3).In terms of trajectories, it reads, cf.Eq. ( 4), Furthermore, from Eq. (B6) we also see that the MGF can be written as where Ĥs is a tilting [51] of the original Hamiltonian, ) This means that in the limit of small temperatures, Eq. (B4), we get that θ(s) = E s − E, where E s is the ground state energy of Ĥs and E that of Ĥs=0 = Ĥ.
For the following we use this to connect quantum phase transitions to transitions in the LD statistics of trajectory observables.

TFIM
We first use the methods of Sec.IV to investigate the trajectory statistics of the TFIM in 2D.The first trajectory observable we consider, cf.Eq. (B2), is the time integral of the two-point correlator, We focus on the low-temperature limit, β ≫ 1, corresponding to the ground state behaviour, where LD theory can be applied.The LD statistics are retrieved from the ground state properties of the tilted Hamiltonian, which is the same as the original Hamiltonian with J → J − s. Figure 7 shows the rate function, −φ(M zz /β), c.f. Eq. (B4), for the 2D TFIM with (a) J/J c = 1, (b) J/J c = 0.5 and (c) J/J c = 1.5.At criticality, we observe a broadening in this distribution, demonstrating the divergence in correlation lengths.In contrast, away from the critical point the distributions become narrower.
The second trajectory quantity we consider is the timeintegral of the order parameter, To simulate dynamics with Eq. (B12), we can use the single-spin update scheme previously described, but we must now consider the trajectories of all other spins when constructing the time-dependent dynamics, due to the coupling introduced via M |z| in Eq. (B12).While this makes the procedure more costly, we are still able to run dynamics for moderate N .We show the rate function in  Here, the broadening at criticality is more pronounced; for comparison, we show the rate function for the longitudinal magnetisation of a single spin (dashed line).This broadening suggests a diverging magnetic susceptibility in the large size limit.

QTPM
We now investigate the effects of the first-order phase transition of the QTPM at the level of the imaginary-time trajectories.We consider the statistics of the trajectory observable  for various J.The rate function at J c is shown in Fig. 8(a): the rate function flattens with increasing system size, indicating the existence of large fluctuations.This behaviour is characteristic of a (dynamical) firstorder phase transition due to the coexistence of two distinct (dynamical) phases.For comparison, we show the LDs of the single-spin problem from Eq. ( 23) (dotted lines), where the value of g is chosen to fix the mean of the distribution.Figure 8(a) shows the LDs of the singlespin problem with ⟨M z ⟩ = 0.42 and 0.8, which approximately match the tails of the distribution for the QTPM.The interpretation is clear: the two coexisting phases are homogeneous phases of distinct M z , and intra-phase fluctuations are uninteresting (thus the modes are well approximated by a single spin); the broadening in φ in the QTPM is due to a Maxwell construct between these modes due to the fact that intermediate values of M z are realised by coexistence, i.e. space-time regions of one phase separated by sharp interfaces from space-time regions of the other phase.In Figs.8(b) and (c) we also show the rate functions away from the transition point for J = 0.8, 1.2, respectively.While there is a slight broadening in the tails of these distributions, they still describe single phases far from coexistence.
It is also possible to reasonably estimate the probability distribution P β (M z ) for finite β through sampling.This is shown in Fig. 9 at the transition point for L = 16, for various inverse temperatures.We use our annealing strategy to sample N traj = 10 6 trajectories for each inverse temperature.The results are used to approximate the probability distribution function using the Gaussian kernel, where M z (ω i ) is the time-integrated magnetisation of the i-th trajectory, and δ = 0.005 is the width of the filter.Notice that at small inverse temperatures the distribution looks approximately Gaussian.However, with increasing β, the distribution becomes bimodal, demonstrating explicitly the coexistence of phases.For the largest inverse temperatures shown here, trajectories with phase coexistence (trajectories with time-averaged M z between the two phases) are improbable.The broadening observed in Fig. 8(a) can only be seen for inverse temperatures much larger than the ones considered here, where phase coexistence within trajectories can be realised by rare transitions between the two phases.

FIG. 1 .
FIG.1.Sampling of a single spin.The magnetisations (a) ⟨ X⟩ and (b) ⟨ Ẑ⟩ for the Hamiltonian Eq. (23), measured from trajectory sampling (data points) and compared to the exact result (solid lines) for various inverse temperatures β.

FIG. 2 .
FIG.2.Single-spin update scheme for the 1D TFIM.The spin l is updated while all the other spins are kept fixed.Only the spins l − 1 and l + 1 (red) contribute to the dynamics of the l-th spin.These spins change at the times τm (dashed grey lines).The update proceeds through the following steps: (a) The trajectory for spin l is discarded.(b) We then sample the configurations σm (blue circles) at the times in which the effective longitudinal field instantaneously changes, as explained.(c) We then sample bridges between each of the σm under a static longitudinal field to construct the full trajectory (blue line).

FIG. 3 .
FIG.3.Ground state properties of the TFIM.Numerical results using the single-spin update for the 1D TFIM (top)system sizes N = 16 to 1024 and the 2D TFIM (bottom) with N = L 2 and L = 4 to 32.We show the average transverse magnetisation, Mx, in panels (a) and (b), the average two-spin correlator, Mzz, in panels (c) and (d) and the average magnetisation, M |z| , in panels (e) and (f).Data points are the results from trajectory sampling.The dashed black line in the top row of panels is from the iMPS, and from the 2D-DMRG with L = 4 for the bottom row.The value of Jc in 1D is Jc = 1 and Jc ≈ 0.32747 in 2D.

FIG. 4 .
FIG. 4. Convergence in the TFIM.The autocorrelation between trajectories Cµ as a function of J/Jc and µ/N ∈ [1, 7] for (a) the 1D TFIM with N = β = 128 and (b) the 2D TFIM with N = L 2 , L = 32 and β = 128.Panels (c) and (d) show the same data but as a function of µ/N for various values of J/Jc.The black dashed lines show exponential decay.

FIG. 5 .
FIG. 5. Elementary transitions in the TPM.A visualisation of the allowed transitions in the TPM.A plaquette of three spins (spins are shown by blue / red arrows), △, can flip.On the level of the individual plaquette, the plaquette can be in four different disconnected sectors, each containing two configurations and behaving like a two-level system.Within each sector, the system can move between the magnetisations M △ ↔ −M △ .The plaquettes in the green (left) box show plaquettes with total magnetisation M △ = ±3.The three orange boxes (right) show plaquettes with total magnetisation M △ = ±1.The flipping of a neighbouring plaquette, △ ′ , causes a single spin within △ to flip, moving △ into a different sector, as demonstrated by the arrows connecting boxes.If it moves between the M △ = |3| and M △ = |1| sectors, then the sign of the magnetisation M △ is preserved; otherwise, the sign is flipped.

Figure 6 (FIG. 6 .
FIG. 6. Ground state phase transition in the QTPM.Results using TPS with exact plaquette updates and simulated annealing for N = L × L and L = {4, 5, 8, 11, 16, 23, 32}.We show (a) the three spin correlator, Mxxx, (b) the longitudinal magnetization, Mz and (c) the magnetic susceptibility, χz(J) = dMz/dJ.The data points show results from trajectory sampling, and for comparison, the dashed lines show results from DMRG for N = 4 × 4. The insets show the behaviour close to the transition point.

FIG. 8 .
FIG. 8. Large deviations in the long imaginary-time trajectories of the QTPM.The long imaginary-time trajectory statistics can be captured through the LDs.We show the rate function, −φ(Mz)/N , for (a) J = 1.0,(b) J = 0.8 and (c) J = 1.2.The data points show results from trajectory sampling for N = L × L with L = {4, 5, 8, 11, 16, 23, 32}, and the dashed line from 2D DMRG for L = 4.For comparison, the dotted lines show the rate functions for the single spin Eq. (23) case, where g is chosen to fit the mean to the peak of the QTPM distributions.

Fig. 7 (
Fig.7(d).Here, the broadening at criticality is more pronounced; for comparison, we show the rate function for the longitudinal magnetisation of a single spin (dashed line).This broadening suggests a diverging magnetic susceptibility in the large size limit.

Figure 8 FIG. 9 .
FIG.9.Trajectory statistics at finite imaginary times in the QTPM.The estimated distribution of the timeintegrated magnetisation for L = 16 for various finite times (inverse temperatures) at J = 1 from trajectory sampling.The curves are calculated using a Gaussian kernel with δ = 0.005.