Prethermalization and dynamic phase transition in an isolated trapped ion spin chain

The dynamic properties of isolated quantum many-body systems have garnered intense interest in recent years [1, 2]. On the theory side, research has been centered on whether and how an isolated quantum system approaches thermal equilibrium. While certain observables are found to relax to equilibrium in some large systems [3–7], it remains unclear under what conditions and on which timescale equilibration occurs in generic systems [8–11]. On the experimental side, recent progress with cold atoms [12–14] and trapped ions [15–20] has made it possible to simulate well controlled simple models, such as the one-dimensional (1D) Bose gas and transverse field Ising model. These quantum systems can be well isolated from the environment and have long coherence times, and their physical properties can be measured at the single-atom level, providing an unprecedented opportunity for studying non-equilibrium dynamics in closed interacting systems.

A particularly intriguing phenomenon in this context is prethermalization [21], which has been shown to emerge in various theoretical setups [22–24], and has been experimentally observed in ultracold atomic gas [14]. Prethermalization is characterized by the establishment of a quasi-stationary state at an intermediate timescale, followed by relaxation to a stationary state at a much longer timescale (thermalization). The physical origin of prethermalization, however, is still elusive, and is speculated to be related to the quasi-integrability of the model [14, 23].

In this paper, we propose a new experimental scheme for observing and studying prethermalization and a related dynamic phase transition in an XY spin model, which can be implemented with the current trapped-ion quantum simulator of [19]. Our model features long-range spin–spin interaction and an inhomogeneous lattice spacing, and unlike many other systems, the prethermalization can occur already for as few as a dozen spins, allowing for its observation in current experimental systems. The prethermalization found in this system results from a non-trivial structure in the energy spectrum, which in turn arises due to a combination of long-range interaction and inhomogeneous lattice spacing. In addition, by tuning the range of interactions with an experimental knob, we find that the dynamic behavior of this system exhibits three different phases: thermalization only, prethermalization followed by thermalization and prethermalization only. The transition between these different phases becomes sharper with an increasing number of spins, hinting at the existence of a dynamic phase transition [25] in the thermodynamic limit.

Our spin model is based on the experimental system of a chain of ions confined in a linear Paul trap (figure 1). Through proper configuration of the Raman beams, the optical dipole force can generate an effective transverse field Ising model [16, 17, 26]

$$\begin{equation} H=\sum_{i<j}^{N}J_{i,j}\sigma_{i}^{x}\sigma_{j}^{x}+B\sum_{i=1}^{N}\sigma_{i}^{z}, \end{equation} \tag{ 1 }$$

where σ_i is the spin-1/2 Pauli matrix for the ith ion qubit. The interaction coefficients J_ij in equation (1) are given by

$$\begin{eqnarray*} &&J_{i,j}=\Omega^{2}\sum_{m=1}^{N}\frac{\eta_{i,m}\eta_{j,m}\omega_{m}}{\mu^{2}-\omega_{m}^{2}}, \end{eqnarray*}$$

where μ is the Raman beatnote frequency, Ω is the effective Rabi frequency, which is assumed to be the same for all of the ions, {ω_m} are the phonon mode frequencies of the ions in the x direction and η_i,m are the Lamb–Dicke parameters measuring the coupling between the ith ion and the mth phonon mode. We are interested in the region where $B\gg \max \{J_{ij}\}$. In this limit, the σ⁺_iσ⁺_j and σ⁻_iσ⁻_j terms in equation (1) will be energetically forbidden, and we end up with the XY Hamiltonian

$$\begin{equation} H\approx H_{XY}=\sum_{i<j}2J_{i,j}(\sigma_{i}^{+}\sigma_{j}^{-}+{\rm h.c.})+B\sum_{i}\sigma_{i}^{z}. \end{equation} \tag{ 2 }$$

A unique feature of the Hamiltonians (1) and (2) realized with the ion system is that the interaction characterized by J_ij is long ranged, and the interaction range is readily tuned by changing the beatnote frequency μ. In particular, in a range of μ, J_ij can be roughly approximated by a power-law decay with J_ij ∼ |i − j|^−α, where α varies from 0 to 3 when we tune μ [19]. In the following analysis, for a given μ, we fit the coefficient J_ij with J_ij ∼ |i − j|^−α and use the fitting parameter α as an indicator of the range of interaction.

Zoom In Zoom Out Reset image size

The 1/r^α Ising model (equation (1)) with no transverse field (essentially classical) can be analytically solved [27, 28]. Its dynamic properties have also been recently studied in [24, 29], where it is found that if α is smaller than the dimensionality of the system, the total interaction energy per site diverges, drastically changing the dynamics after a global quantum quench, with possible emergence of prethermalization. However, the long range XY model (equation (2)) contains non-commuting terms and is generally not exactly solvable, making its dynamic properties largely unknown. Here, we consider a particularly interesting setup that shows rich dynamic properties and is practical for experimental demonstration. We first initialize all of the spins through optical pumping to the spin down state with σ^z_i = −1, which is an eigenstate of H_XY and hence stationary. We then use a focused laser beam to flip the first spin (left end ion) to σ^z_i = 1 (see figure 1). The starting state |ψ(0)〉 = |↑↓↓···↓〉 is no longer an eigenstate of H_XY , and is subject to the dynamics with |ψ(t)〉 = e^{iH_XY t/ℏ}|ψ(0)〉. We consider the time evolution of the local observables 〈σ^z_i(t)〉 and their correlations 〈σ^z_i(t)σ^z_j(t)〉, which can be directly measured in the experiments. For convenience of description of the dynamics, we introduce the operator

$$\begin{eqnarray*} &&C=\sum_{i=1}^{N}f_{i}\frac{\sigma_{i}^{z}+1}{2}, \end{eqnarray*}$$

where the coefficients f_i ≡ (2i − N − 1)/(N − 1) are equally distributed between [ − 1,1] from i = 1 to N. The expectation value of C varies between [ − 1,1], and physically measures the relative position of the spin excitation. It is easy to check that 〈ψ(0)|C|ψ(0)〉 = −1, meaning that the spin excitation is initially at the left edge of the chain. We also note that for any state with spatial inversion symmetry around the center of the chain, 〈C〉 = 0.

To explore the dynamic behavior, we first perform numerical calculations by using exact diagonalization of an N = 16 ion chain, which corresponds to the current size of the experimental trapped ion quantum simulator in [19]. As shown in figure 2, we choose two parameter settings with the corresponding fitting parameter α ≈ 2.6 and 0.52, which represent, respectively, the short-range and the long-range interactions. The distributions of the exact coupling coefficients J_ij are shown in figure 2 for these two cases. Note that our calculation is based on the approximated XY Hamiltonian (equation (2)). However, we have numerically verified that the relative differences of the calculated 〈C〉 and 〈σ_iσ_j〉 obtained by solving the original transverse field Ising Hamiltonian (equation (1)) are negligible at all times (< 10⁻³) when B > 1000J₀. It is worth mentioning that the error caused by the non-$\sum _i \sigma _{i}^{z}$-conserving terms from the finite B field does not accumulate over time: the population of our state outside of the single excitation subspace can be estimated as $\epsilon \sim \sum _{i,j} (J_{ij}/B)^2$, which is independent of time and can be made arbitrarily small by choosing a sufficiently large B. In practice, as J₀ ∼ Ω(η_x/N)(η_xΩ)/(μ − ω_x), and a trapped ion quantum simulator requires η_x ≪ 1 and η_xΩ ≪ μ − ω_x [26], Ω > 1000J₀ is typical for N > 10 ions (see [19] and also example parameters in figure 2). Since the transverse magnetic field is generated by a resonant carrier transition [16, 17], by splitting the Raman laser power evenly between the two sidebands alongwith the carrier transition, B > 1000J₀ can be readily achieved. Moreover, decreases with increasing N, as the J_ij will scale down faster than 1/N for a given laser power [30].

Zoom In Zoom Out Reset image size

For these choices of parameters, the short-time dynamics with t∈[0,2/J₀] for 〈C〉 is shown in figures 3(a) and (b). In the short-range interaction case, we observe that the spin excitation, initially located at the left edge of the chain, travels almost coherently to the other side and oscillates back and forth with relatively small dispersion. In contrast, in the long-range interaction case, the spin excitation diffuses and somehow gets locked before it reaches the middle of the chain (with 〈C〉 ≈ −0.4).

Zoom In Zoom Out Reset image size

To clarify the long-time dynamics, we use the finite-time-averaged quantity $\overline {A(t)}$, defined as $\overline {A(t)}\equiv \frac {1}{t}\int _{0}^{t}\langle A(\tau )\rangle \,{\mathrm { d}}\tau$ for the observable A [31, 32], thereby averaging out temporal fluctuations on short time scales (the following dynamic behaviors are qualitatively the same even without performing any time averaging). The long-time dynamics is shown in figures 3(c) and (d). In the short-range interaction case, the spin excitation position C, as well as the spin correlations σ^z_iσ^z_j, relax to the stationary values at T₀ ≈ 10/J₀ = 500 ms. In the long-range interaction case, the observables first reach quasi-stationary (prethermal) values at a timescale T₀, and further relax to the stationary (thermal) values at a much longer timescale (∼ 10⁴T₀). The emergence of prethermalization in the long-range interaction case is signaled by a nonzero value of $\overline {C}$ at the intermediate timescale T₀. We can use $\overline {C}$ as an order parameter to characterize different dynamic behaviors. By continuously tuning the effective interaction range (indicated by the fitting parameter α) with the beatnote frequency μ, we find that (see figure 4) prethermalization only takes place when α is smaller than a critical value (α_C ≈ 1.3 for N = 16). For larger system size, the prethermalization–thermalization transition still occurs, but α_C becomes smaller and the transition becomes sharper. For the particular case α = 0, the system has uniform coupling and its dynamics can be solved exactly. The exact solution shows that the system stays in the prethermal state forever with $\overline {C}=\frac {2}{N}-1$ [33].

Zoom In Zoom Out Reset image size

We now give a physical explanation of why prethermalization and the variety of observed dynamic behaviors occur in this model. The distinctive short-time dynamics of 〈σ^z_i(t)〉 (figures 3(a) and (b)) can be explained by examining the energy spectrum of H_XY in the single spin excitation subspace (shown in figure 5(a)). In the short-range interaction case, the energy spectrum is close to linear. This is because H_XY can be roughly approximated with only neighboring interaction, and due to inhomogeneous lattice spacing, {J_i,i+1} are similar to those in the 'quantum mirror' model [34, 35], resulting in a near dispersion-free spin wave propagation until nonlinearity sets in. On the other hand, the energy spectrum for the long-range interaction case is highly nonlinear, so the dynamics of the spin excitation is strongly dispersive.

Zoom In Zoom Out Reset image size

The reason for a prethermal stage in the long-time dynamics, however, is much more complicated. Naively, the spin flip–flop matrix J_ij varies smoothly among sites for any α∈(0,3), so the spin excitation should continuously diffuse from one end of the chain to the whole chain, and is not expected to get trapped somewhere in the middle for a long time. The two stage dynamics indicates that there are two different timescales implicit in the Hamiltonian, which is not at all obvious by looking at the J_ij. The time dynamics of any physical observable is simply given by

$$\begin{eqnarray*} &&\langle A(t)\rangle=\sum_{m,n}\rho_{mn}(0)A_{nm}\,{\rm e}^{{\rm i}(E_{m}-E_{n})t/\hbar}, \end{eqnarray*}$$

where ρ_mn(0) is the initial state's density matrix element in energy basis, the different timescales of the dynamics can be understood by computing the eigenenergy differences {E_m − E_n}, as done in figure 5. In the short-range interaction case (figure 5(b)), all {E_m − E_n} are continuously distributed from J₀ to 100J₀, so a single-stage relaxation is expected after T ∼ T₀ = 10/J₀. In the long-range interaction case (figure 5(c)), most {E_m − E_n} still fall into the range of 1–100J₀, but there is a striking separate (i.e. gapped) branch at the much lower rate (∼ 10⁻⁶J₀). This branch results from the near-degenerate pairs ({E_2k − E_2k−1}) of the eigenenergies (figure 5(a)) that make up the first half energy spectrum, and the number of these pairs scale up with system size N. Note that if we put the ions into a ring trap so that the lattice is translationally invariant, we find that there is no separate branch of the {E_m − E_n} (figure 5(d)), and hence there is no prethermalization even with long-range interactions. The appearance of near-degenerate pairs {E_2k − E_2k−1}) in the spectrum of our model seems to arise from the combined effect of long-range interactions and inhomogeneous lattice spacing. Intuitively, from figure 2(b), the inhomogeneity causes the middle part of the chain to have smaller flip–flop interactions J_ij than the ends, potentially leading to a metastable state where the spin excitation is roughly localized in the left half of the chain. However, the existence of two distinct timescales is primarily due to the structure of the energy spectrum, which is not trivially related to inhomogeneity. We are still investigating if there is a wide class of inhomogeneous lattices that can give rise to similar dynamics.

Figure 5(e) shows that the thermal values can be well predicted by the DE, defined as

$$\begin{eqnarray*} &&\rho_{\rm DE}=\rho_{mn}(0)\delta_{mn}. \end{eqnarray*}$$

In the large N limit, under certain conditions implicit in the eigenstate thermalization hypothesis [3–5], the DE prediction will match the canonical ensemble prediction for thermalized states in classical statistical physics.

As prethermalization is due to an emergent small scale of energy differences in the Hamiltonian, to predict the prethermal values here, we can define the PDE:

$$\begin{equation} \rho_{\rm PDE}=\left\{\begin{array}{@{}ll@{}} \rho_{mn}(0)\delta_{mn}, & |\nu_{m}-\nu_{n}|\gtrsim1/T_{0},\\ \rho_{mn}(0), & |\nu_{m}-\nu_{n}|\ll1/T_{0}, \end{array}\right. \end{equation} \tag{ 3 }$$

where {ν_m} are the eigenenergies. We find that the PDE can accurately predict the prethermal values of local observables $\overline {\sigma _{i}^{z}}$, as shown in figure 5(d). Roughly speaking, the prethermalization timescale is determined by the average level spacing (∼ 10/J₀), and the thermalization timescale is determined by the minimum level spacing (∼ 10⁴/J₀ as in figure 5(d)).

The dynamic phase transition can be associated with breaking of the lattice inversion (parity) symmetry, reminiscent of symmetry breaking in the equilibrium phase transitions. Our Hamiltonian H_XY is symmetric under the space inversion around z = 0 (figure 1), but we start from an initial state that does not have this symmetry. The thermal state, with no memory of the initial state, is described by the diagonal ensemble ρ_DE and restores this symmetry as $\overline {C}=0$. However, the prethermal state does not restore the Hamiltonian symmetry due to its non-zero $\overline {C}$ value, which indicates that some 'memory' of the initial state is preserved in the prethermal state. The intermediate timescale T₀ for observation of the prethermal state gives a microscopic interpretation of why this state can break the parity symmetry: one cannot distinguish the near-degenerate pairs of eigenstates in the energy spectrum, so linear combinations within each pair are allowed. Since the two eigenstates of the pair have either even or odd parity, their linear combinations can break the parity symmetry. The dynamic phase diagram shown in figure 4 also has the hint of two non-analytical points: one is where $\overline {C}$ becomes non-zero, representing the appearance of prethermalization, and the other is where $\overline {C}$ approaches $\frac {2}{N}-1$ (α → 0), representing the disappearance of thermalization. Note that the thermodynamic limit of our system is not well defined, as to keep a linear ion chain stable, the length of the chain L cannot be kept proportional to N as N becomes large. Thus, it is not clear if a dynamical phase transition point can be obtained in some generalized thermodynamic limit (we are also limited by a finite number of lattice sites that can be studied using our numerical methods).

The transverse field Ising Hamiltonian equation (1) has already been experimentally simulated in [19] for the N = 16 ions, and highly efficient in situ measurement of spin polarization (σ^z_i) and spin correlation (σ^z_iσ^z_j) has been demonstrated. The XY Hamiltonian (equation (2)) can thus be readily obtained by turning up the effective transverse magnetic field. The non-equilibrium initial state preparation requires a focused laser beam, but is relatively easy due to the large ion spacing near the ends. The laser power and trap frequencies used for generating the interaction pattern J_ij shown in figure 2 are within current experimental reach [17, 19]. The observation of prethermalization and the dynamic phase transition shown in figures 3 and 4 only requires the spin decoherence time to be longer than T₀ = 10/J₀ = 500 ms, and coherence times of up to 2.5 s have been experimentally achieved by using the hyperfine qubit of Yb⁺ ions [36]. However, the second stage of thermalization for the long-range interaction case requires much longer times, and is beyond current experimental reach, similar to the experiment on prethermalization using cold atoms in [14].

In summary, we have proposed a novel scheme to observe a peculiar prethermalization phenomenon and dynamic phase transitions in a trapped ion quantum simulator. The required conditions can be realized by using current experimental technology. We provide an explanation of the mechanism of prethermalization and the dynamic phase transition in our proposed model. We believe that the wealth of interesting dynamic properties for a system with continuously tunable range of interactions is far from exhausted.

We thank C Monroe, A Polkovnikov and A Gorshkov for helpful discussions, and Michael Foss--Feig for improving the writing of this paper. This work was supported by the NBRPC (973 Program) 2011CBA00302, the DARPA OLE program, the IARPA MUSIQC program and the ARO$/$AFOSR MURI program.