Interacting quantum systems cannot generally be efficiently simulated using classical computational techniques, due to the exponential scaling of the information encoded in a quantum state with the size of the system. Accurate modeling of quantum phenomena such as the magnetism of interacting spins, material superconductivity, electronic structure and molecules and their chemistry, Fermi-Hubbard models of electron transport in solids, and interactions within atomic nuclei, are all beyond the reach of classical computation even for just small numbers of interacting degrees of freedom. A quantum simulator exploits the controlled manipulation of a standard quantum system in order to mimick the properties or evolution of a different-perhaps intractable-quantum system (Feynman, 1982;Trabesinger, 2012). A quantum simulator is a restricted quantum computer (Ladd et al., 2010;Nielsen and Chuang, 2000), with operations that may not be universal but instead be tailored to a particular quantum physical model under study. While useful general purpose and universal quantum computers are still a future prospect, special purpose quantum simulations have already been demonstrated and may in fact become the first useful application of quantum computers.
The equivalence of quantum bits (qubits) and their quantum gates to effective spin-1/2 systems and their interactions has guided one of the most important classes of quantum simulations: interacting spin systems and quantum magnetism. The most advanced physical system for quantum bits or effective spins is arguably a collection of trapped atomic ions (Brown et al., 2016;Monroe and Kim, 2013;Steane, 1997;Wineland et al., 1998;Wineland and Blatt, 2008). This is evidenced by the number of controlled and interacting qubits, the quality of quantum gates and interactions, and the high fidelity of quantum state initialization and measurement. Trapped atomic ions, held in a vacuum chamber and confined by electromagnetic fields to be distant from surfaces or solids, laser cooled to be nearly at rest and "wired" together with external laser or microwave fields, offer a very clean quantum system to perform quantum simulations (Blatt and Roos, 2012;Deng et al., 2005;Porras and Cirac, 2004;Taylor and Calarco, 2008).
This review assesses recent progress in the quantum simulation of magnetism and related phenomena using trapped atomic ion crystals. Following an introduction of the mapping of spins to atomic ions, we first cover experimental results on the simulation of magnetic ordering, equilibrium states, and phase transitions in quantum magnetic systems. Then we move to dynamical studies of quantum magnetism, touching on general issues of information and entanglement dynamics, quantum thermalization, inhibitors to thermalization such as many-body localization and prethermalization, "time crystals," dynamical phase transitions, and Hamiltonian engineering and sequencing techniques that translate quantum physics models to effective magnetic spin models. Recent techniques of variational quantum simulations (VQS) and quantum approximate optimization algorithms are also reviewed in the context of trapped ion systems. We conclude with speculations on how these types of simulations with trapped ions may scale in the future and perhaps guide the development of real magnetic material function or more general quantum simulations as special cases of quantum computations.
Trapped ion spin simulators generally exploit the collective motion of the ions to generate many-body spin interactions, with a controlled coupling through collective bosonic harmonic oscollator modes, or phonons. While not a focus of this review, the phonon degree of freedom itself is an interesting medium to host simulations of bosonic models and are complementary to the field of photonic cavity-QED in atomic (Raimond et al., 2001;Ye et al., 2008) and superconducting (Wallraff et al., 2004) systems. There have been many demonstrations of phonon control and phonon-based simulations in trapped ion systems (Ben-Kish et al., 2003;Kienzler et al., 2017;Um et al., 2016), including bosonic interference (Leibfried et al., 2002;Toyoda et al., 2015), blockades (Debnath et al., 2018;Ohira et al., 2020), sampling molecular vibronic spectrum (Shen et al., 2018), quantum walks (Schmitz et al., 2009;Zähringer et al., 2010), phononic "lasing" (Ip et al., 2018;Vahala et al., 2009), quantum thermodynamics (An et al., 2015), and quantum heat engines/refrigerators (Maslennikov et al., 2019;Roßnagel et al., 2016).
We note there are many other quantum systems amenable to quantum spin simulation, including superconducting circuits (Arute et al., 2019;Kandala et al., 2017;Otterbach et al., 2017), neutral atoms in optical lattices (Gross and Bloch, 2017), and solid-state arrays of effective individual spin systems (Choi et al., 2017a,b). Comparisions to these other platforms is made where appropriate, although many of the simulations covered in this review rely on individual quantum state control with long range interactions that are not native to other platforms.
We represent a collection of spins by a crystal of electromagnetically trapped atomic ions, with two electronic energy levels within each ion behaving as an effective spin-1/2 particle. The particular choice of electronic levels depends on the atomic element and also the desired type of control fields used to manipulate and measure the spin state. The most important features of these spin states for executing quantum simulations are (a) the spin levels are long-lived and exhibit excellent coherence properties, (b) the spin states have appropriate strong optical transitions to auxiliary excited states, with selection rules allowing for qubit initialization through optical pumping and qubit detection through fluorescence, and (c) the spins interact through a coherent coupling that can be externally controlled and gated. These requirements restrict the atomic species to a handful of elements and spin states that are encoded in either S 1/2 hyperfine or Zeeman ground states of single outer-electron atoms (e.g., Be + , Mg + , Ca + , Sr + , Ba + , Cd + , Zn + , Hg + , Yb + ) with radiofrequency/microwave frequency splittings or ground and D or F metastable electronic excited states of single or dual outer-electron atoms (e.g., Ca + , Sr + , Ba + , Yb + , B + , Al + , Ga + , In + , Hg + , Tl + , Lu + ) with optical frequency splittings. Some species (e.g., Ba + , Lu + , Yb + ) have sufficiently long D or F metastable excited state lifetimes to host spins in their hyperfine or Zeeman levels with radiofrequency/microwave splittings.
In any of these systems, we label the two relevant spin states as |↓ ≡ |↓ z and |↑ ≡ |↑ z , eigenstates of the Pauli operator σ z separated by energy E ↑ -E ↓ = ω 0 . In the transverse bases of the Bloch sphere, we define by convention the eigenstates of σ x as |↓ x ≡ (|↓ -|↑ )/ √ 2 and |↑ x ≡ (|↓ + |↑ )/ √ 2, and the eigenstates of σ y as |↓ y ≡ (|↓ + i |↑ )/ √ 2 and |↑ y ≡ (i |↓ + |↑ )/ √ 2. We note that in this review, the spin-direction is often not explicitly labeled, in which cases the notation adopts the context of the material.
A typical quantum simulation in the ion trap system is comprised of three sequential steps: initialization, interaction, and measurement, as depicted in Fig. 1. The spins are initialized through an optical pumping process that prepares each spin in a nearly pure quantum state (Happer, 1972). By applying resonant laser radiation that couples the spin states to short-lived excited states, each spin can be initialized with > 99.9% state purity in a few microseconds. This relies on appropriate selection rules for the excited states as well as sufficiently split energy levels of the spin states themselves (Fig. 1a). Laser cooling can prepare the motional states of the ions to near the ground state of harmonic motion (Leibfried et al., 2003), which is important for the control of the spin-spin interactions as detailed below . 11 Each spin can be coherently manipulated by driving the atomic ions with external fields that couple the spin states. This can be accomplished by resonantly driving the spin levels with appropriate radiation at frequency ω 0 , or in Fig. 1b, this is depicted as a two-field Raman process, with a beatnote of two optical fields at ω 0 driving the spin (this will be assumed throughout unless otherwise indicated). This coherent coupling can also drive motional sideband transitions (Leibfried et al., 2003) that couple the spin to the motion of the ion. For multiple ions, this can be used to generate spin-spin couplings mediated by the Coulomb interaction, described in more detail below (Wineland and Blatt, 2008). These external fields provide exquisite control over the effective spin-spin interaction, with the ability to gate the interaction, program different forms of the interaction strength and range, and even reconfigure the interaction graph topology.
At the end of the quantum simulation, the spins are measured by applying resonant laser radiation that couples one of the two spin states to a short-lived excited state through a cycling transition and detecting the resulting fluorescence (Bergquist et al., 1986;Nagourney et al., 1986;Sauter et al., 1986). This is depicted in Fig. 1c, where we take the |↑ or "bright" state as fluorescing and the |↓ or "dark" state as not fluorescing. Even though the photon collection efficiency may be small (typically less than 1%), the effective spin detection efficiency can be well above 99% owing to the low probability of leaving the fluorescence cycle or having the other (dark) state entering the cycle (Acton et al., 2006;Benhelm et al., 2008;Hume et al., 2007;Myerson et al., 2008;Noek et al., 2013). In order to detect the spins in other bases in the Bloch sphere (σ x or σ y ), previous to fluorescence measurement the spins are coherently rotated by polar angle π/2 along the y or x axis of the Bloch sphere, respectively.
Atomic ions can be confined in free space with electromagnetic fields supplied by nearby electrodes. Two types of ion traps are used for quantum simulation experiments: the linear radiofrequency (rf) trap and the Penning trap. The linear rf trap (Fig 2a) (Raizen et al., 1992) juxtaposes a static axial confining potential with a two-dimensional rf quadrupole potential that provides a time-averaged or ponderomotive transverse confinement potential (Dehmelt, 1967;Paul, 1990). The trap anisotropy is typically adjusted so that the static axial confinement is much weaker than the transverse confinement so that laser-cooled ions reside on the axis of the trap where the rf fields are null, resulting in a one-dimensional chain of trapped ion spins. A harmonic axial confinement potential results in an anisotropic linear ion spacing (James, 1998), but they can be made nearly equidistant by applying an appropriate quartic axial confining potential (Lin et al., 2009). The Penning trap (Fig 2b) employs a uniform axial magnetic field with static axial electric field confinement, where the transverse confinement is provided by the E × B drift toward the axis (Bohnet et al., 2016;Brown and Gabrielse, 1986). Here, the trap anisotropy is typically adjusted so that the ions form a two-dimensional crystal perpendicular to and rotating about the axis. Both rf and Penning traps can be modified to support other types of crystals in any number of spatial dimensions, but the quantum simulations reviewed here are either 1-D chains in rf traps or 2-D crystals in Penning traps. However, the dimensionality of the spin-spin interaction graph does not necessarily follow that of the spatial arrangement of spins.
Ions are typically loaded into traps by generating neutral atoms of the desired element and ionizing the atoms once in the trapping volume via electron bombardment or photoionization. Ion trap depths are usually much larger than room temperature, so rare collisions with background gas do not necessarily eject the ion from the trap, but they can temporarily break up the In the case of ground-state (e.g., Zeeman or hyperfine) defined spins separated by frequency ω0, two tones of off-resonant radiation (purple lines) can drive stimulated Raman transitions between the spin states. The two beams have resonant Rabi frequencies g1,2 connecting respective spin states to excited states and are detuned by ∆ γ, and have a difference frequency (beatnote) detuned from the spin resonance by µ. This coherently couples the spin states to create superpositions (µ = 0) and for non-copropagating beams also couples to the motion of the ion crystal (µ = 0). These processes can also be driven directly by radiofrequency or microwave signals, or for optically-defined spin states, a single laser tone. (c) Resonant radiation (blue line) drives the |↑ ↔ |e cycling transition, causes the |↑ state to fluoresce strongly (red line), while the |↓ state is far from resonance and therefore dark. This allows the near-perfect detection of the spin state by the collection of this state-dependent fluorescence depicted by the eyeball.
crystal and scramble the atomic ions in space. Under typical ultra-high-vacuum conditions, these collisional interruptions occur roughly once per hour per ion (Wineland et al., 1998), but cryogenic vacuum chambers can reduce the collision rate by orders of magnitude, where the trapped ions can be undisturbed for weeks or longer between collisions.
When atomic ions are laser-cooled and localized well below their mean spacing, they form a stable crystal, with the Coulomb repulsion balancing the external confinement force. Typical spacings between adjacent ions in trapped ion crystals are ∼ 3 -20 µm, depending on the ion mass, number of ions in the crystal, the characteristic dimensions of the electrodes, and the applied potential values. The equilibrium positions of ions in the crystal can be calculated numerically (James, 1998;Sawyer et al., 2012;Steane, 1997). The motion of the ions away from their equilibrium positions is well-described by harmonic normal modes of oscillation (phonon modes), with frequencies in the range ω m /2π ∼ 0.1 -10 MHz. The thermal motion of laser-cooled ions and also the driven motion by external forces is typically at the 10 -100 nm scale. This is much smaller than the inter-ion spacing, so nonlinearities to the phonon modes (Marquet et al., 2003) can be safely neglected and the harmonic approximation to the phonon modes is justified. Calculations of the phonon mode frequencies and normal mode eigenfunctions follow straighforwardly from the calculated ion spacings (James, 1998;Sawyer et al., 2012;Steane, 1997). For the systems described here, we consider primarily the motion along a single spatial dimension labeled X. We write the X-component of position of the ith ion as Xi = Xi + xi , separating the mean (stationary) position Xi of the ith ion from the small harmonic oscillations described by the quantum position operator xi . The motion of ions in the crystal is tightly coupled by the Coulomb interaction, so it is natural to express the position operator in terms of the m = 1 . . . N normal (phonon) modes, xi = N m=1 b im ξm , where b im is the normal mode transformation matrix, with i b im b in = δ nm and m b im b jm = δ ij . Each phonon mode ξm oscillates at frequency ω m , and can be described as a quantum harmonic oscillator with zero-point spatial spread ξ (0) m = /2M ω m , where M is the mass of a single ion. In the interaction frame for each phonon mode, the position of the ith ion is thus written as
where a † m and a m are bosonic raising and lowering operators for each mode, with [a n , a † m ] = δ nm . In general, the structure of transverse phonon modes of a 1D or 2D ion crystal (motion perpendicular to the axis or plane of the ions) has the center-of-mass (COM) mode as its highest frequency, with the lowest frequency corresponding to zig-zag motion where adjacent ions move in opposite directions, as shown in Fig. 3 for a 1D linear chain of 32 ions and a 2D crystal of about 345 ions. The bandwidth of the transverse modes can be controlled by tuning the spatial anisotropy of the trap.
Effective magnetic fields and spin-spin interactions can be realized by applying external microwave or optical fields to the ions. We consider the case of optical fields, since they can be used to not only provide effective site-dependent magnetic fields for the spins through tight focusing, but the strong dipole forces from laser beams can also drive effective Ising interactions between the spins (Cirac and Zoller, 1995;Milburn et al., 2000;Solano et al., 1999;Sørensen and Mølmer, 1999). Such forces can be applied to pairs of ions in order to execute entangling quantum gates applicable to quantum computing (Wineland and Blatt, 2008). When such forces are instead applied more globally, the resulting interaction network allows the quantum simulation of a wide variety of spin models such as the Ising and Heisenberg spin chains (Deng et al., 2005;Porras and Cirac, 2004;Taylor and Calarco, 2008).
Following Fig. 1b and assuming the spins are encoded in stable (e.g., Zeeman or hyperfine) levels, the ion crystal is addressed with two laser beams detuned from the excited states by much more than the excited state radiative linewidth (∆ γ). By adiabatically eliminating the excited states, this drives coherent Raman transitions between the spin states. Alternatively, optically-defined spin levels can be coupled with a single laser beam (Benhelm et al., 2008), but this is more difficult technically, as the spins acquire an optical phase that requires extreme positional stability of the optical setup. With rf-or microwave-defined spin states, the relevant phase is that of the microwave beat note between the Raman beams, which has much longer wavelength and is easier to control and stabilize.
The two Raman beams are oriented to have a projection δk of their wavevector difference along the X-axis of motion. The Raman beatnote is detuned by frequency ω 0 + µ i from the resonance of spin i with beat note phase δφ i . The resonant Raman Rabi frequency on ion i is Ω i = g i 1 g i 2 /2∆, where g i 1,2 are the direct (single field) Rabi frequencies of the associated transitions through the excited states (see Fig. 1b), proportional to the respective applied optical electric fields. The atom-light interaction Hamiltonian in a frame rotating at the spin resonance frequency ω 0 Ω i (optical rotating wave approximation) takes the form ( = 1) H = 1 2 i Ω i σ i + e i(δk Xi-µit-δφi) + σ i -e -i(δk Xi-µit-δφi) + d i σ i z .
(
The last term is an AC Stark shift of the ith spin by amount d i and arises from differential AC Stark shifts between the two spin levels from the Raman beams. This shift includes the "two-photon" differential Stark shift terms scaled by (g 2 1 + g 2 2 )ω 0 /4∆ 2 and summed over each excited state detuned by ∆, where ω 0 ∆ (see Fig. 1). There are also "four-photon" Stark shifts scaled by Ω 2 i /4µ i and summed over each two-photon Raman resonance detuned by µ, where Ω i µ i . The magnitude of these shifts depends greatly on the atomic energy level structure and light polarization: see Lee et al., 2016 for a discussion of Raman-coupled qubits (Fig 1 ) and Häffner et al., 2003 for direct optically-coupled qubits. These Stark shifts can be either absorbed into the definition of the spin energy levels, or used as an effective axial magnetic field for simulations.
Below, we will assume that the ions are confined to the "Lamb-Dicke limit" (Leibfried et al., 2003;Wineland et al., 1998) where the excursion of ion motion is much less than the associated wavelength of radiation driving transitions: δk x2 1. This is typically a good assumption for trapped ions laser-cooled to near their ground state, with a zero-point spatial spread of all modes ξ (0) m that is typically of order nanometers in experiments.
For a resonant "carrier" interaction (µ = d i = 0) and under the rotating wave approximation (ω m Ω i ), the time dependence of X averages to zero and the Hamiltonian of Eq. ( 2) is just
where the transverse field spin operator is defined by
This Hamiltonian of Eq. ( 3) describes the precession of spin i about an effective transverse magnetic field in the xy plane of the Bloch sphere at an angle
that can be controlled through the phase δφ i . In cases where δk = 0, the Rabi frequency acquires a dependence on the motion of the ions through Debye-Waller factors (Wineland et al., 1998), but these are small in the Lamb-Dicke limit (Leibfried et al., 2003;Sørensen and Mølmer, 2000).
Tuning the Raman laser beat note away from the carrier with µ Ω i generally results in a four-photon AC Stark shift of the spin levels as discussed above, given by the last term in Eq. ( 2) (Häffner et al., 2003;Lee et al., 2016). When each spin is exposed to a unique intensity of light and/or detuning of a single beam, parametrized by d i , this gives rise to a site-dependent effective axial (z) magnetic field.
When the frequency µ i is tuned to the neighborhood of the phonon modes ω m , the spin i is coupled to the ion motion through the spatial variation of the phase factor in Eq. ( 2). This will generate an effective spin-spin interaction between the exposed ions, mediated by the collective transverse vibrations of the chain. For most simulation experiments, the transverse modes of motion are used to mediate the spin-spin interaction because their frequencies are tightly packed and all contribute to the effective spin Hamiltonian, allowing control over the form and range of the interaction, described further below. Transverse modes also oscillate at higher frequencies, leading to better cooling and less sensitivity to external heating and noise (Zhu et al., 2006).
In general, when noncopropagating laser beams form bichromatic beat notes at frequencies ω 0 ± µ i symmetrically detuned from the carrier with µ i ≈ ω m , both upper and lower motion-induced sidebands of the normal modes of motion are driven in the ion crystal, giving rise to a spin-dependent force at frequency µ i (Milburn et al., 2000;Porras and Cirac, 2004;Solano et al., 1999;Sørensen and Mølmer, 1999). Owing to the symmetry of the detuned beatnotes, the four-photon Stark shift is generally negligible. However, when the bichromatic beat notes are asymmetrically detuned from the carrier by ω 0 + µ i+ and ω 0µ i-, the effective spin-dependent force occurs at frequency µ i = (µ i+ + µ i-)/2 and the asymmetry provides a Stark shift that gives rise to an effective uniform axial magnetic field in Eq. ( 2) with
Under the rotating wave approximation (ω 0 µ i Ω i ) with symmetric detuning µ i = µ i+ = µ i-and within the Lamb-Dicke limit, the exponential function in Eq. ( 2) can be expanded to lowest order, resulting in (Zhu et al., 2006)
Here the beatnote detuning on ion i from the mth motional sideband is defined as δ im = µ iω m and the Lamb-Dicke parameter
m describes the coupling between ion i and mode m. The above expression includes two phases: a "spin phase" θ i that determines the angle of the ith spin operator in the XY plane of the Bloch sphere that diagonalizes the spin-dependent force; and a "motion phase" ψ i that determines the phase of the optical forces (but does not play a role in the spin-spin interactions as developed below). These phases depend on the geometry of the bichromatic laser beams and there are two cases written in Eqs. ( 7)-( 8) (Haljan et al., 2005;Lee et al., 2005). When the upper and lower sideband running wave beat notes propagate in the same direction (δk is the same sign for both), this is termed the "phase-sensitive" geometry. On the other hand, when the upper and lower sideband running waves propagate in opposite directions (opposite sign of δk for the two beat notes), this is called "phase-insensitive." The phases for each configuration are written,
Here, δφ i+ and δφ i-are the beat note phases associated with the upper and lower sideband fields, respectively. The additional π/2 phase factors compared with Eq. ( 5) originate from the imaginary linear term in the Lamb-Dicke expansion of e ikxi . The sensitivity of the spin phase to δk Xi has great importance on the practical implementation of spin-spin Hamiltonians. In the phase-sensitive geometry, this dependence may be desired when the phase of other Hamiltonian terms (such as the carrier spin flip operation of Eq. ( 3)) have the same form, as in Eq. ( 5). This allows a phase-sensitive diagnosis of individual spin operations. However, sensitivity to δk Xi over sufficiently long times can lead to decoherence if there are drifts in the relative path length of non-copropagating beams or the ion chain position along the X-direction. The Phase Insensitive configuration is therefore useful for long simulation evolution times.
We note that the phase-insensitive geometry using Raman couplings (Fig. 1b) can remove the spin-phase sensitivity to not only the absolute optical phase of the optical source, but also the relative optical path length difference between the counterpropagating beams, by setting δφ i+ = -δφ i- (Haljan et al., 2005;Lee et al., 2005). This is not possible with direct optical upper and lower sideband transitions on spins with an optical energy splitting, as their optical phases add regardless of the geometry.
For either phase configuration, the evolution operator under this Hamiltonian can be written from the Magnus expansion, which terminates after the first two terms (Zhu et al., 2006),
= exp
The first term of Eq. ( 10) is a spin-phonon coupling with operator
representing spin-dependent coherent displacements (Glauber, 1963;Leibfried et al., 2003) of the mth motional mode through phase space by an amount
The second term of Eq. ( 10) is the key result: a spin-spin interaction between ions i and j with coupling strength
Fractional corrections to this expression arising from higher-order terms in the Lamb-Dicke expansion leading to Eq. ( 6) can be shown to be of order η 2 im η 2 j,m n2 m for each mode (Sørensen and Mølmer, 2000), where nm is the average number of vibrational quanta in mode m.
In this review, we generally consider global spin Hamltonians that are realized by exposing all the ions to the spin-dependent force. However, for the special case of applying a spin-dependent force to just two ions i and j in the chain, which is common for the execution of entangling two-qubit quantum logic gates (Sørensen andMølmer, 1999, 2000), the evolution operator in Eq. ( 10) reduces to
In this expression, the spin projection operators are eigenvectors of σ θi : (Glauber, 1963).
There are two regimes where the collective modes of motion contribute to the spin-spin coupling, taking evolution time τ to be much longer than the ion normal mode oscillation periods (ω m τ 1). In the "resonant" regime (Milburn et al., 2000;Solano et al., 1999;Sørensen and Mølmer, 2000), the optical beatnote detuning µ is close to one or more normal modes and the spins become entangled with the motion through the spin-dependent displacements. However, at certain times of the evolution α i,m (τ ) ≈ 0 for all modes m and the motion nearly decouples from the spin states, which is useful for applying synchronous entangling quantum logic gates between the spins. For closely-spaced modes such as the transverse modes as seen in Fig. 3, resolving individual modes becomes difficult and may require laser-pulse-shaping techniques (Zhu et al., 2006).
For generating pure spin Hamiltonians, we instead operate in the "dispersive" regime (Porras and Cirac, 2004;Sørensen and Mølmer, 1999), where the optical beatnote frequency is far from each normal mode compared to that mode's sideband Rabi frequency (|δ im | η im Ω i ). In this case, the phonons are only virtually excited as the displacements become negligible (|α i,m (t)| 1). The spin-phonon part of the evolution (Eq. ( 11)) is therefore negligible, and the spin-spin interaction evolution (Eq. ( 10)) is dominated by the secular terms of Eqs. ( 14) and ( 15) that are linear in time τ . The final result is an effective fully-connected Ising Hamiltonian,
where the Ising matrix is given by
Here we have used ω rec = (δk) 2 /2M as the recoil frequency associated with the transfer of momentum (δk) to a single ion.
Substituting the exact values for the normal mode matrix b im and assuming that the optical force is detuned at frequencies higher than all phonon modes (δ m > 0), we find that for a uniform Rabi frequency over the ions Ω i = Ω, the Ising matrix is well approximated by a long-range antiferromagnetic (AFM) coupling that fall off as an inverse power law with distance
with nearest-neighbor Ising coupling J 0 . The exponent α that determines the range of the Ising interaction can be set to 0 < α < 3 by simply adjusting the laser detuning µ > ω m (Islam et al., 2011;Porras and Cirac, 2004). The true asymptotic long-chain behavior of a trapped ion chain is more subtle (Nevado and Porras, 2016), but the power law approximation is very good, as shown in Fig. 4, comparing numerically exact couplings with best-fit power laws for various detunings in both a linear and 2D crystal.
When the detuning µ is tuned between the modes of motion, many other patterns of the Ising graph can be realized (Korenblit et al., 2012;Lin et al., 2011). The spin-spin interaction profile J ij can in principle be programmed arbitrarily with sufficient control of the individual spin-phonon couplings on N trapped ions. For example, N unique bichromatic Raman beatnotes applied at frequencies ω 0 ± µ n , with µ n ≈ ω n (n = 1, 2, ..., N ) can be used with local intensity control over each ion to create an arbitrary interaction graph:
Here, Ω i,n is the Rabi frequency corresponding to the n th beatnote on the i th ion. Note that J ij is nonlinear in the O(N 2 ) experimental control parameters Ω i,n and µ n , and hence tuning the quantum simulator requires non-linear optimization methods (Korenblit et al., 2012;Teoh et al., 2020). Alternative approaches to realize a target interaction graph without tuning the full Rabi frequency matrix, Ω i,n include modifying a global Mølmer-Sørensen coupling profile (such as Eq.( 23)) by local spatial control of spins in hybrid analog-digital ways (Hayes et al., 2014;Rajabi et al., 2019). This review will mainly consider Ising interactions in the slow dispersive regime in order to engineer pure spin Hamiltonians given by Eqs. ( 3) and (20) that do not directly involve the bosonic phonon operators. An important class of models in quantum magnetism that will appear throughout this review is the transverse field Ising model, which is one of the simplest physical models that admits a quantum phase transition (Sachdev, 2011), owing to its noncommunting spin operators: 22) in a 2D crystal of 217 ions versus a sampling of the distance dij between ion pairs (circles). The lines are best-fit power law exponents α (lines) for various detunings from the center-of-mass (COM) mode of 795 kHz. Adapted from Britton et al., 2012. In ion trap systems, this model can be generated with a combination of the effective magnetic field in Eq. ( 3) and the Ising interactions in Eq. ( 20) with appropriate settings of the spin phases θ i in Eq. ( 5) and Eqs. ( 7)- (8). When both (noncommuting) terms are simultaneously applied, additional spin-dependent (σ i z ) phonon terms appear in higher orders of the Magnus expansion of Eq. ( 9) (Wang andFreericks, 2010, 2012), which is discussed below in Section II.A.
It is often desired to implement interacting spin models with multiple components of Ising interactions along various axes of the Bloch sphere, with the most general being the anisotropic Heisenberg model involving a sum of all three Ising terms
Subclasses of the Heisenberg model arise naturally in physics and can possess useful symmetries (Sachdev, 2011). For example, the isotropic Heisenberg model (J x ij = J y ij = J z ij ) is relevant to natural 3D magnetic interactions. The XXZ model (J x ij = J y ij ) and the XY model (J x ij = J y ij ; J z ij = 0) result from standard transformations of Fermionic Hubbard models to spin models (Bravyi and Kitaev, 2002;Jordan and Wigner, 1928), which is discussed in Sec. V. The XXZ and XY models also conserve the z-component of the total spin in the system, allowing simplifications to the spin dynamics and the interpretation of measurements.
Ion trap quantum simulators can generate generic Heisenberg models with the primitive of the Ising interaction discussed above, using a variety of extensions. First, the three separate Ising terms in the Hamiltonian can exploit three independent modes (or spatial directions) of motion (Porras and Cirac, 2004), although this may not easily produce the same form or range of Ising interactions for all three axes. Second, the desired interactions can be applied sequentially in a Trotter expansion of the desired Hamiltonian (Johri et al., 2017;Lanyon et al., 2011;Lloyd, 1996;Suzuki, 1985;Trotter, 1959), as will be discussed in section V. Third, a transverse field Ising model (Eq. 25) can be applied with B y J ij so the field overwhelms the Ising interactions. This can be seen by expanding the Ising term as
Here the strong field energetically forbids double spin-flips, or equivalently bestows fast oscillations to the σ i ± σ j ± terms, which average to zero in an effective rotating-wave approximation (Cohen et al., 2015).
A compelling Hamiltonian quantum simulation usually results in some type of nontrivial ground state or dynamics that may elude classical computation. It is therefore important to verify that the desired Hamiltonian is indeed being faithfully implemented by the quantum simulator (Cirac and Zoller, 2012;Hangleiter et al., 2017). For systems that are small enough and tractable for a direct comparison between the simulator's results and a theoretical calculation, this will provide some confidence that the proper simulation has been run. But scaling up the system can introduce additional imperfections that may call into question the accuracy of the applied Hamiltonian.
Two approaches for verifying quantum simulators beyond classically computability are the use of fault-tolerant techniques in the expression of the simulation in terms of discrete error-corrected gates (Gottesman, 1998;Preskill, 1998), and the comparison of the results of multiple quantum simulators built upon different platforms (Leibfried, 2010). Though this list is not comprehensive (see, e.g. (Hangleiter et al., 2017)), we examine these two approaches briefly below.
Universal Hamiltonian digital quantum simulators (DQS) as will be discussed in section V break up the simulation evolution into a series of time steps, and the error introduced by this "Trotter" expansion is bounded and inversely proportional to the number of steps M (Lloyd, 1996). This approach was notably employed for the implementation of various Hamiltonians, including many-body magnetic couplings, in a system of trapped ions (Lanyon et al., 2011). Because DQS relies on discrete quantum gate sets, one possibility for its verification is that this can in principle be accomplished through fault-tolerant error correction on the gates (Terhal, 2015). This increases the number of operations required, and it has been shown that the run time can scale exponentially with the desired precision of the parameter being computed, in which case the resource cost for fault tolerant DQS for parameter estimation becomes similar to that of universal quantum computing (Brown et al., 2006;Clark et al., 2009). Further, simulations with systems that lack universal gate sets or the digital simulation of open systems may render faulttolerance unavailable in a DQS (Hauke et al., 2012). Since the cost of introducing currently known methods for fault-tolerance is too high for precision DQS (Kendon et al., 2010), other methods for verifying non-fault-tolerant-DQS are needed.
Another way one might test quantum simulators is to compare the results of two simulations. This could involve comparing the results of simulations performed on different platforms (Leibfried, 2010), or even comparing the results obtained by using different simulation methods on the same machine. A variant on this second theme is to run a Hamiltonian simulation forward and then backward in time (Cirac and Zoller, 2012), which may reveal flaws that are not undone by the time reversal, such as dissipation. An initial experiment demonstrated this time-reversal technique for a trapped ion quantum simulation by adiabatically ramping from an initial state of high magnetization along y, through a phase transition, and then back again (Islam et al., 2013). Measurements of the magnetization at all three extrema in this time sequence revealed a revival in the magnetization, achieving an average of S y = 68(4)% of the initial value, in agreement with closed-system numerical integration.
Recently, a variational eigensolver approach has been combined with an ion trap quantum simulator to perform variational quantum simulation (VQS) of the lattice Schwinger model that combines some ways to verify some features of the result (Kokail et al., 2019). The VQS uses feedback with a classical computer that translates measurement results from the AQS into expectation values of a software Hamiltonian to find energy eigenstates of the underlying Hamiltonian. Since the conversion between the measurement and its interpretation happens in classical software, it provides a way to perform some verification of the resulting states because both the eigenvalues and their variances are accessible. For instance, the VQS demonstrated in (Kokail et al., 2019) measured the expectation value of the simulated Hamiltonian E = H , as well as the expectation value of (H -E)
2 , which should be zero if the state is an energy eigenstate of H with eigenvalue E. While this does not guarantee that the state found by the VQS is the ground state, this verification can be used to assess the confidence in the state being an eigenstate.
Quantum simulations with trapped ions can be susceptible to unwanted interactions that lead to inaccuracies in the simulation. Many such error sources are common to both simulations and trapped ion quantum computing gates, such as spontaneous emission from the lasers driving spin transitions, and have been examined in detail elsewhere (Ozeri et al., 2007;Wineland et al., 1998). Further, the simulation protocol itself may have known approximations (such as the Trotterization errors and non-adiabatic evolution) that may be rigorously bounded, though their effects may not be fully understood.
The inclusion of a transverse effective magnetic field to the spin-dependent force of Eq. ( 6) includes higher order terms beyond the simple transverse Ising model of Eq. ( 25) (Wang and Freericks, 2012). These additional terms can create substantial spin-motion entanglement that can affect measurements in bases other than the Ising direction. For transverse field strengths that exceed the Ising coupling, the system begins to attain the character of an XY model with both σ i x σ j x and σ i y σ j y couplings, and the strong-field Ising model can break down. The spins in this case do not strictly decouple from the phonons at any point in time, but it has been shown that it can typically be made small for experimentally accessible timescales (Wall et al., 2017).
Protocols that rely on adiabatic ramping through a small gap can be susceptible to the breakdown of the adiabatic approximation in the region where the gap is small. As we discuss below, for a linear ramp of B(t) through a system energy gap of size ∆, the adiabaticity criterion is approximately | Ḃy /∆ 2 | 1. For cases where the gaps are known, the ramp rate can be adaptively matched to the gap to maximize adibaticity, a technique known as local adiabatic evolution (Roland and Cerf, 2002). However, since repeated experiments can be used to gather statistics about the final state, it has been shown that significant non-adiabaticity can be present and still allow the ground state spin configuration to be found due to its statistical prevalence (Richerme et al., 2013b).
State preparation and measurement (SPAM) errors, which are to some degree common to simulators and gate-based quantum computers, are another source of error in lattice spin simulators. Given an uncorrelated, single-shot, single ion SPAM fidelity F , the probability of a single SPAM error is 1 -F N . The current state of the art SPAM fidelity for single qubits is F = 0.99971 (3), which will produce one error on average for single-shot projective measurements of a register with N ≥ ln 2/(1 -F ) ≈ 2400 qubits (Christensen et al., 2020). In some special cases, repeated experiments can be used to mitigate this error through statistical methods (Shen and Duan, 2012). Cross-talk between neighboring ions can also lead to measurement errors, and the ion positions must be calibrated to define a region of interest on the camera for each ion's fluorescence detection. State detection infidelity from neighboring ions is typically no more than a few percent per qubit (Zhang et al., 2017a), although this error can be suppressed to well below 1% (Cetina et al., 2020) and can be made even smaller with ion-shuttling (Kielpinski et al., 2002).
There are various sources of decoherence quantum simulations in the trapped ion platform, such as stray magnetic and electric fields, mode frequency drifts, off-resonant motional excitation and spontaneous emission. While these errors have been analyzed in the context of quantum computing (Wineland et al., 1998), decoherence can set a practical time limit for quantum simulations that may lead to other errors, such as diabatic errors. Since many of these error sources increase with system size, it may be necessary to employ methods of mitigation such as magnetic field shielding (Ruster et al., 2016).
Off-resonant excitation to motional modes is nominally already considered in the spin-motional coupling described by Eq. ( 12). The probability of motion-induced spin-flip errors in a chain of N ions can be estimated by ε = N i,m |α i,m | 2 , summing over all modes. By neglecting the time dependence and assuming equally-contributing modes, this error is expected to scale as ε ∼ N (ηΩ/δ) 2 , where δ is the smallest detuning from any mode. This same approximation was made in the derivation of Eq. ( 22), leading to the scaling J ij ∼ (ηΩ) 2 /δ. Motion induced spin-flip errors can therefore be mitigated by increasing both δ (Eq. ( 6)) and the Raman Rabi frequency Ω = g 1 g 2 /2∆ to keep the same spin-spin interaction. However, increasing the Rabi frequency may increase spontaneous emission errors, which grow linearly with Rabi frequency for large detuning as Γ = γΩ/4∆, where the atomic linewidth γ is defined in section I.A and ∆ γ is assumed. Given a level of motion-induced spin flip error ε, it can be shown that the spontaneous emission error during the interaction time t J0 = 2π/J 0 scales as N/ε (Kim et al., 2010).
For the Ising spin models considered in this review, it is crucial to validate the strength of the Ising coupling matrix J ij of Eqs. ( 20)-( 22) and effective magnetic fields in Eq. ( 3). For small numbers of spins, it is possible to directly extract the Ising couplings and fields by preparing the spins in a σ z eigenstate and subjecting them to the σ θ Ising interactions or field terms. The resulting oscillations in population are given by the energies of the occupied states, so performing a Fourier transform on these oscillations directly provides the energy differences.
An example of directly-measured Ising oscillations and their resultant Fourier transform is shown in Figs. 5a-b for N = 3 spins. The extracted interaction strengths are shown in Fig. 5c. This technique is effective in contexts of both continuous (Kim et al., 2011) and digital (Lanyon et al., 2011) simulations of Ising models. A similar method has been used to extract the strengths of a magnetically induced spin-spin couplings (Khromova et al., 2012;Piltz et al., 2016). Neither technique scales well to more than a few spins, owing to the difficulty in extracting many closely-spaced frequencies in the spin oscillations.
Individual Ising couplings within a given spin chain can be measured using an auxiliary state of the ions, even for large numbers of spins. Because the Ising couplings depend on the vibrational mode spectrum, all ions must be physically present in the trap to obtain a meaningful result, but the spectrum of the population oscillations illustrated in Fig. 5 would be difficult to obtain if all spins participate in the many-body dynamics. An alternative approach is to perform a separate measurement for each individual Ising coupling, by "hiding" all ions except the pair of interest into an auxiliary internal state that does not experience the spin-dependent force giving rise to Ising couplings. In this manner, the frequency with which the ions i and j of interest oscillate between correlated states |00 and |11 for example, allow the determination of the Ising matrix J ij (Jurcevic et al., 2014).
For large collections of spins with long-range interactions, benchmarking of weak Ising interactions can be accomplished by measuring the global precession of the spins. Here, the spins are each prepared in an identical state that is tipped away from the Bloch sphere axis of the Ising interaction and the resulting dynamics of the global Ising interaction can be recorded. Figure 6 shows measurements of this type of benchmarking in a collection of over 200 trapped ion spins (Britton et al., 2012), agreeing well with mean-field theory. For sufficiently strong Ising interactions, as shown in this experiment, the mean-field approximation breaks down, indicating the entanglement between the spins.
Other spectroscopic techniques for probing the energy spectrum of the bare Ising Hamiltonian are also possible (Senko et al., 2014). For instance, in the tranverse Ising model of Eq. ( 25), modulating the effective field B y (t) at a frequency commensurate with an energy difference in the full spin Hamiltonian will drive transitions between the two differing states. Specifically, taking B y (t) = B 0 + B p sin(ω mod t) with B p << J, the frequency ω mod at which such transitions occur are directly related to the Ising couplings J ij (Senko et al., 2014). Figure 7a-b illustrates examples where the Ising matrix is directly measured using this COM )/ω 2 z = 0, -1, and -2.4, respectively. The red squares (J1) and blue circles (J2) are experimentally measured couplings, and the lines are calculated from Eq. ( 22) with no free fit parameters. The particular measurements in (a) and (b) correspond to the scaled laser beat note detuning (µ 2 -ω 2 COM )/ω 2 z = -1.2 indicated by the points highlighted in the violet rectangle in (c). Adapted from Kim et al., 2009 technique to confirm the validity of power-law approximation described in Eq. ( 23) for a handful of spins. The technique of applying a small oscillating term can be generalized to other cases, for example, measuring the critical (minimum) gap between ground and first excited state of the transverse field Ising model, as shown in Fig. 7c.
Finding ground states of a non-trivial Hamiltonian has tremendous importance in various disciplines across condensed matter physics, quantum chemistry and computer science. In condensed matter physics, the rich phenomena of complex quantum systems can be understood by finding the ground states of the corresponding many-body Hamiltonian (Foulkes et al., 2001;Kohn, 1999;Schollwöck, 2005). In quantum chemistry and molecular physics, the central problems is to determine the electronic structure and the ground-state energy of atoms and molecules (Jensen, 1989). In computer science, the ground state of the complex quantum Hamiltonian can encode other computational problems such as satisfiability or optimization (Albash and Lidar, 2018;Kolsgaard et al., 2001;Lloyd, 2008).
The computational tasks of finding the ground state of non-trivial Hamiltonians are classically demanding because of the exponentially increasing Hilbert space of the Hamiltonian. A quantum simulator is expected to provide a solution beyond the limitations of classical computation. Recently various theoretical schemes for the ground-state problem have been proposed and proof-of principle experimental demonstrations have been performed, including adiabatic preparation (Edwards et al., 2010;Friedenauer et al., 2008;Islam et al., 2011Islam et al., , 2013;;Kim et al., 2011Kim et al., , 2010;;Richerme et al., 2013a,b;Schneider et al., 2012), direct cooling by bath-engineering (Barreiro et al., 2011;Lin et al., 2013), and algorithmic cooling schemes (Baugh et al., 2005;Xu et al., 2014;Zhang et al., 2020). For the case of adiabatic method, it has been shown to be closely related to adiabatic quantum computation, which is proved to be equivalent to a universal quantum computer (Albash and Lidar, 2018).
We focus on the adiabatic preparation of the ground state of quantum spin models with trapped atomic ion spins, including a description of the general scheme of the experimental procedure and various adiabatic ramping protocols. Following a discussion of the adiabatic protocol applied to varying numbers of trapped ion spins, we consider how this protocol can be optimized and applied to broader classes of spin models, and briefly discuss the case of a spin-1 system. These quantum spin models clearly show the essences of the adiabatic quantum simulation with wide applications. Moreover, these quantum spin models can describe a large class of many-body quantum physics in condensed matter such as quantum magnetism (Moessner and Ramirez, 2006), spin glasses (Binder and Young, 1986), and spin liquids (Balents, 2010). The solutions of certain spin Hamiltonians are also connected to many other computational problems including optimization problems when the system is extended to 2 dimensions (Albash and Lidar, 2018).
] FIG. 6 Mean-field benchmarking of Ising couplings in a 2D crystal of 206 (10) ions confined in a Penning trap. The spins are all initialized in a separable product state slightly tipped from the axis of a subsequently applied (weak) long-range Ising interaction. The resulting global spin-precession of the ions is measured as a function of the detuning of the optical dipole force laser beams from the axial center-of-mass (COM) mode (see Sec. I and Fig. 3b). Each point is generated by measuring the mean Ising coupling at a particular the laser beam intensity I. The solid line (red) is the prediction of mean-field theory that accounts for couplings to all N transverse modes, with no free parameters.
The lines breadth reflects experimental uncertainty in the initial tipping angle of the spins. The mean-field prediction for the average value of the power-law exponent α from Eq. ( 23), is drawn in green (right axis, linear scale). For stronger Ising interactions, the mean-field approach breaks down. Adapted from Britton et al., 2012.
Adiabatic ground state preparation is analogous to that of adiabatic quantum computation (Albash and Lidar, 2018;Farhi et al., 2000): a quantum system is initialized to the ground state of a trivial Hamiltonian H triv . Next, the Hamiltonian is adiabatically deformed into the Hamiltonian of interest H prob , whose ground state encodes the solution of a problem that has been mapped to this final Hamiltonian. The adiabatic evolution is generated by
where s = s(t) is a time dependent parameter changing from 0 to 1 during the time interval from t = 0 to t = t f . In the context of trapped ion spin models, the trivial spin Hamiltonian can be an effective magnetic field as described by Eq. ( 3) and the Hamiltonian of interest can be a fully-connected transverse Ising model of Eq. ( 25). The determination of ground states of the long range transverse field Ising model cannot always be predicted, even with just a few dozen spins (Sandvik, 2010).
Although the fidelity of remaining in the ground state of Eq. ( 27) can always be improved by evolving more slowly, a practical upper limit on the transition time is enforced by the finite coherence time of the chosen experimental platform. Given a fixed transition time, it is possible to further optimize the preparation fidelity by adjusting the transition rate based on the local energy gap to the nearest excited state (Richerme et al., 2013b). Such "local adiabatic evolution" can be used for improved preparation and determination of many-body ground states in a trapped-ion quantum simulator. Compared with other adiabatic methods, local adiabatic evolution (Roland and Cerf, 2002) yields the highest probability of maintaining the ground state in a system that is made to evolve from an initial Hamiltonian to the Hamiltonian of interest. Compared with optimal control methods (Khaneja et al., 2005;Krotov, 1996), local adiabatic evolution may require knowledge of only the lowest ∼ N eigenstates of the Hamiltonian rather than all 2 N values. These methods have been used in both linear Paul traps (Richerme et al., 2013b) and Penning traps (Safavi-Naini et al., 2018) to demonstrate optimized ground state preparation as well as a method to find the ground state spin ordering, even when the evolution is non-adiabatic.
For example, to find the ground state of a fully-connected Ising Hamiltonian in Eq. ( 25) via an adiabatic protocol, the spins can be initialized to point along the transverse magnetic field direction with B y Max(J ij ). This initial state is, to good approximation, the instantaneous ground state of the full Eq. ( 25). After initialization, the (time-dependent) transverse field B y (t) can then be ramped adiabatically from B y (t = 0) = B 0 to B y (t = t f ) = 0, ensuring that the system remains in its instantaneous ground state during its evolution, as depicted in Fig. 8. At the conclusion of the ramp, the ground state spin ordering of the Ising Hamiltonian [first term in Eq. ( 25)] may be either directly read out or used as a starting point for further experiments. 25) for N = 6, with the ground state energy Eg set to 0, B0 = 5Jmax, and the long-range Jij couplings determined from experimental conditions (see text). Indicated in bold red is the first coupled excited state, the minimum of which determines the critical field Bc and the critical gap ∆c. From Richerme et al., 2013b.
Fig. 9 shows the energy level spectrum for the Hamiltonian in Eq. ( 25) for N = 6 spins. Since the Hamiltonian obeys Z 2 symmetry (Sachdev, 2011) (as well as parity symmetry in the experiments), the ground state |g is coupled to only a subset of the excited energy eigenstates. The first coupled excited state, shown in red in Fig. 9, is the lowest energy excited state |e for which e| σ y |g = 0. This state displays a general property seen in most adiabatic quantum simulations -namely, the existence of a critical gap ∆ c that is central to parameterizing the adiabaticity of a given ramp. Many different ramp profiles allow one to transform from the initial Hamiltonian to the Ising Hamiltonian, each with different implications for adiabaticity and ground state preparation. Three possibilities are discussed below.
Linear Ramps. For a linear ramp, the time-dependent transverse field B y in Eq. ( 25) takes the form B lin y (t) = B 0 (1t/t f ), with a ramp profile shown in Fig. 10a. To determine whether such a ramp is adiabatic or not, it can be compared to the adiabatic criterion (Messiah, 1962)
where Ḃy (t) is the rate at which the transverse field is changed and = Max[ e| dH/dB y |g ] is a number of order unity that parametrizes the coupling strength between the ground state |g and the first coupled excited state |e . Eq. ( 28) highlights that fast ramps and small critical gaps can greatly decrease adiabaticity.
To satisfy the adiabatic criterion, a linear ramp must proceed slowly enough so that the total time t f B 0 /∆ 2 c . For the N = 6 Ising Hamiltonian shown in Fig. 9, B 0 = 3.9 kHz and ∆ c = 0.29 kHz, giving the adiabaticity requirement t f 46 ms. This time is long compared with the typical coherence time of ion trap quantum simulation experiments. It is therefore desirable to seek alternative ways to decrease B(t) more quickly while maintaining adiabaticity.
Exponential Ramps. Decreasing the transverse field exponentially according to B exp y (t) = B 0 exp(-t/τ ), with t f = 6τ , can yield a significantly more adiabatic evolution than linear ramps for the same t f . Figure 9 shows that the instantaneous gap ∆ between the ground and first coupled excited state is large at the beginning of the ramp and small only when B approaches 0. Exponential ramps exploit this gap structure by quickly changing the field at first, then gradually slowing the rate of change as t → t f .Such ramps have been used to produce ground states in several of the previously discussed experiments, such as (Kim et al., 2010), (Islam et al., 2011) and (Islam et al., 2013).
At the critical point of the Hamiltonian shown in Fig. 9, | Ḃexp (t)| = 0.3B 0 /t f . The adiabaticity criterion of Eq. ( 28) then requires t f 14.5 ms, a factor of 3 less time than the requirement found for linear evolution. Note that the adiabaticity gains of exponential ramps can be realized whenever the critical gap occurs towards the end of the ramp (B c /B 0 < τ /t f ), which is generally the case for the transverse Ising Hamiltonian of Eq. (25).
Local Adiabatic Ramps. Local adiabatic ramps seek to keep the adiabaticity fixed at all points along the evolution by adjusting Ḃy (t) based on the instantaneous gap ∆(B y (t)) that varies with the field profile B y (t) (Quan and Zurek, 2010;Roland and Cerf, 2002). If the adiabaticity parameter is defined as then a local adiabatic ramp would follow the profile B y (t) that solves the differential equation 29 with γ fixed. Adiabaticity then requires γ 1. To solve Eq. ( 29), it is necessary to know ∆(t) everywhere along the evolution. This requires knowledge of the first coupled excited state of the N -spin Hamiltonian of Eq. ( 25), which is always the first excited state at B y = 0 and the (N + 1) st excited state at large B y . Determining the local adiabatic evolution profile therefore relies on calculation of only the lowest ∼ N eigenvalues, which is much more computationally approachable than direct diagonalization of a 2 N × 2 N matrix (Lanczos, 1950).
For a local adiabatic ramp, the total evolution time t f may be calculated by integrating Eq. ( 29). Since Ḃy (t) is negative throughout the evolution,
which shows a linear relationship between the total time t f and the adiabaticity parameter γ. Satisfying the adiabaticity condition γ 1 for the Hamiltonian in Fig. 9 implies t f 3.6 ms, a factor of 4 and 12 less time than exponential and linear ramps, respectively. The fact that local adiabatic evolution can lead to faster ramps while satisfying adiabaticity has been well-explored in (Roland and Cerf, 2002), where it was shown that local adiabatic ramps could recover the quadratic speedup of Grover's quantum search algorithm (Nielsen and Chuang, 2000). In contrast, it was found that linear ramps offer no improvement over classical search algorithms (Farhi et al., 2000).
Fig. 10a compares a linear, exponential, and local adiabatic ramp profile for the Hamiltonian shown in Fig. 9. The local adiabatic ramp spends much of its time evolution in the vicinity of the critical point, since the transverse field changes slowly on account of the small instantaneous gap. This is further illustrated in Fig. 10b, which shows that at the critical point, the slope of the local adiabatic ramp is minimized and smaller than slopes of the exponential or linear ramps. As a result, the inverse adiabaticity 1/γ is peaked near the critical point for exponential and linear ramps, greatly increasing the probability of non-adiabatic transitions away from the ground state (see Fig. 10c). By design, the local adiabatic ramp maintains constant adiabaticity for all values of B and does not suffer from large non-adiabaticities near B c .
To determine the effects of a chosen adiabatic ramp protocol, the probability of creating the ground state can be measured following any chosen ramp profile of identical times. (Richerme et al., 2013b) used N = 6 ions and chose the trap voltages and the laser detuning µ to give AFM spin-spin interactions of the form J ij ≈ (0.77 kHz)/|i -j|. These long-range AFM interactions lead to a fully-connected, frustrated system as all couplings cannot be simultaneously satisfied. Nevertheless, the ground state of this system reduces to an equal superposition of the two Néel-ordered AFM states, (|↓↑↓↑↓↑ + |↑↓↑↓↑↓ )/ √ 2. The data in Fig. 11 show how the AFM ground state probability grows during a single 2.4 ms linear, exponential, or local adiabatic ramp. Each data point is the result of 4000 repetitions of the same experiment, with error bars that account for statistical uncertainty as well as estimated drifts in the Ising coupling strengths. In agreement with the arguments above, the data show that local adiabatic ramps prepare the ground state with higher fidelity than exponential or linear ramps. The ground state population grows quickly under local adiabatic evolution since the transverse field B(t) is reduced quickly at first. In contrast, the linear ramp does not approach the paramagnetic to AFM phase transition until ∼ 2 ms, and the AFM probability is suppressed until this time.
The solid lines in Fig. 11 plot the theoretical prediction of the ground state probability with no free parameters. In each case, the Schrödinger equation is numerically integrated using Hamiltonian (25), the desired B(t), and the initial state |ψ(0) = |↓↓↓ . . . y . At the end of the ramp, the overlap between the final state |ψ(t f ) and the AFM ground state (|↓↑↓ . . . + |↑↓↑ . . . )/ √ 2 is calculated to extract the probability of the ground state spin configuration. Effects of decoherence-induced decay in the ground state probability are included by multiplying the calculated probability at time t by exp[-t/t d ], where t d is the measured 1/e coherence time of the spin-spin interactions.
The key to all adiabatic protocols is that the ramp rate must remain slow when compared to the critical energy gap, as shown above in Eq. ( 28). However, determining the scaling of the critical gap with system size is itself a notoriously difficult problem in the general case (Albash and Lidar, 2018). For simple problems, such as the ground state of a Lipkin-Meshkov-Glick model, the gap is known to shrink only polynomially for large-N (Caneva et al., 2008). For more complex problems, such as those in the N P -complexity class, the gap may close exponentially quickly with increasing system size (Das and Chakrabarti, 2008;Morita and Nishimori, 2008). In these cases, or in cases for which the gap scaling is unknown, experimental implementations will require exponentially longer ramp times to ensure that adiabaticity is maintained as the problem size grows.
One obvious question is whether quantum adiabatic protocols are useful for solving complex computation or quantum manybody type problems. There are some reasons to be hopeful. For certain problems, such as the Grover search algorithm, the use of local adiabatic ramp profiles have been shown to provide the same quadratic quantum advantage as found in the circuit model (Roland and Cerf, 2002). For other problems, such as finding the ground state of a transverse-field Ising model, adiabatic quantum protocols can provide a polynomial speedup over simulated annealing (Kadowaki and Nishimori, 1998). Although exponential speedups are more elusive, finding such examples is nevertheless an active area of research.
While adiabatic simulation protocol allows for the preparation of the ground state of non-trivial spin models, maintaining the adiabatic condition (Eq.( 28)) for a large system within the constraint of an experimentally realistic coherence time will be challenging. Alternate protocols have been explored to bypass the strict requirements of adiabaticity, while achieving high ground state probability. For example, a 'Bang-bang' control of the Hamiltonian has been suggested (Balasubramanian et al., 2018;Viola and Lloyd, 1998), where the initial trivial Hamiltonian can be quenched to an intermediate Hamiltonian, followed by a final quench to the problem Hamiltonian. In another approach, a classical-quantum hybrid protocol (the Quantum Approximate Optimization Algorithm, or QAOA (Farhi et al., 2014)) theoretically enables ultra-fast creation of ground states (Ho et al., 2019). Implementations of this method are discussed in detail in section V.B.2. Here, we restrict our discussions to adiabatic simulation protocols.
The adiabatic preparation of the ground state for the transverse Ising model was first demonstrated with two trapped ion spins (Friedenauer et al., 2008;Schneider et al., 2012), followed by experiments with three spins (Edwards et al., 2010;Khromova et al., 2012;Kim et al., 2009Kim et al., , 2011Kim et al., , 2010)). These entry experiments demonstrated the adiabatic evolution from para-magnetic initial state to magnetically ordered ground states and allowed tests of adiabaticity (Edwards et al., 2010) and direct measures of entanglement in the ground state (Kim et al., 2011(Kim et al., , 2010)). The three-spin system moreover supports spin frustration, or a competition between the nearest and next-nearest couplings in the case of antiferromagnetic (AFM) ground states. By tuning the system to have either ferromagnetic (FM) and anti-ferromagnetic (AFM) ground states, two different types of magnetic order were indeed measured, paralleling the two different classes of entanglement known to exist with exactly three spins (Acín et al., 2001). For the case of three Ising spins, the transverse Ising Hamiltonian ( 25) is reduced to
x σ (2) x + σ (2) x σ (3) x ) + J 2 σ (3) x σ (1) x + B y (t)(σ (1) y + σ (2)
where the transverse field and the Ising interaction are chosen to act along the y-axis and x-axis, respectively. This is the simplest Hamiltonian that can exhibit frustration in the ground state due to a compromise between the various Ising couplings.
As seen in Eq. ( 22), the sign and the strength of the Ising couplings J ij can be controlled by the proper choice of the driving field detuning δ m from the motional modes. For three spins, the expected and measured nearest-neighbor (NN) interactions, J 1 ≡ J 1,2 = J 2,3 and the next-nearest-neighbor (NNN) interaction J 2 ≡ J 1,3 are shown in Fig. 5. For certain ranges of the drive field detuning, both NN and NNN couplings have anti-ferromagnetic (AFM) interactions (J 1 , J 2 > 0), and for other domains both show ferromagnetic (FM) interactions (J 1 , J 2 < 0).
Figure 13a shows the time evolution for the Hamiltonian frustrated with nearly uniform AFM couplings and gives almost equal probabilities for the six AFM states |↓↓↑ , |↑↓↓ , |↓↑↑ , |↑↑↓ , |↓↑↓ , and |↑↓↑ (labeled in the x-basis of the Bloch sphere and accounting for 3/4 of all possible spin states) at B y ≈ 0. Because J 2 < 0.8J 1 for this data, a population imbalance also develops between symmetric (|↓↑↓ and |↑↓↑ ) and asymmetric (|↓↓↑ , |↑↓↓ , |↓↑↑ , and |↑↑↓ ) AFM states. Figure 13b shows the evolution to the two ferromagnetic states (|↓↓↓ and |↑↑↑ ) as B y → 0, where all interactions are FM.
The adiabatic evolution of the ground state of Hamiltonian (31) from B y J rms to B y J rms should result in an equal superposition of all ground states and therefore be entangled. For instance, for the isotropic AFM case, the ground state is expected to be |↓↓↑ + |↑↓↓ + |↓↑↑ -|↑↑↓ -|↓↑↓ -|↑↓↑ . For the FM case, the ground state is a GHZ as |↓↓↓ -|↑↑↑ . The entanglement in the system at each point in the adiabatic evolution can be characterized by measuring particular entanglement witness operators (Gühne and Tóth, 2009). When the expectation value of such an operator is negative, this indicates entanglement of a particular type defined by the witness operator. For the AFM (frustrated) case as shown in Fig. 13c, the expectation of the symmetric W state witness, (Gühne and Tóth, 2009). For the FM case as shown in Fig. 13d, the expectation of the symmetric GHZ witness operator
y σ
(2) y σ
(3) y (Gühne and Tóth, 2009;Sackett et al., 2000) is measured, where Î is the identity operator and Ĵi ≡ 1 2 (σ
l ) is proportional to the lth projection of the total effective angular momentum of the three spins. In both cases, as shown in Figs. 13c and(d), entanglement of the corresponding form is clearly observed during the adiabatic evolution.
The ground state in the transverse field Ising model, Eq. ( 25) undergoes a crossover between polarized/paramagnetic and magnetically ordered spin states, as the relative strengths of the transverse field B y and the Ising interactions J ij are varied. For |B y /J ij | 1, the ground state has the spins independently polarized (paramagnetic phase). For |B y /J ij | 1, the ground state is magnetically ordered for |B y /J ij | 1, with ferromagnetic order for J ij < 0 in Eq. ( 25). A second order quantum phase transition is predicted for this model in the thermodynamic limit (Sachdev, 2011) when the magnitude of the transverse field is comparable to the interaction strength.
The presence or absence of a spin-order can be quantified by adopting a suitable order parameter. For example, the average absolute magnetization per site along the Ising direction,
differentiates between a ferromagnetic state and paramagnetic state. Here, P (s) is the probability of finding s spins in the |↑ state
To remove a finite size effect due to the difference between Binomial and Gaussian distributions, a scaled order parameter, mx = (m 0 x,Nm x )/(m 0 x,N -1) can be adopted. Here, m 0
is the average absolute magnetization of the paramagnetic state. This scaled order parameter assumes a value of mx = 1 for the ideal ferromagnetic state while mx ≈ 0 in the paramagnetic state. A finite system does not support a phase transition, but shows a smooth crossover from the paramagnetic to the spin-ordered phases that becomes sharper as the system size is increased. Higher order moments of the distribution of measured spins may be more suitable to extract the phase transition point from experiments performed on finite system sizes. For instance, the fourth magnetization moment is known as the Binder cumulant,
(33)
The Binder cumulant can also be scaled to remove the finite size effect, as before, by defining ḡ = (g 0 Ng)/(g 0 N -1), where g 0 N = 3 -2/N is the Binder cumulant for the paramagnetic phase.
Fig. 14 shows measurements of the mean magnetization and Binder cumulanant in a transverse Ising system ranging from N = 2 to N = 9 atomic ion spins (Islam et al., 2011). The observed sharpening of the crossover from paramagnetic to ferromagnetic spin-order with system size (Fig. 14) is consistent with an onset of the quantum phase transition. In Fig. 14a, theoretical values of both order parameters are shown for up to N = 100 spins in an all-to-all coupled ferromagnetic transverse Ising model, with measurements and comparison to theory in Fig. 14b-d. Both metrics have been scaled to take into account finite size effects (Islam et al., 2011). Ferromagnetic spin order was also observed in adiabatic quantum simulation experiments with up to N = 16 ions by directly measuring a bimodal distribution of magnetization (Islam et al., 2013).
AFM ground states of the transverse field Ising model (Eq. ( 25) with J ij > 0) are more difficult to prepare because the longranged AFM interactions lead to competing pairwise spin order, or frustration, as detailed with N = 3 spins (Kim et al., 2010) in section III.B.1. Intuitively, the longer the range of AFM interactions, the less energy it takes to create spin-flip excitations. Thus the critical field required to destroy the AFM spin order is less than that with a relatively shorter range interaction. The 'critical gap' in the many-body energy spectra also decreases with increasing range of the AFM couplings (Fig. 15). The reduction of the critical gap with increased range of interaction was experimentally probed in a quantum simulation of the transverse field Ising model, Eq.25 with the interaction profile following an approximate AFM power law, Eq. ( 23) for N = 10 spins (Islam et al., 2013). The ratio of the transverse field to the Ising couplings was varied quasi-adiabatically from a high transverse field to a final value of B/J 0 = 0.01. As the interaction range was increased and the critical gap closed, more excitations were created, resulting in a reduction in the ground state order. This was observed through a decrease in the measured structure function,
where the average correlation of spins along the Ising direction x and separated by r sites is
The measured structure function for N = 10 spins (Fig. 15a) shows a clear decline at k = π as the interaction range is made longer, thus quantifying the degradations of nearest-neighbor antiferromagnetic spin order as the ground state gap shrinks.
Interestingly, the ground state spin ordering may still be determined experimentally even when the ramp is non-adiabatic. The key to ground state identification is to examine the probability distribution of all spin configurations at the conclusion of the ramp and select the most prevalent state in the final eigenbasis. Consider an experiment where the spins are initialized into |↓↓↓ . . . y (as usual) and the transverse field B(t) is instantly switched from B = B 0 to B = 0. Measurement along the x-direction would yield an equal superposition of all spin states; in this instance, the ground state is just as probable as any other state. If the transverse field B(t) is instead ramped at a fast but finite rate, the quantum simulation is slightly more adiabatic than the instantaneous case, and the ground state becomes slightly more prevalent than any other state. When B(t) is ramped slowly enough, the ground state population is nearly 100% and dominates over that of any other state.
A close analogy may be drawn with a Landau-Zener process (Zener, 1932) in a two-level system comprised of the ground and first coupled excited states. Adiabatic ramps correspond to half of a Landau-Zener process, in which B(t) starts with B J and ends at B = 0. One can write an analytic expression to calculate the transition probability for this half-Landau-Zener evolution (Damski and Zurek, 2006), which has a maximum value of 0.5 for an instantaneous ramp. Any fast but finite ramp will give a transition probability < 0.5, so the ground state will always be more prevalent than the excited state.
The technique of identifying the most prevalent state as the ground state is subject to some limitations. First, the initial state (before the ramp) should be a uniform superposition of all spin states in the measurement basis -a condition satisfied by preparing the state |↓↓↓ . . . y and measuring along x. If some spin states are more prevalent than the ground state initially, then some non-zero ramp time will be necessary before the ground state probabilities "catch up" and surpass these initially prevalent states. Second, the ramp must not cross any first-order transitions between ordered phases, as non-adiabatic ramps may not allow sufficient evolution time towards the new ground state order. In addition, the initial and final states must share the same FIG. 15 (a) Structure function S(k) for various ranges of AFM interactions, for B/J0 = 0.01 in a system of N = 10 spins. The increased level of frustration for the longer-range interactions reduces the observed antiferromagnetic spin order. The detection errors may be larger than shown here for the longest range of interactions, owing to spatial crosstalk from their closer spacing. (b) Distribution of observed states in the spin system, sorted according to their energy Ei (with E0 denoting the ground state energy) calculated exactly from Eq. ( 25) with B = 0. Data are presented for two ranges (red for α = 1.05 and blue for α = 0.76). The dashed lines indicate the cumulative energy distribution functions for these two ranges. Adapted from (Islam et al., 2013). symmetry properties.
Finally, a good determination of the ground state requires that the difference between the measured ground state probability P g and next excited state probability P e be large compared with the experimental uncertainty, which is fundamentally limited by quantum projection noise ∼ 1/
√ n after n repetitions of the experiment (Itano et al., 1993). This implies that the most prevalent ground state can be determined reliably after repeating the measurement n > (P 2 g + P 2 e )/(P g -P e ) 2 times. Assuming an exponential distribution of populated states during the ramp (as may be expected from Landau-Zener-like transitions), the number of required runs should then scale as n ∼ ( Ē/∆) 2 in the limit Ē ∆, where Ē is the mean energy imparted to the spins during the ramp, and ∆ is the energy splitting between the ground and first coupled excited state.
If the gap shrinks exponentially with the number of spins N (i.e. ∆ ∼ e -N ), ground state identification thus requires an exponential number of measurements n in the simulation. However, in cases where the gap shrinks like a power law (∆ ∼ N -α ), the most prevalent state can be ascertained in a time that scales polynomially with the number of spins. Regardless of the scaling, techniques that improve the ground state probability (such as local adiabatic evolution) can greatly increase the contrast of the most prevalent state and reduce the number of necessary repetitions. Fig. 16 shows a direct identification of the ground state AFM order of N = 14 trapped ion spins by imaging the most prevalent state created after a ramp. Each box in Fig. 16a contains an ion that scatters many photons when in the state |↑ and essentially no photons when in the state |↓ .
demonstrates the resiliency of most-prevalent state selection to ramps that are far from adiabatic. Identification of the ground state is clear, even though the total ground state probability is only ∼ 3%. The requirement of satisfying the adiabatic criterion is replaced only by the requirement that the most prevalent state probabilities are accurately resolvable compared with those of any other states. While the method should remain robust for even larger N , ramps that are more adiabatic (by using longer ramp times or stronger spin-spin couplings) will decrease the number of experimental repetitions needed to clearly resolve the state probabilities.
Adiabatic protocols can also be used to create the ground states of a classical spin model, catalyzed by quantum fluctuations. Consider a system described by the transverse Ising model of Eq. 25 accompanied by a longitudinal field B x :
When the transverse field B y is set equal to 0, and the longitudinal field B x is varied, this Hamiltonian exhibits many distinct ground state phases separated by first-order classical phase transitions. Yet even for just a few spins, the various ground states at different B x are classically inaccessible in a physical system at or near zero temperature due to the absence of thermal fluctuations to drive the phase transitions (Sachdev, 2011). Quantum fluctuations are therefore required to reach these various ground states in a physical system. Such quantum fluctuations can be introduced to the system by applying a transverse magnetic field, which does not commute with the longitudinal-field Ising Hamiltonian. Using N = 6 or N = 10 spins, this technique has been used to experimentally identify the locations of the multiple classical phase transitions and to preferentially populate each of the classical ground states that arise for varying strengths of the longitudinal field (Richerme et al., 2013a). The ground state spin ordering reveals a Wigner-crystal spin structure (Wigner, 1934) that maps to particular energy minimization problems (Katayama and Narihisa, 2001) and shows the first steps of the complete "devil's staircase" (Bak and Bruinsma, 1982) expected to emerge in the N → ∞ limit.
Fig. 17a shows the energy eigenvalues of the Hamiltonian given in Eq. ( 36) with B y = 0 for a system of 6 spins. The ground state passes through three level crossings as B x is increased from 0, indicating three classical first-order phase transitions separating four distinct spin phases. For each B x , there is a critical point at some finite B y characterized by a critical gap ∆ c (inset of Fig. 17b). When B x is near a classical phase transition, the near energy-degeneracy of spin orderings shrinks the critical gap, as shown in Fig. 17b.
Long-range interactions give rise to many more ground state spin phases than does a local Ising model. Consider an N -spin nearest-neighbor AFM model (Ising coupling J) and a ground-state ordering |.. ↓↑↓↑↓↑ .. . An excited state at longitudinal field B x = 0 may have an additional spin polarized along |↓ , either by making a kink of type |.. ↓↑↓↓↑↓ .. or a spin defect of type |.. ↓↑↓↓↓↑ .. . The interaction energy gain of making n kinks is 2nJ, while the field energy loss is 2nB x . At B x /J = 1, multiple energy levels intersect to give a first-order phase transition. Similarly, the energy gain of making n spin defects is 4nJ and the loss is 2nB x , so a second phase transition occurs at B x /J = 2. Only three different ground state spin phases are observable as B x is varied from 0 → ∞, independent of N , and there is a large degeneracy of spin eigenstates at the phase transitions. The presence of long-range interactions lifts this degeneracy and admits
To create these various spin phases, the experiment from (Richerme et al., 2013a) begins by optically pumping the effective spins to the state |↓↓↓ .. z . The spins are then coherently rotated into the equatorial plane of the Bloch sphere so that they point along B = B x x + B y (0)ŷ, with B x varied between different simulations. The Hamiltonian of Eq. ( 36) is then switched on at t = 0 with the chosen value of B x and B y (0) = 5J max . The transverse field (which provides the quantum fluctuations) is ramped down to B y ≈ 0 exponentially with a time constant of 600 µs and a total time of 3 ms, which sacrifices adiabaticity in order to avoid decoherence effects. At t = 3 ms, the Hamiltonian is switched off and the x-component of each spin is measured by applying a global π/2 rotation about the ŷ axis, illuminating the ions with resonant light, and imaging the spindependent fluorescence using an intensified CCD camera. Experiments are repeated 4000 times to determine the probability of each possible spin configuration.
The order parameter of net magnetization along x, M x = N ↑x -N ↓x , can then be investigated as a function of longitudinal field strength. The magnetization of the ground state spin ordering of Eq. ( 36) is expected to yield a staircase with sharp steps at the phase transitions (red line in Fig. 18a) when B y = 0 (Bak and Bruinsma, 1982). The experimental data (blue points in Fig. 18a) show an averaged magnetization with heavily broadened steps due largely to the non-adiabatic exponential ramp of the transverse field. The deviation from sharp staircase-like behavior is predicted by numerical simulations (solid blue line in Fig. 18a) which account for the implemented experimental parameters and ramp profiles. Differences between theory and experiment are largest near the phase transitions, where excitations are easier to make due to the shrinking critical gap (Fig. 17b).
The ground state spin configuration at each value of B x can be extracted by looking at the probability distribution of all spin states and selecting the most prevalent state (inset of Fig. 18a) (Richerme et al., 2013b). The magnetization of the spin states found by this method (black points in Fig. 18a) recover the predicted staircase structure. The steps in the experimental curve agree with the calculated phase transition locations to within experimental error (gray bands in Fig. 18a), which accounts for statistical uncertainty due to quantum projection noise and estimated drifts in the strengths of J ij , B x , and B y . Fig. 18b shows approximately 1000 averaged camera images of the most probable spin configuration observed at each plateau in Fig. 18a. Each box contains an ion that scatters many photons when in the state |↑ and essentially no photons when in the state |↓ along the x-axis of the Bloch sphere. The observed spin orderings in Fig. 18b match the calculated ground states at each magnetization, validating the technique of using quantum fluctuations to preferentially create these classically inaccessible ground states. (For magnetizations M x = 0 and M x = -4, two ground state orderings are observed due to the left-right symmetry of the spin-spin interactions.)
To further illustrate the necessity of using quantum fluctuations to catalyze the magnetic phase transitions, alternate ramp trajectories can be used to reach a final chosen value of B x . Figure 19a shows the ground state phase diagram of the Hamiltonian of Eq. ( 36), with the sharp classical phase transitions visible along the bottom axis (B y /J max = 0). In addition, it shows two possible trajectories through the phase diagram that start in a paramagnetic ground state (which is easy to prepare experimentally) and end at the same value of B x with B y = 0.
The first trajectory, in which B x is fixed and B y is ramped from 5J max to 0, was the one used in Fig. 18 to experimentally verify the locations of the 3 classical phase transitions and to experimentally create the 4 different ground state phases. Along this trajectory, Fig. 19b plots the probability of creating each ground state as a function of B x and find populations of ∼ 40 -80%. A smooth crossover between the four ground state phases was observed, with the classical phase transitions occurring at the crossing points. This arises since distinct spin eigenstates have degenerate energies at the phase transition, causing the critical gap between them to close and allowing quantum fluctuations to populate both states equally (see Fig. 17).
The second trajectory in Fig. 19a is purely classical, with B y set to 0. The spins are initialized into the state |↓↓↓↓↓↓ along x, and B x is ramped from 5J max to its final value at a rate of 5J max /3 ms. Figure 19c shows that in a classical system without thermal or quantum fluctuations, the phase transitions remain undriven and the initial state |↓↓↓↓↓↓ remains dominant for all values of B x . The ground state phases with magnetization 0 and -2 (blue and green in Fig. 19c) are separated from the initial state by several classical phase transitions and have essentially zero probability of being created.
As with the spin-1/2 systems described above, spin-1 systems (spanned by the three basis kets |+ , |0 , and |-) can likewise exhibit a variety of interesting new physics and ground-state phases. As a notable example, Haldane conjectured (Haldane, 1983) that integer-spin Heisenberg chains with nearest-neighbor AFM interactions are gapped, in contrast to gapless half-integer showing Raman beatnotes in relation to Zeeman splittings and motional sidebands for the center-of-mass mode. Level splittings are not to scale. From (Senko et al., 2015) spin chains. This energy gap in integer-spin systems corresponds to short-range exponentially decaying correlation functions, as opposed to long-range power-law decaying correlations in half-integer systems. It was later suggested (den Nijs and Rommelse, 1989) that this Haldane phase of the spin-one chain is governed by a hidden order, which can be characterized by a non-local string order parameter and is consistent with a full breaking of a hidden Z 2 × Z 2 symmetry (Kennedy and Tasaki, 1992). The Haldane phase can also be described by a doubly-degenerate entanglement spectrum (Pollmann et al., 2010), hinting at a topologically protected phase in one-dimension.
Effective spin-1 particles can be represented by three hyperfine levels in the 2 S 1/2 ground manifold of 171 Yb + ions: Spin-1 Ising couplings can be generated analogously to the spin-1/2 case detailed in Section I. Laser beams are applied to the ion chain with a wavevector difference along a principal axis of transverse motion, but here driving stimulated Raman transitions between both the |0 and |states and the |0 and |+ states with balanced Rabi frequencies Ω i on ion i (Kim et al., 2009). To generate spin-1 XY interactions, two beat frequencies are applied at ω -+ µ and ω +µ to these respective transitions, where µω m = δ m is the detuning from the transverse motional mode m sideband, as shown in Fig. 20c. Under the approximations that the beatnotes are far detuned (|δ m | η i,m Ω i ) and that ω ± µ Ω i (the rotating-wave approximation), the resulting interaction Hamiltonian in the Lamb-Dicke regime is (Senko et al., 2015)
where S i ± are the spin-1 raising and lowering operators. The pure "XY" or "flip-flop" spin-spin interaction in the first term of Eq. ( 37) follows the same formula as for generating spin-1/2 Ising interactions in Eq. ( 22) (Kim et al., 2009), As in the case for spin-1/2, when δ m > 0 is larger than transverse mode frequencies, J ij falls off with distance approximately with a power law form from Eq. ( 23), where nearest-neighbor Ising coupling J 0 is typically of order ≈ 1 kHz and α can be tuned between 0 and 3 as discussed above (Islam et al., 2013;Porras and Cirac, 2004).
The second term in Eq. ( 37) represents additional spin-phonon terms and is parameterized by the factor V i,m = (η i,m Ω i ) 2 /(8δ m ). For long-range spin-spin interactions with α 0.5, or for small numbers of ions, the V i,m terms are approximately uniform across the spin chain. In these instances, the V i,m coefficient can be factored out of the sum over ions in Eq. ( 37), leaving only global S i z and (S i z ) 2 terms. For shorter-range interactions or for longer chain lengths, the V i,m terms can be eliminated by adding an additional set of beat frequencies at ω -µ and ω + + µ, which would generate Ising-type interactions between effective spin-1 particles using the Mølmer-Sørensen gate (Mølmer and Sørensen, 1999).
As theoretically proposed in both (Cohen et al., 2015) and (Gong et al., 2016), a range of spin-1 phenomena discussed above can be accessed in trapped ion quantum spin simulators. In (Cohen et al., 2015), for instance, it is shown how to generate the full spin-1 XXZ Hamiltonian
where the S i γ terms are the spin-1 Pauli operators on site i along the γ direction, λ is the ZZ Ising anisotropy, and D is analogous to a magnetic field B term of Eq. ( 3) for spin-1/2 systems. The interacting terms in Eq. ( 39) arise from a generalization of the Mølmer-Sørensen gate (Mølmer and Sørensen, 1999) to spin-1 systems, followed by a transformation to the interaction picture; the on-site D-term can be generated by imposing frequency detunings D on all the previous driving fields.
Generating the ground state of the Haldane phase (Haldane, 1983) can be realized by an adiabatic ramp procedure (Cohen et al., 2015). To begin, the spin-1 system can be initialized into a product state of |0 on each site, which is the trivial ground state when D J. Adiabatically reducing D will then drive the system towards the Haldane phase. As the system size increases and the critical gap between the D and Haldane phase closes, a symmetry-breaking perturbation can be implemented to circumvent the phase transition. For example, adding a site-specific term H pert = -h i (-1) i S i z will break all symmetries of the Haldane phase, allowing for a finite energy gap along the entire ramp path. The ground state can then be characterized using site-specific measurements to determine the spin correlation functions S z i S z j and the string-order correlation S z ij ≡ S z i S z j i<k<j (-1) S z k . Since the interactions J ij in Eq. ( 39) are long-ranged, this can lead to both quantitative and qualitative differences in the phase diagram as compared to the nearest-neighbor XXZ model (Gong et al., 2016). For instance, the positions of the phase boundaries shift for long-ranged AFM interactions, whereas long-ranged FM interactions can destroy the Haldane phase and support a new continuous symmetry-broken phase. Each of these possible phases can be distinguished by comparing measured values of the spin and string-order correlation functions described above.
The first experimental steps towards Haldane physics in an ion-trap quantum simulator implemented the model in Eq. ( 39) with λ = 0 (Senko et al., 2015). To generate the ground states of this effective spin-1 XY model, for 2-and 4-ion spin chains, the spins were initially prepared in the state |00 • • • . This is the approximate ground state of Eq. ( 39) in the presence of a large D-field. This field was then ramped down slowly until D ≈ 0; the resulting state populations, shown in Fig. 21, match reasonably well with the exactly-calculated ground state.
Detection of the spin-1 states was accomplished by imaging the spin-dependent fluorescence (Olmschenk et al., 2007) |± states appear 'bright' during the detection process and are scattered into an incoherent mixture of the |F = 1 states, such a setup does not allow discrimination among all three possible spin states in a single experiment. However, the population of either |+ or |can be measured by repeating the experiment and applying a π rotation to the appropriate |0 ↔ |± transition before the fluorescence imaging. For instance, measuring an ion in the 'dark' state after a π pulse between |0 ↔ |+ indicates that the spin was in the |+ state before detection. This binary discrimination is not a fundamental limit to future experiments, since populations could be "shelved" into atomic states that do not participate in the detection cycle (Christensen et al., 2020).
Measurements of populations in the S z basis necessarily discard phase information about components of the final state. This can be important in many spin models, including the XY model, where such measurements alone cannot discriminate between different eigenstates. For example, the ground state of an XY model with two spin-1 particles is
, differing only by a relative phase. In (Senko et al., 2015), verification of ground state production was accomplished by employing a modified parity entanglement witness procedure (Sackett et al., 2000). First, global π/2 rotations were applied on both the |0 ↔ |+ and |0 ↔ |transitions, with a relative phase ϕ. Then parity Π(ϕ) = 2 n=0 (-1) n P n of the number of spins in state |0 was measured, where P n is the probability of n spins appearing in state |0 . This is expected to result in Π(ϕ) = 3 8 ± 1 2 cos ϕ, where the + andsigns correspond to the ground and highest excited states, respectively. The data shown in Fig. 22 clearly show the phase of the parity oscillations to be consistent with having prepared the 2-spin ground state of the spin-1 XY model.
In contrast to section III dealing with equilibrium properties of quantum spin models, trapped ion simulators are well-suited for studying non-equilibrium phenomena. Non-equilibrium dynamics might even be considered more natural, since the study of equilibrium-like properties requires a specific protocol for preparing the corresponding ground or thermal state, unlike conventional condensed matter materials that are directly cooled through phonon interactions. The simplest non-equilibrium studies, on the other hand, can start with an initial product state and then simply evolve the system under a (time-dependent) Hamiltonian of interest.
Trapped ion quantum simulators allow the study of non-equilibrium dynamics over a broad range of both spatial and temporal resolution. The effective long-range spin-spin interactions described in section I.C.2 can be modulated in time by turning on or off the laser-ion interactions. This allows non-equilibrium states to be prepared via quenches or stroboscopic application of the Ising Hamiltonian while their subsequent dynamics are observed over timescales both shorter and longer than the natural timescale of interactions 1/J ij . Single spin resolution can be achieved even as the system size is scaled to many particles, allowing access to non-trivial observables such as spin-spin correlations and magnetic domain sizes. Since strongly interacting and highly-frustrated Ising spin models are often employed in analytical and numerical studies of non-equilibrium quantum dynamics, the results of trapped ion spin simulations serve as an important benchmark for these theoretical predictions.
Perhaps the most natural non-equilibrium experiments are global and local quenches. In a global quench experiment, a simple initial state evolves under a time-independent Hamiltonian. A global quench originates from situations where the simple initial state can be naturally thought of as the ground state of some simple Hamiltonian, in which case the dynamics ensues when the Hamiltonian is changed (quenched). A local quench allows the comparison of the non-equilibrium dynamics between two initial states that differ by the application of a locally applied unitary operation. A particular example of a local quench is a situation where one of the two initial states is an eigenstate of the Hamiltonian. Short-time spin dynamics (Sec. IV.A) allows the investigation of time scales over which quantum information can propagate across the system. Long-time dynamics, on the other hand, can indicate whether the system eventually approaches some effective steady state and if so, how this steady-state is approached. Effective thermalization is often expected, even in a closed system of spins where a subset of the spin system uses its complement as the bath (D'Alessio et al., 2016)). However, as discussed in Sec. IV.B, disorder can often prevent such thermalization leading to many-body localization. Similarly, as discussed in Sec. IV.C, even in cases where the system eventually thermalizes, it is possible that dynamics dramtically slows and actual thermalization takes a very long-time to occur, in a phenomenon known as prethermalization.
There are many forms of inducing and probing spin dynamics in trapped ion systems. One example is the periodic modulation of a Hamiltonian, which gives rise to stroboscopic Floquet dynamics. In Sec. IV.D, we will discuss such dynamics in trapped ion systems in two complementary contexts. The first will focus on the application of a Hamiltonian and its negative counterpart in order to measure so-called out-of-time-ordered correlation (OTOC) functions. The second will focus on the spontaneous breaking of discrete time translation symmetry, leading to the emergence of time crystalline order.
Many properties of a many-body quantum system depend on how quickly quantum information can propagate in that system. Indeed, the speed of a quantum computation through an array of qubits will be boosted by sending quantum information faster across the array. Similarly, fast information propagation through qubits may allow for faster preparation of highly entangled states of the array. With its tunable approximately-power-law-decaying interactions, a trapped ion chain is an ideal testbed for studying how much such long-range interactions can speed up quantum information propagation relative to nearest-neighbor interactions and for elucidating the implications of such speedups.
To make connection with the native interactions in trapped ion spin crystals, we assume that the interactions fall off with distance as a power-law J 0 /|i -j| α between ions i and j, with 0 ≤ α ≤ 3, as derived in Eq. ( 23) and discussed above. For the purposes of studying the bounds on the speed of information propagation, it is the convenient to express the spin Hamiltonian in the following general form:
where h i is a Hamiltonian acting on spin i and where the two-spin Hamiltonian h ij acting on spins i and j is subject to the bound
Here ||O|| indicates the operator (or spectral) norm of operator O, or the magnitude of its largest absolute eigenvalue. We are interested in studying how quickly quantum information can propagate when the system evolves unitarily under Eq. ( 40).
Various notions of information propagation can be defined, depending on the particular problem (Tran et al., 2020). We first focus on the dynamics following a local-quench (Jurcevic et al., 2014). Let B be a unitary operator acting on a single site, while A is a single-site observable acting on another site a distance r away. Let |ψ be an arbitrary initial state, and let A(t) be the Heisenberg evolution of A under the Hamiltonian H in Eq. ( 40). Then the effect on observable A due to the disturbance B can be defined as the difference between the expectation values of A(t) in the original state |ψ and in the quenched state B |ψ ,
The signal after time t a distance r away is thus bounded by the unequal-time commutator ||[A(t), B]||. Intuitively, the operator A(t) can be thought to originate on a given site at t = 0 (and commuting with operator B on another site at distance r away) and then spreading in time until its support significantly overlaps with the support of B and thus allows for a substantial commutator
Upper bounds on ||[A(t), B]|| subject to the Hamiltonian in Eq. ( 40) are referred to as Lieb-Robinson-type bounds, named after the original work considering nearest-neighbor interactions (α = ∞) (Lieb and Robinson, 1972). The region in the r-t plane outside of which ||[A(t), B]|| must be small is called the causal region, while its boundary is called the (effective) light cone. While the Hamiltonian in Eq. ( 40) is time-independent, Lieb-Robinson bounds also hold for time-dependent h ij subject As the interaction range is increased (Fig. 4b,c), the arrival times of the first maxima in magnetization are seen to appear earlier and earlier, reflecting the ejection of faster and faster quasiparticles from the quench site. Furthermore, the signal decay outside these maxima is very slow: there is an almost instant increase in the magnetization even at large distances (Fig. 4d,top). Clearly we are able to tune our system into a regime where the light-cone picture does not apply and significant amounts of information can propagate directly to distant neighbours. This is consistent with generalized Lieb-Robinson bounds for power laws, which for av1 are trivial, placing no restriction on the speed of information propagation 6-8 .
A quantitative analysis is provided by extracting the maximum quasiparticle group velocity v max g from the data (see Methods and Extended Data Fig. 4). For the shortest-range case, the observed v max g fits well with the nearest-neighbour case (Fig. 4d). As the interaction range is increased, the results are consistent with a divergence of v max g , as recently predicted 12 . Ultimately, the information propagation speed in our system is limited by the propagation of acoustic waves across the ion chain 21 . Note that, despite the faster-moving components in the longer-range data (Fig. 4c), the initial perturbation remains more localized. This is consistent with the predicted flattening of the dispersion relation away from the divergence. For a comparison of data with theory, see Extended Data Fig. 4.
Differences between the observed and ideal quantum dynamics following local quenches largely correspond to imperfect conservation of excitation number. This could be caused by electric field noise leading to heating of the ion's motional state or by unwanted spin-motion entanglement. For global quench dynamics, laser-frequency and magneticfield fluctuations give rise to dephasing.
We have presented a new platform for investigating quantum phenomena-a many-body quantum system in which the states and properties of its quasiparticle excitations can be precisely initialized, controlled and measured. This opens many new paths for experimental investigations, the subjects of which can be broadly split into the following: (1) quantum transport phenomena, concerning how quantum states and entanglement 13 , or excitations 14,27 , propagate across quantum manybody systems; (2) how quantum systems reach equilibrium, including the question of when thermalization 15,28 and localization occur 16 ; (3) entanglement growth and simulation complexity 17 (the interaction range parameter a is known to play a critical role in the growth rate of entanglement and the possibility of simulating the dynamics with conventional computers); and (4) quasiparticle behaviour near phase transitions 1 . For many of these research lines it would be useful, and feasible, to add localized spin excitation absorbers or reflective boundaries, and static or stochastically fluctuating disorder, to our system.
During the final stage of this work, we became aware of complementary recent work investigating global quenches of trapped-ion spin chains 26 .
Ions are held in a linear Paul trap, each encoding a spin-1/2 particle in the electronic states S 1=2 ,m~z1=2i:
;i and D 5=2 ,m~z5=2i: :i. Spins are manipulated with a narrow-linewidth laser at 729 nm (ref. 29). Ions are coherently manipulated with two laser beams intersecting the ion string perpendicularly from opposite directions. The first beam interacts with all the ions with nearly equal strength and is used for carrying out collective spin rotations, as well as implementing effective spin-spin interactions by means of electronic-state-dependent forces 3 . These forces off-resonantly drive the transverse motional modes of the ion string. The interaction range a ð Þ is controlled by how far off-resonant the driving is and the axial trapping confinement. The second beam is strongly focused, steerable, and is used for single-spin rotations. Spatially resolved fluorescence measurements in conjunction with prior single-spin rotations allow us to take single-shot measurements of arbitrary spin correlations.
If our system had only nearest-neighbour interactions, the signal propagation after a local perturbation using s x ' would be bounded by
Þ , where O may be any local operator with norm O j j j j and distance d to the quench site '. As Fig. 4d shows for O~s z i , this bound is only a Methods). Bottom: for a 5 0.75, it does not. e, Maximum group velocity. With increasing a, the measured magnon arrival velocities (red circles) approach the group velocity of the non-renormalized nearest-neighbour model (grey dashdotted line). If renormalized by the algebraic tail, the nearest-neighbour group velocity increases at small a (orange dots), but much less than the increase of the observed magnon velocity. For small a, the measured arrival times are consistent with the divergent behaviour predicted for full power-law interactions (black line).
Macmillan Publishers Limited. All rights reserved ©2014 FIG. 23 (a-c With increasing α, the measured magnon arrival velocities (red circles) approach the group velocity of the non-renormalized nearest-neighbor model (grey dash-dotted line). If renormalized by the algebraic tail, the nearest-neighbor group velocity increases at small α (orange dots), but much less than the increase of the observed magnon velocity. For small α, the measured arrival times are consistent with the divergent behaviour predicted for full power-law interactions (black line). Adapted from Jurcevic et al., 2014. to Eq. ( 41) and for arbitrary time-dependent h i . A growing body of theoretical literature places upper bounds on ||[A(t), B]|| and therefore derives tighter and tighter light cones for different values of α (Chen and Lucas, 2019; Else et al., 2018;Foss-Feig et al., 2015;Gong et al., 2014a;Guo et al., 2019;Hastings and Koma, 2006;Kuwahara and Saito, 2019;Lashkari et al., 2013;Matsuta et al., 2016;Storch et al., 2015;Sweke et al., 2019;Tran et al., 2020Tran et al., , 2019Tran et al., , 2018)). At the same time, a complementary growing body of theoretical literature considers specific Hamiltonians and protocols demonstrating larger and larger causal regions (Bachelard and Kastner, 2013;Buyskikh et al., 2016;Cevolani et al., 2015Cevolani et al., , 2016Cevolani et al., , 2018;;Chen and Zhou, 2018;Chen et al., 2018;Eisert et al., 2013;Eldredge et al., 2017;Frérot et al., 2018;Gong et al., 2014a;Guo et al., 2019;Hauke and Tagliacozzo, 2013;Hazzard et al., 2014Hazzard et al., , 2013;;Jünemann et al., 2013;Kastner, 2011;Kloss and Bar Lev, 2019;Knap et al., 2013;Lepori et al., 2017;Luitz and Bar Lev, 2019;Maghrebi et al., 2016;Nezhadhaghighi and Rajabpour, 2014;Rajabpour and Sotiriadis, 2015;Schachenmayer et al., 2013;Storch et al., 2015;Van Regemortel et al., 2016;Worm et al., 2013). While these upper and lower bounds on information propagation are appearing to converge, there is still no provably tight light cone shape other than at large α, where the light cone is strictly linear (Chen and Lucas, 2019;Kuwahara and Saito, 2019;Tran et al., 2020).
Thanks to Eq. ( 42), Lieb-Robinson bounds on ||[A(t), B]|| directly constrain information propagation after local quenches. A local quench experiment on a trapped-ion chain has realized the XY model of hopping hard-core bosons corresponding to h i = 0 and
in Eq. ( 40) (Jurcevic et al., 2014). Figure 23 shows results for a local quench on 15 spins from the experiment. After flipping up the middle spin (corresponding to B = σ x 8 ) in a chain of down spins, the experiment measures A = σ z i for various i. Since the number of flipped spins is conserved during the time evolution, the evolution in this restricted Hilbert space of spins is well-described using the language of single-flip eigenstates called magnons.
Information propagation can also be studied in global quench experiments, where connected correlations are measured as a simple initial state (such as a product state) evolves under a given Hamiltonian. Linear light-cone-like spreading of such correlations due to a nearest-neighbor Hamiltonian was first measured in neutral atoms confined in an optical lattice (Cheneau et al., 2012). For long-range interacting systems such as trapped-ion spins, the spread of correlations may no longer be confined within a linear light cone. In particular, suppose the system starts in an initial product state |ψ and evolves under the Hamiltonian power-law scaling of the light-cone boundary in the XY model. However, without an exact solution there is no a-priori reason to assume a power-law light-cone edge (used for the fits in Fig. 3); deviations from power-law behaviour might reveal themselves for larger system sizes.
An important observation in Fig. 3j-l is that of faster-than-linear lightcone growth for our shortest-range interaction, with a 5 1.19. Although faster-than-linear growth is expected for a , 1 (see discussion of Ising model), there is no consensus on whether such behaviour is generically expected for a . 1. Our experimental observation has prompted us to numerically check the light-cone shape for a 5 1.19; we find that faster-thanlinear scaling persists in systems of up to 22 spins before our calculations break down (Extended Data Fig. 2).
Whether such scaling continues beyond ,30 spins is a question that, at present, quantum simulators are best positioned to answer. In Figs 2m, n and 3m, n, the excellent agreement between data and theory demonstrates that experiments produce the correct results in a regime still solvable by classical computers. For larger systems, where numerical evolution of the Schro ¨dinger equation fails, the quality of quantum simulations could still be benchmarked against the exact Ising solution of equation ( 4). Finding close agreement in the Ising case would then build confidence in an XY model simulation, which cannot be validated by any other known method.
For the XY model, we additionally study the spatial decay of correlations outside the light-cone boundary. The data (Fig. 4) is well described by fits to exponentially decaying functions. Recent theoretical work 20 predicts an initial decay of spatial correlations bounded by an exponential, followed by a power-law decay; we speculate that much larger in Eqs. ( 40)-(41). At time t = 0, the connected correlation function C ij (t) = O i (t)O j (t) -O i (t) O j (t) (where operator O i acts on site i) vanishes since the first expectation value factorizes. As time goes on, the supports of O i (t) and O j (t) grow, making C ij (t) grow in return. For the case of short-range interactions, C ij (t) is bounded in the r-t plane (where r = |i -j|) by a linear light cone similar to the corresponding light cone for the unequal time commutator ||[A(t), B]|| (Bravyi et al., 2006). For general α, a bound on ||[A(t), B]|| can also be used to derive a bound on C ij (t) (Gong et al., 2014b;Tran et al., 2020), but the relationship between the two light cones is not as trivial as in the nearest-neighbor (α = ∞) case. However, the physical picture is similar: the Lieb-Robinson bound constraints the spread of operators O i (t) and O j (t), so if the support of O i (t) (O j (t)) hasn't spread significantly outside of a ball of radius r/2 around i (j), then O i (t) and O j (t) have approximately disjoint supports, leading to small C ij (t).
Connected correlations following a global quench have been measured in a chain of trapped ions (Richerme et al., 2014). This experiment also studied the XY model corresponding to h i = 0 and 40). Figure 24 shows data on 10 ions following a global quench. Starting with an initial state of all spins pointing down (in the z basis), the experiment measures the time evolution of connected correlations
. The growth of connected correlations following a global quench may also be accompanied by the growth of entanglement, as was observed experimentally in (Friis et al., 2018).
Experiments on ultracold polar molecules (Yan et al., 2013) and defect centers in solid state (Choi et al., 2017a) do not yet have the single-spin resolution necessary for studying the shape of the causal region after local or global quenches in long-rangeinteracting systems. On the other hand, experiments on ultracold neutral atoms interacting via Rydberg-Rydberg interactions (Bernien et al., 2017;Guardado-Sanchez et al., 2018;Lienhard et al., 2018;Zeiher et al., 2017) should be able to access the particularly interesting parameter regime of α = 3 (dipolar interactions) (de Léséleuc et al., 2018) in one, two, and three spatial dimensions.
Many-body localization (MBL) has become one of the most studied nonequilibrium phases of matter, receiving considerable scrutiny in both experiment and theory in the past decade (Abanin et al., 2019;Choi et al., 2016;Nandkishore and Huse, 2015;Oganesyan and Huse, 2007;Schreiber et al., 2015). The localization effect is a generalization of single-particle Anderson localization, which is characterized by a cessation of quasiparticle transport in non-interacting systems subject to a random potential landscape (Anderson, 1958). Surprisingly, in the case of MBL similar insulator-like properties are observed even when particles are strongly interacting (Abanin et al., 2019;Nandkishore and Huse, 2015;Oganesyan and Huse, 2007). When prepared with a quench, the quantum states become highly entangled many-body superpositions of excited eigenstates spanning the entire energy spectrum of the disordered system Hamiltonian. MBL can be distinguished from Anderson localization by the logarithmic growth of entanglement entropy at long-times (Bardarson et al., 2012). The distribution of eigenstates occupied in an MBL phase is decidedly non-thermal and a number of observables have been identified to characterize phase transitions between MBL and thermal states when varying the interaction strength or disorder in the Hamiltonian (Abanin et al., 2019;Luitz et al., 2015;Nandkishore and Huse, 2015;Serbyn et al., 2014)
The disordered potential is implemented with single site resolution across the ion chain such that
where D i is sampled from a uniform distribution, D i ∈ [-W/2, W/2], with width W . For finite system sizes, this Hamiltonian exhibits features consistent with many-body localization and demonstrates a disorder-induced, long-lived memory of the system's initial conditions (Hauke and Heyl, 2015;Maksymov et al., 2017). Understanding the thermodynamic stability of localization with power-law interactions remains an intriguing open question (Burin, 2006(Burin, , 2015;;Nandkishore and Sondhi, 2017;Yao et al., 2014).
In both experiments, the MBL state is created by initially preparing the 10 spin Néel state with staggered order (|ψ 0 = |↓↑↓↑↓↑↓↑↓↑ z ) which is highly excited with respect to the disordered Ising Hamiltonian of Eq. 43. This Hamiltonian is rapidly quenched on and the resulting single spin magnetization dynamics σ z i (t) are measured for times up to t = 10/J 0 . The experiment is repeated under multiple instances of disorder with Stark-shifts (D i ) applied programmatically to each ion using a rastered individual addressing laser (Lee et al., 2016). This individual addressing laser is also used create the initial Néel state using a sequence of controlled spin-flips.
In the absence of disorder these initial spin states will thermalize if the uniform transverse field B is sufficiently large (Deutsch, 1991;Rigol et al., 2008;Srednicki, 1994). In Smith et al., 2016, global rotations are used to prepare eigenstates of both σ x and σ z and measure the resulting single ion magnetization projected into those directions after evolution under H MBL . In the case of a thermalizing system, memory of the initial spin configuration will be lost in all directions of the Bloch sphere, namely σ x i = σ z i = 0 at long times. Above the threshold transverse field (B 4J 0 ) the system rapidly thermalizes to zero magnetization after relatively short timescales (t < 5/J 0 ) (Fig. 25a).
However, with the transverse field held fixed at B = 4J 0 , the data clearly shows that applied disorder localizes the spin chain, retaining memory of the initial Néel state in measurements of the z magnetization σ z i (Fig. 25b). Each measurement of magnetization dynamics for disorder width W is repeated with at least 30 different realizations of disorder, which are subsequently averaged together. This is sufficient to reduce the finite depth disorder sampling error to be of the same order as other noise sources. After some initial decay and oscillations, the magnetization of each spin settles to a steady state value for J 0 t ≥ 5. The degree of localization can be quantified using the normalized Hamming distance (HD)
This observable counts the number of spin flips from the initial state, normalized to the length of the spin chain. At long times, a randomly oriented thermal state shows D = 0.5 while one that remains fully localized has D = 0 (Fig. 25c). The average steady state value D(t) for J 0 t ≥ 5 can serve as an order parameter to display the crossover between the localizing and thermalizing regimes. The most relevant adjustable experimental control parameters for probing the MBL phase diagram are the amplitude of disorder W and the interaction range α. Increasing W pins each spin closer to its initial state and pushes the entire spin chain towards a localized regime (Fig. 25d). Likewise, the localization strengthens as α is increased towards shorter range interactions (Fig. 25e), recovering Anderson localization via a Jordan-Wigner transformation (Jordan and Wigner, 1928) in the α → ∞ limit. Numerical studies have confirmed that full localization occurs within experimentally accessible disorder strengths and interaction ranges (Wu and Das Sarma, 2016).
The slow growth of entanglement entropy (S) has long been understood as a distinguishing feature of localization (Bardarson et al., 2012). A short range interacting MBL state should exhibit a slower entanglement growth rate than interacting quantum states without disorder, where entanglement spreads ballistically. The dynamics of the entanglement entropy are also quite different for non-interacting Anderson localized systems, where the entanglement saturates at short times once the system's dynamics have reached the localization length (Abanin et al., 2019). In a trapped ion quantum simulator with algebraically decaying interactions the entanglement entropy of an MBL state should also grow algebraically, S ∼ t q , but with q < 1 the dynamics are still distinct from those of non-localized or Anderson localized systems (Pino, 2014).
It is generally difficult to measure entanglement entropy in quantum simulators due to the exponential system-size scaling of the number of measurements required for full state tomography. In Smith et al., 2016 the observed slow growth in quantumfisher information (QFI) is used as a proxy for half-chain entanglement entropy, motivated by a similar scaling with disorder and interaction strength as observed in numerical simulations. A more direct measurement is made in Brydges et al., 2019, where they develop a technique to probe the second-order half-chain Rényi entanglement entropy in their 10 ion quantum simulator (S (2) (ρ [1→5] )) using randomized measurements. They find the entanglement growth to be significantly suppressed in the presence of strong disorder, in good agreement with numerical predictions (Fig. 25f).
Interestingly, Anderson localization can be explored using the same long-range disordered Hamiltonian, even though it is not strictly non-interacting, by observing the dynamics of a single spin excitation in the ion chain. For example, in Maier et al., 2019 the transport efficiency of a spin excitation from initial site i = 3 to the target site i = 8 is observed as a function of time in a (Rebentrost et al., 2009), and the quantum Zeno effect are indicated in gray. The data agree well with theoretical simulations of the coin-tossing random process (light bullets) realized in the experiment, while simulations with ideal Markovian white noise (lines) underestimate ENAQT at large γ. Inset (a): Sketch of the transport network. The ions experience a long-range coupling, with darker and thicker connections indicating higher coupling strengths. The green arrows denote the source (3) and the target ( 8) for the excitation in the ion network. Inset (b): Sketch of the ion chain representing interacting spin-1/2 particles as circles, with the spin states denoted by arrows. The ions are subject to random static and dynamic on-site excitation energies, indicated by Bi and Wi(t). 10-ion spin chain (Fig. 26a). The evolution of this single spin excitation can be described by an XY model Hamiltonian
where the disorder field contains single site static (B i ) and time-dependant [W i (t)] components (Fig. 26b). In the absence of disorder, the XY interaction term in this Hamiltonian conserves the total magnetization of the system, allowing a single spin excitation to hop around the chain. The transport efficiency on site i = 8 is then quantified by integrating the instantaneous probability of the excitation appearing on spin 8,
over the full duration of the experiment t max . The transport efficiency is reduced in the presence of strong static disorder (B max > J 0 ), consistent with Anderson localization of the initial spin excitation (Fig. 26 main). However, adding temporal variations in the form of white dephasing noise in W i (t) destroys the localization, a phenomenon known as environment-assisted quantum transport (ENAQT) (Rebentrost et al., 2009). The spectral density of the Markovian-like noise (S(ω) ∝ W 2 max ) in W i (t) determines the rate of dephasing γ = S(ω). This experiment can clearly access regimes where transport is inhibited by either Anderson localization γ < J 0 or the quantum Zeno effect γ > J 0 . In the ENAQT regime, when γ ≈ J 0 , the temporal noise modifies the destructive interference necessary for Anderson localization and transport is revived. Section IV.D further explores experiments studying trapped ion spin dynamics under the influence of both static disorder and a periodically time-varying Hamiltonian.
Many-body localization is a unique case in which a closed quantum system remains non-ergodic and localized even up to infinite times. The trapped ion quantum simulations of Smith et al., 2016 andBrydges et al., 2019 are limited by finite experimental coherence times to only one decade in J 0 t. This makes it difficult for these experiments to quantify how longlived the magnetization or slow the entanglement growth might be. Fortunately, other experiments that have studied MBL using cold neutral atoms can achieve several orders of magnitude longer evolution time relative to their interaction timescale. For example, MBL can be realized using cold fermions in quasi-random 1D optical lattices (Kohlert et al., 2019;Lukin et al., 2019;Lüschen et al., 2017b;Schreiber et al., 2015), verifying MBL-like behavior in a variety of Hubbard Hamiltonians. This system has also been used to confirm the breakdown of MBL in open quantum systems (Bordia et al., 2016;Lüschen et al., 2017a). Moreover, experiments have started to probe whether MBL can exist in systems with dimensionality > 1 (Bordia et al., 2017;Choi et al., 2016;Kondov et al., 2015), where the stability of MBL is still an open question (De Roeck and Imbrie, 2017). Other experimental platforms have used novel metrics to probe many-body localization, including many-body spectroscopy (Roushan et al., 2017), measuring out-of-time order correlators (Wei et al., 2018) (see Sec. IV.D), and performing full state tomography to compute entanglement entropy (Xu et al., 2018).
Hamiltonians that support MBL are believed to be non-ergodic, even after evolution times exponentially long in the system size (Nandkishore and Huse, 2015). There are also systems that are non-ergodic for a shorter amount of time (but often still much longer than the coherence time of typical quantum simulation experiments) before eventually thermalizing. Usually, these systems are not disordered and can be described by models of weakly interacting (quasi-)particles, such as 1D Bose gases (Gring et al., 2012;Kinoshita et al., 2006;Langen et al., 2015). The generic behavior of such a system is called prethermalization, meaning that the system relaxes to a quasi-stationary state different from thermal state before thermalizing eventually. The prethermal quasi-stationary state is usually believed to be described by a generalized Gibbs ensemble (GGE) (Rigol et al., 2007) that corresponds to the model of quasi-particles without the weak interactions. Such state will have a partial memory of the initial state, because the quasi-particle occupation numbers are conserved if interactions are ignored. At sufficiently long time, the weak interactions are expected to break integrability of system and lead to thermalization in the end. This picture of prethermalization has been well studied both in theory (Berges et al., 2004;Manmana et al., 2007;Polkovnikov et al., 2011) and experiment (Gring et al., 2012;Langen et al., 2015).
In a programmable ion-trap quantum simulator, due to long-range spin interactions, new types of prethermalization can occur with prethermal states not described by a standard GGE. An example study was first proposed theoretically (Gong and Duan, 2013) and later demonstrated experimentally (Neyenhuis et al., 2017). The central idea is that with sufficiently long range interactions, a non-disordered and homogenous system can have a strong emergent inhomogeneity due to the open boundary condition of an experimental spin chain. This emergent inhomogeneity can lead to trapping of quasi-particles before the system relaxes to GGE. As both kinetic energy and weak interactions can delocalize trapped quasi-particles, the dynamics of the system can reveal a rich interplay between quantum tunneling and interaction effects, leading to new types of relaxation beyond conventional prethermalization.
The model under study in Refs. (Gong and Duan, 2013;Neyenhuis et al., 2017) is the same transverse field Ising model described by Eq. ( 25). With long-range interactions, H T I is generally nonintegrable, in contrast to the nearest-neighbor case where the 1D model is integrable through a Jordan-Wigner transformation (Sachdev, 2011), so thermalization is anticipated in the long-time limit according to the eigenstate thermalization hypothesis (Rigol et al., 2008). To better understand the dynamics of such a pre-thermal Hamiltonian, we can map each spin excitation along the z direction into a bosonic particle to turn Eq. ( 25) into a bosonic model with two parts: An integrable part composed of noninteracting spin-wave bosons that can be used to construct a GGE, and an integrability-breaking part consisting of interactions among the spin-wave bosons, which is responsible for the thermalization (Neyenhuis et al., 2017). When the initial state has a low spin/ bosonic excitation density and the magnetic field is much larger than the average Ising coupling J ij , the bosonic excitation density will remain low during the dynamics and the interactions among the bosons will remain weak.
The experiment in (Neyenhuis et al., 2017) begins by preparing a single spin excitation on either edge of a 7-ion chain |ψ R = |↓↓↓↓↓↓↑ z or |ψ L = |↑↓↓↓↓↓↓ z . The spins then evolve under Eq. ( 25) and the time evolution of the spin projection in the z-basis is measured. The magnetic field B is at least an order of magnitude larger than J 0 in the experiment so the number of spin excitations along the z direction is approximately conserved in the short time dynamics where the system can be regarded as a single spin excitation. But in the long time dynamics, multiple spin excitations will be created and interacting with each other.
To characterize the dynamics of the spin excitations, we introduce a single observable that measures the relative location of the spin excitation in the chain
where N is the number of ions. The expectation value of C varies between -1 and 1 for a spin excitation on the left and right ends, respectively. The choice of initial states ensures that the initial value of C is either 1 or -1. Due to the spatial inversion symmetry of the underlying Hamiltonian in Eq. ( 25), both the GGE and thermal values of C should be zero. In Fig. 27(c-d) the value of C along with its cumulative time average C are shown for the two initial states with a single spin flips on either end of the spin chain. In the short-range interacting case (α = 1.33), where the system rapidly evolves to a 1 2 3 4 5 6 7 Ion # 1 2 3 4 5 6 7 Ion # 1 2 3 4 5 6 1 2 3 4 5 6 7 Ion # 1 2 3 4 5 6 7 Ion # 1 2 3 4 5 6 prethermal state predicted by the GGE associated (with C = 0) with the integrals of motion corresponding to the momentum space distribution of the single particle representing the spin excitation. The memory of the initial spin excitation location is thus not preserved. However, in the long-range interacting case (α = 0.55), the position of the spin excitation reaches a quasistationary value that retains a memory of the initial state out to the longest experimentally achievable time of 25/J max . This prethermal state differs clearly from both a thermal state and the GGE prediction, which both maintain the left-right spatial symmetry of the system. For evolution times too short to generate more than one spin flip, the dynamics of the Hamiltonian in Eq. ( 25) for the initial states are similar to those of a free-particle in a potential, with the location of the particle representing that of the single spin excitation. For short-range spin interactions, the shape of the potential is approximated by a square well, due to the open boundary condition and no explicit spatial inhomogeneity of the interactions. However, as we increase the range of spin-spin interactions, the shape of the potential distorts from a square well to a double well shaped potential formed by the two hard walls at the ends of the spin chain and the bump at the center of the chain, as shown in Fig. 27a. For a single particle on a lattice with a double well potential, there will be an extensive number of near-degenerate eigenstates that are symmetric and antisymmetric superpositions of wavefunctions in the left and right potential wells. For seven lattice sites, the spectrum of energy differences between all pairs of eigenstates as a function of α is shown in Fig. 27b, together with the overlap of eigenstates with the initial state. For the longest range interaction (α = 0.55), the two lowest energy states are almost degenerate, with an energy difference approximately 1000 times smaller than J max . This stems from the tunneling rate between the two double well, which is exponentially small in the barrier height, resulting in the spin excitation remaining in its initial well until it tunnels across the potential barrier at very long times.
To go beyond the single particle picture above, the experiment in (Neyenhuis et al., 2017) prepares initial states with two spin excitations. In this case, there will be weak interactions between the two particles that represent the spin excitations, similar to the scenario for many-body localization (Smith et al., 2016). Despite the presence of weak interactions, similar prethermal states were found, as shown in the bottom panel of Fig. 27(c-d): Relaxation to GGE is found for shorter-range interactions (α = 1.3) while for longer-range interactions the system clearly does not relax to GGE. Similar results were also found in a spin chain of 22 ions, as shown in Fig. 28 (Neyenhuis et al., 2017). The persistence of the same prethermalization observed with more than a single spin excitation is attributed by the existence of extensive number of nearly degenerate eigenstates for the single particle spectrum in the double well shown in Fig. 28a. Thus an extensive number of spin excitations near one end of the chain will still leading to markedly different thermalization/prethermalization time scales in certain systems (32)(33)(34). We believe that the current experiment, as well as the platform it is built upon, will pave the way to a more complete understanding of the fundamental role long-range interactions play in the quench dynamics and emergent statistical physics of quantum many-body systems.
We generated spin-spin interactions by applying spin-dependent optical dipole forces to ions confined in a three-layer linear Paul trap with a 4.8-MHz radial frequency. Two off-resonant laser beams with a wave vector difference dk along a principal axis of transverse motion globally address the ions and drive stimulated Raman transitions. The two beams contain a pair of beat note frequencies symmetrically detuned from the resonant transition at n 0 = 12.642819 GHz by a frequency m, comparable to the transverse motional mode frequencies. In the Lamb-Dicke regime, this results in the Ising-type Hamiltonian in Eq. 1 (22, 23, 35) with
We initialized ind beam to imprint a a Ramsey or Rabi beam waist of the seven-ion short-ra first optically pump p/2 rotation so th ion addressing bea and then we allow phase compared to global p/2 rotatio method, individua whereas N spin fli We used the R cause site-resolved aration is smaller shift to all of the sp at the hyperfine sp the ions without an single and N spin After quenchi Hamiltonian, we m z direction of the B addresses the cycli fluoresce only if th lected through an N using an intensifie site resolution.
To discriminat respectively), we b of all-bright and a the 2D CCD imag distributions at eac and fluorescence w
We achieved s the experimental d Gaussian distribut varying amplitud compared with a termine whether th tion protocol also misdiagnosing a b Corrected state pr found following th We initialize the spins with a single spin excitation on the left end (middle image). After evolving for 36 J max , the spin excitation is delocalized, but its average position remains on the left half of the chain (bottom image). Effective range of the interaction is a ≈ 0.9. Error bars, 1 SD.
28 Time evolution (light blue) and cumulative time average (orange) of the averaged center of excitation C in a 22-ion chain. The ion spins are initialized with a single spin excitation on the left end (middle image). After evolving in the XY spin Hamiltonian for time t = 36/J0, where J0 is the nearest-neighbor coupling, the spin excitation is delocalized, but its average position remains stuck on the left half of the chain (bottom image), the signature of prethermalization. Effective range of the interaction is α ≈ 0.9. Adapted from Neyenhuis et al., 2017. be localized by the double well before tunneling happens at a later time.
The interplay between single-particle tunneling in the effective double well and particle-particle interactions is always present in this system; even when initialized in a single spin excitation state, the finite transverse field in Eq. ( 25) will create more spin excitations over time. The effect of interactions will thermalize the system, while the effect of tunneling will bring the system to the GGE. Thus, depending on the range of spin-spin interactions, it is possible to either observe the prethermalization to GGE after the prethermalization caused by the trapping in double well potential, or observe thermalization directly after the observed prethermalization. Determining which scenario is relevant would require a much longer coherence time than is possible in current experiments. While improving the experimental coherence time of the spin interactions may be challenging, simulating the long-time dynamics of an non-integrable, long-range interacting spin chain is equally, or even more challenging on a classical computer.
It should be emphasized that the emergence of an effective double-well potential for the spin excitations is a phenomenon unique to an open spin chain with strongly long-range interactions. A spin chain with periodic boundary condition and without spatial inhomogeneity is fully translationally invariant, and translational invariance is not expected to be broken in the long time behavior of the system. It may be surprising that changing boundary conditions of a long-range interacting system can significantly impact its bulk properties. Prethermalization in trapped ion spin crystals is thus a good example of new physics that is possible with programmable quantum simulators, reaching beyond existing condensed matter frameworks, as also demonstrated by other experiments mentioned in this section.
Section IV.B treated static (time-independent) Hamiltonians whose non-equilibrium nature arises from the presence of quenched disorder leading to many-body localization; such localization prevents the system's internal dynamics from thermalizing and leads to certain memory of local initial conditions (Abanin et al., 2019;Nandkishore and Huse, 2015). An alternate setting for exploring non-equilibrium phases is to begin with a time-dependent Hamiltonian whose equations of motion are intrinsically dynamical. This setting is particularly ideal for trapped ion quantum spin simulators, where pulsed control of both the interactions and fields of the transverse field Ising model of Eq. ( 25) naturally allow for the realization of time-dependent Hamiltonians (e.g. see III.A where time-dependent magnetic fields are utilized).
Recently, a tremendous amount of theoretical and experimental work has been devoted to exploring the mildest case of such time-dependence, where the system is governed by a periodic Hamiltonian, H(t + T ) = H(t). Such Floquet systems (Bukov et al., 2015;Floquet, 1883;Kuwahara et al., 2016) are particularly ideal from the perspective of quantum simulation since they do not require the complexities of cooling to the many-body ground state in order to observe novel dynamical phenomena (Bloch et al., 2008). In the case of trapped ions, as we have previously discussed in section II, there exists a natural capability to stroboscopically apply different microscopic Hamiltonians, making this platform an ideal Floquet quantum simulator. Here we focus on two specific examples. First, we describe the implementation of a novel class of measurements termed out-of-time-ordered correlation (OTOC) functions (Larkin and Ovchinnikov, 1969;Maldacena et al., 2016). Such correlators have recently been proposed as powerful diagnostics of quantum chaos and their dynamical behavior remains the subject of intense interest (Landsman et al., 2019;Swingle et al., 2016;Yoshida and Yao, 2019); indeed, the possibility of defining a quantum Lyapunov exponent based on the exponential growth of OTOCs in certain systems has led to a conjectured bound on the rate of thermalization in many-body quantum systems (Maldacena et al., 2016).
The functional form of the OTOC is typically written as the expectation
where W and V are two commuting operators with W (τ ) = e iHτ W e -iHτ the evolved operator W under the Hamiltonian H.
The OTOC compares two quantum states obtained by either (a) applying V, waiting for a time t, and then applying W; or (b) applying W at time t, going back in time to apply V at t = 0, and then letting time resume its forward progression to t. The OTOC F (τ ) thus describes how initially commuting operators W and V fail to commute at later times due to the interactions generated by H, and have a quite different form than the conventional autocorrelation function W † (τ )V (0)W † (τ )V (0) . Operationally, measuring an OTOC requires an intermediate step of time-reversal, representing a major challenge from an experimental implementation perspective (Gärttner et al., 2017;Li et al., 2017).
In the case of ions, it is possible to directly measure an OTOC stroboscopically applying both an interaction Hamiltonian and its negative counterpart, thus effectively reversing time evolution. This was experimentally demonstrated in (Gärttner et al., 2017) by using a two-dimensional array of laser-cooled 9 Be + ions in a Penning trap, summarized in Fig. 29. Here, the spindependent optical dipole force couples to the axial center-of-mass motion of the 2D crystal (see Fig. 2b), resulting in a nearly uniform all-to-all Ising interaction (Eq. ( 23) with α ≈ 0). Such an all-to-all interacting Ising model does not possess the complexity of interactions from a quantum chaos perspective. However, this system is well-suited for OTOC studies, because the sign of the Hamiltonian can be controlled. Since the Ising interaction strength scales as J ∼ 1/δ (Eq. ( 22) with just a single mode contributing and δ m = δ), this allows for the implementation of a Hamiltonian sign reversal by simply changing the sign of the detuning. With the ability to stroboscopically apply both H and -H, a family of different OTOCs can then be measured through the collective magnetization of the system, as shown in Fig. 29 (Gärttner et al., 2017).
The second example of stroboscopic Hamiltonian simulation with trapped ion spin systems is the quantum simulation of Floquet systems, or periodically modulated Hamiltonians (Bukov et al., 2015;Deng et al., 2015;Kuwahara et al., 2016). In a broader context, such time-periodic manipulations have long been used for controlling quantum systems including NMR qubits and atomic ensembles (Goldman et al., 2014;Oka and Kitamura, 2019). However, recent explorations of Floquet systems have stumbled upon an intriguing question beyond the landscape of quantum control; in particular, can Floquet systems host intrinsically new phases of matter that do not have any equilibrium equivalent (Else et al., 2020;Moessner and Sondhi, 2017)?
In the noninteracting (single-particle) case, the question has been affirmitavely answered with the discovery of a host of novel band structures that can only exist in the presence of periodic driving (Cayssol et al., 2013;Kitagawa et al., 2010;Lindner et al., 2011;Rechtsman et al., 2013;Titum et al., 2016). The many-body case is more subtle. On the one hand, one might naturally suspect that new phenomena can in principle arise when the driving frequency is of order the intrinsic energy scales of the system; indeed, this limit is far from the Suzuki-Trotter limit (e.g. as discussed in V.A) where to first order, the effective Hamiltonian describing the Floquet system is simply a sum of its stroboscopic components. On the other hand, it is generally expected that a driven many-body system will absorb energy from the driving field and ultimately heat to infinite temperature (Bukov et al., 2015;Ponte et al., 2015a). However, recent theoretical advances have demonstrated that it is possible to avoid such Floquet heating. One general scheme is to utilize many-body localization as discussed in section IV.B. In principle, a Floquet MBL system (Abanin et al., 2016;Ponte et al., 2015b) can exhibit stable dynamical phases of matter for infinitely long times (Else et al., 2016;Khemani et al., 2016;Potirniche et al., 2017). Interestingly, recent studies suggest an alternative disorderfree approach can also be used to combat Floquet heating, albeit not to infinite times. In particular, for large enough driving frequencies, the system can enter a regime of Floquet-prethermalization (Kuwahara et al., 2016;Mori et al., 2016), where exotic non-equilibrium phases can be observed for exponentially long time-scales (Else et al., 2017;Machado et al., 2020Machado et al., , 2019)). The underlying essence of Floquet-prethermalization is analogous to the discussions of prethermalization in section IV.C. The key difference is that here the lifetime of the quasi-stationary state is controlled by the Floquet drive frequency.
We now turn to recent experiments that demonstrated a Floquet quantum simulation using a one dimensional trapped ion spin chain, depicted in Fig. 30 (Zhang et al., 2017a). In these experiments, a combination of high-precision spatial and temporal control allowed for the implementation of three distinct types of time evolution applied repeatedly in sequence: 1) global spin rotations, 2) long-range Ising interactions, and 3) disordered on-site fields (Fig. 30a).
The stroboscopic combination of these evolutions is the basis for realizing a discrete time crystal (DTC) (Else et al., 2016(Else et al., , 2020;;von Keyserlingk et al., 2016;Khemani et al., 2016;Yao et al., 2017), where a system exhibits a spontaneous breaking of the time-translation symmetry generated by the Floquet evolution. The characteristic signature of a DTC, which is consistent with experimental observations in chains of up to L = 14 ions (Zhang et al., 2017a), is the robust synchronization of oscillations at sub-harmonic frequencies compared to that of the drive, as shown in Fig. 30b-c. Within the decoherence time-scale of the experiments, the observed signatures of DTC order are independent of the initial state. Crucially, the robustness of the subharmonic oscillations depends on the presence of strong interactions in the system; in the absence of interactions, even small perturbations immediately destroy signatures of a time crystal.
In addition to implementations in trapped ion systems, a number of other experimental platforms have also observed signatures of time crystalline order (Choi et al., 2017b;Rovny et al., 2018). Here, we note a particular set of experiments performed using ensembles of nitrogen-vacancy (NV) color centers in diamonds (Choi et al., 2017b). We emphasize this particular platform because it shares a number of similar features with the trapped ion system (i.e. long range interactions and disorder), but also has a number of crucial differences (i.e. three dimensional system with time-dependent disorder). Interestingly, both platforms exhibit similar signatures of interaction stabilized time-translation symmetry breaking, although these signatures appear to be limited to time-scales before local thermalization has fully completed (Else et al., 2020). To this end, such cross-platform verifications are especially valuable once controlled quantum simulators reach a regime where classical computers cannot calculate (Calarco et al., 2018;Leibfried, 2010). In this regard, comparing the results from analog quantum simulators to those from digital quantum computers would also be helpful in cross-checking and assessing validity (Mahadev, 2018).
Having discussed the simulation of non-equilibrium phases in both disordered and periodically driven trapped ion experiments, we now turn to the question of understanding phase transitions in such out-of-equilibrium systems. Novel dynamical phases can emerge after a quantum quench, and the transition between them can be observed experimentally by measuring non-analytic changes in the dynamical response of the many-body spin system (Jurcevic et al., 2017;Zhang et al., 2017b). As also described in sections IV.B and IV.C, out-of-equilibrium systems do not necessarily behave thermodynamically, so it is a fundamental question how to properly establish analogies and differences among thermodynamic equilibrium phases and their dynamical counterparts (Heyl, 2018;Titum et al., 2019), in terms of order parameters (Ajisaka et al., 2014;Žunkovič et al., 2018), scaling and universality (Heyl, 2015), and discrete or continuous symmetry breaking (Huang et al., 2019;Weidinger et al., 2017;Žunkovič et al., 2016).
Dynamical phases can be separated by dynamical quantum phase transitions (DQPT), characterized by a non-analytic response of the physical system as a function of quench parameters. Two types of DQPT signatures have been defined for an interacting spin-1/2 chain ( Žunkovič et al., 2018) governed by the Hamiltonian of Eq. ( 25) with the field along the z-axis:
where the spin-spin interactions J ij and the transverse field B z are generated using the technique explained in section I.C.2. Both types of DQPT have been experimentally observed in a trapped-ion quantum simulator (Jurcevic et al., 2017;Zhang et al., 2017b). The first type of DQPT (type I) is based on the formal analogy between the non-analytic behaviour of the return probability to the initial state |ψ 0 after a quantum quench under the Hamiltonian H, defined as G(t) = ψ 0 | e -iHt |ψ 0 , and the partition function of the corresponding equilibrium system Z = Tr(e -H/k B T ) (Heyl et al., 2013). It is possible to define the complex counterpart of the thermodynamic free energy density f = -N -1 k B T log(Z) using the rate function
. This quantity, in the thermodynamic limit, exhibits dynamical real-time nonanalyticities that play an analogous role as the non-analytic behaviour of the free energy density of a thermodynamic system at equilibrium. It is possible to observe experimentally these nonanalyticities in an interacting spin chain after a quantum quench evolving under the long-range Transverse Field Ising Hamiltonian of Eq. ( 50). This type of DQPT has been observed experimentally with a linear chain of trapped 40 Ca + ion spins (Jurcevic et al., 2017). The spins are initialized in the ground state of the field part of the transverse Ising model, namely |ψ 0 = |↓↓↓ ... ↓ z . Then the transverse field Hamiltonian [Eq. ( 50)] is suddenly switched on (quenched) with B > J 0 , with J 0 being the average nearestneighbour spin-spin coupling. As shown in Fig. 31a, in this regime the rate function λ exhibits pronounced nonanalyticities at the which coincides with Λ for large N, and which we use in our further discussions of the DQPTs.
In Fig. 2(a), we report our first main result, the direct observation of a DQPT through nonanalyticities in the rate function λ. Let us emphasize that the observed DQPTs in λ are neither artificially caused by our definition of PðtÞ nor by the resulting minimum construction. The definition PðtÞ ¼ P ⇒ ðtÞ þ P ⇐ ðtÞ is physically motivated [22,31] by allowing us to connect the DQPTs to other physical quantities, as we demonstrate also below. The rate function λ, on the other hand, provides a tool to quantitatively extract the DQPT already for small systems, resulting in very weak residual finite-size corrections [24], such that we can focus in the following on a single system size. Without the the sharpening of the finite-size crossover for increasing system size, which is much more intricate [24]. While in the following we concentrate our discussion mainly on the firs DQPT, let us emphasize that also the subsequent DQPTs are of the same nature, possibly with the role of P ⇒ ðtÞ and P ⇐ ðtÞ exchanged.
To study the robustness of DQPTs against deformations of the Hamiltonian, we extract the first critical time t c from λðtÞ as a function of the coupling strength J ¼ ðN -1Þ -1 P i>j J ij ; see Fig. 2(b). We find that the temporal nonana lytic behavior is stable over a broad range of J=B and fo different α. For J=B ≪ 1, the critical time τ c -π=4 ∝ ðJ=BÞ 2 exhibits a quadratic dependence on J=B yielding τ c ¼ π=4 for J ¼ 0 where the dynamics becomes equiv alent to Larmor precession of N independent spins. While DQPTs also appear in this apparently simple case, it is important to emphasize that J ¼ 0 represents a singula point in the dynamics due to the absence of nonloca quantum fluctuations as becomes apparent from the entan glement dynamics we discuss later on.
We now present measurements that connect DQPTs to other observables, further corroborating the theory o DQPT as a key framework for understanding quantum many-body dynamics. In Figs. 3(a) and 3(b), we compare λðtÞ and the evolution of the magnetization, M x ðtÞ ¼ hM x ðtÞi with M x ¼ N -1 P i σ x i . The initial state breaks the global Z 2 symmetry σ x i → -σ x i ∀i of the Hamiltonian H The system responds to this symmetry breaking by a repeated crossover between the M x > 0 and M x < 0 sectors, reaching the symmetry-restoring value M x ¼ 0 at specific times. Comparing with λðtÞ, these are tied to the critical times of the DQPT, whose essence is the symmetry restoration in the ground-state manifold.
This connection is tightened by resolving the magneti zation M x ðε; tÞ as a function of energy density ε (see Supplemental Material [24], and Ref. [31]), where ε ¼ E=N and E is the energy measured with the initia Hamiltonian H 0 . The measured data are displayed in Fig. 3(c). The dynamics along ε ¼ 0 (ground-state mani fold) is directly understood from the previous discussion. In large systems, as long as t < t c one has PðtÞ ≈ P ⇒ ðtÞ yielding M x ðε ¼ 0; t < t c Þ ≈ 1. For t > t c , P ⇐ ðtÞ takes over, and M x ðε ¼ 0; tÞ jumps to -1. With increasing energy densities this sudden change smears out. Its influ ence, however, persists up to the system's mean energy density εðtÞ [solid line in Fig. 3(c)], where observables such as M x ðtÞ acquire their dominant contribution [31]. In this way, as sketched in Fig. 1, an extended region of the dynamics is controlled by the DQPT, reminiscent of a quantum critical region at an equilibrium QPT.
As the final result of our work, we now show that DQPTs in the simulated Ising models also control entanglemen production. In this way, we connect entanglement as an important concept for the characterization of equilibrium (a) Measured rate function λðτÞ for three different system sizes at J=B ≈ 0.42, showing a nonanalytical behavior (with τ ¼ tB being the dimensionless time). Dots are experimental data with error bars estimated from quantum projection noise; lines are numerical simulations with experimental parameters. In a lighter black color we have included data for the subdominant contributions at N ¼ 6. Inset: The transition between the normalized ground-state probabilities P ⇒ =P (solid lines) and P ⇐ =P (dashed lines) becomes sharper for larger N. (b) The first-order correction to the critical time τ crit , i.e., the occurrence of the first DQPT, is linear as a function of ðJ=BÞ 2 for small J=B, and approximately independent of interaction range. Error bars are 1σ confidence intervals of the fits on log½P ⇒;⇐ ðτÞ from which we extract τ crit (see Supplemental Material [24]). Inset: DQPT for ðJ=BÞ ¼ 0, 0.392, and 0.734 (light blue, dark blue, and black dots). The grey dashed lines indicate τ crit for ðJ=BÞ ¼ 0.
PRL 119, 080501 ( 2017)
week ending 25 AUGUST 2017 the half-chain entropy SðtÞ measured by quantum tomography (see Supplemental Material [24]). SðtÞ exhibits its strongest growth in the vicinity of a DQPT. While these data are suggestive of entanglement production, SðtÞ is an entanglement measure only for pure states, which does not account for the experimentally inevitable mixing caused by decoherence. Therefore, we additionally measure a mixedstate entanglement witness, the Kitagawa-Ueda spinsqueezing parameter ξ s [33] (see Supplemental Material [24]) signaling entanglement whenever ξ s < 1. As Fig. 4(b) shows, ξ s presents a behavior qualitatively very similar to SðtÞ. Related to common spin-squeezing scenarios [34], the spin squeezing is most effective when the mean spin vector on the Bloch sphere is perpendicular to the direction of the spin-spin interaction. Importantly, this occurs when M x ¼ 0, which we found above to be inherently tied to DQPTs. The presence of the DQPT, moreover, offers a more general interpretation: At exactly t c , the ground-state manifold enters the equal superposition ðj ⇒i þ j⇐iÞ= ffiffi ffi 2 p , a highly entangled Greenberger-Horne-Zeilinger (GHZ) state. Just as for the case of M x , our data suggest that the influence of this state stretches to elevated energy densities, and thus DQPTs critical times t c . This behaviour can be related to other observables, such as the global average magnetization M x = N -1 i σ x i . Since the initial state breaks the Z 2 symmetry of the Hamiltonian (25), the system restores this symmetry during the evolution at the times where the magnetization changes sign, which also corresponds to the critical times in the Loschmidt echo observable, as shown in Fig. 31b.
The second type of DQPT-(type II) has an order parameter defined in terms of long time averaged observables, such as asymptotic late-time steady states of local observables:
where the operator A is the magnetization or higher order correlators between the spins. Here, the DQPT occurs as the ratio B/J 0 is varied and the order parameter changes abruptly from ferromagnetic (B < J 0 ) to paramagnetic order (B > J 0 ). The onset of this non-analytic behaviour can be observed by measuring the late time average values of the two-body correlator
after a quantum quench with Hamiltonian (50). This type of DQPT measurement was observed in a linear chain of trapped 171 Yb + ion spins (Zhang et al., 2017b). Here, the measured late time correlator C 2 exhibits a dip at the critical point that sharpens scaling up the system size N up to 53 171 Yb + ions, as shown in Fig. 32a. This represents one of the largest quantum simulations ever performed on individual spins. Further evidence of the occurrence of the phase transition can be also observed in higher-order correlations, such as the distribution of domain sizes throughout the entire chain, shown in Fig. 32b. The occurrence of the DQPT is observed in the decreased probabilities of observing long strings of aligned ions at the critical point. This is more clearly shown measuring the mean largest domain size as a function of the transverse field strength, for late times and repeated experimental shots, which exhibits a sharp transition across the critical point of the DQPT-type II. Note this measurement is in general an N th order spin correlation function, requiring high fidelity readout of all spins in a single shot, which is an attribute that is unique to trapped ion systems.
The tailored interactions between spins available in trapped ion systems can also be used to realize a digital quantum simulator (DQS), which is a many-body quantum system with enough control to perform a universal set of quantum operations or gates. In other words, it is a quantum computer used to implement Hamiltonian evolution rather than other quantum algorithms (Feynman, 1982). The universality means that these simulations are not limited to the inherent interactions, but any local Hamiltonian can be designed and implemented as a circuit (Lloyd, 1996). Additionally, error bounds and error correction protocols using faulttolerant gate sets will be applicable to large digital simulations (Steane, 1999). Arbitrary Hamiltonians are changed into the spin a b
long-range nature of the Ising interactions 31 .
We further explore many-body dynamical properties of this system by investigating higher-order correlations, which are even harder to calculate classically 25 . Through high-efficiency single-shot state detection of all of the spins, we measure the distribution of domain sizes in the chain directly as a higher-order correlation observable (see Methods). Single-shot images for N = 53 spins are shown in Fig. 4a and are reconstructed from binary thresholding and image convolution of the fluorescence distribution of the ion chain (see Methods).
The occurrence of long domains of correlated spins in the state | ↑ 〉 x (fluorescing spins) signifies the fully polarized initial state, in which the correlations are largely preserved by the interactions. With an increasing transverse field, the absence of spin ordering is reflected by exponentially small probabilities of observing long strings. We plot the domain length statistics at late times in Fig. 4a (see Methods) for three transverse field strengths, / = . . . B J ( 0 1,1 0,1 6) z 0 . The dashed lines in Fig. 4a are fits to exponentials on the histogram of domain sizes. The rare occurrence of especially large domains (see, for example, the red boxes in Fig. 4a) demonstrates the existence of many-body high-order Dashed lines in a and b are calculations using a canonical (thermal) ensemble with an effective temperature corresponding to the initial energy density. For N = 53 spins (d), the correlations are uniformly degraded from residual Stark shifts across the ion chain, so in this case we normalize to the maximum correlation at small field (see Methods). Exact diagonalization for N = 53 spins is not possible, so we instead fit the experimental data to a Lorentzian function with linear background (dashed line). © 2017 Macmillan Publishers Limited, part of Springer Nat
Further signatures of the DPT are observed by measuring the spatially averaged two-spin correlations
From the behaviour of the magnetizations described above, we expect that C 2 → 1 for small B z and C 2 → 1/2 for large B z at long times, because the collective spin precesses around the z axis and C 2 oscillates between 1 and 0. In Fig. 3 we show the cumulative time-averaged correlations.
Near the critical value of B z we observe the emergence of a dip in C 2 , which is a direct signature of the DPT. The sharpening of the dip for larger system sizes is not strong, which might be due to a logarithmic finite-size scaling (see Methods).
For a non-integrable system such as the long-range transverse-field Ising model studied here, it might be conjectured that the spins eventually reach a thermal distribution 30 . However, we find that this is true only for small B z (Fig. 3a,b). We note that the thermal values of the correlator C 2 do not exhibit a dip or signatures of a phase transition with varying / B J z 0 for the system sizes that we are able to model numerically. Interestingly, thermalization appears to break down in this quenched system, which we suspect is a consequence of the inherent long-range nature of the Ising interactions 31 .
We further explore many-body dynamical properties of this system by investigating higher-order correlations, which are even harder to calculate classically 25 . Through high-efficiency single-shot state detection of all of the spins, we measure the distribution of domain sizes in the chain directly as a higher-order correlation observable (see Methods). Single-shot images for N = 53 spins are shown in Fig. 4a and are reconstructed from binary thresholding and image convolution of the fluorescence distribution of the ion chain (see Methods).
The occurrence of long domains of correlated spins in the state | ↑ 〉 x (fluorescing spins) signifies the fully polarized initial state, in which the correlations are largely preserved by the interactions. With an increasing transverse field, the absence of spin ordering is reflected by exponentially small probabilities of observing long strings. We plot the domain length statistics at late times in Fig. 4a (see Methods) for three transverse field strengths, / = . . . B J (0 1, 1 0, 1 6) z 0
. The dashed lines in Fig. 4a are fits to exponentials on the histogram of domain sizes. The rare occurrence of especially large domains (see, for example, the red boxes in Fig. 4a) demonstrates the existence of many-body high-order that C 2 → 1 for small B z and C 2 → 1/2 for large B z at long times, because the collective spin precesses around the z axis and C 2 oscillates between 1 and 0. In Fig. 3 we show the cumulative time-averaged correlations.
Near the critical value of B z we observe the emergence of a dip in C 2 , which is a direct signature of the DPT. The sharpening of the dip for larger system sizes is not strong, which might be due to a logarithmic finite-size scaling (see Methods). For a non-integrable system such as the long-range transverse-field Ising model studied here, it might be conjectured that the spins eventually reach a thermal distribution 30 . However, we find that this is true only for small B z (Fig. 3a,b). We note that the thermal values of the correlator C 2 do not exhibit a dip or signatures of a phase transition with varying / B J z 0 for the system sizes that we are able to model numerically. Interestingly, thermalization appears to break down in this quenched system, which we suspect is a consequence of the inherent long-range nature of the Ising interactions 31 .
We further explore many-body dynamical properties of this system by investigating higher-order correlations, which are even harder to calculate classically 25 . Through high-efficiency single-shot state detection of all of the spins, we measure the distribution of domain sizes in the chain directly as a higher-order correlation observable (see Methods). Single-shot images for N = 53 spins are shown in Fig. 4a and are reconstructed from binary thresholding and image convolution of the fluorescence distribution of the ion chain (see Methods).
The occurrence of long domains of correlated spins in the state | ↑ 〉 x (fluorescing spins) signifies the fully polarized initial state, in which the correlations are largely preserved by the interactions. With an increasing transverse field, the absence of spin ordering is reflected by exponentially small probabilities of observing long strings. We plot the domain length statistics at late times in Fig. 4a (see Methods) for three transverse field strengths, / = . . . B J (0 1, 1 0, 1 6) z 0
. The dashed lines in Fig. 4a are fits to exponentials on the histogram of domain sizes. The rare occurrence of especially large domains (see, for example, the red boxes in Fig. 4a) demonstrates the existence of many-body high-order © 2017 Macmillan Publishers Limited, part of Springer Nat
that C 2 → 1 for small B z and C 2 → 1/2 for large the collective spin precesses around the z axis a 1 and 0. In Fig. 3 we show the cumulative tim Near the critical value of B z we observe the e which is a direct signature of the DPT. The s larger system sizes is not strong, which migh finite-size scaling (see Methods).
For a non-integrable system such as the lon Ising model studied here, it might be conjectu tually reach a thermal distribution 30 . Howeve only for small B z (Fig. 3a,b). We note that th correlator C 2 do not exhibit a dip or signatu with varying / B J z 0 for the system sizes that we ically. Interestingly, thermalization appear quenched system, which we suspect is a cons long-range nature of the Ising interactions 31 .
We further explore many-body dynamical by investigating higher-order correlations, whi culate classically 25 . Through high-efficiency si of all of the spins, we measure the distributio chain directly as a higher-order correlation ob Single-shot images for N = 53 spins are shown structed from binary thresholding and image rescence distribution of the ion chain (see Me
The occurrence of long domains of correla (fluorescing spins) signifies the fully polariz the correlations are largely preserved by th increasing transverse field, the absence of spin exponentially small probabilities of observing domain length statistics at late times in Fig. 4a transverse field strengths, / = . . B J (0 1, 1 0, 1 or qubit basis via the Jordan-Wigner (Jordan and Wigner, 1928) or the Bravyi-Kitaev (Bravyi and Kitaev, 2002) transformation. First-quantization mapping can also been used (Babbush et al., 2017;Bravyi et al., 2017). There are several ways a time evolution or adiabatic ramp has been approximated by discrete operations in a DQS, Trotterized evolution (Suzuki, 1985;Trotter, 1959) or variational approximations (Yuan et al., 2019), in particular the Variational Quantum Eigensolver (VQE) (Peruzzo et al., 2014) and the Quantum Approximate Optimization Algorithm (QAOA) (Farhi et al., 2000). There are other alternatives such as quantum walks or Taylor series expansion. We note that the latter has potential for fault-tolerant devices due to its optimal asymptotic scaling, but has not been demonstrated experimentally (Berry et al., 2015).
In order to decompose the unitary evolution under a given Hamiltonian into discrete operations or gates, Suzuki-Trotter (Suzuki, 1985;Trotter, 1959) formulas are commonly used. These provide an approximate factorization of the unitary timeevolution operator Û . For a time-independent Hamiltonian H containing a sum of terms H k , the first-order Trotter approximation Û = e -i Ĥt/ = lim n→∞ (Π k e -i Ĥk t/n ) n for finite n leads to a sequence of small time steps δt = t/n with local evolution operators Ûk = e -i Ĥδt/ that is repeated n times. These operators can then be deconstructed into the operations from the simulator's universal set. Similarly, a time-dependent Hamiltonian can be simulated by breaking the evolution down into short segments (Poulin et al., 2011). Since the unitaries U k do not commute in general for finite n, the sequence deviates from the target evolution. These Trotter errors are bounded for small step sizes but undergo a dynamical phase transition (Heyl et al., 2019), above which the evolution exhibits many-body quantum chaos (D'Alessio et al., 2016;Sieberer et al., 2019).
This technique was used by (Lanyon et al., 2011) in a trapped-ion experiment to simulate a wide range of different interaction models and interaction graphs of two to six spin-1/2 particles. This spin-model problem matches the system and does not require a transformation. Using the operations available from the toolbox described in section I.C.2, a universal set of operations is realized. The versatility of this approach is shown by simulating the dynamics of two-particles under an XX (Ising), XY, and XYZ-type interaction, and an adiabatic ramp of the interaction strength. Figure 33A shows the dynamics of four spins initialized along z under a long-range (all-to-all) Ising interaction. The oscillation frequencies correspond to energy gaps in the system. Figure 33B presents the collective oscillation of six spins under a σ 1 y σ 2 x σ 3 x σ 4 x σ 5
x σ 6
x interaction, composed digitally of Ising-type interactions and indvidual σ z rotations, periodically producing a six-spin GHZ state.
Interleaving different native operations has also been proposed for creating an XXZ-interaction in a spin-1 system (Cohen et al., 2015). This will allow the study of an integer-spin systems under Heisenberg-type interactions, including the Haldane phase (Haldane, 1983). Adding optical pumping as an incoherent process to the set of operations was used in (Barreiro et al., 2011) to create the ability to engineer open-system dynamics, and show the dissipative preparation of multi-qubit stabilizer states for quantum error correction. Evolution by DQS can also be combined with standard computational gates to measure quantities of interest. A digital simulation of the two-site Fermi-Hubbard model was realized in (Linke et al., 2018). After a trotterized adiabatic ramp composed of two-qubit XX-interactions and single qubit rotations, the second Rényi entropy was measured. In contrast to (Brydges et al., 2019), this is achieved by creating two copies of the state and measuring the expectation value of the SWAP operator following (Horodecki and Ekert, 2002;Johri et al., 2017).
The DQS approach has also been proposed for quantum field theory simulations using trapped ions (Bauls et al., 2019). The Schwinger model, an example of a lattice gauge theory, describes fermions coupled to an electro-magnetic field in one spatial dimension. After mapping the model to qubits with the Jordan-Wigner transformation using open boundary conditions, the Hamiltonian reads (see (Kokail et al., 2019;Muschik et al., 2017) for details):
where n labels the lattice sites. Following a Kogut-Susskind eoncoding, spin-down (spin-up) on an odd (even) lattice site indicates the presence of a positron (electron). The first term represents particle-antiparticle pair creation and annihilation with coupling w. The second is a mass term and the third term represents the energy of the electric field Ln with coupling g. The model reproduces key features of more complex high-energy physics theories such as quantum chromodynamics (Gattringer and Lang, 2010), which are extremely challenging for classical numerical methods.
An experiment realizing the Schwinger model in a trapped-ion spin chain has been implemented in (Martinez et al., 2016), showing very good agreement with the ideal evolution for the real-time dynamics of a 2-site (4-qubit) model. The difficulty in mapping gauge theories to spins lies in the infinite-dimensional Hilbert space of the gauge field operators. Instead of truncating the model to the size of the simulator, the authors analytically integrate out the gauge fields using Gauss' law, allowing the electric field to be expressed in terms of Pauli oprators: Ln = 0 -1 2 Σ N l=n+1 (σ z l + (-1) l ). This creates a spin Hamiltonian with asymmetric long-range couplings (Hamer et al., 1997), which can be directly implemented by a sequence of multi-ion spin-spin operations (Martinez et al., 2016;Muschik et al., 2017). By construction, the scheme constrains the evolution within the subspace allowed by Gauss' law and hence maintains the gauge invariance of the theory (Muschik et al., 2017). These features, as well as the linear scaling of the number of gate operations and qubits with lattice size, make DQS with trapped ions a promising avenue for the simulation of more complex gauge theories. The authors also realize a variational lattice gauge theory simulation (Kokail et al., 2019), which will be discussed below.
Quantum walks are alternative framework to discretize a unitary evolution for a DQSs. Quantum walks, the quantum equivalent of classical random walks (Childs et al., 2003), see a system evolve into on a discrete grid or continuum of spacial locations based on the superposition of a quantum coin. Discrete step quantum walks have been mapped to represent simulations of different physical systems (Chandrashekar et al., 2010;Molfetta and Pérez, 2016), including the evolution of a particle under the Dirac equation (Mallick et al., 2019). The latter was realized in trapped ions by mapping the position space to multi-qubit spin states, where different particle masses were chosen by varying the weight of the quantum coin (Huerta Alderete et al., 2020). Trapped ion experiments have previously realized quantum walks (Schmitz et al., 2009;Zähringer et al., 2010) and separately a Dirac equation (Gerritsma et al., 2010) simulation, using the harmonic oscillator degrees of freedom in the trap rather than mapping to spins. A non-variational digital simulation of spin models was also performed with superconducting qubits (Salathé et al., 2015).
Variational quantum simulation (VQS) implements the versatile classical method of variational simulation (Balian and Veneroni, 1988;Shi et al., 2018;Szabo and Ostlund, 2012) on a quantum computer by employing a quantum-classical hybrid s). The current leading mit both the simulation e thought to be laser ). This is not currently and, once properly adincrease in simulation can be combined with iques for analog simue of systems that can be present work, and cur-rent ion trap development ( 35), digital quantum simulations involving many tens of qubits and hundreds of high-fidelity gates seems feasible in coming years. approach. Hybrid quantum algorithms are a promising way to solve potentially difficult problems on non-fault-tolerant, nearterm quantum systems (Preskill, 2018) by combining classical and quantum resources. They distribute the task between a quantum and a classical computer. The quantum system is running a parameterized sequence of operations to generate a classically intractable state of interest, e.g. following an approximate evolution under a Hamiltonian or a quantum gate sequence. The classical system varies the parameters to minimize a cost function evaluated by measuring the quantum state, a task that is considered relatively easy based on its complexity class (Yuan et al., 2019). Employing trapped ion spin simulators, this concept has been used to train a generative model, a routine from classical machine learning, using an ion trap-based digital quantum circuit with up to 26 parameters (Zhu et al., 2019b). The authors showed that the choice of optimizer is crucial as the parameter space grows. One of the challenges of this approach lies in the cost function landscapes, which can exhibit so-called barren plateaus (McClean et al., 2018), making optimization hard. Finding the classical optimization strategy and making efficient use of quantum resources is an active area of research (Cirstoiu et al., 2019;Sundar et al., 2019;Wecker et al., 2015;Yang et al., 2017).
VQS employs such a hybrid approach to simulating physical models. For solving static problems, the Rayleigh-Ritz variational method is generalized to the quantum regime. The ground state energy of a Hamiltonian Ĥ = λ i ĥi , given as a linear combination of tensor products of local operators ĥi , can be estimated by considering trial wave functions |φ( θ) with parameters θ = (θ 1 , θ 2 , ...). and approximating the ground state energy by an upper bound (Yuan et al., 2019):
E est 0 is obtained by preparing |φ( θ) through a sequence of parameterized quantum operations. Different circuit parameters implement different instances of θ, and φ( a)| Ĥ|φ( a) is obtained by measuring each term φ( a)| ĥi |φ( a) , and calculating the linear combination given by {λ i }. A classical optimization routine, e.g. a gradient descent algorithm, uses this value as a cost function to converge to E est 0 . A clear illustration of this method was presented in (Kokail et al., 2019), where the authors perform a VQS of the Schwinger model described in section V.A above. Starting from |Ψ 0 = | ↑↓ .. ↑↓ , which corresponds to the bare vacuum state, they apply finite-range Ising interactions (see section I.C.2) to all ions, and alternate with individual σz shifts realized with a Stark-shift beam. The set of discrete operations used in this ansatz enforces the symmetries of the problem, namely Gauss' law and gauge invariance. Each operation is parameterized and the cost function is given by Ψ( θ)| ĤS |Ψ( θ) , where ĤS is the Schwinger Hamiltonian from equation 53, and |Ψ( θ) is the trial wave function. The cost function is evaluated by expressing ĤS as a string of Pauli operators and measuring their expectation values by projecting the experimental output state in the appropriate basis. An adaptive classical optimization process based on the Dividing Rectangles (DIRECT) algorithm (Liu et al., 2015) is employed in a quantum-classical feedback loop find the minimum. At the minimum, the classical parameters θ 0 describe a sequence to prepare the approximate ground state wave function, and the cost function gives an upper bound estimate of the ground state energy. The authors use this method with 20 (16) ions to optimize over an up to 15-dimensional parameter space and achieve 85% (90%) ground state fidelity and an estimate of the ground state energy within 2σ uncertainty. One of the challenges in variational quantum algorithms is the verification of convergence since classical simulation cannot provide a measure of success. For eight ions, the authors use the variance of the cost function as a termination criterion by measuring the Hamiltonian in additional sets of bases to determine Ĥ2 S , making the optimization "self-verifying". In a final 8-ion measurement, a phase transition (Byrnes et al., 2002) is observed by varying the mass in the Schwinger Hamiltonian and observing an order parameter (Kokail et al., 2019). The authors also show that the Rényi entropy peaks at the critical point by the method described in (Brydges et al., 2019). To overcome the deep circuits required for digital LGT simulations, an analog approach with tailored spin-spin interactions in an ion chain has recently been proposed (Davoudi et al., 2020).
Different ways to construct circuits or sequences of quantum operations for preparing the trial wave function have been used. On the one hand, there is the so-called hardware-efficient ansatz as used in (Kokail et al., 2019), where the sequence is constructed from available operations without trying to reproduce the precise problem Hamiltonian. It is straightforward to implement but more susceptible to optimization problems like barren plateaus (McClean et al., 2018). The the other hand, the physically-inspired ansatz prescribes a more problem-specific construction for the trial states, for example based on Unitary Coupled Cluster (UCC) theory as mentioned below. The cost function tends to better match the underlying problem, but it can quickly lead to deep circuits as the complexity of the physical system increases (Shehab et al., 2019a).
For the simulation of physical models two flavors of variational algorithms are distinguished, which are both based on the VQS method but used in different domains of application, the variational quantum eigensolver (VQE) (Peruzzo et al., 2014) for the determination of Hamiltonian eigenvalues (Yuan et al., 2019), e.g. in the determination of ground and excited state energies quantum chemistry, and the quantum approximate optimization algorithm (QAOA) (Farhi et al., 2014) for combinatorial optimization such as graph problems and approximate state preparation. The following sections will review results on both from recent work with trapped ions. We note that a probabilistic variant has also been proposed for implementation with trapped-ions (Zhang et al., 2020).
The Variational Quantum Eigensolver (VQE) is a DQS algorithm for the calculation of operator eigenvalues, which has been employed to tackle molecular electronic structure and molecular dynamics problems in quantum chemistry (McClean et al., 2016;Peruzzo et al., 2014;Yung et al., 2014) as an efficient alternative to quantum phase estimation (Aspuru-Guzik et al., 2005). Under the Born-Oppenheimer approximation, which fixes the inter-nuclear distances, the quantum chemistry Hamiltonian can be written in a second quantization formulation using a basis of molecular orbitals. These are formed following the Hartree-Fock method by a linear combination of atomic orbitals (Helgaker et al., 2000). The electronic part of the Hamiltonian then reads:
where the summation is over all of the M molecular basis states. The factors h pq and h pqrs are related to one-electron and two-electron transitions respectively and are calculated numerically (Yung et al., 2014). The fermionic creation and annihilation operators fulfill {â, â † } = 1, which enforces the antisymmetry of the wave function. The Hamiltoninan is mapped to spins (Bravyi and Kitaev, 2002;Jordan and Wigner, 1928) and the expectation values are calculated for trial wave functions of the ground state, generated via an ansatz circuit as described previously. A references state, which represents a classical approximation to the ground state such as a Hartree-Fock solution, usually a product state, can be chosen as the initial state. A physically-inspired ansatz circuit for this type of problem is the Unitary Coupled Cluster (UCC) (Helgaker et al., 2000). The scheme prescribes a series of fermionic operators called cluster operators T = T 1 + T 2 + ... + T N for an N electron system (Peruzzo et al., 2014). These are essentially electron-hopping operators of increasing order that form a unitary to evolve the system into the ground state through the exploration of Hilbert space by going to increasing Hamming distances from the initial state (McClean et al., 2016), ÛUCC = e θ(T -T † ) . Using a pseudo-time evolution of this unitary by first-order Trotterzation is a good approximation of this evolution (Hempel et al., 2018;Yung et al., 2014): ÛUCC Π i e θi(Ti-T † i ) . In the spin basis, this creates circuits of higher-order spin-spin interactions, which can be realized on a quantum computer, particularly in trapped ions via multi-ion Ising-operations, or MS-gates. Classical coupled-cluster calculations can be employed to find reasonable starting values for θ (Romero et al., 2018). The computational complexity of the UCC method scales polynomially with the number of molecular orbitals M . The procedure can be repeated for different nuclear configurations to sample out the energy surface. Alternatively, a quantum phase estimation algorithm can be employed to this end (Yuan et al., 2019).
Below we document recent implementations of VQE with trapped ions, although VQE has also been implemented with other experimental platforms (Arute et al., 2020a;Colless et al., 2018;Dumitrescu et al., 2018;Kandala et al., 2017;O'Malley et al., 2016;Peruzzo et al., 2014). VQE was implemented in trapped 171 Yb + ions (Shen et al., 2017) to simulate the electronic structure of the molecular cation HeH + . Instead of Ising operations, the authors used the four available states in the ground level of a 171 Yb + ion coupled by microwave-and radiofrequency-driven transitions. The experiment matches the ground state energy from exact diagonalization. The authors also look for excited states by measuring the expectation value of a changed target Hamiltonian to ( Ĥ -λ1) 2 . The measurement involved scanning λ and looking for zero-valued expectation values and resolved three out of four excited states.
In (Hempel et al., 2018), the authors mapped two molecular problems to spins realized in chains of trapped 40 Ca + ions, the Hydrogen molecule, H 2 , and Lithium Hydride, LiH. For H 2 , the molecular orbitals are formed by linear combinations of the two 1s orbitals of each atom. The Bravyi-Kitaev (Bravyi and Kitaev, 2002) and Jordan-Wigner (Jordan and Wigner, 1928) encoding are used to map the problem onto two and four qubits, respectively. The former is more efficient since both the Hamiltonian and the UCC ansatz can be reduced to two qubits by appropriate choice of reference state (O'Malley et al., 2016). In both mappings, the UCC operator leads to a single-parameter ansatz circuit, which is implemented using a pair of MS gates based on a technique developed in (Müller et al., 2011). The authors do a full sweep of the parameter as well as a VQE optimization to map out the ground state potential energy using a Nelder-Mead search algorithm, which was found to not converge in some instances due to the presence of noise. For the simulation of LiH, a minimal basis set for the molecular orbitals was chosen. The four electrons are filled into the energy-ordered orbitals, as determined by a Hartree-Fock calculation. The UCC ansatz is truncated to single and double excitations (UCCSD) and an active space of two electrons and three orbitals identified (Kandala et al., 2017). The two core electrons can be considered to not contribute to molecular bonding (Roos et al., 1980). Using the classical theory of configuration interaction of singles and doubles (CISD) (Shavitt, 1984), the active space is further reduced to two singletexcitations. This active space is illustrated in Fig. 34a. The resulting UCC unitary was mapped to spins via Bravyi-Kitaev encoding, which results in a three-qubit operator with two circuit parameters, ÛUCCSD (α, β) = e -iασ x 0 σ y 1 e -iβσ x 0 σ y 2 , which was realized in the circuit shown in Fig. 34c. The authors performed a grid scan of the parameters at different inter-nuclear distances and reconstructed the potential energy curve for the electronic ground state by fitting the minimum in two ways, by Gaussian process regression (GPR) and to 2D-quadratic surface, where the latter yields very good results, reproducing the well depth within statistical uncertainties, see line plots in Fig. 34b. The authors also performed a feedback VQE using a classical optimizer that incorporates elements of simulated annealing and matches the results of the grid scan, see data points in Fig. 34b.
In (Nam et al., 2020) the authors estimated the ground state energy of the water molecule (H 2 O). The UCC ansatz circuits are optimized by taking advantage of an active subspace and by sorting the UCC operators into bosonic and non-bosonic excitatios from the Hartree-Fock ground state. Bosonic terms refer to excitations of two electrons into a molecular orbital such that their spins pair up. These excitations can be represented by a single-qubit operation. Non-bosonic excitations are implemented via the Jordan-Wigner transformation. The resulting circuits are further optimized by term ordering (Hastings et al., 2014). The authors present the cost of adding an increasing number of ansatz terms to the circuit in terms of qubits and MS gates, see fig. 35a. Their results show that chemical accuracy can be reached with about 15 terms, which corresponds to 11 qubits and 140 MS-gates. Experimentally, the authors realize circuits with up to three terms and three qubits. The measurement results agree within errors with classical calculations, see fig. 35b. Finally, the authors proposed the implementation of this simulation problem with a growing number of terms as a benchmark for quantum computing platforms (Nam et al., 2020).
VQE was used in (Shehab et al., 2019b) to solve a nuclear physics problem, finding the binding energy of the Deuteron nucleus. The nuclear interaction is modeled in pionless effective field theory (EFT) following (Dumitrescu et al., 2018). The Hamiltonian is expressed in the N -oscillator basis,
where the operators â † n and ân create and annihilate, respectively, a deuteron in the harmonic-oscillator s-wave state. Using finite squeezed state that behaves as a thermal state for one mode when traced over the other mode, and is of great interest in a number of areas of physics. These states provide a way to prepare thermal (Gibbs) states of a many-body Hamiltonian in a quantum simulator, which underpin phenomena like high temperature superconductivity (Lee et al., 2006) and quark confinement in quantum chromodynamics (Gross et al., 1981). TFD states also play a key role in the holographic correspondence where they represent so-called traversable wormholes (Gao et al., 2017;Maldacena et al., 2017) and a number of approaches for their preparation have been proposed (Cottrell et al., 2019;Maldacena and Qi, 2018;Martyn and Swingle, 2019;Wu and Hsieh, 2019). A TFD state corresponding to inverse temperature β is defined on a joint system of two identical Hilbert spaces A and B:
where Z(β) is a normalization factor and |n A and |n B are the eigenstates in space A and B with corresponding energy E n .
In (Zhu et al., 2019a), the authors generated a TFD state similar in form to Eq. ( 58 ρ A is maximally mixed. The evolution alternates between inter-system coupling H AB = i X i,A X i,B + Z i,A Z i,B and the intra-system Hamiltonians H A + H B , where H A and H B are identical transverse-field Ising Hamiltonians of the subsystems, H XX +gH Z ≡ L i=1 X i X i+1 +g L i=1 Z i . This evolution was approximated by variational application of these Hamiltonians. The sequence was executed with theoretically optimal parameters for p = 1, and found to reproduce the expected inter-and intrasystem correlators well. The authors also use a related approach (Ho et al., 2019) to directly prepare the zero-temperature ground state of the quantum critical transverse field Ising model with seven trapped ion spins using quantum-classical feedback, starting from a random point for p = 1 and from the ideal parameters for p = 2.
The QAOA approach was used in an ion trap system to approximately generate the critical ground state energy of the longrange anti-ferromagnetic transverse field Ising model (Eq. ( 25)) (Pagano et al., 2019), where the spin-spin interaction J ij is defined in Eqs. (22,23) and B is the transverse field. Here we set H A to be the Ising coupling term and H B to be the field term in Eq. ( 25). QAOA was employed with up to 40 qubits in this system, using both grid search and closed loop approaches. In the grid search approach, the whole parameter space ( β, γ) was explored experimentally to find the optimum; in the closed-loop approach the analog trapped-ion quantum simulator was interfaced with a greedy gradient-descent algorithm to optimize the measured energy (Fig. 36a) starting from an initial guess. In the p = 1 QAOA, the optimization trajectory in the experiment can be visualized on the theoretical performance surface as shown in Fig. 36b.
In (Shehab et al., 2019a) the authors mapped the MaxCUT problem for a specific 5-node graph to a five-qubit trapped ion processor and solved it approximately using an optimized QAOA protocol, that reduces circuit width and depth by splitting it into circuits for each Hamiltonian term or correlator at the cost of a larger number of measurements. For the five-node Dragon graph, the authors reported a 36% improvement in the accuracy of the solution over the standard circuit in their system.
Finally, we note that VQS can also be used to study dynamics by employing generalizations of time-dependent variational principles (Broeckhove et al., 1988) to simulate real (Heya et al., 2019;Li and Benjamin, 2017) and imaginary (Chen et al., 2019;Jones et al., 2019;McArdle et al., 2019a) time evolution of pure quantum systems, as well as mixed states (Yuan et al., 2019).
Quantum interacting spin models are among the simplest many-body quantum systems with nontrivial features that can elude classical computational approaches. Trapped atomic ion spins offer the ability to implement and control quantum spin models with tunability of the interaction form and range. Much of the research in this field has concentrated on studies of the long-range transverse Ising model, which can feature frustrated ground states with associated degeneracies and entanglement in the ground state. The many types of phase transitions and dynamical processes in this system form a fruitful test-bed for studying quantum non-equilibrium processes, in many cases challenging classical computational power even for small numbers of spins. There are many extensions in this physical system to simulating more complex spin models with trapped ions such as Heisenberg couplings (Graß and Lewenstein, 2014), higher-dimensional spin models, and interactions involving three or more spins. These future directions may allow the quantum simulation of more exotic spin phases such as spin liquids (Balents, 2010), or topological orders in spin systems such as the Haldane chain (Haldane, 1983) or the Kitaev lattice (Kitaev, 2006).
There are a number of technical challenges to overcome in the path to scaling up trapped-ion systems to ∼ 100 qubits and beyond. However none of these obstacles are fundamental and many groups have shown effective mitigation techniques in this growing area of research. With a single linear chain of ions, it should be possible to extend the number of ions to hundreds using anharmonic trapping potentials that generate uniform spacing and thus avoid structural phase transitions for larger ion crystals (Lin et al., 2009;Pagano et al., 2018). Two technical challenges include the increased sensitivity of equally-spaced ion crystals to background fluctuating electric fields (Cetina et al., 2020) and the shorter lifetime of the ion chain due to collisions with the background gas. Solutions to these challenges include the use of cryogenic trapped-ion systems, that significantly lower the vacuum pressure so that storage times of hours or longer can be achieved with large ion crystals (Pagano et al., 2018;Poitzsch et al., 1996), and the use of sympathetic cooling with multiple species ions in the system (Chou et al., 2010;Larson et al., 1986;Lin et al., 2009;Pino et al., 2020), allowing the background heating to be quenched without disturbing the spins in the simulation.
With a 1D ion chain, it is possible to simulate any frustrated spin graph with multiple engineered laser pulses (Davoudi et al., 2020;Korenblit et al., 2012). However, it is possible to further increase the number of ions in a quantum simulator by using higher-dimensional ion crystals (Wang et al., 2015). The Penning trap geometry is amenable to 2D and 3D crystal structures, as discussed in this review. In a rf trap, 2D or 3D ion crystals are necessarily accompanied with significant rf micro-motion for ions not positioned at a rf null (Dehmelt, 1967), which is limited to a point or a line in space. Although micro-motion will introduce a number of subtleties for laser cooling, gate operations, and ion detection, it is a coherent and well-controlled motion whose effects can be mitigated through proper gate design and the engineering of spin interactions (Shen et al., 2014;Wang et al., 2015). Alternatively, an array of rf traps can be engineered to create multidimensional designer ion crystals (Cirac and Zoller, 2000;Kumph et al., 2016), although this may require very small distances between the ions and the electrodes to generate sufficiently strong spin interactions. Apart from offering an avenue to increase the ion number, 2D and 3D ion crystals provide a natural platform to realize rich frustrated spin Hamiltonians in a more complicated geometry that can match the native dimension of the underlying target spin Hamiltonian.
Finally, we note that well-known modular approaches to scaling trapped ion systems can be directly applied to quantum spin simulations. For instance, large numbers of ions can be grouped in spatially-separated modules, where each module contains many trapped ion spins with Coulomb-mediated spin-spin couplings discussed in section I. Here, the system is scaled by shuttling subsets of ions between modules to extend the size of the system (Kielpinski et al., 2002;Pino et al., 2020). Because the idle spin states are nearly perfectly decoupled from ion motion between modules, the scaling procedure is largely a systems engineering task. Ultimately, separated modules of trapped ion crystals can also be connected via heralded photonic interconnects, even between separated trap structures, for a remote modular scalable architecture (Duan and Monroe, 2010;Li and Benjamin, 2012;Monroe et al., 2014).
Useful quantum simulations should have unambiguous benchmarks to indicate performance beyond what is possible with classical computational simulation. As opposed to recent simulations of random quantum circuit evolution (Arute et al., 2019), we desire a benchmark that also indicates the usefulness of the underlying problem. When we use a quantum simulator to probe the ground-state of a quantum many-body Hamiltonian for instance, a possible heuristic benchmark could be the ground-state variational energy of the Hamiltonian. Although in general there is no quantum algorithm to guarantee that we can successfully find the genuine ground-state of a many-body Hamiltonian, the associated variationtal energy of the approximate ground state realized by a quantum simulator can always be efficiently measured by the experimental setup. If such a measured energy is lower than the energy obtained by any classical simulation even with the most powerful computers, it is an indication that the quantum simulator finds a useful solution that better approximates the ground state of the target Hamiltonian compared with any known classical method.
Of course, there is a close relationship between spin simulations and quantum computations with qubits, and the underlying mechanism behind the Ising couplings in trapped ion spin simulations is exactly that used for discrete quantum gates between trapped ion qubits (Sørensen and Mølmer, 2000), which are sometimes called Ising gates. Quantum simulations in this sense can be considered as a special case of a quantum computation, and it should be expected that as trapped ion quantum computers scale in the future (Brown et al., 2016;Monroe and Kim, 2013;Pino et al., 2020;Wright et al., 2019), so will the reach of trapped ion quantum spin simulators.