Importance Sampling Scheme for the Stochastic Simulation of Quantum Spin Dynamics

The numerical simulation of dynamical phenomena in interacting quantum systems is a notoriously hard problem. Although a number of promising numerical methods exist, they often have limited applicability due to the growth of entanglement or the presence of the so-called sign problem. In this work, we develop an importance sampling scheme for the simulation of quantum spin dynamics, building on a recent approach mapping quantum spin systems to classical stochastic processes. The importance sampling scheme is based on identifying the classical trajectory that yields the largest contribution to a given quantum observable. An exact transformation is then carried out to preferentially sample trajectories that are close to the dominant one. We demonstrate that this approach is capable of reducing the temporal growth of fluctuations in the stochastic quantities, thus extending the range of accessible times and system sizes compared to direct sampling. We discuss advantages and limitations of the proposed approach, outlining directions for further developments.


I. INTRODUCTION
Experimental breakthroughs in the simulation of isolated many-body quantum systems [1,2] have led to great theoretical interest in their far-from-equilibrium dynamics [3]. Concepts such as the thermalization of isolated quantum systems [4,5], or the absence thereof [6,7], and the discovery of novel non-equilibrium phenomena [8][9][10], have been the subject of intense experimental and theoretical exploration. Great progress has been achieved for one-dimensional systems, which include analytically solvable integrable models [11] and are often amenable to efficient numerical treatment via tensor-network based approaches [12][13][14]. However, the limitations of existing techniques call for the development of additional analytical and numerical tools to describe non-equilibrium quantum dynamics. This is particularly important in higher-dimensional settings, where no exact solutions are generally available and the applicability of tensor network methods is limited [15,16]. A number of directions are currently being explored, including linked cluster expansions [17] and neural network approaches [18][19][20].
An alternative technique, recently applied to manybody quantum spin systems, consists in exactly mapping unitary quantum dynamics to an ensemble of classical stochastic processes [21][22][23][24][25]. This approach is based on the disentanglement formalism [21][22][23]26], which provides an exact functional integral representation of the time-evolution operator. In this formalism, interactions are decoupled by means of Hubbard-Stratonovich (HS) transformations [27,28], and the resulting single-spin dynamics is then parameterized in terms of a set of classical disentangling variables, defined by a Lie-algebraic transformation [21-23, 26, 29, 30]. Quantum expectation values are then obtained as classical averages over ensembles of stochastic trajectories of the disentangling variables. Early investigations have shown that the stochastic approach is immediately applicable to higherdimensional settings, but its practical performance is limited by the exponential growth of fluctuations in the stochastic quantities as a function of time [24]. In recent work, the disentanglement formalism was applied in imaginary time, providing an analytical and numerical framework to study the ground states of quantum spin systems [31]. In this context, the identification of the saddle point trajectory, which provides the dominant contribution to observables, was used to perform an exact measure transformation. This resulted in an importance sampling scheme which greatly improves the performance of the numerical stochastic approach.
In this manuscript, we generalize the importance sampling scheme to real-time evolution. We begin by briefly recapping the main aspects of the disentanglement formalism in Section II. The real-time importance sampling scheme in then introduced in Section III, emphasizing similarities and differences with the imaginary-time case. The approach is first applied to local observables in Section IV, comparing it to direct sampling and discussing the role of fluctuations. Return probabilities are then considered in Section V in the context of dynamical quantum phase transitions [8,32]. We conclude in Section VI, summarizing our findings and discussing directions for further research.

II. DISENTANGLEMENT FORMALISM
The dynamics of a quantum state |ψ 0 under the action of a HamiltonianĤ is encoded in the time-evolution op- , where we set = 1 and Ì denotes time ordering: |ψ(t) =Û (t)|ψ 0 . We consider a generic quadratic spin Hamiltonian where the indices j, k correspond to lattice sites and a, b, run over the generators of SU (2). The external fields h a i and the interaction matrix J ab jk can in general be timedependent and no specific boundary conditions are assumed. The constant J is an overall coupling strength.
The time-evolution operator corresponding to the Hamiltonian (1) is represented by a matrix whose size grows exponentially with the system size N , forbidding its direct evaluation for large systems. However, this issue can be circumvented by means of the disentanglement formalism [21-24, 26, 31]. In this approach,Û is exactly written as a functional integral featuring matrices whose size is determined by the local Hilbert space dimension, e.g. they are 2 × 2 matrices for spin-1/2 systems. Below we recap the main features of this approach; additional details can be found in Ref. [24]. For simplicity, we consider systems initialized in a product state |ψ 0 = ⊗ i |ψ 0 i ; more general initial conditions can also be treated within the same formalism, e.g. by first performing imaginary-time evolution from a product state and subsequently evolving in real time. The time-evolution operatorÛ admits the exact functional integral representation [21][22][23] where the noise action S 0 is given by and the disentangling variables ξ a i (t) satisfy the equations of motion [22] with Φ a i ≡ h a j + Jϕ a j and initial conditions ξ a i = 0. The system of equations (4) encodes the details of the system at hand. The fields ϕ a i are in general complex valued. The noise action (3) can be diagonalized by a linear transformation ij and (J −1 ) ab ij are treated as matrices by grouping the (a, i) and (b, j) indices. This yields [22][23][24] The real-valuedness of the fields φ a j is required for convergence of the integral (2) and can be viewed as defining the appropriate integration lines for the complex-valued fields ϕ a j [24]. Due to the Gaussian action (5), the functional integral in Eq. (2) can be seen as an average over realizations of Gaussian white noise variables φ [22], Eqs (4) can then be interpreted as stochastic differential equations (SDEs) [21][22][23]. Eq. (6) makes it possible to establish an exact map between quantum expectation values and classical averages. Time-evolved expectation values are given by O(t) = Û † (t)ÔÛ (t) , where . . . denotes the expectation value with respect to a chosen initial state. Quantum observables can thus be cast in functional integral form by expressing their time dependence in terms of time-evolution operatorsÛ (t) and substituting the representation (6) [23,24]. In order to disentangle two time-evolution operators, as in the case of local expectation values, one can introduce independent HS fields φ x = {φ a x,i } and the associated disentangling variables ξ x = {ξ a x,i }, with x ∈ {f, b} denoting forwards and backwards evolution. Eq. (6) features a product of single-site operators, whose action on product states can be straightforwardly evaluated; this yields a classi- Local observables correspond to classical functions of the form where F ½ (t) φ gives the norm of the state andF O is a product of a finite number of terms, each featuring the disentangling variables relative to a single site [31]. For instance, for spin-1/2 systems the spin operators are represented by the Pauli matrices,Ŝ a i = σ a i /2; inserting this in Eq. (6) one readily finds that for an initial ⊗ i |↓ i state the normalization function is given by while the on-site longitudinal magnetization M z Expressions such as (9) are easily obtained for any physical observable and take the same form for different models, in real or imaginary time; see Refs [24,31]. Different initial conditions correspond to different F O [23,24], or, alternatively, can be encoded in the initial conditions of the disentangling variables [25,31]. In general, any time-evolving quantity can be expressed within the disentanglement approach. However, the approach is best suited to quantities that are readily expressed in terms of time-evolution operators, since these are the objects that are replaced by their disentangled counterparts (2). This is the case of the local observables discussed above.
Quantum expectation values O(t) can be numerically computed by averaging the classical functions F O over realization of the stochastic processes φ(t). As demonstrated in Refs [24,25], the numerical evaluation of such averages is stymied by the exponential growth of fluctuations in F O with time and the system size. In the imaginary-time case, it was recently shown that the growth of fluctuations can be greatly suppressed by applying an importance sampling method whereby, when randomly generating trajectories, strongly-contributing ones are sampled preferentially [31]. In contrast, when sampling according to the original measure (5) one mostly draws trajectories that are nearly non-interacting and give little contribution. Below we generalize the importance sampling approach to real-time evolution, showing that it leads to a significant reduction in fluctuations compared to the direct sampling approach of Refs [23][24][25].

III. IMPORTANCE SAMPLING
For analytical computations, it is convenient to work with the ϕ fields. For a given observable O, we seek to identify the saddle point (SP) trajectory ϕ SP yielding the largest contribution to the corresponding functional integral. This is obtained by extremizing the effective The action S O features the fields ξ a i via the function F O ; these are themselves functionals of ϕ, such that Euler-Lagrange equations cannot be derived. Instead, the SP equation is obtained by direct extremization of S O : This condition yields a functional integral equation for ϕ SP , which can be solved recursively. One can then perform an exact measure transformation such that trajectories around the saddle point configuration are sampled preferentially, as illustrated in Fig. 1. This is carried out by performing the change of variables φ(t) → φ(t) + φ SP (t) in the functional integral (10) [31], where φ SP is readily obtained from ϕ SP as indicated in Section II. This transformation does not constitute a saddle point approximation: a change of measure does not truncate fluctuations, so that the resulting expressions are still formally exact. Due to the Gaussianity of the action (5), the importance sampling method then amounts to numerically sampling a modified functional x,i , and the index x runs over all sets of disentangling fields, e.g. forwards and backwards for local observables.
Eq. (12) can be evaluated numerically in order to compute quantum expectation values, yielding an importance sampling scheme for the stochastic approach. In practice, one generates an ensemble of stochastic trajectories ξ a i whose time evolution is determined via the SDEs (4) with the modified field φ(t) → φ(t) + φ SP (t). To the best of our current knowledge, the SDEs (4) can only be Direct Sampling Importance Sampling FIG. 1. Illustration of the importance sampling approach. When directly sampling according to the measure (5), one predominantly generates trajectories around φ = 0 (dashed black line), which might however carry a small contribution to a given observable. In contrast, in the importance sampling approach one preferentially generates trajectories in the vicinity of the saddle point trajectory (dashed blue line), which carries the largest contribution. The solid black and blue lines respectively represent stochastic trajectories generated according to direct and importance sampling.
solved analytically in certain special cases, such as noninteracting systems or Hamiltonians made up of commuting terms [24]. Thus, a discrete-time numerical method is generally needed in order to solve (4). Different numerical integration schemes have been previously used to this end, including the Euler-Maruyama [23,24,33] and the stochastic Heun [25,33,34] schemes. Here, unless otherwise stated, we use the explicit strong order-1 scheme [33,35] with time-step ∆t = 0.01, which was found to perform comparatively well. The numerical time evolution of the disentangling variables ξ a i is known to give rise to divergences whereby |ξ + i (t)| → ∞ at finite t [23]; this issue can be avoided by means of a suitable reparameterization of the disentangling variables [25,36]. We explicitly normalize observables by F ½ φ , which improves the final accuracy [25]. For results obtained as classical averages of stochastic quantities, we estimate error bars as σ/ √ n B , where σ is the standard deviation over n B = 5 batches of independent simulations unless otherwise specified.
For definiteness, we will illustrate the importance sampling approach by considering the quantum Ising model in D spatial dimensions. For a system with N = N 1 × · · · × N D sites, this is given bŷ where we use D-dimensional spatial indices i = (i 1 , · · · , i D ) with i k ∈ (1, · · · , N k ) and ij denotes pairs of nearest neighbors 1 . We consider periodic boundary conditions and ferromagnetic (FM) interactions, J > 0. When D = 1, the model (13) reduces to the quantum Ising chain, and for h = 0 it can be solved exactly in terms of free fermions, harboring a quantum phase transition (QPT) at Γ c = J/2 in the present units [37]. The Hamiltonian (13) is encoded in the SDEs (4) with [22,23]. In our numerical results below we set J = 1 and consider D ∈ {1, 2}, illustrating the applicability of the importance sampling method to higher-dimensional systems.

IV. LOCAL OBSERVABLES
In general, the SP equation (11) must be solved numerically, and the resulting SP field configuration depends on the chosen end time t f , i.e. ϕ SP ≡ ϕ SP (t|t f ). However, as further discussed below, for local observables of the form O = Û †ÔÛ the functional equation (11) can be reduced to a differential equation. The effective action for an observable of this form is given by Eq. (14) features forwards and backwards fields ϕ x = {ϕ a x,j } with x ∈ {f, b}, introduced to decouple the two time-evolution operators in O, and the respective noise actions S 0 [ϕ x ] and disentangling variables ξ a x,j , given by (3) and (4). Eq. (14) is general to spin-1/2 systems in any dimension. It is apparent that only the last term of Eq. (14) is observable-dependent. Let us begin by considering the caseF O = 1, corresponding to the normalization function F ½ . Extremization of (14) with respect to ϕ a f,i (t ′ ) leads to the SP equation where we used δξ a j /δϕ b k ∝ δ jk ; explicit expressions for the functional derivatives are given in Appendix A. Symmetry between the forwards and backwards fields implies ϕ a f,i | SP = ϕ a b,i | SP = ϕ a SP,i and similarly ξ a f,i | SP = ξ a b,i | SP ≡ ξ a SP,i . For ground state expectation values, SP equations analogous to (15) can be reduced to algebraic ones by considering the infinite imaginary-time limit [31]; this is however not possible in the present context of real-time evolution. However, it can be shown that the solution of Eq. (15) satisfies ∂ t f ϕ a SP,i (t|t f ) = 0; see Appendix B. As a consequence, the end time t f in Eq. (15) can be chosen freely so as to simplify the computation of the SP configuration. It is convenient to set t f = t ′ ; the SP equation (15) then readily yields For this model, if any of the dimensions N k of the system is a multiple of 4, the corresponding interaction matrix needs to be regularized by including a diagonal term, which does not affect the resulting dynamics; see Ref. [24].
where the vectors v a k ≡ with a ∈ (+, z, −), feature the normalized expectation values of the spin operatorsŜ a k under the dynamics induced by the SP field. Eq. (16) is thus equivalent to a mean field condition, ϕ a SP,j (t) = bk J ab jk Ŝ b k (t) | SP , whereby the effective field acting on each spin is produced by the magnetization of its neighbors. The above steps should in principle be repeated for each different observable, adding the corresponding term − logF O to the action. However, for translationally invariant observables it can be shown that theF O -dependent term in the SP equation becomes negligible in the thermodynamic limit; see Appendix C. This makes it possible to use the SP configuration given by (16) to perform importance sampling for local observables given a sufficiently large system. These findings generalize the results of Ref. [31], where it was shown that in the limit of infinite imaginary time the dominant contribution to ground-state expectation values corresponds to the mean-field ground state. The physical interpretation is analogous: the optimal approximation to the full dynamics of a quantum system within the manifold of single-spin trajectories, in the spirit of the time-dependent variational principle (TDVP) [38], is given by mean field. In contrast to TDVP, however, here we do not restrict ourselves to the optimal trajectory, but perform a sum over trajectories: this restores entanglement, and the resulting time evolution is formally exact. The set of coupled equations (16) can be solved numerically together with (4); the solution matches the direct recursive solution of (15), but is much more efficient.
Having obtained the SP configuration from Eq. (16), we can perform importance sampling according to Eq. (12) to compute local observables. To illustrate the difference in performance between direct and importance sampling, in Fig. 2 we consider results obtained using the same numerical solution scheme, discretization time step, and number of simulations. We consider the longitudinal magnetization of a 3 × 3 quantum Ising model initialized in the symmetry-broken ferromagnetic ground state for Γ = 0, h < 0, |⇓ ≡ ⊗ i |↓ i , and evolved with Γ = h = 2J, comparing our results to exact diagonalization (ED) performed using the QuSpin package [39]. Generating each data set required ≈ 12 minutes on a laptop, using 2 Intel Core i5 processors with a clock speed of 2.9 GHz. The results obtained using the importance sampling approach are in much better agreement with ED compared to direct sampling. The difference between the two approaches is also reflected in the behavior of fluctuations around mean values. Due to the strong fluctuations in the direct sampling results, dividing the data set into small batches leads to an underestimation of fluctuations. Instead, in Fig. 2 we estimate fluctuations as the standard error σ/ √ N over the full data set, where the standard deviation σ over N independent simulations is obtained from the standard deviations of the numerator and denominator by using uncertainty propagation [31]. The bars clearly show that the importance sampling scheme results in a significant mitigation of fluctuations. To investigate this quantitatively, in the inset of Fig. 2(b) we show the variance σ 2 F of the stochastic function F ½ yielding the normalization, whose behavior is representative of all local observables due to Eq. (7) [31]. Beyond a transient regime, the time evolution of this quantity is well-approximated by exponential growth, σ 2 F (t) ∼ α exp(βN t) with α ≈ 10 −3 , β ≈ 1. In contrast, for direct sampling one has a comparable β ≈ 1 but a much larger prefactor α ≈ 10. Thus, the importance sampling approach does not eliminate the exponential growth of fluctuations, but it can suppress by several orders of magnitude the associated prefactor. This allows importance sampling to access larger systems and later times than it was previously possible using the direct approach [23][24][25]. As an example, in Fig. 3 we consider 5×5 and 13 × 13 systems initialized in the |⇓ state and timeevolved using the Hamiltonian (13) with Γ = h = J/4. The suppression of fluctuations is increasingly effective as the transverse field Γ is reduced and the classical Γ = 0 limit is approached, as previously reported for imaginarytime evolution [31]. This can again be understood in light of the physical interpretation of the importance sampling approach, whereby the sampling accounts for the presence of entanglement on top of the optimal mean-field trajectory. As a consequence, although the above derivation did not assume a particular regime, the approach can be expected to be most effective in regimes where mean field theory would provide a reasonably good approximation to the true quantum dynamics.

V. LOSCHMIDT AMPLITUDE
The importance sampling approach is not restricted to local observables and can be applied to compute global quantities. As an illustration, we consider the Loschmidt amplitude A(t) ≡ ψ(0)|ψ(t) , where |A(t)| 2 gives the probability for the system to return to its initial state following unitary evolution for a time t. Since A(t) is exponentially suppressed as N → ∞, one typically considers the rate function λ(t) ≡ − log |A(t)| 2 /N , also known as the fidelity density, which has a well-defined thermodynamic limit [32]. This quantity has recently received great theoretical [8] and experimental [40] interest in the context of dynamical quantum phase transitions (DQPTs), a proposed generalization of equilibrium QPTs whereby λ(t) becomes non-analytic as a function of time [8,32]. Notably, DQPTs can only occur in the thermodynamic limit, but signatures of their presence can be observed in large but finite systems [23,24,32,40]. Furthermore, these phenomena typically occur at early times, making them promising candidates for currently available experimental platforms.
To illustrate how the importance sampling method can be applied to reveal DQPTs, we consider the 2D Ising model (13) with h = 0. This model has a QPT at Γ ≈ 1.523 J [41]. We initialize the system in the symmetric FM ground state |ψ FM ≡ (|⇑ + |⇓ )/ √ 2, where |⇑ ≡ ⊗ i |↑ i , and consider a quench deep into the paramagnetic phase; the return probability is then given by whereÛ =Û (t) and we used | ⇑|Û |⇑ | = | ⇓|Û |⇓ | and | ⇑|Û |⇓ | = | ⇓|Û |⇑ |. DQPT in this quantity can be understood as arising from the crossing of the contributions coming from |A ud | 2 and |A dd | 2 . In the stochastic approach, the amplitudes are obtained as The SP equation for each amplitude can be obtained as in (11) and solved numerically to find the SP field. In contrast to local observables, here the SP field is not equal to the mean field and the whole configuration depends on the chosen end time t f , ϕ SP ≡ ϕ SP (t|t f ). One can then perform importance sampling separately for the two contributions in (18), using the SP field obtained numerically for each. Fig. 4 shows the results obtained in this way for systems of increasing size. In the main panel we compare our results to ED for a 5 × 5 The main panel shows that the results obtained from the importance sampling (IS) approach (dots) are in good agreement with ED (full line) for a 5 × 5 system. As the system size is increased, the peak in the rate function sharpens, eventually becoming non-analytic in the thermodynamic limit; this is demonstrated in the inset, where we show results for 7 × 7 and 15 × 15 systems for which ED cannot be performed. The numerical results were obtained from 5 × 10 4 , 10 5 and 10 6 independent simulations, in order of increasing system size.
The error bars are nearly invisible on the scale of the plot.
system time-evolved with Γ = 8J, finding good agreement. Here we use the stochastic Heun scheme applied in Ref. [25], setting ∆t = 10 −4 . Notably, the number of simulations required by the importance sampling scheme to accurately reproduce the ED result is more than two orders of magnitudes smaller in comparison to direct sampling using the measure (5) [25]. The inset shows that the peak becomes increasingly sharp as the system size is increased to 7 × 7 and 15 × 15. It is known that non-analytic points can also occur in the time evolution of the term |A dd | 2 , which corresponds to the return probability for a quench from the symmetry-broken initial state |⇓ , and the relative rate function λ dd = − lim N →∞ log |A dd | 2 /N [32]. In contrast to the previously considered case of a quench from |ψ FM , such DQPTs cannot be straightforwardly understood as arising from a sum of competing contributions. However, insights about the origin of DQPTs in the amplitude |A dd | can be obtained from the solution of the SP equation itself. We illustrate this for the 1D Ising chain, where the exact location of the non-analytic point can be computed analytically or determined to arbitrary numerical precision using infinite time-evolving block decimation (iTEBD) [12].
In Fig. 5 we consider the solution ϕ SP (t|t f ) of Eq. (11) as a function of the stopping time t f for a system initialized as |⇓ and evolved with Γ = 8J. For this quench, a DQPT occurs at t * ≈ 0.397. It can be seen that the entire SP configuration evolves as a function of t f , with an abrupt change occurring at t f ≈ t * whereby the real and imaginary parts of the SP field change sign. As shown in the inset, the number of iterations N SP required to recursively solve Eq. (11) sharply increases in a region near t f ≈ t * . In particular, we find that the recursive solution does not converge within 10 3 iterations for 0.404 < t f < 0.409. Such behavior points to an instability of the recursive solution in this region, which could be due to the coexistence of different minima yielding a comparable action; the leading and subleading contributions would then switch at a critical time t * SP ≈ 0.406. These observations suggest that DQPTs in |A dd | 2 are associated with an abrupt change of the product-state configuration carrying the largest contribution to the quantum dynamics, consistently with recent findings [42]. This immediately generalizes the phenomenology of ground-state QPTs in the disentanglement formalism, where QPTs are associated with an abrupt change of the dominant trajectory as a function of a control parameter [31]. Similar switching behavior near DQPTs was previously observed by considering generalized expectation values in fermionic systems [43]. In the present context, this phenomenology would not be reproduced if the SP field corresponded to the mean field, since it is closely tied to the end-time dependence of Eq. (11). The full quantum dynamics, obtained from sampling, includes the additional effects of entanglement on top of the optimal SP productstate dynamics. This potentially explains the small difference between t * SP and the DQPT critical time t * : the product-state approximation provided by the SP field can only determine the DQPT location approximately, as a degree of entanglement is always needed to fully capture these phenomena [42].

VI. CONCLUSIONS
In this manuscript, we have introduced an importance sampling scheme for the real-time dynamics of many-body quantum spin systems. The importance sampling method is based on the disentanglement approach, whereby unitary quantum dynamics is exactly mapped to an ensemble of classical stochastic processes. Quantum expectation values are then obtained as averages over stochastic trajectories, which can be generated numerically. We have shown that the dominant contribution to a given observable is given by a saddle point trajectory, which can be obtained by extremizing an appropriate effective action. By preferentially sampling trajectories close to the saddle point trajectory, it is possible to significantly improve the performance of the method compared to direct sampling, as we have demonstrated for both local observables and return probabilities. This improvement in performance is due to a suppression in the strength of fluctuations in the stochastic quantities, which determine the numerical efficiency of the method. However, fluctuations are still found to grow exponentially with time and the system size; while the accessible parameter regions can be extended by using importance sampling, large systems and late times remain challenging to capture. Further progress will thus be needed in order for the method to access interesting regimes of higher-dimensional quantum dynamics; several directions for development can be envisaged, including cluster approaches and the development of approximate schemes that truncate fluctuations. Since the importance sampling scheme completely eliminates fluctuations in the limit of a Hamiltonian made up of commuting terms, it would also be interesting to explore its application to higher-spin systems, where the performance of the approach might benefit from the proximity to a classical limit. Furthermore, the broader framework introduced in this work might prove useful beyond its application to importance sampling. In the context of imaginary-time evolution, it has recently been shown that corrections to the mean-field estimate for ground state expectation values can be analytically obtained order-by-order by viewing the disentanglement approach as a field theory [31]. A similar development for real-time evolution would make it possible to systematically include the effect of entanglement on top of the leading-order mean-field dynamics.
Acknowledgments During the preparation of this manuscript we became aware of the work [44], in which an importance sampling scheme is developed by considering the Hermiticity of the effective Hamiltonian governing the stochastic evolution; this also leads to an improvement in the accessible timescale for a given number of simulations.
solution for ϕ a SP,k obtained below can be self-consistently checked to satisfy these conditions. The initial equality is then verified if ∂ t f ϕ a SP,i (t ′ ) = 0 ∀ t ′ = t f , so that L also vanishes; a t f -independent solution is thus consistent with Eq. (15). We therefore arbitrarily choose t f in (15); the simplest choice is given by t f = t ′ , which yields: bk [J −1 ] ab jk ϕ b SP,k = − 1 2 δ az − 2δ a− ξ + SP,j + ξ + * SP,j 1 + ξ + SP,j ξ + * SP,j δ a+ + δ az ξ + SP,j − δ a− (ξ + SP,j ) 2 , where the explicit time-dependence has been suppressed, since all variables are evaluated at the same time t, and we used Eqs (A2) and (A4). Considering the different cases a = (+, z, −) readily reproduces Eq. (16), which is consistent with the conditions ϕ + SP,k = (ϕ − SP,k ) * , ϕ z SP,k ∈ Ê. This can be checked to be a solution of (15) by substitution. For the quantum Ising model, a direct numerical solution of Eq. (15) matches the solution (16).