Dynamical quantum correlations of Ising models on an arbitrary lattice and their resilience to decoherence

Ising models, and the physical systems described by them, play a central role in generating entangled states for use in quantum metrology and quantum information. In particular, ultracold atomic gases, trapped ion systems, and Rydberg atoms realize long-ranged Ising models, which even in the absence of a transverse field can give rise to highly non-classical dynamics and long-range quantum correlations. In the first part of this paper, we present a detailed theoretical framework for studying the dynamics of such systems driven (at time t=0) into arbitrary unentangled non-equilibrium states, thus greatly extending and unifying the work of Ref. [1]. Specifically, we derive exact expressions for closed-time-path ordered correlation functions, and use these to study experimentally relevant observables, e.g. Bloch vector and spin-squeezing dynamics. In the second part, these correlation functions are then used to derive closed-form expressions for the dynamics of arbitrary spin-spin correlation functions in the presence of both T_1 (spontaneous spin relaxation/excitation) and T_2 (dephasing) type decoherence processes. Even though the decoherence is local, our solution reveals that the competition between Ising dynamics and T_1 decoherence gives rise to an emergent non-local dephasing effect, thereby drastically amplifying the degradation of quantum correlations. In addition to identifying the mechanism of this deleterious effect, our solution points toward a scheme to eliminate it via measurement-based coherent feedback.

Abstract. Ising models, and the physical systems described by them, play a central role in generating entangled states for use in quantum metrology and quantum information. In particular, ultracold atomic gases, trapped ion systems, and Rydberg atoms realize long-ranged Ising models, which even in the absence of a transverse field can give rise to highly non-classical dynamics and long-range quantum correlations. In the first part of this paper, we present a detailed theoretical framework for studying the dynamics of such systems driven (at time t = 0) into arbitrary unentangled non-equilibrium states, thus greatly extending and unifying the work of Ref. [1]. Specifically, we derive exact expressions for closed-timepath ordered correlation functions, and use these to study experimentally relevant observables, e.g. Bloch vector and spin-squeezing dynamics. In the second part, these correlation functions are then used to derive closed-form expressions for the dynamics of arbitrary spin-spin correlation functions in the presence of both T 1 (spontaneous spin relaxation/excitation) and T 2 (dephasing) type decoherence processes. Even though the decoherence is local, our solution reveals that the competition between Ising dynamics and T 1 decoherence gives rise to an emergent non-local dephasing effect, thereby drastically amplifying the degradation of quantum correlations. In addition to identifying the mechanism of this deleterious effect, our solution points toward a scheme to eliminate it via measurement-based coherent feedback.

Introduction
Interacting spin models provide a remarkably accurate description of a diverse set of physical systems, ranging from quantum magnetic materials [2,3,4] to quantum dots [5], nitrogen vacancy centers [6], superconducting qubit arrays [7], ultracold atomic gases [8], and trapped ions [9,10]. Despite being relatively simple, and often admitting accurate theoretical descriptions, they support a variety of complex equilibrium properties found in real materials, e.g. emergent spatial ordering [2], quantum criticality [11], and nontrivial topological phases [12,13]. While the equilibrium physics of the simplest quantum spin models is, with many notable exceptions, fairly well understood, the study of driven, dissipative, and otherwise non-equilibrium behavior is comparatively full of open questions: Under what circumstances can equilibrium correlations survive coupling to a noisy environment [14,15,16]? To what extent do the concepts of criticality and universality extend to dynamics and non-equilibrium steady-states [17,18,19,20]? Can interesting quantum-correlated states be stabilized by (rather than degraded by) decoherence [21,22,23,24,25,26,27,28,29]? In recent years, it has become increasingly apparent that non-equilibrium dynamics is ideally suited to investigation by quantum CONTENTS 3 simulation [30], making such questions especially timely and important. Moreover, there are many examples where interesting non-equilibrium states of matter are more readily achievable than low temperature equilibrium states in ultracold neutral gases [31], polar molecules [32], and trapped ions [33].
With these motivations in mind, in this manuscript we develop a general formalism for calculating unequal-time correlation functions of arbitrary-range Ising models driven far out of equilibrium at time t = 0, thus establishing a comprehensive toolbox for the description of non-equilibrium dynamics in a simple context. In addition to providing a tractable example of quantum many-body spin dynamics, the Ising model is realized to a good approximation in a variety of experimentally relevant systems. And, despite its simplicity, Ising spin dynamics is known to be useful for the production of entangled states with applications in quantum information and precision metrology [34]. Our results constitute a unified approach to describing experiments aimed at producing such states [35,36,37,38,39,33], and facilitate a quantitative treatment of a variety of unavoidable experimental complications, e.g. long-range (but not infinite-range) interactions and initial-state imperfections.
Ultracold atomic systems are also well suited to the controlled inclusion of dissipation, prompting a number of theoretical proposals to exploit dissipation for the creation of interesting quantum states [40,19,21,22]-remarkably, such ideas are already coming to experimental fruition [24,41]. Verifying that these experimental systems behave in the expected manner in the presence of dissipation, however, is extremely challenging, in large part due to the numerical complexity of simulating dynamics in open quantum systems and the scarcity of exact solutions. The Ising model, especially as implemented in trapped ion experiments [39,42,33], poses a unique opportunity to study the effects of dissipation in a controlled and, as we will show, theoretically tractable setting. In the absence of dissipation, an important issue in the Ising model is whether ground state correlations survive the application of an equilibrium coherent drive that does not commute with the interactions-i.e. a transverse field. In the dissipative Ising model an analogous question can be posed: How does the system respond to being driven incoherently by processes that do not commute with the Ising interaction? Our formalism for the calculation of unequal-time correlation functions allows us to definitively answer this question.
Quite surprisingly, non-equilibrium dynamics in the Ising model remains solvable in the presence of non-commuting dissipation [1], even for completely arbitrary spatial dependence of the Ising couplings (and therefore in any dimension). This manuscript substantially extends the groundwork laid in Ref. [1], where the quantum trajectories technique was used to obtain a closed-form solution of non-equilibrium dynamics for a special class of initial states. The present work not only provides a more direct, unified, and comprehensive exposition of the relevant theory, but also generalizes those results to include a broader class of initial conditions and observables, and applies the solutions to a number of experimentally relevant problems (most of which had previously been explored only by numerical or approximate techniques, if at all). We also develop a CONTENTS 4 clear physical picture of the interplay between coherent interactions and spontaneous spin flips, which reveals that T 1 decoherence is much more detrimental to entanglement generation than might be naively expected. However, our solution also points toward a measurement-based feedback scheme that can mitigate its detrimental effects.
The organization of the manuscript is as follows. In Sec. 2 we consider the coherent (Hamiltonian) far-from-equilibrium dynamics of an Ising model with arbitrary spinspin couplings. Our results comprise a unified framework for calculating unequal-time correlation functions starting from arbitrary unentangled pure states. As special cases, these results will be applied to calculating Bloch-vector dynamics, arbitrary equal time spin-spin correlation functions, and two-time dynamical correlation functions. These results are substantially more general than any already available in the literature [43,44,34,45,1], and will help quantify the quantum-enhanced precision in metrology experiments using trapped ions and ultracold neutral atomic gases. In Sec. 3 we consider the effect of Markovian decoherence on this dynamics, incorporating dephasing, spontaneous excitation, and spontaneous relaxation. Because the excitation and relaxation processes do not commute with the Ising dynamics, including them is especially nontrivial: we work in the interaction picture of the Ising Hamiltonian and incorporate them as time-dependent perturbations. Terms in the perturbative expansion are evaluated using the tools laid out in Sec. 2, and by summing the perturbation theory to all orders we obtain an exact (closed-form) description of the dissipative dynamics of arbitrary two-point correlation functions. This is in stark contrast to the behavior of a coherently driven Ising model, where such a perturbative expansion cannot in general be resummed. An interesting feature revealed by our exact solution is that the spin dynamics undergoes an oscillatory-to-damped transition at a critical dissipation strength, which-in the absence of a coherent drive-cannot occur at the single-particle or mean-field level. This feature, along with the more general structure of our solution, is demonstrated by solving for spin dynamics in a nearest-neighbor Ising model. We conclude this section by casting our solution in terms of a clear physical picture, in which T 1 decoherence (spontaneous excitation/relaxation), through its interplay with the Ising dynamics, gives rise to an emergent non-local dephasing process. Section 4 applies the solution to calculating experimentally relevant observables in a dissipative version of the one-axis twisting model. We show that the emergent dephasing discussed in Sec. 3 severely diminishes the precision enhancement achievable compared to that obtained in the absence of decoherence. However, we also show that, in special cases, this degradation can be prevented by a measurement-based feedback mechanism. In Sec. 5 we summarize our results and pose a number of unanswered questions that would be interesting to address in future work.

Coherent dynamics in Ising models
Our goal is to develop a unified strategy for describing the dynamics of a collection of spin-1/2 particles interacting via Ising couplings and initially (at time t = 0) driven far CONTENTS 5 out of equilibrium. In the absence of a magnetic field, the most general form for an Ising model is ‡ whereσ z j are z Pauli matrices and the indices j, k label lattice sites located at spatial positions r j . The coupling constants J jk are left completely arbitrary, and hence there is not necessarily any notion of dimensionality. In many physical realizations of this Hamiltonian, such as trapped ions, neutral atoms, or Rydberg atoms, the couplings exhibit a roughly power-law spatial dependence, Because there is no transverse field (a term ∝ h jσ x j ), the eigenstates of H can always be chosen to be simultaneous eigenstates of all theσ z j (with eigenvalues σ z j = ±1). As a result, the partition function Z(β) = tr e −βH (with β the inverse temperature) is identical to that of a classical Ising model, and the equivalence of all equilibrium properties follows. This is the sense in which the Ising model without a transverse field is often said to be "classical" (even though it is a quantum Hamiltonian acting on vectors in a Hilbert space). In passing we note that classical Ising models can, of course, be highly nontrivial: for example, disordered or frustrated couplings give rise to classical glassiness [46,47,48]. Out of equilibrium, however, this notion of classicality is inapplicable. While a thermal density matrix ρ(β) = Z −1 e −βH commutes with H, the density matrix describing some non-equilibrium initial conditions will not in general commute with the Hamiltonian, and nontrivial dynamics will ensue. This dynamicswhich has no direct analogue in the classical Ising model-is generically characterized by the growth of entanglement, leading in some special cases to spin-states with applications in quantum information and precision metrology [34].
Everywhere in this manuscript, we assume the system starts in a pure state that is a direct product between the various spins [ Fig. 1(a)]. § The most general such state can be specified by choosing spherical angles θ j and φ j describing the orientation of the spin at each site j [ Fig. 1(b)]. Defining and states |σ j that are eigenstates ofσ z j with eigenvalues σ j = ±1, such a state can be written For uniform θ j = θ and φ j = φ, the state |Ψ(0) is frequently encountered in experiments ‡ Note that many references studying long-ranged Ising models-i.e. those for which j J ij is an extensive quantity-often normalize the interaction by dividing by the number of spins N . We drop this constant here to avoid cluttering the notation. § Unentangled but mixed initial density matrices can be easily accounted for by suitable averaging of the expressions given over the initial ensemble.  Figure 1. The only restriction on the initial state is that it be unentangled (i.e. a product state over the various sites in the spin model). For example, one can imagine (a) an initial state with some slow spatial variations in the spin angles due to inhomogeneities in the pulse strength of a Ramsey-type experiment. The notation used to characterize the state of any one spin is shown on the Bloch sphere in (b).
implementing Ramsey spectroscopy [49,33,32], and spatially varying angles could be used, for example, to describe the effects of defects or excitation inhomogeneities in such experiments [50,51,52].
Essentially all properties of the non-equilibrium dynamics are contained in unequaltime correlation functions of the spin operatorsσ ± j andσ z j (these subsume, of course, the time evolution of all equal-time correlation functions). We focus first on the case where only operatorsσ ± j occur Here a, b = ±, and the time dependence of the operators is given in the Heisenberg picture of Hσ a j (t) = e itHσa j e −itH .
The time-ordering operator T C orders all operators along a closed-time-path C shown in Fig. 2, with times t occurring on the forward path and times t * occurring along the backwards path. This closed-time path ordering occurs naturally, for instance, in any perturbative treatment of additional non-commuting terms in the Hamiltonian. In Sec. 3 we encounter this situation when treating a dissipative coupling to an environment, but the same structure occurs for coherent couplings, e.g. a transverse field h jσ x j . Our goal in what follows is to obtain analytic expressions for such correlation functions in full generality, and then apply them to calculating a variety of experimentally relevant quantities. In order to describe G concisely, it is useful to define a variable α j on each site such that α j = 1 if there are no occurrences of the operatorsσ a j in G, and α j = 0 otherwise. Now we recognize that if an operatorσ a=± j occurs in G one or more times, the operatorσ z j (appearing in the time evolution operator) is forced to take on a well defined value σ z j (t) at all points in time (see Fig. 2). As a result, we can rewrite the correlation function G as Graphical representation of a sample correlation function G = . Here α j = α k = 0, α l =j,k = 1, and the time-dependent functions σ z j (t) and σ z k (t) are shown in the bottom two panels. where (t = 0 * marking the end of the backwards trajectory, Fig. 2),f is the complex conjugate of f , C is a time integral that runs along the closed-time path, and The first thing to notice is that, since C dt = 0, the final time-independent term in K vanishes. Since this is the only non separable term in K, the remaining time evolution is straightforward to compute. It is helpful to define the following parameters that depend on the functions σ z j (t) in terms of which The time-ordered correlation functions of Eq. (7) can now be compactly written = e −iϑ with (g − j will be useful momentarily). The equivalence between Eqs. (13) and (14) can be understood by explicitly comparing the expressions inside the product for the two CONTENTS 8 possible situations α j = 0, 1: If α j = 0, then ϕ j = 0 and the expectation value is unity, whereas when α j = 1 we find p j = 0 and the expectation value gives g + (ϕ j ).
The insertion of an operatorσ z j (t) inside a correlation function G, which we denote by writing G → G z j , is relatively straightforward. If α j = 0, then clearly the substitution σ z j → σ z j (t) does the trick. If α j = 1,σ z j can be inserted by recognizing that the variable ϕ j couples toσ z j as a source term, and thus the insertion ofσ z j (t) is equivalent to applying i ∂ ∂ϕ j to G. Both possibilities are captured by writing which, using ∂g + (ϕ)/∂ϕ = −ig − (ϕ), can be simplified as Notice that if all operators occur at the same time t, e.g. when calculating equaltime correlation functions, then ϑ = 0 and ϕ j = k J jk (1 − α k )α j (±2t) (with the ± depending on whetherσ ± k is applied to the spin on site k).

Bloch vector dynamics
It is now straightforward to calculate the dynamics of the Bloch vectors Becauseσ z j commutes with the Ising interaction the z component of spin is time independent, and given by S z j = 1 2 g − j (0) = 1 2 cos θ j . The transverse spin components S x j (t) and S y j (t) can be obtained from the real and imaginary parts, respectively, of σ + j (t) . A straightforward application of Eq. (13) gives For the special case where all spins point along θ = π/2 at t = 0 there is pure decay of the Bloch vector without any rotation. In Fig. (3) we show the projection into the xy plane of the Bloch vector S 0 (t) (where j = 0 labels the central site of a 55-site triangular lattice) for an initial state in which all spins point in a single direction lying outside of the xy plane (θ = π/4, φ = 0). The Bloch vector spirals inwards: The precession can be understood as a mean-field effect [33], with the spin rotating due to the average magnetization of the other spins, while the decay is due to the development of quantum correlations. Note that, in a finite system, the length S(t) = |S(t)| of the total Bloch vector S(t) = j S j (t) decays even at the mean-field level for any ζ = 0. This decay, however, is due to the existence of a spatially inhomogeneous mean-field, and cannot be associated with the development of correlations. Figure 3. Trajectories of the Bloch vector S 0 projected into the xy plane, for allto-all (left) and dipolar (right) couplings. Here j = 0 labels the central site (green dot in left pannel) of a 55-site triangular lattice, and all spins are initialized at {θ j , φ j } = {π/4, 0}. In both cases we choose a nearest neighbor coupling J, and scale the time by J tot = j =0 J/(r j − r 0 ) ζ . This rescaling of time is used so that different range interactions give rise to comparable precession rates. Note that at mean-field level the trajectory would close on itself. The inward spiral indicates the growth of quantum correlations and resultant decay of the spin length, and-in these rescaled time units-is more significant for shorter-range interactions.

Equal time correlation functions
Spin-spin correlation functions can be calculated just as easily from Eqs. (13) and (16). All two-point correlation functions can be calculated from the four quantities and their complex conjugates. Since the Hamiltonian commutes with allσ z j , the first one is given trivially by C zz jk = g − j (0)g − k (0) = cos θ j cos θ k . The second one can be obtained from Eqs. (16) and (19) as and the third and fourth are These correlation functions can be used, for example, to calculate the time dependence of the spin squeezing parameter Here ∆S min (t) is defined to be the minimum uncertainty along a direction perpendicular to S(t) . The squeezing parameter determines the phase sensitivity in a suitably performed Ramsey experiment, which is enhanced over the standard quantum limit whenever ξ < 1 [34]. For ζ = 0 (infinite-range interactions) the calculation was first performed in [34]. However, in many experimentally relevant situations the interactions have some finite range and the maximum achievable squeezing is diminished (Fig. 4).

Unequal-time correlation functions
It is also possible to calculate correlation functions involving the application of spin operators at different times, which describe the propagation in time of a perturbation to the system. As an example, we can easily calculate dynamical response functions of the form These can be combined to calculate dynamical response functions involving arbitrary Pauli matrices, some examples of which are shown in Fig. 5. As we will see in Sec. 3, such dynamical correlation functions allow us to calculate the effect of spontaneous relaxation and excitation on the dynamics of the system.
If we choose our x-axis to be in the direction of S(t), and defineŜ ψ = 1 for a 100 site 1D chain with ζ = 1 and nearest-neighbor coupling J. In (a) we plot S yy i,i+r (0, t) for i = 50 and r = {2, 3, 4, 5} (from top to bottom). In (b) we plot S yy i,i+1 (t, t + δt) (again for i = 50) as a function of t and δt. In both plots the initial state consists of all spins pointing in the +x direction (θ = π/2 and φ = 0).
For instance, for all spins initially polarized along the x axis, the effect of the sudden relaxation of a spin on site k at time t s on the time dependence ofσ x j (for j = k) is given by Note that this is the same time evolution we would obtain if no spontaneous relaxation had occurred, the k'th spin were simply absent, and the j'th spin were coupled to a longitudinal magnetic field of strength 2J jk (2t s /t − 1). The reason for this behavior is straightforward when considering the time evolution in the Schrödinger picture. The application ofσ − k to the wave function at time t s not only forces the k'th spin to point down between times t s and t, but also destroys the piece of the initial wave function having weight into states with σ k = −1. Hence it is as if the k'th spin pointed up for a time t up = t s , and down for a time t down = t − t s , thus contributing an inhomogeneous longitudinal magnetic field of strength 2J jk (t up − t down )/t = 2J jk (2t s /t − 1) (the factor of 2 arises because the J jk couple to Pauli matrices rather than the spin-1 2 matrices).

Inclusion of dissipation
In the previous sections we have treated our Ising spins as a closed system. That is, we neglected any coupling that might exist between our system and the outside world, and thus initially pure states remained pure throughout the dynamics. In any physical realization of the Ising model, this is clearly an idealization; decoherence occurs and often must be accounted for. For example, Rydberg atoms suffer from spontaneous emission [53], while ions couple strongly to fluctuating classical (electric and magnetic) fields and can decohere through off-resonant light scattering from the spin-dependent optical dipole forces used to engineer the Ising interactions (this off-resonant light scattering can produce spontaneous excitation/relaxation and dephasing) [54]. One way to envision dynamics in an open system is by considering the probabilistic occurrence of sudden perturbations of the system-quantum jumps-due to the system-environment coupling [55]. As suggested in Sec. 2.3, and as will be explained in detail below, the strategy we have developed for computing unequal-time correlation functions is well suited to describing such effects.

Description of the problem
Given a density matrix describing a system coupled to a reservoir, it is always possible to express the expectation value of a system operatorÂ in terms of the system reduced density matrix ρ = tr R [ ] as Â = tr S [ρÂ]. Here tr R(S) denotes a trace over the reservoir (system) degrees of freedom-we will drop these subscripts from now on, since all future instances of tr refer to a trace over system degrees of freedom only. In this language, the effect of a finite system-reservoir coupling is that an initially pure system density matrix ρ(0) ≡ |Ψ(0) Ψ(0)| will evolve into a mixed state (reflecting entanglement between the system and reservoir degrees of freedom). When the system-environment coupling is weak (i.e. small compared to the inverse of relevant system time scales) and the reservoir correlation time is small, the Born-Markov approximation is justified and the reduced system density matrix obeys a Markovian master equation of Lindblad form [56]. We choose a very general master equation appropriate for describing the various types of decoherence relevant to trapped ions [54], Rydberg atoms [53], and condensed matter systems such as quantum dots [5] and nitrogen vacancy centers [6]: where The first term involving a commutator describes coherent evolution due to the Ising interaction, and the various terms having subscripts "ud", "du", and "el" correspond respectively to spontaneous relaxation, spontaneous excitation, and dephasing ¶.
Equation (28) has the formal solution ρ(t) = U (t)ρ(0), with The exponential of super-operators is meant to be understood via its series expansion, in which the multiplication of two Lindblad super-operators implies composition Our goal in what follows is to compute the time dependence of an arbitrary operatorÂ at time t, given in the Schrödinger picture by

Dephasing (T 2 decoherence)
An immediate simplification follows from the observation that That the last two commutators vanish is less obvious than the first, but physically it has a very clear meaning: Spontaneous relaxation/excitation on a site j causes the j'th spin to have a well defined value of σ z j , and thus to be unentangled with the rest of the system. Since the dephasing jump operatorσ z j changes the relative phase between the states |σ z j = ±1 , whether spontaneous relaxation/excitation occurs before or after a dephasing event only affects the sign of the overall wave function, which is irrelevant. As a result, we can write and the time dependence of an arbitrary observable A(t) = Tr ρ(t)Â can be written + Note that application of a superoperator does not commute with operator products , and we use the convention that a superoperator should act on the operator immediately to its right. The effect of the time evolution due to L el can be understood by considering its effect on the Pauli operators: In light of equations (37) and (38), we are free to ignore the dephasing terms in the master equation at the expense of attaching a factor of e −Γ el t/2 to every operator σ x,y j occurring inside an expectation value: where ρ(t) on the right-hand side evolves under the master equation without the dephasing term. + Note that we apply the operator e −tL el toÂ rather than ρ(0). This is justified by the identity tr L el (Ô 1 )Ô 2 = tr Ô 1 L el (Ô 2 ) , true for arbitrary operatorsÔ 1 andÔ 2 , which holds because the jump operatorsσ z j are hermitian.

Spontaneous relaxation and excitation (T 1 decoherence)
Because the effects of dephasing are fully included by Eq. (39), the remaining problem is to compute the time dependence of operators whose expectation values are taken in a density matrix evolving simultaneously under H , L ud , and L du : Formally the challenge of including the effects of spontaneous relaxation and excitation is related to the nontrivial commutation relation Physically, the obstacle is that spontaneous relaxation and excitation change the value of σ z for the spin which they affect, and this change feeds back on the system through the Ising couplings. From the results on coherent dynamics presented in Sec. 2, we know that time evolution under H alone is tractable. This suggests that we attempt to solve Eq. (40) by rewriting where R contains all and only terms that do not commute with H , and then doing perturbation theory in R. The above separation is accomplished by defining an effective (non-Hermitian) Hamiltonian and its corresponding superoperator and the recycling term in terms of which the time evolution operator is U (t) = e −t(iH eff −R) . In Eq. (43) we have defined γ = 1 4 (Γ ud − Γ du ) and Γ r = Γ ud + Γ du . Defining U 0 (t) = e −itH ef f , we can now expand the time evolution operator as a power series in R in order to obtain the time-dependent expectation value A(t): This expansion is the underlying object being evaluated when Monte Carlo wave function methods [55] (quantum trajectories) are used to approximate the density matrix. In Appendix A, we show in detail how the series in Eq. (45) leads to an expression for A(t) in terms of the closed-time path ordered correlation functions obtained in Sec. 2, and the summation of that series is carried out in Appendix B. Here we will simply summarize the calculation, and explain in physical terms the essential structure of the Hamiltonian and decoherence that allows the result to be cast in closed form. We begin by noting that when writing A(t) as a sum over closed-time path ordered correlation functions, each operator inserted along the forward leg of the time-contour is accompanied by its hermitian conjugate appearing at the same time along the backward part of the time contour. If we had explicitly included an environment and attempted to trace over it (rather than starting with a Markovian master equation), this feature of the problem would emerge as a direct consequence of the Markov approximation. As a result, it suffices to describe any term in the series expansion by specifying the occurrence of operators on the forward time contour. To facilitate this description, we introduce notation describing the occurrence of operators belonging to a particular site j (see Fig. 6 for a summary). We take R ± j to be the number of times the operatorσ ± j occurs along the forward time path, {t j 1 , . . . , t j R j } to be the set of times at which jump operators are applied to site j, R j = R + j + R − j , κ j = ±1 depending on whether the operator at the latest time along the forward path isσ ± j , and We will also use bold symbols R, κ, and τ to specify the complete set of these variables on all lattice sites. Note that specifying R j and κ j determines both R + j and R − j , since two consecutive (in the closed-time-path-ordering) applications of an operatorσ ± j gives zero. Therefore, we will only include the R j and κ j as explicit arguments in the correlation functions below. These variables are sufficient to determine the value of any term in the series expansion of A(t), so we do not need to keep track of the individual times at which each jump operator is applied.  Figure 6.
A few examples of how spin-raising and spin-lowering operators belonging to the j'th site may occur along the forward time evolution, and the notation used to characterize these occurrences.
All nonzero terms in Eq. (45) are captured by summing over the R j and κ j , and integrating over the times t j 1 , . . . , t j R j , denoted IfÂ can be written as a product of operatorsσ b j j (b j = ±) on sites j contained in a set η, then A(t) can be compactly expressed (see Appendix A) as The prefactor is closely related to the probability that a series of jumps described by the variables R, κ, and τ has occurred. Defining the symbol s as the vector of quantities s j = ϕ j /t−2iγ, the correlation functions are of the general form presented in Sec. 2. Careful bookkeeping reveals that the variables ϕ and ϑ defined in Sec. 2 are given by Note that the R j dependence in G j is hidden in the implicit dependence of p j and g + j on α j (which for j / ∈ η satisfies α j = δ R j ,0 ). When evaluating the sums and integrals in Eq. (48), one must keep in mind that, whenever j ∈ η, terms with R j = 0 vanish because they contain the consecutive application of either σ + j or σ − j . Physically, the vanishing of such terms reflects the lack of coherence for any spin that has undergone even a single spontaneous spin flip.
The factorization of P into functions P j of local site variables {R j , κ j , τ j } is a direct consequence of the single particle nature of the anti-Hermitian part of H eff . Physically, this factorization occurs because the dissipation we are considering is uncorrelated from site to site (in contrast to the collective relaxation processes that arise, e.g., in the context of Dicke superradience [57]). The factorization of the correlation function G into functions G j of local site variables {R j , τ j } is a more surprising result, and depends crucially on the occurrence of jump operators at the same times on the forward and backward time evolution (Appendix A). The γ appearing in the argument of g + j (s j t = ϕ j − 2iγt) affects the value of any term in Eq. (48) in which no jump operators are applied to site j (such that α j = 1 and therefore g + j (s j t) = 0). It arises from the term −iγσ z in H eff [Eq. (43)], and decreases the expectation value ofσ z j for γ > 0 (when spontaneous relaxation outweighs spontaneous excitation). This effect is often referred to as null measurement state reduction: Gaining knowledge that the spin on site j has not spontaneously flipped affects the expected value when measuringσ z j .
Because G and P both factorize, we need only to evaluate the quantities in terms of which The explicit dependence on s j is included to remind the reader that, after the sums and integrals (over R j , κ j , and t j 1 . . . t j R l ) have been carried out, s j is the only site-dependent quantity on which Φ j (s j , t) depends. We evaluate these sums and integrals in Appendix B, obtaining where r = Γ ud Γ du . IfÂ also contains an operatorσ z l (t) (with l / ∈ η), we denote its expectation value by A z l (t). The insertion ofσ z l must be dealt with at the point of Eq. (54). Keeping in mind the discussion surrounding Eq. (16), and remembering that s l = ϕ l /t − 2iγ, we must replace Φ l (s l , t) with Therefore we have and in Appendix B we find Ψ l (s l , t) = e −Γr/2 cos t s 2 j − r + is l + 2γ − Γ r 2 cos θ l t sinc t s 2 j − r .

A simple application: under-damped to over-damped transitions
These equations reveal that correlation functions will generally undergo a qualitative transition in dynamics-from over-damped to oscillatory-whenever the condition s 2 l = r is satisfied. This behavior is the most clearly manifest when the couplings J ij have a simple structure, such as nearest neighbor or all-to-all. For instance, for nearestneighbor coupling in 1D, assuming {Γ el , Γ ud , Γ du } = {0, Γ, Γ}, and choosing the initial state to point along the x−axis, we find (ignoring boundary effects) which becomes critically damped at Γ c = 2J (see Fig. 7). It is interesting to note that this solution is similar in structure to the damping of a classical harmonic oscillator or a coherently driven two-level system (c.f. the weak-coupling limit of the spin boson problem [58,14]). It is important to contrast this behavior with that of a single spin coupled to a Markovian bath, where the decoherence we consider only causes a dampedto-oscillatory transition in the presence of a transverse magnetic field; the Hamiltonian dynamics must be able to restore coherence in the basis for which the environment induces a measurement. In the present case, there is no transverse field, and it is not a priori obvious that such behavior should emerge. In fact, a simple mean-field estimate of the dynamics fails to capture the oscillatory-to-damped transition. Using a sitefactorized ansatz for the density matrix, ρ = j ρ j , it is straightforward to see that [33] S x MF (t) = N 2 e −Γt cos(4J cos(θ)t).
Thus mean-field theory, which assumes an unentangled density matrix, always predicts under-damped dynamics; the transition to over-damped behavior captured by the exact solution depends crucially on the competition between decoherence and entanglement.

Qualitative insights into decoherence in interacting many-body systems
In addition to providing an efficient way to compute arbitrary observables for an open many-body system, the above calculation provides significant insight into the interplay beween decoherence and interactions. To keep the notation as simple as possible, we will focus on the case of nearest-neighbor coupling on a lattice with coordination number z, and calculate the time dependence of the transverse spin length of a single spin s x (z, t) = σ x j (for an infinite system it is independent of j) for a state in which all spins are initially polarized along the x-axis. In the absence of decoherence but in the presence of a longitudinal magnetic field of strength h, application of results in Sec. 2 gives In the absence of a longitudinal field but including equal rates of spontaneous relaxation/excitation (Γ ud,du ≡ Γ, γ = 0), we find instead [from Eqs. (55) and (56)] However, it is instructive to temporarily hold off evaluating the sums and integrals implicit in the Φ j , and work directly with Eqs. (54) and (55): In Eq. 64 we have divided the product over lattice sites into three parts: the j'th site, the nearest neighbor sites (product labeled by NN) and the rest of the lattice (product labeled by D, reminding us that these sites are Disconnected from site j). First, we observe that Φ j (s j , t) = 1 2 e −Γrt/2 . Next we evaluate D l Φ(s l = 0, t) = 1 [the s l = 0 because these sites are disconnected from the j'th site, and hence J jl = 0, see Eq. (51)], which follows from a sum rule l P l (R l , κ l , τ l )G l (0, R l , τ l ) = 1 (derivable from tr [ρ(t)] = 1). It remains only to evaluate where we've adopted an abbreviated notation NN = NN j j for the sums and integrals over the nearest-neighbor sites. Utilizing * where R = NN j R j , τ = NN j τ j , and h = Jτ /t, and noting that s x h (z − R, t) only depends on τ and R (and not any other combinations of the local site variables), we can then write Here is obtained by carrying out the sums and integrals while holding R and τ fixed (the factor of t/J arises when changing variables from τ to h in the remaining integral). Equation (67) has a very suggestive form. The factor of e −Γrt/2 out front is the single-particle contribution of T 1 decoherence, and would be present in the absence of interactions; it reflects the probability that the j'th spin has not spontaneously flipped before the time t (a prerequisite for having any coherence along x). As the z nearest neighbors evolve in time, they can undergo spin flips which cause them to fluctuate in time between pointing along +z and −z [ Fig. 8(a)]. Even once flipped, they influence (via the Ising couplings) the j'th spin in a manner formally equivalent to a longitudinal magnetic field of strength h = J(τ /t) [ Fig. 8(b)]. The quantity P (R, h) describes the probability that R nearest neighbors have spontaneously flipped, collectively contributing an effective magnetic field of strength h to the time evolution of the j'th spin. The sum over R and integral over h then average the resulting dynamics for the j'th spin, s x h (z − R, t), over the possible behaviors of its neighbors. For a given R, the integral over h reduces the Bloch vector length by an amount depending on the width in h of the distribution P (R, h). Physically, this integral captures the phase diffusion of the j'th spin due to the stochastic temporal fluctuation of its immediate environment (neighboring spins, as shown in Fig. 8). Thus we see very clearly that the interplay between spontaneous spin flips and coherent interactions leads to an emergent source of dephasing: Flipping spins act as fluctuating magnetic fields (mediated by the Ising interactions) on other spins, even if these latter spins have not been directly affected by decoherence.

One-axis twisting in an open system
The expressions in Eqs. (56) and (59) furnish a complete description of correlation functions, and in special cases afford descriptions of common experimental observables and entanglement witnesses. As a concrete example, we will use these expressions to study the development and loss of entanglement in an open-system version of the oneaxis twisting (OAT) model. It is important to keep in mind, however, that most of the following results can be generalized to take into account arbitrary Ising couplings. Defining the collective spin operatorŜ z = 1 In (a) we show the optimal squeezing (optimized over angle ψ) as a function of time.
In (b) we plot the normalized variance, evaluated at time t s , as a function of the squeezing angle ψ. Note that at the angle of maximal squeezing (ψ ≈ 0), dephasing is much less detrimental than spontaneous relaxation/excitation, whereas both have a similar effect at the angle of maximal antisqueezing (ψ ≈ π/2). In all plots, we choose Γ = 1/t s , with t s = N −2/3 (3 1/6 /2J) being the time of optimal squeezing in the absence of decoherence.
which is the ζ = 0 limit of our more general Ising model [Eq. (1)]. For an initial state polarized along the x-axis (θ = π/2, φ = 0), it is well known [34] that the OAT Hamiltonian generates spin squeezed states at short times [the squeezing is optimal at t s = N −2/3 (3 1/6 /2J)], a fact that has been exploited in a number of beautiful experiments [35,36]. In principle (i.e. in the absence of any decoherence or other imperfections), these spin squeezed states allow for precision metrology with a phase sensitivity that scales as N −5/6 , thus beating the N −1/2 scaling of the standard quantum limit. At time t * = π/4J, the OAT Hamiltonian gives rise to a GHZ (or Schrödinger cat) state [59], which in principle affords Heisenberg limited (∼ N −1 ) sensitivity in phase estimation. In the following subsections, we use the results of Sec. 3 to extend calculate spin squeezing and characterize the metrological utility of (and GHZ-type entanglement of) the state at t*, which would be the GHZ state in the absence of decoherence.

Spin Squeezing
Given the results in Sec. 3, analytic calculation of the squeezing parameter in the presence of arbitrary decoherence rates Γ el , Γ ud , and Γ du is now straightforward. As can be seen in Fig. 9(a), the effect of T 2 decoherence is much less severe than that of T 1 decoherence. One reason for this behavior is that the minimum variance ∆S 2 min occurs at an angle ψ (in the yz plane) only slightly deviating from the z axis. Therefore dephasing, which can be thought of as random rotations of the individual spins around the z-axis, does not introduce much noise in the squeezed quadrature (see Fig. 9b). To the contrary, spontaneous relaxation/excitation processes introduce noise directly into the squeezed quadrature.

Macroscopic superposition states
In the absence of decoherence, the OAT Hamiltonian is known to give rise to N -spin GHZ states at a time t * = π/4J, where | ⇑ x (| ⇓ x ) denotes the state where all spins point along the positive (negative) x-axis [59]. These entangled states afford Heisenberg-limited phase sensitivity [60], and are a resource for certain types of fault-tolerant quantum computation ( [37] and references therein). However, they are also a canonical example of a fragile quantum state, and their usefulness is easily destroyed by decoherence [61]. The effect of dephasing on the production of GHZ state via one-axis twisting is well understood [61,62]. With the results of Sec. 3, however, we can easily calculate the effects of dephasing and spontaneous relaxation/excitation on the production of a GHZ state by one-axis twisting in a unified way. In this section we explicitly compare the effects of dephasing to the those of pure spontaneous relaxation ({Γ el , Γ ud , Γ du } = {0, Γ, 0}). We first characterize the GHZ state by its phase coherence, obtained from the expectation value C = tr ρ(t)Ĉ of the operator Strictly speaking the form given only applies when the particle number N is even. For odd N the GHZ state created looks similar, but is composed of states polarized along the ±y direction.
The quantity C characterizes the extent to which the superposition between the macroscopically distinct states | ⇑ x and | ⇓ x is quantum mechanical (rather than a classical mixture). Formally, C serves as a witness to N -particle entanglement of the GHZ type, with entanglement guarantied whenever |C| > 1/4 is satisfied [63,37]. Application of the results in Sec. 3 yields and in the absence of decoherence one finds |C(t * )| = 1/2. As can be seen in Fig. 10, the effect of spontaneous relaxation on the coherence C is comparable to (but worse than) the effect of dephasing. Both types of decoherence cause a loss of phase coherence when Γ ∼ J/N , with the factor of N responsible for the fragility of a GHZ state composed of a large number of spins. We can also directly calculate the metrological usefulness of a GHZ state prepared in the presence of spontaneous relaxation. The favorable sensitivity of the state ρ * ≡ ρ(t * ) to rotations by angle Ω around the x-axis can be understood as the strong dependence of the expectation value of the parity operatorπ = jσ where ∆P (Ω) = 1 − P (Ω) 2 (taking into consideration thatπ 2 = 1) is the uncertainty of the operatorπ calculated in the state ρ * . The approximation in Eq. (75) is simply that P (0) ≈ 0 and therefore ∆P (0) ≈ 1. This can be checked explicitly by looking at the large N limit of In the absence of decoherence, M = N and the Heisenberg limit of phase sensitivity is obtained. By generalizing Eq. (58) to the case whereσ z j is inserted on N − 1 of the sites, we find that in the presence of decoherence This result is plotted in Fig. 11 for different types of decoherence. The enhancement in M survives T 1 decoherence only if N Γ r t * 1 is satisfied (the scaling by N is shown in the inset). This result should be contrasted with the effect of T 2 decoherence ({Γ el , Γ ud , Γ du } = {Γ, 0, 0}). In this case Ψ(2J, t * ) = 1, yielding and hence the precision enhancement decays on a timescale that is independent of N (consistent with results in Ref. [62]). In contrast, the entanglement witness C decays at an N -enhanced rate for either type of decoherence.

Removing the effects of decoherence via measurement-base feedback
In Sec. 4.2, we showed that spontaneous relaxation significantly degrades the precision enhancement of a GHZ state unless Γ J/N is satisfied. In this section, we will show that a time-resolved record of spontaneous relaxation events provides sufficient information to restore the phase enhancement under the much less stringent constraint Γ J. In particular, we are imagining a situation where spontaneous spin flips are accompanied by the real spontaneous emission of a photon, such that they can be measured by photodetection.
For reasons that will become clear in what follows, we take our initial state to have all spins pointing at an arbitrary angle θ (rather than θ = π/2, as assumed in Sec. 4.2). Our goal is to evaluate the expectation value of the operator jσ y j l =jσ z l at time t * , which was accomplished above by appealing directly to Eqs. (59) and (56). Pursuing a strategy similar to that employed in Sec. 3.5, we first observe that n spins initialized at θ = π/2, evolving in the absence of decoherence but in the presence of a longitudinal field of strength h, would yield the phase-sensitivity enhancement † † M (n, h) = n Im e −2iht * (i sin 2Jt * ) n−1 (79) = |sin(2ht * + π(n − 1)/2)| .
Now we calculate M in the presence of decoherence, but we hold off evaluating the multiple sums and integrals that yield the functions Ψ l (s l , t) in closed form, and instead obtain M by working directly with Eq. (48) First, we evaluate G(s, R, τ ), obtaining which only depends on the R j and τ j only through their sums R = j R j and h = (J/t) j τ j . With the judicious choice θ = π − 2 tan −1 e 2γt * , we can rewrite and thus we obtain The initial value of θ, which places the initial spins slightly above the xy plane, was carefully chosen so that g − (2Jt − 2iγ) ∝ sin(2Jt). This finite tipping angle is required to precisely cancel the null measurement effect, which causes the z-projection of each spin to change in time due to the lack of emission of a photon, and thus ensures that at time t * the unflipped spins are brought down into the xy plane. Finally, we can write where P (R, h) is obtained by carrying out in Eq. (84) with R and τ = ht/J held fixed. As in section 3.5, we can interpret P (R, h) as the probability to have R flipped spins contributing an effective magnetic field h to the dynamics of the remaining spins. For each particular value of R and h, the function M (n, h) (where n = N − R) yields the precision enhancement of a GHZ state produced from n spins in a longitudinal field of strength h.
If an experiment can record the times (t 1 , . . . , t R ) at which photons are emitted, thus gaining access to both R and then that experiment produces the conditional density matrix ρ(n, h), corresponding to a GHZ state of n spins produced by one-axis twisting in the presence of the longitudinal field h. The effect of this field is simply to rotate the system around the z-axis by a (shot-to-shot random) angle δ = (ht * +π(n−1)/2). However, rotation by a random angle only causes decoherence if the value of that angle is not known. Because the experiment measures δ indirectly via the photon emission record, it is possible to remove the effect of h in any experimental shot by applying the rotation operator R(δ) = exp (−iS z δ) to the conditional density matrix ρ(n, h). In this way we create which is GHZ state of the form in Eq. (70) containing n spins, and thus obtain a precision enhancement of n. Because the expected value of n will decay only at the bare rate Γ, there is no longer an enhancement by N in the decay of precision. This measurementbased coherent feedback can also be applied in the context of spin-squeezing, where once again it can vastly improve the metrological usefulness of a state generated by one-axis twisting in the presence of spontaneous relaxation. Before concluding, we note that this feedback strategy could, in principle, be applied to situations where both spontaneous relaxation and spontaneous excitation are present, and even when the coupling constants J ij are not uniform. However, in the former case it is necessary to independently record the photon emission record corresponding to excitation and relaxation processes, and in both cases it is necessary to obtain siteresolved (in addition to time-resolved) information about the photon emissions.

Conclusions and future directions
In this paper we have presented a comprehensive theoretical toolbox for understanding far-from equilibrium dynamics in Ising models both with and without decoherence. The underlying objects of interest are unequal-time correlation functions, which are then used to compute spin squeezing, dynamical response functions, entanglement witnesses, and the effects of dephasing, spontaneous excitation, and spontaneous relaxation on the system dynamics. We believe these tools will be of fundamental importance in understanding and optimizing a diverse array of systems in which entanglement is engineered by Ising interactions. In particular, these tools enable the quantification of detrimental effects due to system-environment coupling, even when the coupling does not commute with the Ising interactions.
The ability to compute dynamics in any dimension and in the presence of noncommuting noise is a particularly surprising result; it is well known that the inclusion of non-commuting but coherent linear couplings admits solutions only in highly specialized geometries, such as 1D nearest neighbor chains. The key structures that allow our solution for the open system to proceed are (1) the statistical independence of decoherence processes on different sites and (2) a symmetry between the forward and backward time evolution along a closed-time path, i.e. the Markov approximation. It would be interesting to understand to what extent this simplification generalizes to other models where the incorporation of decoherence would-at first sight-appear to be intractable.
We gratefully acknowledge Salvatore Manmana, Michael Kastner, Dominic Meiser, Minghui Xu, and Murray Holland for helpful discussions. This work was supported by NIST, the NSF (PIF and PFC grants), AFOSR and ARO individual investigator awards, and the ARO with funding from the DARPA-OLE program. K. R. A. H. and M. F.-F. thank the NRC for support. This manuscript is the contribution of NIST and is not subject to U.S. copyright.

Appendix A.
The main goal in this Appendix is to explicitly cast the series expansion for arbitrary observables in terms of the time-ordered correlation functions encountered in Sec. 2 of the text, thus bridging the gap between Eqs. (45) and (48) of the text. The summation of the series is carried out later in Appendix B.
Our starting point is the series expansion for the time-evolution superoperator This leads immediately to the expression for A(t) given in the manuscript [Eq. (45), reproduced here for convenience]: In order to simplify notation in the following equations, we define time dependent jump operators in the Heisenberg picture of the effective Hamiltoniañ  we can now express A(t) as tr Â e −itH effJ (j n , a n , t n ) . . .J (j 1 , a 1 , t 1 )ρ(0)J † (j 1 , a 1 , t 1 ) . . .J † (j n , a n , t n )e itH † eff = J † (j 1 , a 1 , t 1 ) . . .J † (j n , a n , t n )Ã(t)J (j n , a n , t n ) . . .J (j 1 , a 1 , t 1 ) . (A.5) In the above expression, the time dependence of the operatorÃ is defined as which is distinct from the time dependence assigned to the operatorsσ a j in defining theJ (j, a, t) (of course this distinction vanishes when considering time evolution under a Hermitian Hamiltonian). In the final line the trace has been removed (upon rearrangement, it becomes a completeness identity), and the expectation value is in the initial pure state |ψ(0) .
The correlation functions in Eq. (A.5) are explicitly closed-time path ordered, with the time ordering enforced in the limits of integration. The remaining task is to explicitly separate all time-dependence of the correlation functions in Eq. (A.5) due to the anti-Hermitian part of H eff , so that we can directly employ the results from Sec. 2. Because the anti-Hermitian part of the effective Hamiltonian commutes with H, this is relatively straightforward. As in the manuscript, we restrict ourselves at this point to considering operatorsÂ that can be written as products of spin-lowering and spin-raising operatorŝ σ b j j (b j = ±1) on sites j ∈ η. We can then writẽ , we are essentially ready to read off the results of a given term from the expressions for correlation functions in Sec. 2. However, anticipating that the terms in the series will be site-factorizable, we pause here to introduce some notation describing the occurrence of jump operators belonging to a particular site. Because jump operators always occur at the same time along the forward and backward evolution (and each operator on the backward path is the complex conjugate of a corresponding operator on the forward path), we only need to describe the jump operators applied during the forward time evolution. First, we define R ± j to be the total number of jump operators of typeσ ± j applied to site j along the forward time evolution, and R j = R + j + R − j . We also define {t j 1 , . . . , t j R j } to be the set of times at which those jump operators are applied to the site j, and to be the total amount of time the j'th spin has spent pointing up minus the time it has spend pointing down, again during the forward evolution only (note that C σ z j (t) = 0 ∀j = η). Finally we use the bold symbols R, κ, τ to represent vectors of the quantities R j , κ j , and τ j , respectively. Now we can write A(t) as A(t) = P(R, κ, τ ) σ −a 1 j 1 (t 1 ) . . .σ −an jn (t n )Â(t)σ an jn (t n ) . . .σ a 1 j 1 (t 1 ) exp where P(R, κ, τ ) = j P j (R j , κ j , τ j ) (A.14) The exponential of operatorsσ z j Eq. (A.13) leads to the so-called null measurement state reduction: it causes the j th spin to drift out of the xy plane in the event that it has not spontaneously relaxed (α j = 1).
The expectation value in Eq. (A.13) is now precisely of the form given in Eq. (14), with the exception that we must map ϕ j → s j t ≡ ϕ j − 2iγt in order to account for the term in square brackets. Thus we obtain G(s, R, τ ) ≡ σ −a 1 j 1 (t 1 ) . . .σ −an jn (t n )Â(t)σ an jn (t n ) . . .σ where s is a vector of quantities s j = ϕ j /t − 2iγ. Careful bookkeeping reveals that which allows us to factorize G(s, R, τ ) = j G j (s j , R j , τ j ) (A.21) = j e −iϑ j (τ j ) p j + g + j (s j t) , (A. 22) thus completing the derivation of Eqs. (48)(49)(50) in the manuscript.

Appendix B.
In Sec. 3 we encountered the functions Φ j (s j , t) = j P j (R j , κ j , τ j )G j (s j , R j , τ j ) (B.1) which we now show how to evaluate. When j ∈ η, only the R j = 0 term in j survives, and we obtain Φ j (s j , t) = e −Γrt/2 p j . (B. 3) The function Ψ j (s j , t) is only defined for j / ∈ η, so this case does not apply to it. We next consider the case j / ∈ η, and drop the site index j (the derivation does not depend on the specific site j, though the answer will depend on s, which can be reindexed at the end). We begin with Φ(s, t), which can be simplified as   We will carry out the integrals first, for which it is convenient to treat the cases R even and R odd separately. Remembering that τ = (1 − α) Here the j µ are spherical Bessel functions, and the integral has been carried out by changing variables to allow evaluation of all but one integral while τ is held fixed, and then evaluating the remaining integral over τ (this is where the Bessel functions In terms of these, we finally obtain