Holographic quantum algorithms for simulating correlated spin systems

We present a suite of"holographic"quantum algorithms for efficient ground-state preparation and dynamical evolution of correlated spin-systems, which require far-fewer qubits than the number of spins being simulated. The algorithms exploit the equivalence between matrix-product states (MPS) and quantum channels, along with partial measurement and qubit re-use, in order to simulate a $D$-dimensional spin system using only a ($D$-1)-dimensional subset of qubits along with an ancillary qubit register whose size scales logarithmically in the amount of entanglement present in the simulated state. Ground states can either be directly prepared from a known MPS representation, or obtained via a holographic variational quantum eigensolver (holoVQE). Dynamics of MPS under local Hamiltonians for time $t$ can also be simulated with an additional (multiplicative) ${\rm poly}(t)$ overhead in qubit resources. These techniques open the door to efficient quantum simulation of MPS with exponentially large bond-dimension, including ground-states of 2D and 3D systems, or thermalizing dynamics with rapid entanglement growth. As a demonstration of the potential resource savings, we implement a holoVQE simulation of the antiferromagnetic Heisenberg chain on a trapped-ion quantum computer, achieving within $10(3)\%$ of the exact ground-state energy of an infinite chain using only a pair of qubits.

We present a suite of "holographic" quantum algorithms for efficient ground-state preparation and dynamical evolution of correlated spin-systems, which require far-fewer qubits than the number of spins being simulated. The algorithms exploit the equivalence between matrix-product states (MPS) and quantum channels, along with partial measurement and qubit re-use, in order to simulate a D-dimensional spin system using only a (D-1)-dimensional subset of qubits along with an ancillary qubit register whose size scales logarithmically in the amount of entanglement present in the simulated state. Ground states can either be directly prepared from a known MPS representation, or obtained via a holographic variational quantum eigensolver (holoVQE). Dynamics of MPS under local Hamiltonians for time t can also be simulated with an additional (multiplicative) poly(t) overhead in qubit resources. These techniques open the door to efficient quantum simulation of MPS with exponentially large bond-dimension, including ground-states of 2D and 3D systems, or thermalizing dynamics with rapid entanglement growth. As a demonstration of the potential resource savings, we implement a holoVQE simulation of the antiferromagnetic Heisenberg chain on a trapped-ion quantum computer, achieving within 10(3)% of the exact ground-state energy of an infinite chain using only a pair of qubits.

I. INTRODUCTION
One of the most promising near-term applications of quantum computers is the simulation of correlated quantum systems in which entanglement plays a crucial role, for which accurate classical simulations are often intractable. Examples include predicting low-temperature properties of correlated materials [1], calculating reaction rates or photoabsorption spectra of large molecules [2], and simulating lattice gauge theories of particle physics [3]. While these applications have generated considerable excitement, it is far from clear how large and how accurate a quantum computer will need to be in order to address classically hard questions of practical scientific and technological relevance. Many problems of interest require extracting information about systems in the thermodynamic limit, which often requires finite-size scaling to be performed on simulation results obtained from systems with hundreds (if not many thousands) of spins. At present, there are no circuit-model quantum computers that can directly simulate spin systems of these sizes.
However, it is well known that system size alone does not determine the classical hardness of simulating a quantum system. While the system size N determines the Hilbert space dimension (D ∼ e N ), which in turn sets the classical complexity of simulating the system's wave function by brute force, the Hilbert space actually explored by physical systems is highly structured, enabling efficient parameterizations of physical wave functions. The study of tensor networks over the last few decades has brought this point into sharp focus: Tensor network simulations generally require resources that scale no worse than algebraically with the system size, and only suffer an exponential scaling with respect to the amount of entanglement, quantified as the bipartite entanglement entropy. This realization, and the scaling laws connecting entanglement entropy to equilibrium and nonequilibrium phases of matter, has made it possible to judge by inspection of general properties of a model-e.g. whether it is in equilibrium (and if so if it's at zero or finite temperature), its geometry, its spacial dimension, or its topological properties-whether it can be simulated efficiently on a classical computer or truly requires quantum resources to simulate. Situations in which the latter case is realized are compelling examples of hard (and practically relevant) problems for which quantum computing could provide a significant near-term benefit.
In general, the existence of a simple tensor-network representation for a state does not guarantee that properties of that state can be calculated efficiently, because the network may be difficult to contract. Important examples include when the size (bond-dimension) of the tensors needs to be extremely large to achieve a good approximation, as happens generically for matrix-product state (MPS) simulations of higher-dimensional systems or long-time evolution, or because the tensor network topology does not permit an efficient contraction [4]. In the last few years, several proposals have pointed out that near-term quantum computers may be capable of carrying out tensor-network calculations that are beyond the reach of classical computers [5][6][7][8], with a key insight that the size of that quantum computer can be far smaller than the physical system described by the tensor network. Very recently, these ideas have been exploited to provide variational energy estimates for the 2D Heisenberg model by simulating small quantum circuits [9].
In this manuscript we present a toolbox for constructing and time-evolving high-bond dimension MPS states on a quantum computer. We refer to these techniques as "holographic" [6,7] because they enable simulation of a D-dimensional system using only a (D-1)-dimensional cross-section's worth of qubits by simulating the transfer matrix for the MPS as a quantum channel [10]. The channel effectively moves along the MPS by one unit of distance per channel iteration. Operationally, a purified version of this channel is implemented via unitary operations between a cross-section of spins (physical qubits) and an ancillary quantum memory (bond qubits), followed by partial measurement of the physical qubits. Each iteration of the quantum channel moves one step along the stacking direction of the cross-sections, with physical qubits reset between iterations and reused without duplication, thereby trading spatial resources (qubit number) for time resources (circuit depth). We present a detailed description of these techniques, and benchmark them via classical simulations of algorithm performance on solvable spin-chains.
Next, we show that this representation can actually be made to work by implementing it on a Honeywell trapped-ion quantum computer, and using it to estimate the ground-state energy of the anti-ferromagnetic Heisenberg chain. A crucial technical ingredient required to perform holographic simulations is the ability to selectively measure and initialize a subset of qubits mid-circuit, without affecting the remaining qubits. The quantum charge-coupled device (QCCD) architecture [11], in which individual qubits can be dynamically positioned far from other qubits during the execution of a circuit, enables individual addressing (including gates, measurement, and state preparation) with extremely low cross talk, and is therefore very well suited for these types of algorithms.
Finally, we extend these ideas to the simulation of quench dynamics starting from a holographically generated MPS. Naively, one would expect the simulation of N initially correlated qubits evolving for a time t under a local Hamiltonian to require a circuit of width w = N and depth d ∼ poly(t). A constructive algorithm achieving such scaling for k-local Hamiltonians was found more than 20 years ago [12], and in the years since the dependence of the circuit depth d on t, N , and the error tolerance have all improved (see [13,14] and references therein). It might seem that the circuit width requirement w ∼ N is fundamental, or even tautological. While w = N is indeed required in certain worst-case scenarios, in this paper we explore how far fewer than N qubits suffice in many cases of practical interest. Consider for example a 1D system initially in a state with maximum bipartite entanglement entropy S: Time evolution of this state by a local Hamiltonian can be implemented by a circuit of width poly(t) + S and depth N × poly(t). In other words, the required number of qubits is determined by the evolution time (with a modest constant offset to accommodate the initial entanglement of the state), while the physical size of the system being simulated can be accommodated by increasing the depth of the circuit. Moreover, if the initial state has a finite correlation length ξ, the N → ∞ limit can be well approximated by a circuit of depth ∼ ξ × poly(t). This reshuffling of resources from circuit width to circuit depth is what we mean by the term "holographic".

II. HOLOGRAPHIC SIMULATION OF MATRIX PRODUCT STATES
The basis for our quantum simulation algorithms will be the MPS representation of quantum states, which provides an efficient compressed approximation of quantum states with less-than maximal entanglement. We will construct methods for simulating ground-state properties and quench dynamics of local lattice models, defined on a Hilbert space H that decomposes into a tensor product of "sites" as H = ⊗ i=1 H i . Here the sites are arranged along a 1D line with open boundary conditions and have finite local-dimension |H i | = Q, and by local we mean that interactions act on at most k sequential sites. Note that this class of systems includes not only 1D spin-chains, where each site simply represents a single spin, but also D-dimensional models which can be sliced into a 1D stack of (D−1)-dimensional cross-sections (e.g. for a D-dimensional cube = N 1/D ).
MPS describe quantum states on such 1D stacks using far-fewer parameters: ∼ O( Qχ 2 ), with χ the bond dimension-than the worst-case Q required to specify a generic state, and describe a very special class of states with entanglement entropy across any cut bounded above by log χ. For ground states of 1D or quasi-1D systems, MPS states can be implemented using classical resources that scale at worst polynomially in the desired accuracy and system size, and provide efficient classical algorithms for simulating low-dimensional ground-states and shorttime dynamics. However, classical MPS methods fail for 2D and 3D systems, and for longer-time dynamics where substantial amounts of entanglement have been generated. Nevertheless, many of these systems have far less than the maximal amount of entanglement obtained by random states, and the MPS description provides a dramatic compression. Our goal is to devise an efficient method to prepare and time evolve an MPS-representable state on a quantum computer with an economical use of qubit resources, which gives access to MPS with classically inaccessible bond-dimension, while still leveraging the economical MPS representation for states with lessthan maximal entanglement.
To set the stage for these algorithms, we first briefly review a few key properties of MPS that are essential to understanding how they can be represented holographically on a small quantum computer.
For each site j, V σj is a set of Q square matrices [15], or equivalently a rank-3 tensor with "physical" index σ j = 1 . . . Q and "bond" indices α, β = 1 . . . χ, where χ is referred to as the bond dimension. The χ-dimensional vectors L and R specify the left and right boundary conditions, respectively. The standard graphical representation of an MPS is shown in Fig. 1. An individual tensor is drawn as a box with a leg for each index, as in Fig. 1a. Joined legs imply contraction of the associated tensor indices, such that Fig. 1b gives the wave function components of |Ψ [the contractions are implied as matrix-matrix or matrix-vector multiplication in Eq. (1)].
One can imagine creating an -site MPS as a physical state of a quantum computer by letting an ancilla register (containing "bond qubits") interact unitarily and sequentially with physical registers [10], each representing one site of the MPS, as in Fig. 2(a-c). In these circuit diagrams and elsewhere, open circles denote initialization of a qubit (or register of qubits) to the |0 state. To implement this construction with unitary circuit elements, one must exploit the gauge redundancy of the MPS description [16] to place the MPS in right canonical form The relationship between RCF and unitary embedding will be discussed further in the next section, but for now we simply want to emphasize that the aforementioned MPS gauge redundancy ensures that we can, without loss of generality, assume that Eq. (2) holds for the state in Eq. (1). Note that this definition of RCF is slightly non-standard: the boundary tensors are typically ordinarily also canonical, however we choose this version of RCF because it simplifies much of what follows. While this choice imposes some limitations on how faithfully the right-boundary condition can be imposed in a quantum circuit, bulk physics will be unaffected Note also that for a normalized state |Ψ , the imposition of RCF implies that in general both L and R cannot be simultaneously normalized. Without loss of generality we take L to be normalized but not necessarily R.
Correlation functions of local operators, Ô iÔj , can be computed by contracting the tensor network depicted in Fig. 1c. An efficient method to contract such a network for a long chain is to first contract the physical legs of the tensors on each site to form transfer matrices on each site k = i, j (withV the complex conjugate of V ), and begin multiplying these transfer matrices from left to right. When site i is encountered, we apply the modified transfer matrix and similarly for site j. If we interpret the bond vector space as the Hilbert space of some quantum system, then the transfer matrix is a linear superoperator acting on bond-space density matrices ρ βδ (β, δ ∈ {1, 2, . . . χ}) via the mapping Together with the RCF conditions in Eq. (2), Eq. (5) establishes E as a quantum channel (trace-preserving completely positive map) on the bond space [16], with the MPS matrices V σ as the Krauss operators of the channel. In this language, the contraction depicted in Fig. 1c FIG. 3. A circuit implementing holographic state preparation and correlation function measurement of an MPS. In this simple example we use one system qubit (enabling simulation of an infinite 1D spin-1 /2 chain) and n b bond qubits for an MPS with bond-dimension χ = 2 n b . Open circles denote qubit initialization to the |0 state, and the measurements are implied to be in the eigenbasis of the operatorsÔn andÔn+r. Averaging over the outcome of this circuit (with the two measurement results multiplied together) produces the correlation function Ô nÔn+r .
can be expressed as the overlap of a "time-evolved" initial bond-space density matrix ρ i = |L⟫⟪L| (with ket | . . .⟫ indicating a state in the bond Hilbert space) with the final un-normalized state |R⟫, Indeed, the earliest comprehensive treatment of MPS in the literature (in which they were called finitelycorrelated states) defined them in terms of quantum channels [17].

B. Holographic MPS generation
Equation 6 demonstrates that correlation functions of an MPS can be encoded in the dynamics of a quantum system with size independent of ; the spatial structure of the physical Hilbert space has been converted into a (discrete) time direction of the bond Hilbert space. Thus one can simulate the state of a D-dimensional system using a system of dimension D − 1, inspiring the moniker "holographic" [5]. However, it is important to keep in mind that this dynamics is not unitary. The holographic algorithms described here can be viewed as explicit purifications of this non-unitary dynamics in the form of quantum circuits. Alternatively, one can understand these algorithms starting from the known representation of MPS as quantum circuits, previously described in Sec. II A and illustrated in Fig. 2. From this perspective the dimensional reduction can be understood by looking at the causal structure of that circuit and recognizing that the physical register corresponding to site j can be measured and reset before the bond qubit register interacts with the physical register of site j + 1, implying that only a single site worth of physical qubits is required to implement the entire circuit [8]. Despite the constant erasure of information in the physical qubits, long-range spatial correlations in the system are retained as memory in the bond qubits.
Generic versions of such holographic algorithms for 2D systems were previously outlined in Ref. [5], though with-out explicit discussion of connections to the MPS formalism, and the known representation of MPS as quantum channels. Subsequent work [6,7] revealed an intriguing element of noise resilience in holographic simulation techniques. Due to the repeated partial-measurement and reset in the holographic technique, errors do not propagate indefinitely as for purely unitary circuits. Consequently, a finite density of errors produces a finite imprecision on the measured correlation function, in contrast to a purely unitary circuit, for which a single error can spread and contaminate all outputs.
In what follows, we unify these perspectives with the framework of MPS, and develop concrete variational ground-state preparation and quantum dynamics simulation techniques using this framework. We begin with a detailed description of the holographic MPS preparation/measurement protocol, summarized in Fig. 3. This protocol utilizes a register of n b = log 2 χ "bond" qubits (representing the χ-dimensional bond Hilbert space) initialized in state |L⟫, and a register of n p = log 2 Q "physical"-qubits (representing the Q-dimensional physical Hilbert space of a single lattice site) prepared in a fixed reference state |0 . The channel E is realized by applying unitary gates between the physical qubits and bond qubits, and then tracing out (i.e. discarding) the physical qubit. Such a unitary purification of the MPSchannel can always be constructed via the Stinespring dilation. Specifically, one can embed each MPS tensor (V σ ) α,β as the columns of a unitary matrix U α,σ;β,σ with fixed index σ = 0, i.e. (V σ ) α,β = σ|⟪α|U |0 |β⟫. Since, in right canonical form, V forms an isometry from C χ → C χ ⊗ C Q , the columns with σ = 0 form an orthonormal set, which can always be completed into a full orthonormal basis for C χ ⊗ C Q to obtain U .
Any correlation function, C = Ψ| ⊗ iÔi |Ψ , of the corresponding MPS can be sampled by iterating this quantum channel to step through the sites of the chain from 1 to in the following sequence of steps: 1. State prep: Start with C = 1. Prepare the bondqubit register in a given state |L⟫, which sets the left-boundary condition for the MPS. Then, start-ing with site i = 1: 2. While i ≤ : Iteratively apply the quantum channel for the MPS matrix to step along the chain from site 1 toward site by the following steps: (a) Prepare the physical qubit register in a reference state, |0 . (b) Act on the physical qubits and bond qubits with a unitary circuit U , which is a purification of the MPS on-site tensor for site i. (c) Measure the physical qubits in the eigenbasis of operatorÔ i (which may be the identity operator on most sites, for which the measurement is unnecessary). Denote the eigenvalue ofÔ i corresponding to the measurement outcome by λ i , and multiply 3. Measure the bond-qubits in a basis containing |R⟫.
If we post-select on outcomes for which the bond qubits are found to be in state |R⟫, then the above algorithm samples the correlator C in the state |Ψ , and the expectation value can be estimated by averaging over sufficiently many repetitions to achieve a desired statistical precision [18]. Note that fixing the right-boundary condition through post selection incurs a multiplicative overhead ∼ χ, which will be very large in cases where quantum computation is required. However, since we are primarily interested in bulk properties, and since correlations decay exponentially in distance from the boundary, it is generally sufficient to skip this post-selection, which provides a weighted average of the correlation function over the right-boundary condition. If we are interested in boundary effects, for example when examining impurity models or boundary conformal field theories, one can study the left boundary (where the boundary condition is set by the state-preparation for the bond-register).

C. Holographic entanglement measurements
Measures of entanglement provide detailed insights into quantum many-body systems beyond what can be drawn from local correlation functions, revealing nonlocal correlations, topological and symmetry-protected topological orders [19][20][21] , and diagnosing thermalization, scrambling, and many-body localization [22]. Methods to measure Rényi entropies for a subsystem A, with integer index n, have been developed based on creating replica copies of the system and measuring operators that cyclically permute the quantum states of various copies [23], or by examining crosscorrelations in randomized measurements to virtually implement the desired replicas [24]. These methods can be directly adapted to holographically represented states by performing the desired measurements on the physical qubit registers, as described in the previous section. Such measurements would enable, for example, the estimation of free energy from holographically generated thermofield double states [25,26]. If one is interested in the entanglement entropy of a bipartition of the chain, it can be directly obtained from measurements of the bond-qubit register. Namely, since the holographic simulation method recreates an MPS in right canonical form, the entanglement spectrum for the physical qubits bipartitioned by cutting between sites j and j −1 is precisely equal to the spectrum of the density matrix for the bond-qubits after j iterations of the holographic simulation algorithm. In this case, the replica-SWAP or randomized-measurement techniques can be applied directly to the bond qubits, with a number of measurements that grows with χ, but not the interval size, offering a potential savings in measurement complexity. For replica-SWAP based measurements of entanglement entropy, this holographic method provides a potentially huge savings in qubits required, as one needs only to replicate the bond-qubits, without replicating an extensive number of physical qubits for every site in the chain.

D. Expressivity of holographic MPS
While a unitary circuit representing the (purified) quantum channel of any MPS is formally guaranteed to exist, the crux for practical use of holographic simulation techniques will be constructing effective methods for implementing channels for physically relevant systems using low-depth circuits. Namely, arbitrary unitary synthesis from a local gate set generically requires gate counts that scale exponentially with qubit number, and is clearly not a viable technique for large bond dimension. However, physical systems with local Hamiltonians are far from "generic", and have considerable structure that could be exploited for efficient simulation. This observation poses the following basic question: What class of quantum states can be efficiently holographically represented on a quantum computer with exponentially large χ ∼ 2 n b , but using low-depth [i.e. with poly(n b ) gates] quantum circuits? Though we cannot conclusively answer this general question, in the following we develop holographic algorithms for simulating non-equilibrium dynamics starting from correlated ground states, and provide numerical evidence demonstrating that low-depth circuits may suffice in many situations of practical interest.

III. HOLOGRAPHIC VARIATIONAL QUANTUM EIGENSOLVER (holoVQE)
A central task for quantum materials simulation is to accurately approximate the ground-state correlations of local Hamiltonians H = j h j , where each term h j acts on sites within a distance at most k from site j. Hybrid classical/quantum variational algorithms, like the variational quantum eigensolver (VQE) [27], offer promising methodologies for attacking this problem on moderate scale quantum computers. In VQE, one prepares a trial wave-function |ψ(θ) on a quantum computer by evolving a fixed initial state with a quantum circuit composed of gates parameterized by rotation angles θ ∈ R p (p being the number of variational parameters). The expectation value of the energy, E(θ) = j ψ(θ)|h j |ψ(θ) , is subsequently estimated by measuring each individual term in the sum to the desired precision. Then, a classical computer updates the parameters θ to lower the variational energy in order to find the best approximation of the true ground state within the family of states, |ψ(θ) .
The holographic representation of MPS on a quantum computer naturally suggests a holographic extension of VQE (holoVQE), which uses the MPS representation described in the previous section with the unitary U V represented by a parameterized circuit. The expectation value of energy can be computed by measuring each term in the Hamiltonian using the above-described procedure for measuring correlation functions. Then, the variational ansatz can be optimized by using a classical algorithm to minimize E(θ).
The implementation of holoVQE is simplified for crystalline materials with translation invariant Hamiltonians, where there are only a finite number of terms in the Hamiltonian that must be independently measured in the thermodynamic limit. The holographic method produces an MPS with open boundary conditions which is not translation invariant. However, in an MPS the boundary's influence decays exponentially with the correlation length ξ, so one can simply measure the distinct terms in the Hamiltonian a distance r > ∼ ξ from the boundary. In the holographic correspondence, we move a large distance into the spatial bulk by iterating the MPS quantum channel for a long time to burn in its steady state, as in Fig. 3.
In the following sections, we demonstrate a simple application of the holoVQE technique for approximating the ground-state energy of an XXZ chain using only a single ancillary bond qubit, and then physically implement this technique on a trapped ion quantum computer to analyze the SU(2)-symmetric Heisenberg point. We show that symmetry principles can be incorporated into the variational circuit ansatz to reduce the number of variational parameters and simplify the optimization. Next, we consider a model with lower symmetry: the transverse-field Ising model (TFIM). We apply a more generic circuit ansatz, and analyze the scaling of algorithm performance for ground-state preparation and reconstruction of critical correlation functions with increasing number of bond qubits. For small circuit sizes, we are able to reproduce the optimal MPS ansatz at bond dimension χ = 2 n b , obtaining relative accuracy on the (infinite size) critical ground-state energy below 10 −4 with only a few qubits.

A. Role of symmetries
The crux of effectively implementing a variational procedure is constructing a good variational ansatz. For present purposes, this entails identifying a parameterized unitary circuit family that effectively implements the purified transfer matrix of the desired MPS approximation to the ground state. Consideration of symmetries can guide the design of these circuits. Symmetries of the model can be strictly enforced on the variational states by restricting attention to symmetry preserving circuit families. Here, by symmetrypreserving circuit, we mean that we choose a particular linear representation of the symmetry action on the bond qubits (possibly including projective representations [28] when dealing with potential symmetry-protected topological states), and ensure that parameterizations of the variational circuit preserve the total symmetry quantum numbers of the physical and bond qubits together.
We note that it is not always desirable to explicitly enforce all symmetries of the model. For example, we may wish to assess whether or not the ground state spontaneously breaks symmetries, in which case one could build an ansatz around various possible symmetry-broken configurations. Moreover, it is often the case that, for a fixed bond dimension, the lowest energy state may not preserve the full set of symmetries possessed by the ground state. For example, when χ = 1 (no bond qubits) an MPS simply corresponds to a mean-field (best productstate) ansatz, for which energy minimization often yields symmetry-broken solutions. The examples below are indicative of these various possibilities.

B. XXZ chain
The XXZ spin chain is a canonical model of strongly correlated 1D systems, describing both one-dimensional quantum magnetism and superfluidity (by mapping the spins to hard-core bosons). For nearest-neighbor interactions the Hamiltonian is where (without loss of generality) we take J > 0. The model has a global U(1) symmetry, enlarged to a full SU(2) symmetry at the Heisenberg points ∆ = ±1. For ∆ < −1 [∆ > 1] the spectrum is gapped and the symmetry-broken ground state is ferromagnetically [antiferromagnetically] ordered, while for |∆| ≤ 1 the model is gapless and has no long-range spin order. For all values of ∆ the model is exactly solvable by Bethe ansatz [29].
In the antiferromagnetic phase (∆ ≥ 1) the meanfield solution is antiferromagnetically ordered, spontaneously breaking discrete-translational symmetry (and SU(2) symmetry for ∆ = 1). Since this state is consistent with the known value of S z = 0 for the true ground state, we can build a χ = 2 MPS by introducing a single bond qubit, and allowing it to interact with the system qubit via unitaries that conserve total S z (note that by breaking discrete translational symmetry, this χ = 2 MPS achieves the same energy as a χ = 4 translationally invariant MPS). Choosing we can now use a holographic representation of the MPS to measure the energy for a given choice of parameters (θ, φ), and minimize using a classical feedback loop. Using gradient descent and simulating 256 shots per energy measurement, we obtain the results shown for ∆ ≥ 1 in Fig. 4.
Additional care must be taken when computing the energy for 0 < ∆ < 1. In this case, the mean-field ground state breaks the U(1) symmetry of the model by spontaneously aligning (antiferromagnetically) along some direction in the XY plane. Since this state does not live in the correct symmetry sector with respect to U(1), it is not sufficient to restrict our attention to circuits that conserve total S z . In these cases, we find that the idea MPS at χ = 2 can be obtained using the three-parameter

C. Trapped ion implementation: Heisenberg chain
At the antiferromagnetic Heisenberg point, we implement the holoVQE procedure experimentally using Honeywell's QCCD trapped-ion quantum computer described in Ref. [11]. We utilize a subset of the 5 designated "gate zones" (orange/blue in Fig. 5a), which suffices to run two parallel instances of the holographic state preparation protocol with a single bond qubit for each [30]. At the Heisenberg point, it can be shown that the two-parameter ansatz described above for ∆ ≥ 1 is actually unnecessarily flexible, and it suffices to restrict our attention to φ = 0. Thus the entangling unitary between physical-and bond-qubit is The native two-qubit gate for our architecture (the Mølmer Sørensen gate [31]) is local-unitary-equivalent to a controlled-Z (cZ) gate, at least two of which are necessary to synthesize G θ for arbitrary θ (a minimal decomposition is shown in Fig. 5b). The holographic MPS circuit is then built by alternating applications of G θ andG θ = X p G θ (this alternation corresponds to starting with a classical antiferromagnet, as discussed above), with reinitialization of the physical qubit in between each entangler (see Fig. 5c). For the small bond dimension accessed in this example, we find that a "burn in" distance of 4 lattice sites is sufficient to approximate bulk expectation values to well within shot noise. Figure 5d shows the results of holoVQE. Starting with a randomly chosen parameter θ, we use gradient descent with derivatives estimated from finite differences of the measured energies E(θ), increasing the shot count for each energy measurement from 500 to 2000 as gradient descent proceeds. Error bars are 2σ confidence intervals obtained from a nonparametric bootstrap resampling of the data. Averaging over the final four data points (taken at the highest shot counts) we obtain an estimate of E = −1.59(5)J for the per-site ground-state energy. For comparison, the meanfield ground state-which lowers the energy as much as possible without entanglement-achieves E = −J (blue dot-dashed line in Fig. 5d). The optimal MPS with bond dimension χ = 2 n b = 2 achieves E ≈ −1.712J, which the experimental results would converge to in the absence of noise or other imperfections (purple dashed line in Fig. 5d), while the exact ground-state (χ = ∞) has E = J(1 − 4 log 2) ≈ −1.773J (from Bethe ansatz, black solid line in Fig. 5d). Note that this measured value provides a proper variational upper bound for the infinite chain, despite being obtained from a small quantum circuit. To highlight the resource savings of holoVQE, we note that achieving comparable accuracies using brute-force simulation of an L-site chain would require L = 6 (rather than 2) qubits to sufficiently suppress finite-size effects. If the circuit infidelities were reduced by either improving the gate fidelities or using error mitigation techniques, the minimum achievable energy with 2 qubits for holoVQE, E min ≈ −1.712J, would require 10 (perfect) qubits to achieve by brute force. In the previous examples, we used just a single bond qubit corresponding to MPS with bond-dimension two. These results already demonstrate the dramatic compression of resources enabled by the holographic simulation method in achieving reasonable accuracy on an infinite, critical spin-chain using only a pair of qubits. However, turning holoVQE into a useful algorithm requires a method to systematically improve the accuracy of the simulations. We now explore the performance of holoVQE upon including additional bond qubits, focusing on the task of ground-state energy estimation for the 1D transverse-field Ising model (TFIM), with Hamiltonian The TFIM exhibits a ground-state phase transition from an ordered (h < J) to a disordered (h > J) phase, both of which are gapped and can be well described by MPS of fixed (system-size independent) bond dimension. These phases are separated by a self-dual critical point at J = h, described by a conformal field theory (CFT) with central  (4) is decomposed into three native two-qubit gates and 8 single-qubit gates, with a total of 15 real variational parameters (as in Ref. [33]). (b) Energy obtained using "star" circuits with an increasing number of bond qubits. The agreement with exact MPS energy minimization at bond dimension χ = 2 n b is excellent, although we were not able to obtain reliable energies from this ansatz for n b ≥ 3.
charge c = 1/2, whose non-constant entanglement scaling requires a bond dimension that grows with system size as χ > ∼ L c/3 = L 1/6 to achieve asymptotically accurate correlations. Nevertheless, it turns out that the ground-state energy and moderate-range spin-correlation functions of this model can be captured with fairly high accuracy using modest bond dimension MPS, even at the critical point [32].
To explore the efficacy of holoVQE for this paradigmatic toy model, we numerically simulate the holoVQE procedure at the critical point (h = J) for a sequence of variationally parameterized circuits with a variable number of bond qubits. For n b bond-qubits, we construct a "star" circuit that involves only two-qubit gates that sequentially entangle the physical qubit and the j th bondqubit for j = 1, . . . n b , allowing each individual gate to be an arbitrary ∈ SU(4) two-qubit unitary, as in Fig. 6a. The primary challenge in these calculations is to reliably find the global minimum of a constrained non-linear optimization problem; we were only able to find reliable results for n b = 0, 1, 2, for which simulated annealing worked well. We note that classical MPS calculations also suffer from this challenge, which is typically overcome by breaking translational invariance in order to make the problem linear (at the cost of greatly expanding the parameter space). It is clear that scaling holoVQE to large circuits (and therefore large effective bond dimension) will require significant further development along these lines.
The results obtained by brute-force global optimization are shown in Fig. 6b, along with those obtained by an unconstrained MPS optimization using bond-dimension χ = 2 n b . Surprisingly, this simple circuit design finds the best possible MPS even for n b = 2, for which the parametrization is not exhaustive of all n b + 1 = 3 qubit unitaries. Because global optimization strategies did not yield reliable improvements for n b ≥ 3, we do not know if this feature is generic or restricted to small circuits.
Part of the challenge in achieving further improvements for the TFIM is the extremely rapid convergence of variational energy with bond-dimension to the exact ground-state energy. We note that, in an actual quantum computation, resolving very small energy differences will become impractical due to large statistical sampling overhead. This issue is especially pronounced in the TFIM, likely due to its integrability and small central charge, c = 1 2 (the smallest of any minimal model). In particular, assuming a near-optimal variational ansatz, the effective correlation length scales with the number of bond qubits like ξ ∼ n (3/c) log 2 b (at criticality), so that smaller c yields larger correlation length, and more rapid convergence of finite-range correlations with the addition of bond qubits.
Since the Hamiltonian is comprised of nearest-neighbor interactions, the holoVQE procedure only requires measurement of nearest-neighbor correlation functions. Once the circuit is optimized one can freely use the holographic state preparation subroutine to extract any desired observables in the variational solution. For example, in Fig. 7 we show the critical (J = h) transverse and longitudinal connected spin-spin correlation functions: where we've assumed translational invariance. The points are obtained from the holoVQE method, while the solid black lines are exact results from fermionization. The holoVQE results show clear signs of the universal scaling behavior for the Ising transition, C X (r) ∼ r −1 , and C Z (r) ∼ r −1/4 over moderate length scales (as much as 20 sites for the transverse correlations), despite using no more than three qubits.

IV. DYNAMICS
Calculating the dynamical properties of interacting quantum systems is an essential challenge for practical applications such as predicting chemical kinetics, computing non-equilibrium electronic and optical properties of quantum materials and devices, and analyzing NMR spectra [34]. Quantum dynamics also underpins foundational scientific questions ranging from the nature of thermalization and quantum chaos, to properties of quarkgluon plasmas in heavy-ion collisions, to understanding cosmological scenarios for defect production. Despite its fundamental importance, simulating quantum dynamics remains among the most challenging tasks for classical computers, generically requiring exponential classical resources even in low-dimensional systems where groundstates can be efficiently simulated. This area is therefore a promising candidate for achieving a practical quantum advantage on near-term quantum hardware.
In the following, we develop a quantum simulation algorithm that incorporates the holographic representation of MPS initial states described above to simulate timeevolution of an initial state under a quantum quench: Here, the initial state |ψ(0) is represented by an MPS (potentially with interesting correlations and entanglement), and H(t) = α h α is any geometrically-local time-dependent Hamiltonian, where each term h α acts on at most k adjacent physical sites. We dub this technique holographic quantum dynamics simulation (holo-QUADS). HoloQUADS enables a simulation of arbitrary time-ordered correlation functions using only ∼ poly(t) log Q + log(χ) qubits, which is independent of [35]. By comparison, the classical resources required to simulate time-evolution from MPS initial states generically scale exponentially with t due to rapid entanglement growth, even in 1D systems [36].

A. Holographic quantum dynamics simulation (holoQUADS)
If the dynamics we care about is naturally generated by some circuit of depth r, we can proceed immediately with a holographic circuit construction as detailed below. If we are concerned with continuous time evolution under a Hamiltonian H(t), the first step is to approximate the resulting unitary by a circuit consisting of r layers. There are many ways this approximation can be accomplished (see Refs. [13,14] for a helpful review of some state-ofthe-art Hamiltonian simulation techniques), with different techniques having different scalings of r with time, qubit number, and desired simulation accuracy. For our purposes, it suffices to note that algorithms exist for which, at fixed error, r scales no worse than ∼ ε t 1+ε for any ε > 0 [13,37]. We consider this scaling to be essentially linear in t and independent of . Taking a more pragmatic approach, we note that simple product formulas based on Trotterization generally require a number of layers that is far smaller than most rigorous bounds suggest, and-at least for local observables like 2-body correlation functions-produce accurate results using a modest (and system-size independent) number of time steps.
Next, we seek a holographic description of the state resulting from applying this discretized evolution to the initial MPS. To do so, it is useful to adopt space-timeinspired terminology to discuss different geometrical regions of the circuit. Denote the layer of the circuit by a discrete time index τ , and the position along the spin chain by x. Each wire in the circuit has an implied directionality, as shown by the arrows in Fig. 8. We define the past light cone of the point (x, τ ) to be the set of all points (x , τ ) from which one can arrive at (x, τ ) by flowing along the circuit in a forward direction, exiting gates along any outgoing wire. Inspection of Fig. 8 shows that measurements on the first k (k = 2 in this example) sites of the chain depend on the past light cone of the k th qubit (gray-shaded-region of Fig. 8). The circuit in this region can be implemented using only the physical qubits for the first r + k − 1 sites (with r scaling polynomially in t for fixed error) along with the log 2 (χ) bond qubits: First implement the unitary circuits to prepare the MPS state within the gray-shaded region at τ = 0, and then apply the layers of gates from τ = 0 to τ = r  Fig. 8 by attaching wires at exiting the top of the circuit to those at the bottom wherever the respective exiting qubits are reset and reused at the bottom of the circuit. The space direction wraps diagonally around the cylinder indefinitely (the site of the physical system at position x is labeled sx), and the circumference of the cylinder, which determines the required number of qubits, is determined by the evolution time. A time-ordered correlation function Ô (x1, t1)Ô(x2, t2) is obtained by measuring the corresponding operators at appropriate places in the circuit, as shown here (green circles representing measurements) for the example x1 = 2, x2 = 2r + 1.
that fall within this region. The first k physical qubits at τ = r can be measured in any desired basis, and then reset and reused to represent the next k physical qubits at τ = 0. These can be initialized into the correct state for the initial MPS by a horizontal (left-to-right) sequence of interactions with the bond qubits, and then propagated diagonally to τ = r by acting with the remaining gates lying within the past causal cone of the 2k th qubit at τ = r, labeled as the diagonal "slice 2" in Fig. 8.
Repeating this process, one implements the full timeevolution circuit from left to right by sequentially implementing left-facing diagonal slices of the circuit. Effectively the top of the circuit is being sheared and reattached to the bottom, such that the actual circuit fits naturally on the geometry of a cylinder as drawn in Fig. 9. By measuring the physical qubits at desired space-time points, one can reconstruct the time-ordered correlation functions of any local operators (relevant for dynamical response both near and far from equilibrium). Further, out-of-time ordered correlators that provide insight into thermalization, scrambling, and many-body quantum chaos can be simulated by including intervals of reverse-time evolution.
One can alternatively view this construction as a generalization of the holographic-MPS representation obtained by slicing the time-evolved-MPS circuit into diagonal slices whose boundaries are left-future-null-trajectories, and considering all qubit lines entering the slice from the left as "bond-qubits", and those exiting the slice vertically as physical qubits. With this interpretation, one can measure entanglement Renyi entropies of the physical chain by measuring the corresponding quantities of the bond-qubits as described above for general holographic MPS.

B. Comparison to classical methods
We now contrast this method with classical timeevolution techniques for simulating quantum dynamics. In 1D systems the leading classical method for simulating time-evolution of an MPS is time-evolving block decimation (TEBD). TEBD works by converting an infinitesimal time evolution step e −iH∆t ≈ (1 − iH∆t) into a matrix-product operator (MPO), applying the MPO to the initial MPS state, and reinterpreting the result as an MPS with bond space given by the tensor product of the MPO and MPS bond spaces.
At each stage, the tensors of the MPS are compressed (if possible) by discarding subleading singular values below the target accuracy. The compression step is effective if little entanglement is generated during the Trotterstep compared to the maximum possible, but does not save classical resources in cases where significant entanglement is generated during each Trotter-step, e.g. as in a strong quench with a thermalizing Hamiltonian. In fact, in extreme examples with maximally entangling dynamics, such as simulating stroboscopic dynamics of random circuits [38,39], no compression whatsoever can be achieved. In contrast, holoQUADS exhibits polynomial scaling of the required qubit resources with evolution time regardless of the amount of entanglement generated per Trotter step, and will exhibit a maximal advantage in cases where low-rank classical compression is ineffective, or in higher dimensions where bond dimension can be prohibitively high for classical simulation from the outset.
As an aside, we note that applying an MPO to an MPS (which forms the basis of many classical methods) cannot be directly implemented holographically by unitary circuits plus measurement, since application of an MPO does not generally preserve the right canonical form of an MPS. We leave as an open question for future work whether holoQUADS can be generalized to incorporate different, non-MPO-based quantum analogs of such compression schemes that operate efficiently on exponentially large bond spaces [40][41][42] to further save on qubit resources.

V. DISCUSSION
These holographic methods will provide a quantum advantage for situations where classical MPS techniques are intractable due to prohibitively high bond dimen-sion. Physically relevant examples include ground-state preparation of systems in dimensions D > 1, and timeevolution under thermalizing Hamiltonians. State-of-theart classical techniques run out of steam for 2D spinsystems of widths of around 10 spins (less for gapless systems), or 1D time evolution with thermalizing dynamics over a few tens of interaction times. Holographic quantum algorithms might provide quantum advantage on these tasks with as few as 30-40 qubits, which are employed to directly tackle the difficult, highly-entangled quantum aspects of these problems, rather than spending these resources to capture large sections of Hilbert space that are not accessed in physically relevant systems.
For higher-dimensional simulations, the classical simulation cost for MPS techniques grows exponentially in the system width, even for area-law entangled states. One could employ a holoVQE method analogous to higher dimensional DMRG, obtained by slicing the system into a 1D stack of (D−1)-dimensional cross-sectional slices, and specifying a circuit architecture to implement a quantum channel connecting one slice to another. The exponentially large bond dimension of the resulting MPS could be captured with polynomially many bond qubits. Importantly, for physically relevant low-energy states, MPS provide an exponential compression over a full wave-function description even in D > 1, so that the holographic representation affords substantial savings in qubit resources.
A more intrinsically higher-D generalization of these holographic methods would be to implement holographic simulation of isometric tensor networks [43]. These recently constructed higher-D generalizations of MPS capture a broad range of correlated states including (nonchiral) long-range entangled topological orders [44], and can be implemented straightforwardly as unitary circuits between physical qubits (now for each site in the lattice rather than each cross-sectional slice) and ancillary bond qubits. The isometric tensor networks invoke an explicit ordering of operations, and can be recast as a quantum channel that can be implemented holographically by resetting and reusing physical qubits that have already completed their participation in the circuit. The advantage of this technique over the usual boustrophedonic sweeping for 2D DMRG techniques is that it imposes a natural geometrically local 2D structure onto the physical and bond qubits.
It would also be desirable to extend these techniques to treat fermionic systems, using fermionic MPS representations [45], in order to simulate electronic materials. Other topics for future study include identifying effective heuristics for iteratively increasing bond-dimension to improve the holoVQE accuracy, and developing a formal and systematic understanding to the class of states that can be efficiently represented holographically using large qubit numbers, but reasonable circuit depths, to implement each MPS tensor.