Fluctuations and Non-Hermiticity in the Stochastic Approach to Quantum Spins

We investigate the non-equilibrium dynamics of isolated quantum spin systems via an exact mapping to classical stochastic differential equations. We show that one can address significantly larger system sizes than recently obtained, including two-dimensional systems with up to 49 spins. We demonstrate that the results for physical observables are in excellent agreement with exact results and alternative numerical techniques where available. We further develop a hybrid stochastic approach involving matrix product states. In the presence of finite numerical sampling, we show that the non-Hermitian character of the stochastic representation leads to the growth of the norm of the time-evolving quantum state and to departures for physical observables at late times. We demonstrate approaches that correct for this and discuss the prospects for further development.

We investigate the non-equilibrium dynamics of isolated quantum spin systems via an exact mapping to classical stochastic differential equations. We show that one can address significantly larger system sizes than recently obtained, including two-dimensional systems with up to 49 spins. We demonstrate that the results for physical observables are in excellent agreement with exact results and alternative numerical techniques where available. We further develop a hybrid stochastic approach involving matrix product states. In the presence of finite numerical sampling, we show that the non-Hermitian character of the stochastic representation leads to the growth of the norm of the time-evolving quantum state and to departures for physical observables at late times. We demonstrate approaches that correct for this and discuss the prospects for further development.
Experimental progress on cold atomic gases and trapped ions has led to pristine realizations of isolated quantum spin systems in and out of equilibrium [1][2][3][4]. This has stimulated intense theoretical activity to expose the unitary dynamics of paradigmatic spin Hamiltonians, with a view towards extracting universal results [5][6][7]. Much of the attention has focused upon one-dimensional spin models due to the availability of analytical [8][9][10][11][12][13] and numerical [14,15] techniques. This has yielded fundamental insights into the nature of thermalization [16][17][18] and to the development of new techniques [9,19,20]. The prediction of dynamical quantum phase transitions (DQPTs) occuring in the time-domain [21] has been confirmed by experiment on Ising Hamiltonians realized with trapped ions [22]. This opens the door to time-resolved dynamics in tunable quantum spin systems, allowing direct comparison between theory and experiment.
A recent theoretical approach to non-equilibrium quantum spin systems permits an exact mapping to classical stochastic differential equations (SDEs) [23][24][25][26][27]. The time-evolution of quantum observables is encoded by classical averages over independent realizations of the stochastic process. The method is therefore inherently parallelizable and can be implemented by numerically sampling the SDEs [26,27]. The stochastic approach is rather general, since it applies to both integrable and non-integrable Hamiltonians, including those in higher dimensions. The stochastic framework also reveals deep connections between classical and quantum dynamics, as recently illustrated in the context of DQPTs [26,27].
In this work, we show that the stochastic approach to quantum spin systems can address significantly larger system sizes than previously possible [26,27]. This is obtained through the use of a Heun integration scheme and the elimination of divergent stochastic trajectories. We show that the results obtained for the onedimensional (1D) quantum Ising model are in very good agreement with those obtained from free fermions [21] and via matrix product operator (MPO) methods [28].
We also provide results for the two-dimensional (2D) quantum Ising model with up to 49 spins. We relate the growth of stochastic fluctuations at late times to the non-Hermiticity of the effective stochastic Hamiltonian. Due to the impact of finite numerical sampling, this leads to an increase in the norm of the time-evolving quantum state and to departures for observables at late times. This can be partially corrected by rescaling by the norm. We show that a hybrid numerical scheme combining SDEs with matrix product states can reduce the number of noise variables, thereby extending the simulation time. We conclude and provide directions for research.
Stochastic Approach.-Here, we briefly recall the principal steps in the stochastic approach to quantum spin systems [23][24][25] following the notations in [25,26]. We begin with a generic Heisenberg Hamiltonian where i and j indicate lattice sites, J ab ij are exchange interactions, and h a i are magnetic fields. The spin operators,Ŝ a i , obey the canonical commutation relations, [Ŝ a i ,Ŝ b j ] = i abc δ ijŜ c i , with = 1 and a, b, c ∈ {x, y, z}. The corresponding time-evolution operator between times t i and t f is of the formÛ (t f , t i ) = The key idea is that the exchange interactions can be decoupled using a Hubbard-Stratonovich transformation, which introduces fluctuating stochastic fields ϕ [23][24][25]. The time-evolution operator can be expressed aŝ with the Gaussian noise measure [23,25]. This formulation describes N decoupled spins evolving under effective "magnetic" fields.
It is inherently arXiv:1912.10551v1 [cond-mat.str-el] 22 Dec 2019 non-Hermitian due to the factor of 1/ √ i, and the presence of complex fields ϕ [26,27]. Focusing upon the case where a = b in (1), we can diagonalize the N × N matrix (J aa ) −1 , for a given a. Explicitly, we may write (J aa ) −1 = V a D a (V a ) −1 where D a is a diagonal matrix and V a is an eigenvector matrix, where a labels the matrices and not their components. It is also convenient to introduce the white noises [25] It follows from (2) that the resulting dynamics are described by a local stochastic Hamiltonian,Ĥ s .. denote the average with respect to the Gaussian measure, the time-evolution operator can be expressed as an average over stochastic evolution operators. Explicitly, ; here we writeÛ (t) ≡Û (t, 0) for brevity. SinceÛ s j (t) has an exponent that is linear in theŜ a j with complex coefficients,Û s j (t) is an element of SL(2, C). As such it can be re-expressed as a product of exponentials without time-ordering, via a disentanglement transformation [25,29,30]. Using the Gauss parametrization [25,31] where ξ ±,z j ∈ C are referred to as disentangling variables [25,26]. The disentanglement is achieved independently on each site by solving the Schrödinger equation This yields [25] iξ . The equations (5) are SDEs for the ξ-variables due to the stochastic fields Φ a j (t) [25]. Quantum observables, Ô (t) are calculated as classical averages over functions f (ξ) of the ξ-variables, via Ô = f (ξ) φ . In order to evolve the ξ-variables forward in time, we solve the SDEs (5) using a stochastic Heun predictor-corrector method in the Stratonovich formalism [32,33]. We find that this is capable of maintaining accuracy with larger time-steps than the Euler-Maruyama scheme used previously [26,27], thereby reducing the computational cost.
Parametrization.-To gain some intuition into the dynamics of the SDEs (5), it is instructive to consider the parametrization (15) in more detail. The stochastic evolution operatorÛ s j (t) has a particularly simple form when FIG. 1. Projection of the Bloch sphere for an un-normalized quantum spin onto the complex plane, parametrized by ξ + . The point P on the unit sphere is projected onto the point P N via the North pole, N. The point ξ + = 0 corresponds to spindown |↓ , whilst |ξ + | → ∞ corresponds to spin-up |↑ . Potential divergences associated with |ξ + | → ∞ can be avoided via a two-patch parametrization: the upper (lower) hemisphere is parametrized by projection from the South (North) pole. A mapping between the two patches is performed at the equator.
acting on spin-down states [26,27], due to the explicit form of the parametrization (15): where the variable ξ − j (t) drops out. Any stochastic state |ψ s (t) = j |ψ s j (t) can be parametrized in this way by introducing a preparation stage in which the initial state, |ψ s j (0) , is obtained as a rotation from a spin-down state |↓ ; see Supplementary Material.
For a normalized spin state, spin-down |↓ corresponds to ξ + j (t) = 0 and spin-up |↑ corresponds to |ξ + j (t)| → ∞; this is a stereographic projection of the Bloch sphere, via the North pole, as shown in Fig. 1. The complex parameter ξ z j (t) determines the amplitude and phase of the spin state. Divergences in the SDEs (5) [26,27] corresponding to |ξ + j (t)| → ∞ can be avoided by a two-patch parametrization of the Bloch sphere by projecting from the South pole for states in the upper hemisphere. This can be implemented by the change of variables whenever the spins cross the equator. The corresponding SDEs for the new coordinates are given in the Supplementary Material. This approach for avoiding divergences in SDEs has also been used in [34]. We will use this two-patch approach throughout the manuscript. For simplicity, we focus on the nearest neighbor spin-1/2 ferromagnetic quantum Ising model where we impose periodic boundary conditions. Throughout this paper we set J = 1 in the simulations.  [21] in the thermodynamic limit (dashed line). It is readily seen that the rounding of the Loschmidt peak is a finite-size effect. In both figures we use a time-step of dt = 10 −2 , except in the vicinity of the peaks, where dt = 10 −3 is used. The results are obtained by averaging over 10 7 stochastic samples. (b) 2D case for a 5 × 5 lattice using 1.5 × 10 7 samples. The results are in agreement with QuSpin [35] (solid line). The inset shows results for a 7 × 7 lattice, which cannot be obtained using QuSpin. Convergence is checked by changing the number of samples.
Loschmidt Amplitude.-As discussed in [26,27], one of the simplest quantities to examine in the stochastic approach is the Loschmidt amplitude, The corresponding rate function λ L (t) ≡ − 1 N ln |A(t)| 2 plays a similar role to the equilibrium free energy density: as N → ∞ it exhibits nonanalytic peaks at DQPTs [21]. We first consider the onedimensional case. In order to compare to results obtained in the thermodynamic limit [21] it is convenient to evolve from the ground state |ψ(0) = 1 √ 2 (|⇓ + |⇑ ) at Γ = 0.
Here |⇑ and |⇓ correspond to the states with all the spins pointing up and down respectively. Time-evolving these separately using the SDEs (5) one obtains: where ψ(0)|Û s (t)|⇓ = 1 , and ψ(0)|Û s (t)|⇑ is obtained by ξ a →ξ a . The SDEs are solved with the initial conditions ξ a (0) = 0 andξ a (0) = 0 respectively. The results for λ L (t) corresponding to timeevolution with Γ = 8J are shown in Fig. 2(a), for a 1D system with N = 50 spins. The results go beyond what is achievable using Exact Diagonalization (ED) and are in excellent agreement with those obtained via the MPO W I method [28], implemented using ITensor [36]. Deviations are observed for t 1/J as the stochastic fluctuations become harder to sample. The inset shows the first Loschmidt peak for N = 75, which again demonstrates excellent agreement. For comparison we display exact results obtained in the thermodynamic limit, N → ∞ [21]. Although the finite-size effects are stronger in the vicinity of the peak, the SDE and MPO results remain coincident for all of the system sizes considered. Results for the same quench on a 5×5 lattice in 2D are shown in Fig.  2(b) (dots), and are verified against those obtained using QuSpin's time-evolution solver [35] (solid line). The inset shows the first Loschmidt peak for a 7 × 7 lattice, which goes beyond what we can readily verify using other techniques. We check for convergence near the peak by doubling the number of samples, N , and noting that the results change by less than 0.5% in this region.
Growth of Fluctuations.-To quantify the role of stochastic fluctuations it is instructive to consider the spectrum of an effective Hamiltonian,Ĥ eff (t), defined bŷ (10) in analogy to Floquet systems [37]. SinceÛ s (t) is nonunitary, the eigenvalues ofĤ eff , ε = ε R + iε I , are generically complex. The spectrum ofĤ eff (t) can be calculated directly from (10) by time-evolving the SDEs to the time of interest. This can also be obtained by noting that U s (t) can be calculated directly as a product of random matrices,Û s (t) =Û s (t, t − δ)Û s (t − δ, t − 2δ)...Û s (δ, 0), by time-slicing into small intervals of size δ, without the disentangling transformation. In Fig. 3 we show the time-evolution of the eigenvalue distribution ofĤ eff , for 100 stochastic samples with Γ = 8J and N = 10 in 1D. It can be seen that the distribution of ε R is uniform at late times, whilst that of ε I is well approximated by a normal distribution. In the Supplementary Material we show that the variance of the distribution of ε I , denoted by σ 2 (t) = ε 2 I − ε I 2 , exhibits damped oscillations as a function of time, with extrema that occur in proximity to those in the time-dependent magnetization. In general, the presence of the positive imaginary eigenvalues results in the growth of the norm of individual stochastic states over time. Due to the effect of finite numerical sampling, this leads to the growth of the norm of the overall quantum state, and to departures for physical observables. As we will see below, this can be partially compensated by rescaling by the norm.
Magnetization Dynamics.-As discussed in [26,27], time-dependent physical observables can be obtained by using two Hubbard-Stratonovich transformations to decouple the forwards and backwards evolution operators: where φ andφ are independent noise variables. In this representation, the local magnetization is given by [26] Ŝ z j = − where χ ≡ i ξ z i ,χ ≡ iξ z i , and we implicitly take the real part; in general, observables have imaginary parts which vanish in the limit of infinite sampling [19,34,38]. In Fig. 4(a) we show results for the time-dependent magnetization M(t) = 1 N j Ŝ z j (t) , following a quantum quench from the fully-polarized initial state |⇓ to Γ = 8J for a 1D system with N = 25 spins. The results are in good agreement with MPO calculations until times t 1/J when stochastic fluctuations become large. In Fig. 4(b) we show results for the norm of the timeevolving state as computed from the SDEs: where again, we take the real part. It is readily seen that the norm departs from unity once the stochastic fluctuations become significant. In Fig. 4(c) we show results for the re-scaled magnetization M res (t) = M(t)/|ψ(t)| 2 which provides much better agreement with the MPO results until later times. Fluctuations in M res (t) still occur however, especially when the norm of the state is close to zero. This clearly highlights the importance of normalization and stochastic sampling in the computation of observables. Matrix Product States.-Another approach to reducing fluctuations in observables is to decompose the timeevolving state into a matrix product state (MPS): where A σi are matrices with physical spin indices, σ i , and auxiliary indices d i ∈ 1, . . . , D i , where D i are the bond dimensions; see [39] for an introduction. This reduces the number of noise variables required since only a single Hubbard-Stratonovich transformation is needed for time-evolution; see Supplementary Material. One may also calculate the norm of the state using MPS techniques, thereby eliminating fluctuations from the stochastic sampling of |ψ(t)| 2 . In Fig 4(d) we show the results of the hybrid stochastic-MPS approach for a quantum quench in a 1D system with N = 50 spins. The results are in excellent agreement with the MPO approach, in spite of doubling the system size and reducing the number of stochastic samples. A notable disadvantage of this hybrid approach is that one must store the MPS state in memory at the expense of the stochastic parallelization. Nonetheless, the marriage of these approaches may be useful for future developments.
Conclusions.-In this work we have demonstrated that the stochastic approach to non-equilibrium quantum spin systems can address significantly larger systems than recently obtained, in both one and two dimensions. We have shown that the non-Hermitian character of the representation leads to a growth of the norm of |ψ(t) , due to the effect of finite numerical sampling. However, this can be compensated for by rescaling by the norm. We have shown that the approach can be combined with a decomposition in terms of matrix product states, for the calculation of time-dependent observables. There are many directions for research, including extensions to larger system sizes and later times, particularly in higher dimensions where few techniques are available.
Acknowledgements.-We acknowledge helpful conversations with S. De Nicola, B. Doyon, D. O'Dell and S. Wüster. We also thank J. Morley for assistance with the implementation of MPS algorithms. SEB is supported by the EPSRC CDT in Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES) via grant number EP/L015854/1. We are grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC (EP/P020194/1). The MPO calculations were performed using the ITensor Library [36]. MJB acknowledges the support of the ICTS (Bengaluru) during the program on Non-Hermitian Physics PHHQP XVIII. AGG acknowledges EPSRC grant EP/P013449/1.

SUPPLEMENTARY MATERIAL I. PARAMETRIZATION
As discussed in the main text, we may eliminate divergent trajectories from the SDEs (5) by a suitable parametrization of the stochastic time-evolution operator,Û s j (t). Adopting the Gauss parametrization of SL(2,C) [25] U s For spin-1/2 systems this can be represented in matrix form aŝ where ξ ±,z j ∈ C, and the initial conditions ξ ±,z j (0) = 0 ensure thatÛ s j (0) = 1 is the identity operator. In general, this corresponds to a representation of the group SL(2,C) with three complex parameters. The action ofÛ s j (t) on a generic initial state |ψ j (0) = a|↓ + b|↑ , with a, b ∈ C yields Although (17) is formally exact, the parametrization contains some redundancy: an arbitrary un-normalized spin state can be represented by four parameters, including the overall phase. To see that a reduction is possible it is instructive to consider an initial spin-down state |↓ corresponding to a = 1 and b = 0. This yields where the complex parameter ξ − j (t) has dropped out. In this case, the divergence in the SDE (5a) corresponding to |ξ + j | → ∞ is associated with an inability to parametrize the spin-up state |↑ , using a projective representation of the Bloch sphere. As discussed in the main text and in Section II below, this can be avoided by using a two-patch parametrization. Although these considerations apply only for an initially spin-down state, more general initial states can always be prepared by rotation from this state. Explicitly, we may introduce a state-preparation protocol starting at t = −δ, with δ > 0, and evolving deterministically until t = 0. In this approach, the time-evolution operator takes the formÛ s and the coefficients α a j specify the initial conditions according to |ψ s j (0) =Û s j (0, −δ)|↓ . In practice, this is equivalent to setting non-trivial initial conditions for the ξ-variables and evolving under the SDEs (5). For example, the initial state |ψ j (0) = 1 √ 2 |↓ + |↑ corresponds to ξ + j (0) = 1 and ξ z j (0) = ln 2, as follows directly from (18). In this approach the trivial initial conditions correspond to t = −δ, so that ξ + j (−δ) = ξ − j (−δ) = ξ z j (−δ) = 0 andÛ s j (−δ, −δ) = 1. In general, the time-evolution of an arbitrary product state is given by where the initial conditions ξ a j (0) specify the initial spin-orientation at each site. A generic superposition can be obtained by summing over (20) with the appropriate initial conditions. Stochastic expressions for physical observables are readily obtained from the projective representation (18). For example, the Loschmidt amplitude to remain in the spin-down state is given by φ , in agreement with (9) and [26]. In a similar way, the quantum expectation value of the spin operatorŜ j is given by where n j (t) = 1 2 corresponds to the position of a spin on the Bloch sphere. The factor of is the norm of the state |ψ s i (t) . In writing (21), (22) and (23), it is implicit that ξ a * j is an independent variable from ξ a j , which we denote byξ a * j = ξ a * j in the main text. The result (21) is readily generalized to multipoint correlation functions. For example, The expressions for entangled states can be obtained by averaging over the initial conditions for ξ(0) andξ(0).

II. ELIMINATING DIVERGENT TRAJECTORIES
As discussed above and in the main text, the SDE (5a) exhibits divergences corresponding to |ξ + j | → ∞. These can be avoided by two-patch parametrization of the Bloch sphere. To this end, it is convenient to define new variablesξ a via the generalization of equation (18), where the roles of |↓ and |↑ are interchanged: Equating the coefficients of (18) and (25), one obtains the identifications This coordinate system is related to the original Gauss parametrization (15) by swapping the pole of projection from the North to the South pole. Performing this change of variables, the SDEs (5a) and (5b) become A convenient place to perform the change of variables (26) is when the spins cross the equator, since the magnitude of |ξ + j (t)| = |ξ + j (t)| = 1 at this point, and the numerical error associated with simulating the nonlinear term in (5a) is minimized. This approach is also used in [34]. In practice, the initial state will determine which parametrization is initialized on each site. Denoting one may track the dynamics of this single variable over all times. In Fig. 5 we plot the time-evolution of the singlepatch variable |ξ + 1 (t)| following a quantum quench in the 1D quantum Ising model. For the chosen parameters and the specific noise realization, it can be seen that this quantity diverges at time t div = 6.25. As a result,ξ + 1 (t) overflows due to the ξ + 2 1 (t) term in (5a), and numerical integration fails for this trajectory. In contrast, the two-patch variable |ζ + 1 (t)| remains finite, and we can evolve beyond t div . This enables us to retain all the stochastic trajectories when computing the time-dependent magnetization, in contrast to previous work [26,27].
To analyze the spectrum ofÛ s (t, 0), or equivalentlyĤ eff = i t lnÛ s (t, 0) as illustrated in Fig. 3, the parameter ξ − j is also required; this dropped out in (18). Under the transformation (ξ + j , ξ z j ) → (ξ + j ,ξ z j ) the SDE (5c) becomes In order to ensure that (29) is well-behaved asξ + j → 0 it is convenient to make the change of variables The resulting SDE forξ − j is given by which mirrors (5c) up to a sign change and Φ − j → Φ + j . The maps (26) and (30) can now be conducted simultaneously to avoid the divergence associated with ξ + j → ∞, allowingÛ s (t, 0) to be calculated at much later times. Eventually this strategy will break down if (5c) or (31) cannot be integrated due to e ξ z j → ∞ or eξ z j → ∞ respectively. However, ξ − j andξ − j are not required for the time-dependent magnetization, this is not a limitation. The spectrum ofĤ eff (10) is calculated directly from trajectories by mapping between the ξ a j andξ a j variables in accordance with the prescription (28). As discussed in the main text, damped oscillations of the variance of the imaginary eigenvalues ofĤ eff occur as a function of time; see Fig. 6.

III. HYBRID TECHNIQUE WITH MATRIX PRODUCT STATES
As discussed in the main text, the time-evolving quantum state evaluated within the stochastic approach can be represented as a matrix product state. This halves the number of noise variables required since only a single Hubbard-Stratonovich transformation is needed to evaluate |ψ(t) =Û (t)|ψ(0) . However, this comes at the cost of storing the state in memory, so the method is no longer fully parallelizable. Here we demonstrate how this representation is obtained. The quantum state is first written as the sample average of the stochastic state: where r = 1, ..., N is the sample index. A matrix product state (MPS) is given by forms a product state that can be identified as one term in the sum (32). We can identify a trivial, inefficient MPS representation of (32) by taking a diagonal form for the MPS tensors. For example where the factor of 1/N in (32) can be absorbed into one of the A σi matrices. The state |ψ can be compressed to a lower bond-order, full rank MPS by using a sequence of singular value decompositions (SVDs). In practice, we may perform this procedure in batches. The MPS tensors in each batch can be further combined and compressed by first collecting them into a block diagonal MPS tensor given by where the batch index runs from 1 to k. A subsequent sequence of SVDs would again lead to an MPS with reduced bond order and full rank. For example, the results in Fig. 4 were obtained by dividing the 50, 000 trajectories into 50 batches of size 1000. After the initial compression, the batches were combined in pairs according to (35) and compressed again. This was repeated two more times in groups of five batches. In practice, if the bond dimension after a compression exceeds the nominal value of 20 we truncate it back to this value. Since the mapping to MPS must be carried out at each time where observables are calculated, it is more efficient to only carry it out at the times of interest. For all the MPO W I [28] simulations carried out in the main text we use a maximum bond dimension of 50, a minimum singular value cut-off of 10 −14 , and a time-step of dt = 0.001.

IV. SCALING
In order to quantify the scaling properties of the real-time stochastic approach, we consider the time-scale over which the simulations are accurate as a function of the system size and the number of samples. Since the stochastic approach only produces perfectly normalized quantum states in the limit N → ∞, deviations of the norm from unity can provide an estimate of the simulation's convergence and the time-scale over which the method can be trusted. As illustrated in Fig. 3, rescaling by the norm can lead to good approximations for physical observables, even if the norm deviates significantly from unity. In view of this, we define the breakdown time, t b , as the earliest time for which a 10% error is observed in the norm. In Fig. 7(a) we show t b as a function of the inverse system size N −1 , for quenches in the 1D quantum Ising model from the fully polarized initial state |⇓ to Γ = 8J, for a fixed number of samples N = 10 6 . The data are well approximated by the linear relation Jt b = 16.94N −1 + 0.12. (36) That is to say, for a fixed number of samples the breakdown time scales with the inverse of the system size. This is consistent with the results of [27]. In Fig. 7(b) we also show t b as a function of the number of samples, for the same quench, but for a fixed system size, with N = 7 spins. The data are compatible with the relation Jt b = 0.22 ln N −0.6, i.e. the number of samples required to reach a given time t b therefore scales exponentially, N ∝ e αJt b , where α ≈ 4.5, in this example. This is consistent with the scaling of fluctuations analyzed in [27] using different diagnostics.