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

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–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–8], rigorous proofs of equilibration have been obtained for generic Hamiltonians [9–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 long-range 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 ⁹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

$$\begin{equation} H=-\sum_{i < j}J_{i,j}\sigma_i^z\sigma_j^z-{\boldsymbol{B}}_\mu\sum_{i}{\boldsymbol{\sigma}}_i. \end{equation} \tag{ 1 }$$

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 ([19, equation (4)]). 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–22] and, along similar lines, benchmarking of quantum magnets emulated by means of ultracold molecules is possible [23].

Zoom In Zoom Out Reset image size

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 this 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 time-evolved 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 [9–14, 24, 25]. 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 this 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.

For the theoretical analysis of relaxation times and the decay of correlations, we study a model that is more general than the Hamiltonian (1) in some respects, and more restricted in others. We generalize 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

$$\begin{equation} H_\ell=-\sum_{i < j}J_{i,j}\sigma_i^z\sigma_j^z-B\sum_{i}\sigma_i^z, \end{equation} \tag{ 2 }$$

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 non-negative 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 of whether the corresponding energies are situated in the low- or high-energy regime.

To study the relaxation to equilibrium in this model, our aim is to compute the time evolution of spin–spin correlation functions

$$\begin{equation} \langle\sigma_i^x\sigma_j^x\rangle(t)={\rm Tr} \left({\rm e}^{{\rm i} H_\ell t} \sigma_i^x\sigma_j^x \,{\rm e}^{-{\rm i} H_\ell t}\rho_0\right). \end{equation} \tag{ 3 }$$

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 ρ₀ that are diagonal matrices in the σ^x tensor-product eigenbasis,

$$\begin{equation} \rho_0=\frac{1}{2^N}\left({\mathbbm{1}}+\sum_{i}\sigma_i^x\left(s_i^x + \sum_{j>i}\sigma_j^x\left(s_{ij}^{xx} + \sum_{k>j}\sigma_k^x\left(s_{ijk}^{xxx} + \sum_{l>k}\cdots\right)\right)\right)\right), \end{equation} \tag{ 4 }$$

where ${\mathbbm {1}}$ denotes the identity operator on the tensor product Hilbert space $\mathcal {H}={\mathbb {C}}^2\otimes \cdots \otimes {\mathbb {C}}^2$. The indices i, j, k, l in (4) are summed over the lattice sites. The set of all ρ₀ as defined in (4) is a very large class of initial states, parameterized by 2^N − 1 real continuous parameters $s_i^x=\langle \sigma _i^x \rangle (0)$, $s_{ij}^{xx}=\langle \sigma _i^x\sigma _j^x\rangle (0),\ldots \;$ The reader may wish to compare the size of this class with the much studied quantum quenches, a rather restrictive class of initial states parameterized 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 ρ₀. This extension essentially leaves the results of this 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–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

$$\begin{equation} \langle\sigma_i^x\sigma_j^y\rangle (t)=\langle\sigma_i^x\sigma_j^z\rangle(t)=\left\langle\sigma_i^z\sigma_j^z\right\rangle (t)=0, \end{equation} \tag{ 5a }$$

$$\begin{equation} \langle\sigma_i^y\sigma_j^z\rangle (t)=\langle\sigma_i^x\rangle (0)\sin(2tJ_{i,j})\prod_{k\neq i,j}\cos(2tJ_{i,k}),\end{equation} \tag{ 5b }$$

$$\begin{equation} \langle\sigma_i^x\sigma_j^x\rangle (t)=P_{i,j}^- + P_{i,j}^+,\quad\langle\sigma_i^y\sigma_j^y\rangle (t)=P_{i,j}^- - P_{i,j}^+ \end{equation} \tag{ 5c }$$

$$\begin{equation} P_{i,j}^\pm=\frac{1}{2}\langle\sigma_i^x\sigma_j^x\rangle (0)\prod_{k\neq i,j}\cos\left[2\left(J_{k,i}\pm J_{k,j}\right)t\right]. \end{equation} \tag{ 5d }$$

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 ρ₀ 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 ${\mathbb {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{\in }{\mathbb {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 〈σ^x_iσ^x_j〉(t)/〈σ^x_iσ^x_j〉(0), 〈σ^y_iσ^y_j〉(t)/〈σ^x_iσ^x_j〉(0), 〈σ^y_iσ^z_j〉(t)/〈σ^x_i〉(0), 〈σ^x_i〉(t)/〈σ^x_i〉(0), 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 B.1. 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

$$\begin{equation} \langle\sigma_i^x\sigma_j^x\rangle (t),\quad \left\langle\sigma_i^y\sigma_j^y\right\rangle (t)\leqslant\bar{P}_{i,j}^- + \bar{P}_{i,j}^+ \end{equation} \tag{ 6a }$$

with

$$\begin{equation} \bar{P}_{i,j}^\pm=\frac{1}{2}\langle\sigma_i^x\sigma_j^x\rangle (0)\exp[-C_{i,j}^\pm(\alpha)N^{q^\pm(\alpha)}t^{p^\pm(\alpha)}] \end{equation} \tag{ 6b }$$

and

$$\begin{equation} q^+(\alpha)=\max\{0,1-2\alpha/d\}, \end{equation} \tag{ 6c }$$

$$\begin{equation} q^-(\alpha)=\max\{0,1-2(1+\alpha)/d\},\end{equation} \tag{ 6d }$$

$$\begin{equation} p^+(\alpha)=\min\{2,d/\alpha\},\end{equation} \tag{ 6e }$$

$$\begin{equation} p^-(\alpha)=\min\{2,d/(1+\alpha)\}. \end{equation} \tag{ 6f }$$

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 two-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 the $\bar {P}_{i,j}^-$-term. For d ⩽ 2, relaxation to zero is therefore governed by a term of the form $\exp \!\left [-C_{i,j}^-(\alpha )t^{d/(1+\alpha )}\right ]$, 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

$$\begin{eqnarray} C_{i,j}^+(\alpha) =\left\{\begin{array}{@{}ll} \displaystyle\frac{(8J)^2\,3^{\alpha-1}}{\pi^2(1-\alpha)} & {\rm for} \ \alpha<1,\\ \displaystyle\frac{1}{\alpha-1}\left(\frac{8J}{\pi}\right)^{2/\alpha} & {\rm for} \ \alpha>1,\end{array}\right.\end{eqnarray} \tag{ 7a }$$

$$\begin{equation} C_{i,j}^-(\alpha) = \frac{D_{i,j}^2}{4\alpha}\left(\frac{4\alpha J}{\pi}\right)^{2/(1+\alpha)}. \end{equation} \tag{ 7b }$$

$\langle \sigma _i^x\sigma _j^x \rangle$

$\langle \sigma _i^x\sigma _j^x \rangle$

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) on $\bar {P}_{i,j}^-$, one can define an upper bound on the relaxation time by requiring the exponent in (6b) to be, say, minus unity,

$$\begin{equation} C_{i,j}^-(\alpha)t^{d/(1+\alpha)}=1. \end{equation} \tag{ 8 }$$

Inserting (7b) and solving for t, an upper bound is obtained on the relaxation time τ for a triangular lattice in the thermodynamic limit

$$\begin{equation} \tau(\alpha)\leqslant\left(\frac{2\sqrt{\alpha}}{D_{i,j}}\right)^{1+\alpha}\frac{\pi}{4\alpha J}. \end{equation} \tag{ 9 }$$

In figure 3 (right) this bound on the relaxation time is plotted versus the exponent α. Finite-system relaxation times as read off from figure 2 for α = 1/2 and 3/2 are also shown for comparison.

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 $\langle \sigma _i^x\sigma _j^x \rangle$ and $\langle \sigma _i^y\sigma _j^y \rangle$ display a Gaussian decay on a fast timescale τ₁ to one half of the $\langle \sigma _i^x\sigma _j^x\rangle$ initial value, before finally relaxing to their vanishing equilibrium value on a much longer timescale τ₂. τ₁ 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 τ₁/τ₂ 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 τ₂ diverges, whereas the fast timescale τ₁ 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 τ₂ 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

$$\begin{equation} \gamma_{i_1,\ldots,i_n}(t)={\rm Tr}(\rho_{i_1,\ldots,i_n}^2(t)), \end{equation} \tag{ 10 }$$

where

$$\begin{equation} \rho_{i_1,\ldots,i_n}={\rm Tr}_{\Lambda\setminus\{i_1,\ldots,i_n\}}\rho \end{equation} \tag{ 11 }$$

is the n-spin reduced density operator, as obtained by tracing the density operator ρ over all sites of the lattice Λ except $i_1,\ldots ,i_n$. Knowledge of the one-spin expectation values s^a_i = 〈σ^a_i〉 with a∈{x,y,z} allows for the reconstruction of

$$\begin{equation} \rho_i=\frac{1}{2}\left(\!\!\begin{array}{cc}1+s_i^z & s_i^x-{\rm i} s_i^y\\ s_i^x+{\rm i} s_i^y & 1-s_i^z\end{array}\!\!\right) \end{equation} \tag{ 12 }$$

and additional knowledge of the spin–spin correlations s^ab_ij with a∈{x,y,z} facilitates the reconstruction of

$$\begin{eqnarray} \fl \rho_{i,j}=\frac{1}{4}\!\left(\!\!{\begin{array}{cccc} 1+s_i^z+s_j^z+s_{ij}^{zz} & s_j^x-{\rm i} s_j^y+s_{ij}^{zx}-{\rm i} s_{ij}^{zy} & s_i^x-{\rm i} s_i^y+s_{ij}^{xz}-{\rm i} s_{ij}^{yz} & s_{ij}^{xx}-s_{ij}^{yy}-{\rm i} s_{ij}^{xy}-{\rm i} s_{ij}^{yx}\\ & 1+s_i^z-s_j^z-s_{ij}^{zz} & s_{ij}^{xx}+s_{ij}^{yy}+{\rm i} s_{ij}^{xy}-{\rm i} s_{ij}^{yx} & s_i^x-{\rm i} s_i^y-s_{ij}^{xz}+{\rm i} s_{ij}^{yz}\\ & & 1-s_i^z+s_j^z-s_{ij}^{zz} & s_j^x-{\rm i} s_j^y-s_{ij}^{zx}+{\rm i} s_{ij}^{zy}\\ {\rm h.\, c.} & & & 1-s_i^z-s_j^z+s_{ij}^{zz} \end{array} }\!\!\right)\!. \end{eqnarray} \tag{ 13 }$$

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 one- and 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 two-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:

$$\begin{equation} |(\rho_{ij})_{kk}| = 1/4 \quad\forall k\in\{1,2,3,4\}, \end{equation} \tag{ 14a }$$

$$\begin{equation}|(\rho_{ij})_{14}| = |(\rho_{ij})_{41}| = s_{ij}^{xx}-s_{ij}^{yy}|/4=P_{i,j}^+/2, \end{equation} \tag{ 14b }$$

$$\begin{equation} |(\rho_{ij})_{23}| =|(\rho_{ij})_{32}| = |s_{ij}^{xx}+s_{ij}^{yy}|/4=P_{i,j}^-/2, \end{equation} \tag{ 14c }$$

$$\begin{eqnarray} |(\rho_{ij})_{12}| & = & |\left(\rho_{ij}\right)_{13}| = |(\rho_{ij})_{21}| = |(\rho_{ij})_{24}| = |(\rho_{ij})_{31}| = |(\rho_{ij})_{34}|\nonumber\\ & = & |(\rho_{ij})_{42}| = |(\rho_{ij})_{43}| = |s_i^x+{\rm i} s_{ij}^{yz}|/4. \end{eqnarray} \tag{ 14d }$$

$(\rho _{ij})_{23}$

$\left (\rho _{ij}\right )_{32}$

$(\rho _{ij})_{kk}$

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 such as 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 timescales P^±_i,j discussed in the previous sections. Preparation of diagonal states in the σ^x tensor product eigenbasis is straightforward 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 timescales P^±_i,j (see equation (5d)) that are short compared with other relaxation timescales have been demonstrated with 10–20 ions in linear radio-frequency traps [19–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 timescales. However, with straightforward 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° (see  [15, figure 1]) 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 centre 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₂ ≳ 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 $\left |\uparrow \right \rangle _{z}$ state and the total number of spins in the $\left |\downarrow \right \rangle _{z}$ state). Shot-to-shot fluctuations in these measurements are sensitive to the second-order moment $\left \langle J_{z}^{2}\right \rangle$ of the total z-component of the spin $J_{z}\equiv \sum _{i=1}^{N}\sigma _{i}^{z}/2$. 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

$$\begin{equation}\fl \langle J_{z}^{2}\rangle _{\theta}=\sum_{i,j=1}^{N}\frac{1}{4}[\sin^{2}\theta\langle \sigma_{i}^{y}\sigma_{j}^{y}\rangle -\sin\theta\cos\theta\langle \sigma_{i}^{y}\sigma_{j}^{z}+\sigma_{i}^{z}\sigma_{j}^{y}\rangle +\cos^{2}\theta\langle \sigma_{i}^{z}\sigma_{j}^{z}\rangle]. \end{equation} \tag{ 15 }$$

These so-called spin-squeezing measurements [37] are sensitive to pairwise correlations summed over all the spins in the ensemble.

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 non-equivalent 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 are not presently available.

BCS and JJB acknowledge support from the DARPA OLE program. MK acknowledges support by the Incentive Funding for Rated Researchers programme of the National Research Foundation of South Africa. Papers with contributions from the US National Institute of Standards and Technology are not subject to US copyright. During the preparation of this paper we were informed of related work on spin–spin correlation functions of long-range Ising models [42]. In contrast to our work, which focuses on large time asymptotic behaviour and relaxation to equilibrium, Foss-Feig et al [42] investigate the effect of decoherence on the spin–spin correlation functions.

Starting point of the calculation is the general expression of the expectation value

$$\begin{equation} \langle\sigma_i^\pm\sigma_j^\pm\rangle (t)={\rm Tr} \left({\rm e}^{{\rm i} H_\ell t} \sigma_i^\pm\sigma_j^\pm\, {\rm e}^{-{\rm i} H_\ell t}\rho_0\right) \end{equation} \tag{ A.1 }$$

with respect to the initial density operator ρ₀, where we assume ρ₀ 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

$$\begin{equation} \exp\!\left(-{\rm i} H_\ell t\right)=\prod_{k < l}\exp\!\left({\rm i} J_{k,l}\sigma_k^z\sigma_l^z t\right)\prod_m\exp({\rm i} B\sigma_m^z t) \end{equation} \tag{ A.2 }$$

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

$$\begin{eqnarray}&&\fl \prod_{k\neq i,j}\exp(-{\rm i} J_{k,i}\sigma_k^z\sigma_i^z t)\exp\!\left(-{\rm i} J_{k,j}\sigma_k^z\sigma_j^z t\right)\exp\!\left(-{\rm i} B\sigma_i^z t\right)\exp\!\left(-{\rm i} B\sigma_j^z t\right)\sigma_i^\pm\sigma_j^\pm\nonumber\\ &&\times\exp({\rm i} B\sigma_j^z t)\exp\!\left({\rm i} B\sigma_i^z t\right)\prod_{l\neq i,j}\exp\!\left({\rm i} J_{l,j}\sigma_l^z\sigma_j^z t\right)\exp\!\left({\rm i} J_{l,i}\sigma_l^z\sigma_i^z t\right). \end{eqnarray} \tag{ A.3 }$$

Making use of $[\sigma ^z,\sigma ^\pm ]=\pm 2\sigma ^\pm$, the time evolution due to the magnetic field B simplifies to

$$\begin{equation} \fl \exp(-{\rm i} B\sigma_i^z t)\exp\!\left(-{\rm i} B\sigma_j^z t\right)\sigma_i^\pm\sigma_j^\pm\exp\!\left({\rm i} B\sigma_j^z t\right)\exp\!\left({\rm i} B\sigma_i^z t\right)=\sigma_i^\pm\sigma_j^\pm\exp\!\left(\mp4{\rm i} Bt\right). \end{equation} \tag{ 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

$$\begin{eqnarray}&&\fl \exp\!\left(-{\rm i} J_{k,i}\sigma_k^z\sigma_i^z t\right)\exp\!\left(-{\rm i} J_{k,j}\sigma_k^z\sigma_j^z t\right)\sigma_i^\pm\sigma_j^\pm \exp\!\left({\rm i} J_{k,j}\sigma_k^z\sigma_j^z t\right)\exp\!\left({\rm i} J_{k,i}\sigma_k^z\sigma_i^z t\right)\nonumber\\ &&=\sigma_i^\pm\sigma_j^\pm\cos\left[2t(J_{i,k}+J_{j,k})\right]\mp {\rm i}\sigma_i^\pm\sigma_j^\pm\sigma_k^z\sin\left[2t(J_{i,k}+J_{j,k})\right]. \end{eqnarray} \tag{ A.5 }$$

Since the initial state ρ₀ is assumed to be diagonal in the σ_x tensor product eigenbasis, only diagonal elements of the operator e^iH_ℓtσ^±_iσ^±_je^−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

$$\begin{equation} \langle\sigma_i^\pm\sigma_j^\pm\rangle (t)=\langle\sigma_i^\pm\sigma_j^\pm\rangle (0)\exp\!\left(\mp4{\rm i} Bt\right) \prod_{k\neq i,j}\cos\left[2t(J_{i,k}+J_{j,k})\right], \end{equation} \tag{ A.6 }$$

where $\langle \sigma _i^\pm \sigma _j^\pm \rangle (0)={\rm Tr}\left (\sigma _i^\pm \sigma _j^\pm \rho _0\right )$. A similar calculation yields

$$\begin{equation} \langle\sigma_i^\pm\sigma_j^\mp\rangle (t)=\left\langle\sigma_i^\pm\sigma_j^\mp\right\rangle (0)\prod_{k\neq i,j}\cos\left[2t(J_{i,k}- J_{j,k})\right] \end{equation} \tag{ A.7 }$$

a notable difference to (A.6) being that the presence of an external magnetic field B has no effect on this expectation value. From

$$\begin{equation} 4\sigma_i^\pm\sigma_j^\pm=\sigma_i^x\sigma_j^x-\sigma_i^y\sigma_j^y\pm{\rm i}\sigma_i^x\sigma_j^y\pm{\rm i}\sigma_i^y\sigma_j^x, \end{equation} \tag{ A.8 }$$

$$\begin{equation} 4\sigma_i^\pm\sigma_j^\mp=\sigma_i^x\sigma_j^x+\sigma_i^y\sigma_j^y\mp{\rm i}\sigma_i^x\sigma_j^y\pm{\rm i}\sigma_i^y\sigma_j^x \end{equation} \tag{ A.9 }$$

we obtain

$$\begin{equation} \sigma_i^x\sigma_j^x=\phantom{\mp}2\Re\left(\sigma_i^\pm\sigma_j^\mp+\sigma_i^\pm\sigma_j^\pm\right),\end{equation} \tag{ A.10 }$$

$$\begin{equation} \sigma_i^y\sigma_j^y=\phantom{\mp}2\Re\left(\sigma_i^\pm\sigma_j^\mp-\sigma_i^\pm\sigma_j^\pm\right),\end{equation} \tag{ A.11 }$$

$$\begin{equation} \sigma_i^x\sigma_j^y=\mp2\Im\left(\sigma_i^\pm\sigma_j^\mp-\sigma_i^\pm\sigma_j^\pm\right),\end{equation} \tag{ A.12 }$$

$$\begin{equation} \sigma_i^y\sigma_j^x=\pm2\Im\left(\sigma_i^\pm\sigma_j^\mp+\sigma_i^\pm\sigma_j^\pm\right) \end{equation} \tag{ A.13 }$$

and furthermore, using the properties of the initial density operator ρ₀,

$$\begin{equation} 4\langle\sigma_i^\pm\sigma_j^\pm\rangle (0) = 4\langle\sigma_i^\pm\sigma_j^\mp\rangle (0) = \langle\sigma_i^x\sigma_j^x\rangle (0). \end{equation} \tag{ A.14 }$$

Inserting (A.6), (A.7) and (A.14) into (A.10)–(A.13), we arrive at the final expressions

$$\begin{equation} \left\langle\sigma_i^x\sigma_j^x\right\rangle (t)=P_{i,j}^- + \cos(4Bt)P_{i,j}^+,\end{equation} \tag{ A.15 }$$

$$\begin{equation} \langle\sigma_i^y\sigma_j^y\rangle (t)=P_{i,j}^- - \cos(4Bt)P_{i,j}^+,\end{equation} \tag{ A.16 }$$

$$\begin{equation} \langle\sigma_i^x\sigma_j^y\rangle (t)= -\sin(4Bt)P_{i,j}^+ \end{equation} \tag{ A.17 }$$

with P^±_i,j as defined in (5d). A similar calculation yields

$$\begin{equation} \langle\sigma_i^x\sigma_j^z\rangle (t)=\sin(2Bt)P_{i,j}^z,\quad \langle\sigma_i^y\sigma_j^z\rangle (t)=\cos(2Bt)P_{i,j}^z \end{equation} \tag{ A.18 }$$

with

$$\begin{equation} P_{i,j}^z=-\langle\sigma_i^x\rangle (0)\sin(2tJ_{i,j})\prod_{k\neq i,j}\cos(2tJ_{i,k}). \end{equation} \tag{ A.19 }$$

Finally, from symmetry considerations we obtain

$$\begin{equation} \langle\sigma_i^z\sigma_j^z\rangle (t)=0. \end{equation} \tag{ A.20 }$$

A generalization of these calculations to n-spin correlation functions is straightforward, yielding for example

$$\begin{equation} \fl \langle\sigma_i^\pm\sigma_j^\pm\sigma_k^\pm \rangle (t) =\langle\sigma_i^\pm\sigma_j^\pm\sigma_k^\pm\rangle (0)\exp\!\left(\mp6{\rm i} Bt\right) \prod_{l\neq i,j,k}\cos\left[2t(J_{i,l}+J_{j,l}+J_{k,l})\right] \end{equation} \tag{ A.21 }$$

and similar expressions for other correlation functions.

In the thermodynamic limit of infinite system size, a bound on the product

$$\begin{equation} \mathcal{P}_{i,j}^\pm=\prod_{\scriptstyle k=1\atop \scriptstyle k\neq i,j}^N\left|\cos\left[2(J_{k,i}\pm J_{k,j})t\right]\right| \end{equation} \tag{ B.1 }$$

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

$$\begin{equation} |2(J_{k,i}\pm J_{k,j})t|<\pi/2 \end{equation} \tag{ B.2 }$$

for all k∈G\g^±_i,j(t). Since |cos x| ⩽ 1, we have

$$\begin{equation} \mathcal{P}_{i,j}^\pm\leqslant\prod_{k\notin g_{i,j}^\pm(t)}\left|\cos\left[2(J_{k,i}\pm J_{k,j})t\right]\right|. \end{equation} \tag{ B.3 }$$

Next, using the inequality

$$\begin{equation} |\cos x|\leqslant1-\left(\frac{2x}{\pi}\right)^2\leqslant\exp\!\left[-\left(\frac{2x}{\pi}\right)^2\right]\quad\mbox{for} \;x < \pi/2 \end{equation} \tag{ B.4 }$$

we obtain

$$\begin{equation} \mathcal{P}_{i,j}^\pm\leqslant\exp\!\left[-\frac{16t^2}{\pi^2}\sum_{k\notin g_{i,j}^\pm(t)}\left(J_{k,i}\pm J_{k,j}\right)^2\right]. \end{equation} \tag{ B.5 }$$

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.

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)

$$\begin{equation} \sum_{k\in G\setminus g_{i,j}^\pm(t)}(J_{k,i}\pm J_{k,j})^2{\geqslant}\sum_{k\in h_{i,j} \cap (G\setminus g^{\pm}_{i,j}(t))}(J_{k,i}\pm J_{k,j})^2. \end{equation} \tag{ B.6 }$$

Each individual summand $(J_{k,i}\pm J_{k,j})^2$ can be minorized by $(J_{K,i}\pm 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

$$\begin{equation} \sum_{k\in G\setminus g_{i,j}^\pm(t)}(J_{k,i}\pm J_{k,j})^2\geqslant2J^2\sum_{r=R_0^\pm(t)}^R r(D_{i,K(r)}^{-\alpha}\pm D_{j,K(r)}^{-\alpha})^2. \end{equation} \tag{ 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^±₀(t) has to be chosen large enough. An asymptotic analysis of (B.2) shows that this condition can be implemented by choosing

$$\begin{equation} R_0^+(t)>\left(\frac{8Jt}{\pi}\right)^{1/\alpha},\quad R_0^-(t)\sim\left(\frac{4\alpha Jt}{\pi}\right)^{1/(1+\alpha)}, \end{equation} \tag{ B.8 }$$

where the condition on R⁺₀(t) works always, the one for R⁻₀(t) only for sufficiently large t. Inserting the distances

$$\begin{equation} D_{i,K(r)}=r,\quad D_{j,K(r)}=\sqrt{\delta^2+r\delta+r^2} \end{equation} \tag{ B.9 }$$

and bounding the sum by an integral, we obtain

$$\begin{equation} \mathcal{P}_{i,j}^\pm\leqslant\exp\left[-\frac{32J^2t^2}{\pi^2}\int\nolimits_{R_0^\pm(t)}^R{\rm d} r\, r \left(\frac{1}{r^\alpha}\pm\frac{1}{\sqrt{\delta^2+r\delta+r^2}^\alpha}\right)^2\right] \end{equation} \tag{ B.10 }$$

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

$$\begin{equation} r \left(\frac{1}{r^\alpha}+\frac{1}{\sqrt{\delta^2+r\delta+r^2}^\alpha}\right)^2 \sim 4r^{1-2\alpha}, \end{equation} \tag{ B.11 }$$

$$\begin{equation} r \left(\frac{1}{r^\alpha}-\frac{1}{\sqrt{\delta^2+r\delta+r^2}^\alpha}\right)^2 \sim \left(\frac{\alpha\delta}{2}\right)^2 r^{-1-2\alpha}. \end{equation} \tag{ B.12 }$$

Inserting these expressions into (B.10), the integral can be solved by elementary means. Making use of the asymptotic equality R² ∼ N/3 and performing the limit N → ∞, we obtain the final results

$$\begin{equation} \mathcal{P}_{i,j}^+\leqslant \left\{\begin{array}{@{}ll} \exp\!\left[-\frac{64J^2t^2}{\pi^2(1-\alpha)}\left(\frac{N}{3}\right)^{1-\alpha}\right] {\rm for} \ \alpha < 1,\\ \displaystyle\exp\!\left[-\frac{1}{\alpha-1}\left(\frac{8Jt}{\pi}\right)^{2/\alpha}\right] {\rm for} \ \alpha>1\end{array}\right. \end{equation} \tag{ B.13 }$$

and

$$\begin{equation} \mathcal{P}_{i,j}^-\leqslant\exp\!\left[-\frac{\delta^2}{4\alpha}\left(\frac{4\alpha Jt}{\pi}\right)^{2/(1+\alpha)}\right] \end{equation} \tag{ 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 B.1.

Zoom In Zoom Out Reset image size