Trajectory-Resolved Weiss Fields for Quantum Spin Dynamics

We explore the dynamics of quantum spin systems in two and three dimensions using an exact mapping to classical stochastic processes. In recent work we explored the effectiveness of sampling around the mean field evolution as determined by a stochastically averaged Weiss field. Here, we show that this approach can be significantly extended by sampling around the instantaneous Weiss field associated with each stochastic trajectory taken separately. This trajectory-resolved approach incorporates sample to sample fluctuations and allows for longer simulation times. We demonstrate the utility of this approach for quenches in the two-dimensional and three-dimensional quantum Ising model. We show that the method is particularly advantageous in situations where the average Weiss-field vanishes, but the trajectory-resolved Weiss fields are non-zero. We discuss the connection to the gauge-P phase space approach, where the trajectory-resolved Weiss field can be interpreted as a gauge degree of freedom.


I. INTRODUCTION
Quantum spin systems play a prominent role in the field of many-body physics with diverse applications ranging from quantum magnetism to non-equilibrium dynamics [1][2][3].They also play a crucial role in the development of theoretical techniques ranging from methods of integrability [4][5][6][7][8] to numerical algorithms [9][10][11][12][13][14].A notable challenge is the theoretical description of two and three-dimensional quantum spin systems, where the lack of integrability and the dimension of the Hilbert space stymies progress.The growth of quantum entanglement in real-time dynamics also impedes the description of non-equilibrium phenomena beyond short timescales; this is particularly severe in higher dimensions, as discussed in Refs.[15][16][17][18][19].In addition, tensor network representations are computationally less tractable than in one-dimension, with the number of network contractions scaling exponentially with the system size [20]; this significantly reduces the accessible time-scales.Recent progress has been made using machine learning techniques [14,19,21], and via semi-classical approaches based on the truncated Wigner approximation [22][23][24][25].
Recently, a stochastic approach to quantum spin systems has been developed, based on a Hubbard-Stratonovich decoupling of the exchange interactions [26][27][28][29][30][31][32][33][34].This approach provides an exact reformulation of the quantum dynamics in terms of classical stochastic processes.Quantum expectation values are computed by averaging over independent stochastic trajectories, and the method can be applied in arbitrary dimensions.In recent work we showed that one could extend the accessible simulation times in this approach through the use of an effective Weiss-field [33].Similar conclusions have also been inferred using saddle point techniques [32,34].These approaches allow one to expand around the stochastically averaged time-evolution in order to reduce the need for large stochastic fluctuations.We further discussed [33] the connection to the gauge-P phase space formulation [35], and the possibility to make efficient choices of gauge.This complements a large body of phase space techniques which have emerged in recent years [22][23][24][25][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52].
In this work, we show that the stochastic approach can be significantly enhanced by allowing the Weiss field to be determined on a trajectory by trajectory basis.This allows one to generalize the expansion around a single mean-field trajectory and incorporate sample to sample fluctuations.This trajectory-resolved approach is particularly advantageous in situations where the stochastically averaged Weiss field vanishes.We show that the approach can lead to longer simulation times for a range of quantum quenches in both two and three dimensions.In particular, we demonstrate an exponential improvement in the sampling efficiency.This establishes the technique as a viable tool for the simulation of non-equilibrium quantum spin systems beyond one-dimension.
The layout of this paper is as follows.In Section II we provide an overview of the stochastic approach.In Sections III and IV we discuss the use of trajectory-averaged and trajectory-resolved Weiss fields respectively.In Section V we use these Weiss fields to simulate quantum quenches in the two and three-dimensional quantum Ising model.In Sections VI and VII we discuss the growth of stochastic fluctuations under time-evolution.We conclude in Section VIII and provide an appendix on the associated links [33] to the gauge-P phase space formalism [35,39,41,43,[53][54][55]. We also provide details of our numerical simulations.

II. STOCHASTIC FORMALISM
In this section we recall the principal features of the stochastic approach to quantum spin systems [26][27][28][29][30][31][32][33][34].The method is applicable to generic quadratic spin arXiv:2209.09409v3[cond-mat.str-el]19 Feb 2024 Hamiltonians in arbitrary dimensions: where J ab jk is the exchange interaction between sites j and k and h a j is an applied magnetic field.The spin operators, Ŝa j , obey the canonical commutation relations [ Ŝa j , Ŝb k ] = iϵ abc δ jk Ŝc k , where a, b ∈ {x, y, z} label the spin components, ϵ abc is the antisymmetric symbol and we set ℏ = 1.
The dynamics of the quantum spin system is encoded in the time-evolution operator Û (t) = Te −i t 0 Ĥ(t)dt , where T denotes time-ordering.Decoupling the interactions by means of a Hubbard-Stratonovich transformation [56] over auxiliary fields φ j one obtains where Φ a j = 1 √ i φ a j + h a j ∈ C and the integration is performed over all paths with Dφ = ja Dφ a j .The parameter Φ a j plays the role of an effective, complex magnetic field.In writing Eq. ( 2) we define the so-called noise action [30,[32][33][34], which allows one to regard the fields φ as correlated random noises with the Gaussian measure Dφ e −S [φ] .The time-evolution operator (2) can be recast as where Ĥs ≡ − ja Φ a j Ŝa j is referred to as the stochastic Hamiltonian and ⟨...⟩ φ denotes averaging over the Gaussian noise variables.Eq. ( 4) motivates the introduction of the stochastic evolution operator Û s (t) = Te −i t 0 Ĥs (t ′ )dt ′ and the stochastic state |ψ s (t)⟩ = Û s (t)|ψ s (0)⟩.As can be seen from Eq. ( 3), the exchange interactions enter the representation via the correlations of the noise.For both analytical and numerical calculations it is convenient to convert these noises to Gaussian white noise by diagonalizing the noise action (3) [28][29][30][31][32][33][34].Specifically, we decompose the original fields φ a j in terms of white noise variables, ϕ a k , as φ a j = kb O ab jk ϕ b k , where O T J −1 O = 1 and the bold symbols indicate matrices.The time-evolution operator (4) can be further simplified by means of a so-called disentangling transformation [28].This eliminates the time-ordering operation by introducing a new set of variables ξ a j : The variables ξ satisfy the stochastic differential equations (SDEs) [28][29][30][31][32][33] − i ξ+ where Φ ± j = 1 2 (Φ x j ∓ iΦ y j ) and ξ a j (0) = 0. Time-evolution can therefore be achieved by solving these SDEs numerically.In order to evaluate generic local observables ⟨ Ô (t)⟩ = ⟨ψ(0)| Û † Ô Û |ψ(0)⟩, one may decouple the forwards and backwards time-evolution operators independently [29].Expectation values therefore reduce to averages over these classical stochastic variables: where φ and φ correspond to the forwards and backwards evolution respectively as indicated by the subscripts on the evolution operators.Employing the disentangling transformation given in (5) one may recast this in the form ⟨ Ô(t)⟩ = ⟨f (ξ, ξ)⟩ φ, φ, where the function f (ξ, ξ) depends on the operator Ô. Evaluating this expression involves solving the SDEs (6) and computing the classical average over realizations of the stochastic process.
In general, the variables ξ grow without bound as they approach coordinate singularities, leading to a failure of numerical integration schemes.This can be seen by considering the action of Û s (t) on an initial down-state |↓⟩ : It is evident that the up-state |↑⟩ is associated with the diverging quantity |ξ + j (t)| → ∞.As discussed in Ref. [31], these singularities can be eliminated by a suitable parameterization of the Bloch sphere.Specifically, one introduces a second coordinate patch ξ for the Bloch sphere, which instead has the coordinate singularity at the down-state |↓⟩.Singularities can therefore be avoided by mapping between the coordinate patches whenever the evolution crosses the equator of the Bloch sphere, associated with |ξ + j | = 1.Note that by using the parametrization (8), the variable ξ − in (6c) is not required.Since any spin state can be obtained as a rotation from a downstate |↓⟩, the parametrization (8) can be used generically, provided one introduces a state preparation protocol.For example, one may rotate the down-state |↓⟩ to the initial state |ψ j (0)⟩ = Û s j (0, −δ)|↓⟩, over some arbitrary timeinterval δ before t = 0.In practice, this evolution need not be computed since we are only interested in calculating the time-evolution from t = 0; the initial conditions on ξ are constrained by the initial state according to Û s j (ξ j (0))|↓⟩ = |ψ j (0)⟩.Entangled initial states can be treated by introducing a probability distribution over these initial conditions [31].Further extensions also exist for combining the SDEs (6) with matrix product states [31].Although this results in improvements in 1D, we do ).We consider a quantum quench in the two-dimensional quantum Ising model from the disordered to the ordered phase.We examine a 3 × 3 array of spins which are initialized in the state |+⟩ which is fully-polarized along the positive x-direction, and time-evolve with Γ/J = 0.1.The spins evolve rapidly over the Bloch sphere and gather at the poles, as highlighted by the bright regions.It can be seen that the trajectory-averaged position ⟨n z ⟩φ (dashed line) doesn't approximate the trajectory dynamics.This motivates the use of trajectory-resolved Weiss fields in the stochastic approach.The data correspond to all of the spins with N = 1000 stochastic samples.
not pursue this here; for 2D and 3D systems the number of tensor network contractions scales exponentially with the system size [20] and no advantage is gained.Instead, we turn our attention to the use of Weiss fields.

III. TRAJECTORY-AVERAGED WEISS FIELD
In recent work in both Euclidean [32] and real-time evolution [33,34] it has been noted that the sampling efficiency can be improved by shifting the Hubbard-Stratonovich fields by a constant at each time slice t: φ(t) → φ(t) + ∆φ(t).As discussed in Ref. [33] it is convenient to parameterize this as ∆φ a j (t) = √ i kb J ab jk m b k (t).This introduces a term of the form kb J ab jk m b k (t) Ŝa j into the stochastic Hamiltonian, which can be interpreted as a Weiss field: where Î is the identity operator.In writing (9) we have included Î-dependent terms in the definition of the stochastic Hamiltonian, to ensure that the Gaussian measure is still in the original form (3). The Weiss field allows one to sample the noise fluctuations around a single deterministic 'mean-field' trajectory determined by m a j (t).This can be selected to reduce stochastic fluctuations.In previous work [33], we set this equal to the trajectory-averaged Weiss field m a j = ⟨n a j ⟩ φ , where is the expectation value of the a-component of the spin.In writing (10), we include the normalization factor |ψ s (t)| 2 since the stochastic state |ψ s (t)⟩ is unnormalized.Eq. ( 10) can be calculated self-consistently from a relatively small number of trajectories [33]; we use four iterations of O( 103 ) samples.The stochastic sampling is then carried out around the deterministic trajectory determined by the first two terms in Eq. ( 9).The trajectories are re-weighted via the non-Hermitian term in (9).This procedure therefore corresponds to a form of importance sampling [32,57].As discussed above, without loss of generality we may consider the SDEs starting from an initial spin-down state: where the Î-dependent terms in (9) enter into the evolution of ξ z j [33] and the effective magnetic field becomes It can be seen that this consists of the applied magnetic fields h a j , the Hubbard-Stratonovich field φ a j , and the Weiss field contribution.At this stage, we note that the O(m 2 ) terms in ( 9) and ( 12) may be safely ignored since they result in a deterministic phase for the stochastic state |ψ s (t)⟩.
While the trajectory-averaged Weiss field has been shown to improve simulation times for a range of quenches [33], it decays to zero over time; the trajectories spread out over the Bloch sphere.The situation can be summarized by considering a quench of the 2D quantum Ising model with nearest neighbor interactions where we set J = 1 and use periodic boundary conditions.In Fig. 1 we show the distribution of n z j over time following a quench from the disordered to the ordered phase.Over time the trajectories spread out over the Bloch sphere, accumulating at either pole.However, the mean behavior ⟨n z j ⟩ φ remains zero and fails to approximate the trajectories.Large non-Hermitian fluctuations are therefore required, since the noise must drive the state towards the poles.In Section IV we address this issue by allowing the Weiss field to vary on a trajectory by trajectory basis.In Appendices A and C we give two independent derivations of the approach.

IV. TRAJECTORY-RESOLVED WEISS FIELD
In order to treat problems where the trajectoryaveraged Weiss field vanishes, we consider the instantaneous Weiss field for each trajectory taken separately.We restrict our attention to the quantum Ising model but discuss the application to more general models of the form (1) in the Appendices.Specifically, we set m z j (t) = n z j (t) where n z j (t) is the instantaneous value of the z-component of the spin for a single trajectory.In this approach, each trajectory develops a different Weiss field due to the effect of the accumulated noise.This may be regarded as a form of feedback [50][51][52], in which quantum fluctuations recailibrate the reference trajectory.The trajectory-resolved Weiss fields do not necessarily decay to zero at long times.In Section VII, we show that the use of these Weiss fields allows samples to contribute more equally to quantum expectation values.The trajectory dependence of the Weiss field also changes the effect of 1 2 ij J ij m i m j Î in (9) on the stochastic state |ψ s (t)⟩; it goes from a removable deterministic phase to a stochastic phase that determines how trajectories 'interfere'.
For simulations it is convenient to introduce a small time-delay δ into the trajectory-resolved Weiss field so that m z j (t) = n z j (t − δ).This enhances the numerical stability without approximation; see Appendix B. Throughout this work we solve the SDEs ( 11)-( 12) using the Stratonovich-Heun predictor-corrector scheme.In simulations with a trajectory-resolved Weiss field we set δ = ∆, where ∆ is the time-step, unless stated otherwise.In general, we find that a time-step of ∆ = 0.01 is sufficient for most of our simulations, unless the observable in question involves very small quantities, such as the Loschmidt amplitude.The use of a different time-step is highlighted in each instance.

V. SIMULATIONS
In this section we demonstrate the improvements for numerical simulations when using trajectory-resolved Weiss fields.For simplicity, we focus on quantum quenches of the nearest neighbor quantum Ising model ( 14) in both two and three dimensions, although improvements can also be seen in one-dimension.In Fig. 2 we show results for the Loschmidt rate function following a quantum quench in two dimensions from an initial state with Γ = 0 to Γ = 8J.Explicitly, we consider the initial state |ψ(0)⟩ = 1 √ 2 |⇓⟩ + |⇑⟩ , which is the superposition of the symmetry broken groundstates: 2 |⇓⟩ + |⇑⟩ corresponding to the superposition of degenerate ground states and time-evolved with Γ/J = 8, thereby quenching across the quantum critical point at Γ/J ≈ 1.52 [58].The results obtained using a trajectory-resolved Weiss field (orange) improve upon those obtained by the trajectory-averaged Weiss field (blue).For comparison, the results obtained using QuSpin's ODE Solver [59] on a 5 × 5 lattice are shown (dotted) in both panels; comparison results for a 7×7 lattice are currently beyond reach.The SDE results are obtained from 5 batches of N = 10 5 stochastic samples with a time-step ∆ = 0.001.The latter is chosen to accurately resolve the Loschmidt peaks.The faded lines indicate the standard error of the mean.
In Figs. 2 (a) and 2 (b) we show the results for a 5 × 5 and a 7 × 7 lattice respectively.It can be seen from Fig. 2 (a) that the trajectory-resolved approach allows simulations to be carried out for longer time durations, before departures arise.For comparison, we show results obtained by QuSpin's ODE solver [59] for a 5 × 5 lattice.The results are in excellent agreement until t ∼ 3/J, after which the stochastic fluctuations are not well sampled.Similar behavior is evident in Fig. 2 (b), although we can only compare to the QuSpin results for the smaller system size, as 7 × 7 is not possible at present.
In Fig. 3 we show similar results for a 3D lattice with 125 sites on a 5 × 5 × 5 grid.This is a far more challenging problem due to the increase in dimensonality.The results for the trajectory-resolved approach show clear Loschmidt peaks, which extend beyond those obtained by the trajectory-averaged approach.This complements earlier investigations [60] which were unable to resolve the sharp non-analyticities due to significant finite-size effects.
Having established the utility of trajectory-resolved 2 |⇓⟩ + |⇑⟩ corresponding to the superposition of degenerate ground states and time-evolved with Γ/J = 8, thereby quenching across the quantum critical point at Γ/J ≈ 2.58 [58].The results obtained using a trajectoryresolved Weiss field (orange) improve upon those obtained by the trajectory-averaged Weiss field (blue).The SDE results are obtained using 5 batches of N = 10 5 stochastic samples and a time-step of ∆ = 0.001.The latter is chosen to resolve the Loschmidt peaks.The faded lines indicate the standard error of the mean.
Weiss fields, we now examine the differences in performance for quenches that have significant trajectoryaveraged Weiss fields, and quenches that don't.In Fig. 4 (a) we show results for the transverse magnetization M x = 1 N j ⟨ Ŝx j ⟩ following a quantum quench within the ordered phase of the 2D quantum Ising model ( 14) on a 7 × 7 square lattice with 49 sites.We start from the fully-polarized state |⇓⟩ with all spins down, and timeevolve using the Hamiltonian with Γ/J = 0.2.As can be seen in Fig. 4 (c), this quench is associated with a nonzero trajectory-averaged Weiss field over the duration of the simulation, corresponding to a non-vanishing mean field in the initial state.It is evident from panel (a) that the trajectory-resolved Weiss fields perform better than the trajectory-averaged ones, although the relative gains from their use is modest in this case.This can also be seen from the behavior of the norm of the quantum state, |ψ(t)| 2 , which remains close to unity when the stochastic fluctuations are adequately sampled [31,33].In the Appendices, we provide data for a quench in 3D within the ordered phase.The results demonstrate similar performance improvements when using the trajectory-resolved Weiss field.
The difference between the two approaches is noticebly greater for quantum quenches without a significant trajectory-averaged Weiss field.This is illustrated in Fig. 4 (b) for a quantum quench in the 2D quantum Ising model from the disordered phase to the ordered phase.Specifically, the system is initialized in the fully-polarized state in the x-direction, |+⟩ = resolved Weiss field is approximately double that of the trajectory-averaged approach.For this particular quench the stochastic trajectories rapidly spread out over the Bloch-sphere and the trajectory-averaged Weiss field remains close to zero throughout the evolution.As shown in Fig. 4 (c), the use of a trajectory-averaged Weiss field offers little advantage over the case with m z j = 0.In contrast to the situation in equilibrium, where the utility of mean-fields increases with dimensionality, their use for dynamics is more subtle.In particular, the utility of Weiss fields can depend on the details of the quantum quench, and not just the dimensionality.The utility of the trajectory-resolved approach can also be seen in the scaling of the accessible timescale for numerical simulations with the number of samples required.Following Refs.[31,33] we define the breakdown time t b as the time when the normalization |ψ(t)| 2 differs from unity by 10%.As discussed in Refs.[31][32][33], the number of samples N scales exponentially with the breakdown time, corresponding to the onset of strong fluctuations.Specifically, N ∼ ce αt b where the growth exponent α depends on the details of the quench.In Fig. 5 we compare the scaling of the trajectory-averaged approach with the trajectory-resolved approach, for a quench in the twodimensional quantum Ising model with a 3 × 3 array of spins.It can be seen that the growth-exponent α is significantly reduced for the trajectory-resolved case, in comparison with the trajectory-averaged case.In Sections VI and VII we show that this is related to a reduction in the fluctuations of the normalization of the stochastic state.

VI. STOCHASTIC STATE NORMALIZATION
In this section we discuss how the normalization of the stochastic state determines the sampling efficiency of the approach.To see this we note that an arbitrary normalized state |ψ(t)⟩ can be expressed as the average over normalized stochastic states, ||ψ s (t)⟩, according to for a single trajectory following a quench of the 2D quantum Ising model with a 3 × 3 array of spins.The initial state |ψ(0)⟩ = |+⟩ is evolved with Γ/J = 0.3.We compare a direct numerical evaluation of |ψ s (t)| 2 (solid line) against the theoretical prediction (18) (dashed line).For the latter we use the value of n z j (t) obtained from the numerical procedure as an input.The data is obtained using the Stratonovich-Heun scheme [61,62] without time-delay and a time-step of ∆ = 0.001.The growth of the stochastic state normalization is required in order to maintain the overall normalization of the quantum state.
where W (t) = |⟨ψ s (t)|ψ s (t)⟩| 1/2 is the norm of the stochastic state.It can be seen that W (t) corresponds to the weight of each sample in the ensemble.A large spread directly inhibits sampling.As we show in Appendix A 3, the norm of the stochastic state grows monotonically in time.For the quantum Ising model this is given by where ν j = k |O zz jk | 2 ; see Appendix A 3. We verify (18) in Fig. 6 for a single stochastic trajectory following a quench in the 2D quantum Ising model.It can be seen from ( 18) that the normalization is controlled by deviations from the fully-polarized spin state.The growth of the norm with time is required in order to maintain the overall normalization of the quantum state, |ψ(t)| 2 = 1.
To see this we note that where r and r ′ are independent sample indices.If |ψ s r (t)⟩ was normalized then the overlaps in (19) would be less than or equal to unity.As such, |ψ(t)| 2 < 1.It follows that the normalization of |ψ s r (t)⟩ must grow with time in order to maintain the condition that |ψ(t)| 2 = 1.This reflects the independent decouplings of the forwards and backward time-evolution in the stochastic approach.The variance grows linearly following an initial transient.The growth rate for the trajectory-resolved case is lower than the other two cases, which again aids sampling.The variance for the trajectory-averaged case is only comparable to the trajectoryresolved case at very short timescales, when the state is wellapproximated by fluctuations around a product state.Inset: time-evolution of the mean of the distribution ⟨γ⟩ ϕ .The linear growth rate is similar for all three cases.

VII. GROWTH OF FLUCTUATIONS
Having demonstrated that the norm of the stochastic state grows with time, we now examine its distribution.Given the exponential scaling of W with time it is convenient to consider the distribution of γ = ln W .In Fig. 7 (a) we show the distribution P (γ) at a fixed time for different values of the Weiss field.It can be seen that P (γ) is normally distributed, as indicated by the solid lines.The use of a trajectory-resolved Weiss field results in a narrower distribution than the other cases.In par-ticular, it leads to a reduction in the extremal values of W = e γ which contribute the most to ensemble averages.This leads to an improvement in the sampling efficiency.In Fig. 7 (b) we show the growth of the variance var(γ) as a function of time.It can be seen that the growth rate is reduced when using the trajectory-resolved Weiss field, even at late times.In contrast, the growth rate for the trajectory-averaged Weiss field eventually follows the case without a Weiss field; this was also observed in Ref. [32], by sampling around a saddle-point trajectory in Euclidean time.In the inset of Fig. 7 (b) we show that the mean of the distribution ⟨γ⟩ ϕ changes very little with the choice of Weiss field.The improvements in the scaling are therefore attributed to the reduction of the variance with a trajectory-resolved Weiss field.Although the focus of this work is on 2D and 3D systems, similar improvements can also be seen in 1D, as shown in Appendix D.
In closing this section, we briefly comment on the role of fluctuations on the success of the trajectory-resolved Weiss field.Over timescales in which the dynamics can be approximated by fluctuations around a well-chosen product state trajectory, both the trajectory-averaged and the trajectory-resolved Weiss fields can efficiently encode the evolution.However, a key advantage of the trajectory-resolved approach, is that it can remain efficient beyond this timescale.In the trajectory-resolved approach, an entangled superposition such as a triplet state 1 √ 2 (|↑↓⟩ + |↓↑⟩) can be obtained from two separate trajectories, |↑↓⟩ and |↓↑⟩, with each encoding its own Weiss field.In contrast, there is no single product state that approximates this triplet state.The trajectoryresolved approach therefore has a notable advantage for simulating quantum dynamics.

VIII. CONCLUSION
In this work we have investigated the real-time dynamics of quantum spin systems in two and three-dimensions by means of an exact stochastic approach.We have shown that the use of a trajectory-resolved Weiss field can significantly extend the accessible simulation times, with an exponential improvement in the sampling efficiency.We have illustrated the utility of this approach for exploring dynamical quantum phase transitions in two and three dimensions, although the applicability is broader.Our results address a critical shortage of exact simulation techniques for non-equilibrium quantum spin systems in two and three dimensions.It would be interesting to see if trajectory-resolved Weiss fields can be used in other contexts, for example in situations which traditionally involve expanding around a single mean-field configuration.for a single trajectory using the Ito Euler-Maruyama (EM), Stratonovich-Heun (SH) and delayed Stratonovich-Heun (SHδ) schemes.We consider the set-up as (a) but with Γ/J = 0.1 and a small time-step of ∆ = 10 −5 .The results for the different integration schemes are almost indistinguishable, which demonstrates that the time-delayed result converges to the undelayed result in the limit of a small delay δ.In this simulation we use complex noises residing on each bond; see [33] for a discussion.
numerical stability at large time-steps, as afforded by other Stratonovich schemes.For times t < δ we set m z (t − δ) equal to the initial magnetization.
To perform the numerical integration of ( 11) and ( 12) for the quantum Ising model in the presence of a delayed trajectory-resolved Weiss field we employ the predictorcorrector scheme used in Ref. [64].For clarity, we restate the SDEs from the main text in the presence of a trajectory-resolved Weiss field with delay:

FIG. 1 .
FIG.1.Time-evolution of the distribution of stochastic trajectories for the z-component of spin P (n z ).We consider a quantum quench in the two-dimensional quantum Ising model from the disordered to the ordered phase.We examine a 3 × 3 array of spins which are initialized in the state |+⟩ which is fully-polarized along the positive x-direction, and time-evolve with Γ/J = 0.1.The spins evolve rapidly over the Bloch sphere and gather at the poles, as highlighted by the bright regions.It can be seen that the trajectory-averaged position ⟨n z ⟩φ (dashed line) doesn't approximate the trajectory dynamics.This motivates the use of trajectory-resolved Weiss fields in the stochastic approach.The data correspond to all of the spins with N = 1000 stochastic samples.

N j=1 1 √ 2 (FIG. 4 .
FIG. 4. Time-evolution of the transverse magnetization M x (t) following quenches of the 2D quantum Ising model for a 7 × 7 lattice with 49 spins.(a) The spins are initialized in the fully-polarized state along the z-axis |⇓⟩ and time-evolved with Γ/J = 0.2.The results obtained with the trajectoryresolved Weiss field (orange) and the trajectory-averaged Weiss field (blue) are in very good agreement.The normalization |ψ(t)| 2 indicates that the former reaches slightly longer timescales.The success of the latter is consistent with the presence of a non-vanishing trajectory-averaged Weiss field as shown in panel (c).(b) The spins are prepared in the fully-polarized state |+⟩ along the x-axis and time-evolved with Γ/J = 0.1.In this case, the trajectory-resolved Weiss field performs much better than the trajectory-averaged Weiss field.This is consistent with a vanishing trajectory-averaged Weiss field as shown in panel (c).The results are in very good agreement until the norm of the quantum state |ψ(t)| 2 departs from unity.The results in (a) correspond to 5 batches of N = 2 × 10 5 stochastic samples while the results in (b) correspond to 5 batches of N = 6 × 10 5 stochastic samples, all obtained with a time-step ∆ = 0.01.The faded lines indicate the standard error of the mean.

FIG. 5 .
FIG. 5. Exponential scaling of the required number of samples N versus the breakdown time t b of the simulations: N ∼ ce αt b .The data correspond to quantum quenches in the 2D quantum Ising model for a 3 × 3 array of spins.The spins are initialized in the fully-polarized state |+⟩ along the x-axis and are time-evolved with Γ/J = 0.1.The trajectory-resolved Weiss field (circles) is shown to require less samples than the trajectory-averaged Weiss field (diamonds) to reach a given time.This is confirmed by the coefficient α of the linear fit, which suggests a reduction in the exponent by a factor of approximately 2.7 for this simulation.Each of the data points correspond to the mean of 10 batches of simulations of N runs with the standard error indicated by a bar.

FIG. 6 .
FIG.6.Growth of the stochastic state normalization |ψ s (t)| 2 for a single trajectory following a quench of the 2D quantum Ising model with a 3 × 3 array of spins.The initial state |ψ(0)⟩ = |+⟩ is evolved with Γ/J = 0.3.We compare a direct numerical evaluation of |ψ s (t)| 2 (solid line) against the theoretical prediction (18) (dashed line).For the latter we use the value of n z j (t) obtained from the numerical procedure as an input.The data is obtained using the Stratonovich-Heun scheme[61,62] without time-delay and a time-step of ∆ = 0.001.The growth of the stochastic state normalization is required in order to maintain the overall normalization of the quantum state.

FIG. 7 .
FIG. 7. (a)Distribution P (γ) at a fixed time Jt = 20 following a quench in the two-dimensional quantum Ising model for a 3 × 3 lattice of spins.We start in the fully-polarized state |⇓⟩ and quench to Γ/J = 0.7, taking N = 5 × 10 4 stochastic samples.The use of a trajectory-resolved Weiss field (orange) results in a narrower distribution than the trajectoryaveraged case (blue) and the case without a Weiss field (red); this makes the dynamics easier to sample.The distributions are approximately normal, as indicated by the solid line fits.(b) Time-evolution of the variance var(γ).The variance grows linearly following an initial transient.The growth rate for the trajectory-resolved case is lower than the other two cases, which again aids sampling.The variance for the trajectory-averaged case is only comparable to the trajectoryresolved case at very short timescales, when the state is wellapproximated by fluctuations around a product state.Inset: time-evolution of the mean of the distribution ⟨γ⟩ ϕ .The linear growth rate is similar for all three cases.

FIG. 8 .
FIG. 8. (a) Transverse magnetization M x (t) following a quench of the quantum Ising model from N i |+⟩i to Γ/J = 0.3 with N = 3.The data correspond to a trajectory-resolved Weiss-field m z i = n z i (t−δ) with a large delay δ = 0.2 (orange).The results are in excellent agreement with exact diagonalization (ED) until the breakdown time t b associated with finite sampling.We use 10 batches of N = 10 5 samples, and a time-step of ∆ = 0.01.The solid line corresponds to the mean and the light orange lines indicate the standard error.(b) Time-evolution of the real parts of ξ + 1 and ξ z 1

FIG. 10 .
FIG. 10.(a) P (γ) at a fixed time Jt = 20 following a quench in the one-dimensional quantum Ising model with N = 9.We start in the fully-polarized state |⇓⟩ and quench to Γ/J = 0.3, taking N = 5 × 10 4 stochastic samples.The use of a trajectory-resolved Weiss field (orange) results in a narrower distribution than the trajectory-averaged case (blue) and the case without a Weiss field (red).The distributions are approximately normal, as indicated by the solid line fits.(b) Time-evolution of the variance var(γ).The variance grows linearly following an initial transient.The growth rate for the trajectory-resolved case is lower than the other two cases.(c) Time-evolution of the mean of the distribution ⟨γ⟩ ϕ .The linear growth rate is similar for all three cases.
IX. ACKNOWLEDGEMENTS SEB acknowledges support from the EPSRC CDT in Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES) via grant number EP/L015854/1, as well as the support of the Young Scientist Training program at the Asia Pacific Center for Theoretical Physics (APCTP).MJB acknowledges support of the London Mathematical Laboratory.AGG acknowledges EPSRC grant EP/S005021/1.We are grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EP-SRC (EP/P020194/1 and EP/T022213/1).For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.The data supporting this article is openly available from the King's College London research data repository, KORDS, at https://doi.org/10.18742/24438967.