Relaxation timescales and decay of correlations in a long-range interacting quantum simulator

We study the time evolution of correlation functions in long-range interacting quantum Ising models. For a large class of initial conditions, exact analytic results are obtained in arbitrary lattice dimension, both for ferromagnetic and antiferromagnetic coupling, and hence also in the presence of geometric frustration. In contrast to the nearest-neighbour case, we find that correlations decay like stretched or compressed exponentials in time. Provided the long-range character of the interactions is sufficiently strong, pronounced prethermalization plateaus are observed and relaxation timescales are widely separated. Specializing to a triangular lattice in two spatial dimensions, we propose to utilize these results for benchmarking of a recently developed ion-trap based quantum simulator.


Introduction
When a physical system is coupled to a heat bath, one expects to observe thermalization to an equilibrium state whose temperature is determined by the bath properties. For an isolated many-body system, i.e. in the absence of a heat bath, the situation is less clear, although some kind of relaxation to equilibrium may be expected for sufficiently large generic systems and suitable observables. Recent progress in experiments with cold atoms and ions [1,2,3] has stimulated intense theoretical interest in equilibration and thermalization behaviour of isolated many-body quantum systems [4,5]. General mechanisms leading to thermalization have been proposed [6,7,8], rigorous proofs of equilibration have been obtained for generic Hamiltonians [9,10,11,12,13,14], and analytic as well as numeric model studies have been reported (see [5] for a list of references). Much less is known about the timescales on which relaxation to equilibrium takes place. In this paper we study the time evolution of correlation functions in isolated long-range interacting quantum Ising systems on arbitrary lattices. We derive exact analytic results that further our understanding of different relaxation timescales and general non-equilibrium properties of long-range interacting systems. In particular, we find that for sufficiently long-range interactions, a second relaxation timescale occurs, widely separated from and substantially slower than the timescale on which single spin observables approach equilibrium. As a consequence, pronounced prethermalization plateaus are observed.
Apart from broadening the theoretical understanding, our results may also prove beneficial for the interpretation of data from cold atom or ion experiments where longrange Ising interactions can occur. For example, our results can be used to benchmark a recently developed ion-trap based quantum simulator consisting of several hundred beryllium ions stored in a Penning trap and confined to a single plane [15]. Due to their mutual electrostatic repulsion, the ions arrange into a two-dimensional Coulomb crystal on a triangular lattice (see figure 1). The valence electron spins of the 9 Be + ions are the relevant degrees of freedom used for quantum simulation. Effective interactions between these spins are induced by transverse motional modes of the Coulomb crystal. Under the assumption of small, coherent ion displacements, it was shown in [16,17] that the resulting interactions are described by the Ising Hamiltonian Here i and j label the N ions on the triangular lattice, σ i = (σ x i , σ y i , σ z i ) denotes the vector of Pauli matrices for ion i, the J i,j are coupling coefficients, and B µ is an effective magnetic field generated by externally applied microwave radiation. Unlike in the conventional Ising model [18], spin-spin coupling is not restricted to nearest neighbours on the lattice, but extends over all pairs of ions. The coupling coefficients J i,j can be expressed in terms of the transverse phonon eigenfunctions of the lattice and a few other experimental parameters (equation (4) of [19]). A numerical evaluation of that expression for the given lattice geometry shows that the approximation J i,j ∝ D −α i,j holds to a very good degree [15]. Here D i,j denotes the Euclidean distance between sites i and j on the lattice, and α is an exponent that can in principle be tuned within the range 0 α 3 by varying the difference in frequency of two off-resonant lasers used in the experiment. The absolute values and even the signs of the coefficients J i,j can also be tuned, allowing for the investigation of ferromagnetic as well as anti-ferromagnetic couplings. In the latter case, due to geometric frustration on the triangular lattice, spin liquids and other exotic quantum phases may possibly occur. The results can also be applied to the linear radio-frequency ion trap quantum simulators [20,21,22] and, along similar lines, benchmarking of quantum magnets emulated by means of ultracold molecules is possible [23].

Equilibration and thermalization of isolated quantum systems.
Equilibration and thermalization are two related, but not equivalent, notions concerning the long-time evolution of a dynamical system. Quoting Linden et al. [10], 'a system equilibrates if its state evolves toward some particular state [. . . ] and remains in that state (or close to it) for almost all times'. These authors then continue their description of equilibration by pointing out that the corresponding equilibrium state need not necessarily be a Gibbs state or any other special state, and it even may depend on the initial state of the system. In contrast, thermalization is a stronger notion. It requires that the equilibrium state attained does not depend on the details of the initial state, but at most on a few relevant parameters (like the total energy or the bath temperature), and that the equilibrium state is a Gibbs state (microcanonical, canonical, or possibly a generalization thereof). The calculations reported in the present paper are valid only for a certain class of initial states as specified in section 4, and we are therefore in no position to make claims about initial state independence. For this reason, we will speak only about equilibration in the following.
It is important to note that the definition of equilibration does not require the timeevolved state (or density operator) to converge to the equilibrium state in the long-time limit: It is perfectly admissible to have fluctuations around the equilibrium state that do not fade away in the long-time limit, as long as such fluctuations are either sufficiently small, or large but very rare. In any quantum system on a finite-dimensional Hilbert space, time evolution is periodic or quasi-periodic. As a consequence, the aforementioned large but rare fluctuations will inevitably occur in the form of recurrences (or Loschmidt echos). We will observe such quasi-periodic behaviour for finite-system correlation functions of the long-range Ising model in section 4.
From a technical point of view, there are (at least) two different ways in which equilibration in the presence of fluctuations can be studied. One possible approach is probabilistic and has become known under the name of typicality [24,9,10,11,12,25,13,14]. This approach amounts to studying the probability in time that the time-evolved state differs (in operator norm) from some equilibrium state by more than some small ǫ > 0. If this probability is less than some small δ(ǫ) > 0, the system is considered as equilibrated. A second approach of dealing with fluctuations consists in taking a suitable infinite-system limit, thereby pushing the recurrence times to infinity and simultaneously making the amplitude of small fluctuations vanishingly small. We follow this second strategy in the present work.
Finally, a comment is in order on the dichotomy of closed versus open systems. We consider in the following a closed (isolated) spin system, in the sense that no thermal bath or other reservoir is coupled to the system. However, for such a closed system, we do not study the relaxation to equilibrium in the above described sense of density operators being close to the equilibrium state. Instead, we restrict our attention to the long-time behaviour of spin-spin correlation functions between lattice sites i and j. This can be viewed as studying relaxation to equilibrium in a closed system, but only for a restricted class of observables. Alternatively, since such correlation functions are fully determined by the reduced density operator of sites i and j, we can consider these two sites as an open system, coupled to a bath consisting of all the spins on the remaining lattice sites. We will take advantage of this open-system point of view when investigating dephasing and purities in section 6.

Long-range quantum Ising model.
For the theoretical analysis of relaxation times and the decay of correlations, we study a model which is more general than the Hamiltonian (1) in some respects, and more restricted in others: We generalism to lattices of arbitrary spatial dimension d and arbitrary lattice structure, but come back to the two-dimensional triangular lattice later in this paper. For our analytic calculations, the effective magnetic field B µ in (1) is required to point in the z-direction, yielding a Hamiltonian of the form where the index ℓ indicates the presence of a longitudinal field B. At first sight it may look as if the resulting problem is purely classical, as all terms in H ℓ commute with each other, and equilibrium properties such as the partition function are indeed identical to those of the corresponding classical Ising model. For non-equilibrium calculations, however, quantumness enters through observables and/or initial density operators that do not commute with the Hamiltonian (2), and indeed one can prove that entanglement and other genuine quantum properties are generated under the time evolution of (2). Observables consisting only of σ z i operators commute with the Hamiltonian and therefore show no relaxation behaviour, but other observables do, as we will confirm in the following. In spite of its simplicity, the model (2) exhibits non-trivial and remarkably rich non-equilibrium behaviour.
In equilibrium, the long-range interacting Ising model (2) with power law decaying interaction strength J i,j ∝ D −α i,j is known to undergo a transition from a ferromagnetic phase at low temperature to a paramagnetic phase at high temperature. In spatial dimension d 2 such a transition occurs for all nonnegative values of the exponent α, whereas for d = 1 the transition is present only for α 2 [26]. Remarkably and surprisingly, the long-time dynamics investigated in this paper turns out to be independent on whether the corresponding energies are situated in the low-or highenergy regime.

Time evolution of correlation functions.
To study the relaxation to equilibrium in this model, our aim is to compute the time evolution of spin-spin correlation functions Here and in the following, we set = 1 for convenience. For analytical calculations, a significant simplification is achieved by restricting the initial states to density operators ρ 0 that are diagonal matrices in the σ x tensor-product eigenbasis, where 1 denotes the identity operator on the tensor product Hilbert space H = C 2 ⊗ · · · ⊗ C 2 . The indices i, j, k, l in (4) are summed over the lattice sites. The set of all ρ 0 as defined in (4) is is a very large class of initial states, parametrized by 2 N − 1 real continuous parameters s . . The reader may wish to compare the size of this class to the much studied quantum quenches, a rather restrictive class of initial states parametrized by only two parameters (namely the quench parameter before and after the quench). We could, as a matter of fact, extend the class of initial states (4) even further, allowing also σ y -type contributions to ρ 0 . This extension essentially leaves the results of the present paper unaltered, but the presentation becomes somewhat more involved. We decided to discuss this generalization in a future paper devoted to the more technical aspects of long-range quantum spin models [27].
The idea of considering initial states of the form (4) has been used in [28,29,30] for the computation of expectation values of one-spin observables. In this paper, these restrictions are relaxed and the exact analytic calculations are extended to correlation functions for arbitrary finite system sizes. Details of the calculation are reported in Appendix A. For vanishing longitudinal magnetic field B = 0, the final results are with These expressions are valid for arbitrary coupling constants J i,j , and therefore apply in arbitrary spatial dimension and on arbitrary lattices. Besides an overall factor σ x i (0) or σ x i σ x j (0), equations (5a)-(5d) do not depend on the particular choice of the initial state, but apply to all ρ 0 from the large class of initial states as given in (4). Moreover, the equations can easily be evaluated on a personal computer for millions of spins. The presence of a longitudinal magnetic field B = 0 modifies equations (5a)-(5d) by imposing an additional spin precession at an angular frequency proportional to B (see Appendix A). For the relaxation to equilibrium we are interested in, such an oscillatory spin precession is irrelevant and we will therefore restrict the discussion to the B = 0 case in the following. A generalization of the calculation in Appendix A to arbitrary n-spin correlation functions is also feasible by similar methods, but is not reported here. Since any operator on C 2 N can be expanded in terms of products of Pauli operators, such results for correlation functions of arbitrary order permit, in principle, the computation of the time evolution of expectation values of arbitrary operators. This does not add much to the discussion of the physical phenomena we are interested in here, and we will therefore restrict the discussion to two-spin correlations. As an example, the time evolution (5a)-(5d) of spin-spin correlation functions is illustrated for system parameters similar to those of the ion trap experiments of Britton et al. [15]. We consider hexagonal patches of triangular lattices (figure 1), and couplings J i,j = JD −α i,j proportional to the αth power of the inverse Euclidean distance of sites i and j on that lattice, with J ∈ R being a coupling constant. We have also performed numerical calculations for other two-dimensional lattice structures and geometries, and the behaviour we found is qualitatively the same and also quantitatively very similar.
The plots in figure 2 show the normalized expectation values σ , as given by (5a)-(5d). Normalized in this way, the plots are valid and numerically exact for all initial conditions of the form (4), and the same is true for figures 3, 4, and B1. The plotted curves suggest that correlations decay to their microcanonical equilibrium values σ a i σ b j mc = 0 where a, b ∈ {x, y, z}, but this decay is only apparent: Since the expressions in (5b)-(5d) consist of products of N trigonometric functions, the evolution in time of spin-spin correlations is quasi-periodic for any finite number N of spins. This is exactly the situation described in section 2 for general quantum systems on a finite-dimensional Hilbert space. Accordingly, recurrences will occur and repeatedly bring the system arbitrarily close to its initial state. These recurrences occur on a timescale that is exponentially large in N and, already for moderate system sizes, they do not show up on the timescales plotted in figure 2. (For example, for a lattice with side length L = 4 and α = 3/2, the first significant recurrence occurs at a time three orders of magnitude larger than the relaxation time). In the large-N limit, this recurrence time is expected to be pushed to infinity. We will confirm this expectation in section 5 by proving, in the thermodynamic limit, an upper bound on spin-spin correlation functions which converges to zero for large times t.

Upper bound on correlations in the thermodynamic limit.
With (5b)-(5d) as a starting point, it is possible to derive upper bounds on the time evolution of spin-spin correlations in the thermodynamic limit for different lattice structures (see Appendix B for the derivation in the case of a triangular lattice). All these bounds are of the form and The constants C ± i,j in (6b) depend on the lattice structure, the exponent α, and the position of i and j on the lattice. Other 2-spin correlation functions have similar properties, and details of the derivation are provided in Appendix B. All the bounds in the above equations decay to zero for large t. As can be read off from these equations, the slowest decay is always due to theP − i,j -term. For d 2, relaxation to zero is therefore governed by a term of the form exp −C − i,j (α)t d/(1+α) , implying a compressed exponential decay for α < d − 1, and stretched exponential decay otherwise.
For the example of a triangular lattice in two dimensions, we obtain With i, j as in figure 1 we have D i,j = 2, and setting α = 3/2 we obtain C + i,j ≈ 11.04 and C − i,j ≈ 1.119. We inserted these values into the bound (6a) and compared the outcome to exact finite-system results. Figure 3 (left) shows the correlation function σ x i σ x j and its upper bound (6a) vs. the rescaled time t 2/(1+α) . The rescaling of time is chosen such that the upper bound (6a) is a linearly decaying function. The plot reveals that, over a range of more than hundred orders of magnitude, the exact finite-system result for σ x i σ x j likewise decays with a linear trend, superimposed by fluctuations. This confirms that the upper bound (6a) indeed correctly captures the stretched or compressed exponential decay of the correlation function, albeit with prefactors C ± i,j that overestimate the actual behaviour.
The stretched or compressed exponential relaxation of the Ising model with power law interactions discussed above is different from the exponential decay of correlations known to occur in the nearest-neighbour Ising model [31]. In the limit α → ∞ of the power law interactions one would usually expect to approach the nearest-neighbour-interacting model and recover exponentially decaying correlations, but this does not seem to be the case. This puzzle is resolved by observing that, in the limit α → ∞, the coefficient C − i,j in (7b) also diverges, implying that the bound (6a) is not meaningful in this limit.
On the basis of the bound (6b) onP − i,j , one can define an upper bound on the relaxation time by requiring the exponent in (6b) to be, say, minus unity, Inserting (7b) and solving for t, an upper bound is obtained on the relaxation time τ for a triangular lattice in the thermodynamic limit, In figure 3 (right) this bound on the relaxation time is plotted vs. the exponent α.
Finite-system relaxation times as read off from figure 2 for α = 1/2 and α = 3/2 are also shown for comparison.

Prethermalization.
Apart from the overall relaxation times, there are other striking α-dependent aspects of the relaxation dynamics: The plot in the right panel of figure 2, representative for exponents α d/2, shows a simple relaxation to equilibrium with a single relevant timescale. The plot in the left panel of figure 2, representative for exponents 0 α d/2, differs considerably: In a first step, the correlation functions σ x i σ x j and σ y i σ y j display a Gaussian decay on a fast timescale τ 1 to one half of the σ x i σ x j initial value, before finally relaxing to their vanishing equilibrium value on a much longer timescale τ 2 . τ 1 is roughly the same as the timescale on which one-spin observables (green lines in figure 2) relax, while correlations are not yet equilibrated. Then the system remains prethermalized [32] for a relatively long period of time, a behaviour observed in the left panel of figure 2 as a conspicuous plateau. Different from earlier observations of prethermalization in quantum dynamics, the ratio τ 1 /τ 2 goes to zero in the large-N limit, i.e. separation of timescales become more pronounced with increasing system size. Such a behaviour of N-dependent relaxation timescales had previously been observed in classical long-range interacting systems, either for mean-field-type spin models [33], or for self-gravitating systems in the astrophysical context [34]. Besides the large-N limit, also the α → 0 limit leads to a more pronounced separation of timescales in the long-range Ising model. In the α → 0 limit, the slow relaxation timescale τ 2 diverges, while the fast timescale τ 1 remains finite for finite N. This fact can be explained by the presence of N(N − 1) additional symmetries in the α = 0 case, where all operators σ x i σ x j + σ y i σ y j commute with the Hamiltonian (2). Although these symmetries are absent for α 0, the slow timescale τ 2 can be seen as a remnant of the weakly destroyed α = 0 symmetry. While two-spin correlations exhibit relaxation in a two-step process for α < d/2, we have checked that no third timescale emerges when the calculations are generalized to three-spin correlations.
Multiple timescales of relaxation may emerge for various reasons. One possible scenario that has been observed previously, both in theory and experiment [35], is that dephasing is responsible for prethermalization, while a collisional mechanism causes relaxation on a slower timescale. In order to test which of these mechanism is at work in the long-range Ising model, we have computed the time evolution of the n-spin purity where ρ i 1 ,...,in = Tr Λ\{i 1 ,...,in} ρ is the n-spin reduced density operator, as obtained by tracing the density operator ρ over all sites of the lattice Λ except i 1 , . . . , i n . Knowledge of the one-spin expectation values s a i = σ a i with a ∈ {x, y, z} allows for the reconstruction of and additional knowledge of the spin-spin correlations s ab ij with a ∈ {x, y, z} facilitates the reconstruction of Inserting our results for the time evolution of spin expectation values of the long-range Ising model into these expressions, we readily obtain the time evolution of the oneand two-spin purities (10). As can be seen in figure 4 (left), both relaxation steps of spin-spin correlations we observed in figure 2 turn out to be associated with a drop in the purity γ ij , which indicates that both relaxation steps are caused by dephasing.
To demonstrate the dephasing yet more clearly, we also studied directly the modulus of the off-diagonal elements of the 2-spin density operator ρ ij in the σ z tensor product eigenbasis. For vanishing external magnetic field B = 0 and the class of initial conditions (4), we have s y i = s z i = s xy ij = s xz ij = s yx ij = s zx ij = s zz ij = 0 for all times t. With this simplification and the additional symmetries s a i = s a j and s ab ij = s ba ij , only four of the 16 matrix elements of ρ ij have different moduli, We have plotted these quantities in figure 4 (right) for several sizes of triangular lattices. All nondiagonal matrix elements decay to zero, but they do so on two different timescales: The matrix elements (ρ ij ) 23 and (ρ ij ) 32 (yellow in figure 4 (right)) are given by P − i,j /2 and hence decay on the longer of the two timescales we had previously observed. All other non-diagonal matrix elements decay on faster, N-dependent timescales. This demonstrates that both steps of the two-step relaxation process we observed for spinspin correlation functions in figure 2 are caused by dephasing. All the diagonal elements (ρ ij ) kk (black dashed line in figure 4 (right)) are constant in time, and therefore no redistribution of mode occupation numbers takes place. Such a redistribution would signal collisional relaxation due to inelastic scattering processes between these modes. Since the Ising model (2) has only commuting terms, no inelastic collisional mechanism is available.

Application to trapped-ion quantum simulators.
Long-range Ising interactions between effective two-level systems or spins can be implemented in crystalline arrays of trapped ions. Qualitatively, spin-dependent forces are used to modify the Coulomb interaction energy of the trapped ions and generate an effective Ising interaction. With small numbers of ions in linear radio-frequency traps, long-range Ising interactions like those in equation (1) are the basis of quantum gates, and have been implemented with high fidelity (>98%). These techniques have been recently extended to single-plane triangular lattices of several hundred ions stored in a Penning trap. To date, the Ising interactions implemented in the Penning ion trap simulator have been benchmarked for short timescales through comparison with mean field theory [15], a classical limit of the Hamiltonian (1) where quantum fluctuations and correlations are ignored. The exact results on correlation functions reported here will enable a much higher level of benchmarking, and a determination of the timescales over which the expected emulation of quantum effects of an Ising model is indeed realized.
Application of the results derived here for the ion trap simulators requires initial states that are diagonal in the σ x tensor product eigenbasis, and assurance that the trapped ion simulator is an isolated system for the quantum many-body equilibration time scales P ± i,j discussed in the previous sections. Preparation of diagonal states in the σ x tensor product eigenbasis is straight forward in ion traps. Specifically, optical pumping and coherent control techniques can initialize each spin to point along the x-axis with high fidelity. In the Penning ion trap simulator this initialization can be done with a fidelity of greater than 99% [36]. Equilibration time scales P ± i,j (see equation (5d)) that are short compared with other relaxation time scales have been demonstrated with 10-20 ions in linear radio-frequency traps [19,20,21,22]. In the current Penning trapped-ion simulator [15], spontaneous emission from the off-resonant laser beams used to engineer the long-range Ising interaction sets an experimental relaxation timescale that is comparable to the quantum many-body equilibration time scales. However, with straight-forward improvements to the setup described in [15], the quantum many-body relaxation timescales can be short compared to the spontaneous emission timescale. Specifically, with N = 217 ions and an optical dipole force generated from two 10 mW beams crossing with an angular separation of 35 o (see figure 1 of [15]) and frequency difference tuned to generate an α = 1/2 power law interaction, we calculate the two relaxation timescales P ± i,j in (5d) to be ≈ 30 µs and ≈ 430 µs. Here i and j are chosen to be two sites on opposite sides of the center site as shown in figure 1. The spontaneous emission time for this configuration is (1/Γ) ∼ 4 ms. All of these timescales are short compared to the T 2 50 ms coherence time of the Be + valence electron spin qubits.
To compare with the exact results reported here requires an experimental measurement of the spin-spin correlations after the Ising interaction has been applied for a variable period. This is readily accomplished with trapped ions by using spin-dependent resonance fluorescence. Britton et al [15] use spin-dependent resonance fluorescence to measure spin orientation in the σ z basis. With resolved imaging of the ions array (see experimental image in figure 1), the spin orientation of each ion can be detected and arbitrary pair correlation functions calculated. However, more simply, the global fluorescence collected from all the ions in the array can be used to measure the global spin state of the system (i.e., the total number of spins in the |↑ z state and the total number of spins in the |↓ z state). Shot-to-shot fluctuations in these measurements are sensitive to the second-order moment J 2 z of the total z-component of the spin Of particular interest are measurements of the second-order moments in directions perpendicular to the mean composite spin vector. For example, with all spins initially pointing along the x-axis, fluctuations in measurements of J z after rotation about the x-axis by an angle θ are sensitive to These so-called spin-squeezing measurements [37] are sensitive to pairwise correlations summed over all the spins in the ensemble.

Outlook and summary.
Once such a quantum simulator is benchmarked, it may in turn be used for the testing of approximate calculations of the quantum dynamics of the long-range Ising model in the presence of non-longitudinal fields (resulting in non-commuting terms in the Hamiltonian (1)). Moreover, there are certain phenomena which are believed to be peculiar to long-range interacting systems and which might receive their first experimental confirmation in the ion trap setup. One example is the threshold at α = d/2 below which a second timescale of relaxation emerges. There is evidence that not only the existence of this dynamical long-range threshold, but also its numerical value d/2, is universal for classical as well as quantum mechanical long-range interacting systems [38]. Other long-range peculiarities include nonequivalent equilibrium statistical ensembles [39,40], a phenomenon that has been known in the astrophysical context for decades [41] but has not seen experimental verification in a tabletop experiment. In summary, under convenient restrictions on the initial conditions, we have obtained exact analytic results for spin-spin correlations of long-range quantum Ising models in arbitrary spatial dimensions and for various lattice structures. Prethermalization, widely separated timescales, and other phenomena that further our understanding of the timescales governing the relaxation to equilibrium were discussed. Comparison of our analytic results to experimental measurements of spin-spin correlations will enable benchmarking of the trapped-ion quantum simulator of [15] in a regime where quantum effects are paramount. Subsequently, the benchmarked quantum simulator may be employed to investigate relaxation timescales in parameter regimes and for initial states where analytic results presently are not available.
Appendix A. Exact calculation of spin-spin correlations for finite N.
Starting point of the calculation is the general expression of the expectation value with respect to the initial density operator ρ 0 , where we assume ρ 0 to be diagonal in the σ x tensor product eigenbasis. Spin ladder operators are defined as σ ± = (σ x ± iσ y )/2. Since all terms in the Hamiltonian H ℓ commute, we can factorize the time evolution operator, and similarly for the Hermitian conjugate. All factors in (A.2) that contain neither σ z i nor σ z j commute with σ ± i σ ± j . To compute the time evolution in (A.1), we therefore have to deal with the expression Making use of [σ z , σ ± ] = ±2σ ± , the time evolution due to the magnetic field B simplifies to . (A.4) Picking one lattice site k = i, j, the time evolution of σ ± i σ ± j due to the interaction with the spin at k can be written as Since the initial state ρ 0 is assumed to be diagonal in the σ x tensor product eigenbasis, only diagonal elements of the operator e iH ℓ t σ ± i σ ± j e −iH ℓ t in the same basis contribute to the trace in (A.1). For this reason we can drop the second term on the right-hand side of (A.5), as it is proportional to σ z k . Inserting (A.4) and (A.5) into (A.1), we obtain where σ ± i σ ± j (0) = Tr σ ± i σ ± j ρ 0 . A similar calculation yields a notable difference to (A.6) being that the presence of an external magnetic field B has no effect on this expectation value. From we obtain 13) and furthermore, using the properties of the initial density operator ρ 0 , Inserting (A.6), (A.7) and (A.14) into (A.10)-(A.13), we arrive at the final expressions with P ± i,j as defined in (5d). A similar calculation yields Finally, from symmetry considerations we obtain A generalization of these calculations to n-spin correlation functions is straightforward, yielding for example and similar expressions for other correlation functions.
Appendix B. Upper bounds on correlations in the thermodynamic limit.
In the thermodynamic limit of infinite system size, a bound on the product in equation (5d) with J i,j = D −α i,j and α 0 is derived. For any given t, one can find a compact region g ± i,j (t) (containing i and j) of the infinite triangular lattice G such that |2(J k,i ± J k,j )t| < π/2 (B.2) for all k ∈ G \ g ± i,j (t). Since | cos x| 1, we have Next, using the inequality we obtain For simplicity, we consider lattice sites i and j symmetrically arranged to the right and left of lattice site 0, as illustrated below for the example of distance δ = D i,j = 2, but other arrangements can be treated in a similar way.
K K (4) Disregarding all sites that lie outside the grey shaded area h i,j in the above illustration, we obtain a lower bound on the sum in (B.5), Each individual summand (J k,i ± J k,j ) 2 can be minorized by (J K,i ± J K,j ) 2 , where K is determined from k by going vertically upwards along the dashed (red) line in the above illustration until a boundary point of the shaded region is reached. To a boundary point labelled by K(r) there correspond r + 1 points in the sum, and this allows us to bound the sum by (B.7) To exclude lattice sites in the region g ± i,j (t), i.e. in order to satisfy the inequality (B.2) used earlier, R ± 0 (t) has to be chosen large enough. An asymptotic analysis of (B.2) shows that this condition can be implemented by choosing where the condition on R + 0 (t) works always, the one for R − 0 (t) only for sufficiently large t. Inserting the distances D i,K(r) = r, D j,K(r) = √ δ 2 + rδ + r 2 (B.9) and bounding the sum by an integral, we obtain i, j Figure B1. Logarithmic plot of the products P + i,j (top) and P − i,j (bottom) as a function of rescaled time t p ± (α) , evaluated for a hexagonal patch of a triangular lattice with side length L = 16 for α = 3/4 (left) and α = 5/2 (right). The exponents p + (α) = min{2, 2/α} and p − (α) = 2/(1 + α) of the time rescaling are those predicted by the bounds (B.13) and (B.14) for the asymptotic behaviour at large L and t. The linearly decaying trend in the plot, superimposed by fluctuations, confirms over a range of more than hundred orders of magnitude that the asymptotic behaviour of P ± i,j is indeed as given by the bounds. The straight lines in the plots are the bounds (B.13) and (B.14), indicating that the numerical constants in the exponents of (B.13) and (B.14) are, as expected, underestimated.
valid in the limit of large R. Since we are interested in the limit of large R and t, we can expand the integrand to leading order in 1/r, obtaining Inserting these expressions into (B.10), the integral can be solved by elementary means. Making use of the asymptotic equality R 2 ∼ N/3 and performing the limit N → ∞, we obtain the final results , (B.14) valid in the limit of large R and t. A comparison of these bounds with an exact evaluation of (B.1) for finite lattices is shown, over more than hundred orders of magnitude, in figure B1.