Probing phases of quantum matter with an ion-trap tensor-network quantum eigensolver

Tensor-Network (TN) states are efficient parametric representations of ground states of local quantum Hamiltonians extensively used in numerical simulations. Here we encode a TN ansatz state directly into a quantum simulator, which can potentially offer an exponential advantage over purely numerical simulation. In particular, we demonstrate the optimization of a quantum-encoded TN ansatz state using a variational quantum eigensolver on an ion-trap quantum computer by preparing the ground states of the extended Su-Schrieffer-Heeger model. The generated states are characterized by estimating the topological invariants, verifying their topological order. Our TN encoding as a trapped ion circuit employs only single-site addressing optical pulses - the native operations naturally available on the platform. We reduce nearest-neighbor crosstalk by selecting different magnetic sublevels with well-separated transition frequencies to encode even and odd qubits.


I. INTRODUCTION
Quantum phases of matter that do not have a classical counterpart, known as exotic phases, are pivotal in the development of novel quantum materials and devices [1,2]. Recently, topologically ordered phases have drawn much interest as potential fault-tolerant components for quantum technologies [3][4][5]. Yet, our complete understanding of such collective quantum phenomena remains incomplete [6]. While analytical treatment supports only few (integrable) examples [7][8][9], numerical studies are limited to small systems due to the exponential growth of the state-space in many-body quantum problems.
Numerical approximations are thus used at intermediate and large systems. At these sizes, tensor-network (TN) methods provide a successful entanglement-based approach [10][11][12]. TN ansatz states can efficiently represent ground states of bulk-gapped Hamiltonians in one and higher dimensions [13][14][15], including topological insulators [16,17]. In one spatial dimension, TNs − and Matrix Product States (MPS) as special case − are efficient when evaluating observables. However, in higher dimensions, numerical resources for exact evaluations often increase exponentially in the system size. Various efficient numerical approximate techniques are known [18,19], although the precision of these approaches is often not tunable.
Quantum simulators [20] provide an alternative pathway towards understanding many-body quantum phenomena. Recent progress in quantum hardware [21,22] * These authors contributed equally to this work. enabled state-of-the-art experiments that are at the cusp of outperforming classical computers [4,[23][24][25][26][27][28][29][30][31][32][33][34]. However, the current Noisy Intermediate-Scale Quantum (NISQ) devices [35] are still fairly limited, either by size, controllability or noise. Quantum-classical algorithms, such as Variational Quantum Eigensolvers (VQE) [36][37][38][39][40], are protocols designed for NISQ devices [41] and, in the near-term, have the potential to outperform either purely quantum or classical approaches. Such hybrid methods aim to solve problems where implementing a given numerical task with specialized quantum resources can provide an advantage over the classical hardware. In this context, NISQ processors can be tailored to represent and sample TN states [42][43][44][45][46][47][48], e.g. by encoding specific classes of variational TNs as programmable quantum circuits. Said circuits thus allow the implementation of an efficient variational ansatz, used for instance in VQE [38,39], to prepare ground states of local gapped Hamiltonians on a quantum simulator. The optimization of the circuit parameters is performed via classical computation. Meanwhile, observables (including the energy cost function) are sampled directly on the quantum device, potentially providing an advantage over classical approaches.
In the present work, we utilize trapped ion quantum resources, sketched in Fig. 1(a), to implement a tensor-network-based VQE (TN-VQE). Its purpose is to prepare a many-body entangled ground state on the ion qubits. The ions are sequentially interacting with a common motional (phonon) mode, which acts as an entanglement carrier [49][50][51] or quantum data bus (QDB), as highlighted in Fig. 1(b). In our approach, the interactions are optimized by a classical algorithm with the ultimate goal to realize a target many-body quantum state with the phonon mode being disentangled from the ions at the end of the circuit [46].
As a relevant target problem, we use our TN-VQE to demonstrate topologically ordered ground-state phases in the non-integrable extended Su-Schrieffer-Heeger model (eSSH) [52,53], shown in Fig. 2(a). The target states allow an accurate, efficient representation as MPS since they are ground states of 1D gapped Hamiltonians [11,12,54]. We experimentally approximate the ground states of various model phases by running the TN-VQE algorithm. Finally, we measure Many-Body Topological Invariants [53] (MBTI) in the prepared states to detect phases and identify their topological order [55,56], see Fig. 2

II. TENSOR-NETWORK VARIATIONAL QUANTUM EIGENSOLVER WITH TRAPPED IONS
In this section, we briefly review the concept of the TN-VQE approach in trapped ions previously proposed in Ref. [46]. Further details on our experimental implementation are provided in Sec. IV. We consider a chain of N ions in a linear Paul trap as sketched in Fig.1(a).
For each site k, the qubit (|0 k , |1 k ) is encoded in a pair of electronic levels of the respective ion. The electronic levels can be controllably coupled to collective vibrational (phonon) modes of the ion chain. The entanglement among the ions is distributed by bosonic excitations of the phonon modes, which act as 'quantum data buses' (QDBs). In typical setups, a common approach to create entanglement relies on the simultaneous off-resonant coupling of multiple target ions to the QDBs via a bichromatic field [50,51]. In this approach, phonons are excited only virtually, which ensures that the QDBs are not entangled with the qubits. However, the underlying process is only at second order in the Lamb-Dike parameter η, which is typically small η 1 (Lamb-Dicke regime [57]). Here we consider elementary quantum processing resources that couple qubits and a single phonon mode as a first-order process in η, thus actively creating entanglement between ions and phonons. In practice, we employ controlled quantum dynamics generated by the Anti-Jaynes-Cummings Hamiltonian at any single site k with Rabi frequency Ω, in the Lamb-Dicke regime.
Here a † σ + k excites the qubit |0 k −→ σ + |0 k = |1 k and simultaneously creates a phonon |n −→ a † |n = √ n + 1 |n + 1 , while aσ − k does the opposite. For optical qubits, such dynamics are realized by driving with blue-detuned local laser pulses in resonance with the motional sideband of a selected phonon mode [58] as sketched in Fig. 1(a). Here, we consider the axial center-of-mass (COM) mode due to the near homogeneous coupling to all ions in the chain and the large frequency gap to higher order modes [59]. Originally, these sideband operations were proposed and used for implementing the C-NOT gate in trapped ions [49,60,61]. However, the accurate implementation requires precise calibration of the experimental setup as well as cooling of the phonon modes to zero temperature, as the effective Rabi frequency of the sideband operations heavily depends on the phonon population as ∼ √ n. Miscalibrations and any finite phonon population lead to imperfect gates and, as a result, residual entanglement between the QDB and the qubits at the end of a circuit. Recently, it was proposed [46] to tackle this problem by employing adaptive feedback-loop strategies, such as the Variational Quantum Eigensolver (VQE) [36][37][38][39][40]. VQE is a NISQ protocol which aims to prepare ground states of interacting Hamiltonians on programmable quantum hardware, while mitigating the imperfections. VQE can address Hamiltonians inaccessible for a given quantum platform and provide the best possible outcome for the available resources. However, identifying the best programmable resources in a specific quantum platform, for a given task, determines the efficiency of VQE. For ion traps, it was shown [46] that VQE with sideband operations can suppress the detrimental impact of the finite motional temperatures. The goal is to prepare the pure ground state of a given Hamiltonian H targ using the phonon mode as a data bus, which is required to be disentangled (only) at the end of the circuit. To this end, sideband operations are arranged to build a variational circuit U (θ) = l exp(−iθ l H k(l) ), where the l-th operation acts on the k(l)-th ion, and the parameters θ ≡ {θ l } are controlled by the durations of the laser pulses. For simplicity, let us consider the ideal case where the ions are initialized in some accessible pure state |ψ in , and the QDB is prepared at zero temperature |0 COM , with a|0 COM = 0. Unlike VQE in closed systems, the output variational states of the qubits are not restricted to pure states and are generally mixed due to residual ion-phonon entanglement at the end of the circuit. Nevertheless, the procedure to optimize the variational parameters is analogous to standard VQE. Namely, we experimentally measure the energy of the output state for the target model H targ. , by evaluating all contributing observables. As usual, several experimental shots are required to evaluate the energy expectation value within a desired errorbar [26]. The measured energy then acts as a cost function of a numerical variational optimizer, which iteratively proposes new parameter sets θ to find a global minimum in the landscape of the variational ansatz manifold. The optimization algorithm we employ is discussed in detail in Sec. V. Upon convergence to the ground state energy, for some optimal parameters θ opt. , the prepared state ρ out (θ opt. ) approaches the manifold of (nearly) degenerate ground states. If the ground state is unique, then ρ out (θ) = |ψ(θ opt. ) ψ(θ opt. )|, thus qubits are disentangled from the QDB at the end of the preparation. At this optimal scenario, our circuit ansatz guarantees [42] that the optimized output state can be represented as as sketched graphically in Fig. 1(c) [46]. Here, j k ∈ {0, 1} are indices of the canonical basis of the k-th qubit and A 1 are D k ×D k+1 parametric matrices as indicated in Fig. 1(c), with D 1 = D N +1 = 1 under open boundary conditions. States (4) are MPS [42,62], and the bond dimension D = max k D k is the refinement parameter of the ansatz. The entanglement entropy of the real-space bipartitions in MPS is bounded as s ≤ log D, matching the area law of entanglement for the ground states of gapped 1D Hamiltonians [63]. Indeed, it was shown that MPS provides an efficient and faithful ansatz to approximate such states [11,12,54]. In our circuit, D is controlled by the lattice size m of the sideband operation blocks (colored boxes in Fig. 1), such that D 2 m−1 .
in the bulk. In our circuit this can be achieved by repeatedly using the same subset of variational parameters along the circuit, with period : In the language of Fig. 1, boxes of the same colors use the same parameters. The optimal size of the edge blocks with the unique parameters can be identified variationally in the experiment, and it is expected to grow solely with the quantum correlation length of the target state. Consequently, the number of independent scalar parameters in the circuit ansatz does not ultimately grow for increasing system size.

III. TARGET MODEL
We use TN VQE in an ion trap to study phases in the interacting extension of the Su-Schrieffer-Heeger model [64,65], a one-dimensional spin-chain Hamiltonian capturing the transport properties of polymer molecules [66]. Open boundary conditions are required to exhibit topological order, and the Hamiltonian for N sites reads in dimensionless units, where σ µ k are the Pauli operators acting on qubit k. The coupling t − controls the 'staggerization' of the interaction strength, separating even-odd from odd-even pairs as sketched in Fig. 2(a). Additionally, δ defines the anisotropy of the XXZ-type interaction, and the 'standard' SSH model [52] coincides with the case δ = 0. For δ = {0, 1}, the model is not integrable, thus requiring numerical techniques to study its behavior. In the antiferromagnetic regime (δ ≥ 0), the model exhibits three different energy-gapped phases at zero temperature [53], as indicated in Fig. 2(b): (i) a nondegenerate trivial dimer phase, (ii) a four-fold (quasi-) degenerate Symmetry-Protected Topological (SPT) dimer phase, with soft excitations localized at the edges, (iii) a spontaneously symmetry-broken, Ising antiferromagnetic phase emergent at large δ.
The transition between symmetry protected topological phases can be detected using many-body topological invariants [55,56] associated with the symmetries protecting the corresponding phases of the 1D system. In case of the eSSH model, these symmetries [67,68] are the dihedral group of π-rotations about two orthogonal axes, reflection symmetry with respect to the center of the bond, or time-reversal symmetry. In present work, we consider the partial reflection and the partial timereversal MBTIs [69] corresponding to the two last sym- metries. Each of of the MBTIs can be used as a separate phase detector. The partial reflection MBTI reads [53,69] Here, ρ I , ρ I1 , and ρ I2 are reduced states of subsystems I = I 1 ∪ I 2 , where I includes an even number n < N of sites in the middle, and I 1 and I 2 are their left and right subsystems. The parity operatorR I =R † I =R −1 I reflects the 1D lattice around its center, practically swapping each symmetric pair of qubits. Analogously, the time-reversal MBTI [53,69] is where T 1 indicates the partial transpose operation on the partition I 1 , and u T = i∈I1 σ y i , so that the anti-unitary mapping ρ I → u T ρ T1 I u † T completely inverts the spin (thus time-reversal) of each qubit inside I. In the thermodynamic limit, with N n → ∞ (both N and n being multiples of four),Z R andZ T approach discrete values in the gapped phases, as demonstrated in Fig. 2(b) abs (c). Thus they are actual invariants under any fluctuation which does not disrupt phase properties. Precisely, they acquire the valueZ R =Z T = 1 in the trivial phase, −1 in the topological phase, and 0 in the symmetry-broken phase. The subsystem size n required to achieve the convergence depends on the correlation length in the system [53]. Thus, in the finite-size system, the measured MBTI become smoothed in proximity of phase boundaries, as shown in Fig. 2(c) forZ R , where we take N = 8 and n = 4.
Far from the phase boundaries, the eSSH ground states are gapped phases, satisfying the entanglement area law [14,15], and thus allowing an efficient MPS representation that can be realized with TN-VQE. Notably, blue-detuned sideband operations of Eq. (1) are a sufficient resource to prepare ground states of any real, zmagnetization conserving Hamiltonian [46]. Moreover, these operations take into account specific symmetry properties of H eSSH , eventually simplifying the optimization problem within TN-VQE. Precisely, due to its symmetries, H eSSH always exhibits at least one ground state with zero z-magnetization and all real amplitudes in the canonical basis. Accordingly, blue-detuned sideband operations realize real-valued variational unitaries U (θ) since H k are fully imaginary. Furthermore, they protect an 'extended' magnetization symmetry 1/2 k σ z k − a † a. Thus, if the initial qubit state has well-defined z-magnetization, so does the output qubit state once the COM mode has been variationally disentangled [46]. Thus, by protecting these two symmetries (magnetization, complex conjugation), which are present in the target model, we substantially simplify the variational optimization problem [26]. Conversely, we can not protect the SPTO-symmetries (reflection, time reversal, etc.): doing so would prevent us from establishing topological order from a trivial input state.

IV. EXPERIMENTAL IMPLEMENTATION
We implement TN-VQE on an ion trap quantum simulator, using laser-cooled 40 Ca + ions in vacuum [70]. All experimental results presented here are carried out on linear strings of eight ions. We encode the qubits in the optical transitions from the electronic ground state 4S1 /2 to the 3D5 /2 manifold, such that |0 = |4S1 /2 and |1 = |3D5 /2 . The 3D5 /2 state is long-lived with a life time of T 1 = 1.15 s. Initially, the system is prepared in the electronic and motional ground state by a sequence of Doppler-, polarization-gradient, and sideband-cooling and optical pumping [70]; we achieve a mean phonon numbern ≤ 0.05. The quantum state of each qubit is manipulated by a sequence of laser pulses. For each circuit we aim to implement the sideband operations Eq. (1) on a single qubit at a time, requiring precise control over the position of the laser beam [71,72]. In our setup, we address the ions individually via a high numerical aperture objective aligned at an angle of 22.5 • with respect to the ion string axis. Such single qubit operations will inevitably introduce cross-talk-errors in the non-addressed ions. Due to the geometry of the laser beam, a finite residual field will overlap with -most prominently -immediate neighbors, inducing excess rotations and phase errors [73]. Comparing the Rabi frequency Ω j of a target qubit on site j with its neighbors j ± 1 we measure the ratio Ωj±1 /Ωj to be as large as 5.9(5)% in the center of the ion string, which corresponds to a relative laser intensity of 0.35%. We seek to minimize the cross-talk errors arising during state preparation and measurement via two schemes. First, resonant errors on the |S ←→ |D carrier transition are suppressed by implementing a sequence of qubit rotations, decoupling the spectators from the target ion j. Here, any rotation R of the qubit by an angle θ on the Bloch sphere is decomposed into a set of single-qubit rotations where X j rotates the state vector about the x-axis and Z j about the z-axis. We implement X j by resonantly driving the carrier transition, while Z j is realized by a far-detuned laser pulse, inducing an AC Stark shift, in turn rotating the state of the qubit in the equatorial plane of the Bloch sphere. Since the AC Stark shift is proportional to the intensity of the laser rather than the amplitude as in the case of the resonant transition, the spectator ions become quadratically decoupled from the operation on the target qubit. Thus, resonant cross-talk can be attributed to the X operations in the sequence in Eq. (8) alone and as both of them are subjected to the same systematic errors, any resonant error is cancelled out by the opposing phase angles. With this method we manage to suppress resonant cross-talk to well below 1%.
However, this scheme is not applicable when driving sideband operations. First, the Stark shift experienced by the individual states depends on the phonon occupation n [74]. Second, the effective Rabi frequency scales as √ n, yielding a different evolution for each state. Here we implement a novel scheme using the internal structure of 40 Ca. We encode the qubit in different, well separated transitions in the 4S1 /2 and 3D5 /2 manifolds -specifically, we select two transitions with magnetic quantum numbers m = ± 1 /2 as shown in Fig. 3. Thus, for a given site k the encoding is defined as |0 = |4S, m = (−1) k 1 /2 and |1 = |3D, m = (−1) k 1 /2 , respectively. This method prevents any unwanted sideband excitation on neighboring ions. However, it introduces additional efforts during state preparation. For each ion required to be initialized in the |4S, m = + 1 /2 ground state two additional laser pulses would be needed. Conveniently, as our circuits start out from a Néel state |0101 . . . , we initialize |3D, m = + 1 /2 on all even sites where only a single pulse on the transition |4S, m = − 1 /2 ←→ |3D, m = + 1 /2 is required for each respective site. We quantify the benefits of this scheme via state tomography of a circuit with M = 14 sideband operations on four qubits using parameters optimized in numerical simulation, given in App. B, and compare the data with the state obtained by the circuit simulation. We achieve a fidelity of F = 89.7(11)% compared to 82.5(12)% in the case with all qubits encoded in the |4S, m = − 1 /2 ←→ |3D, m = − 1 /2 transition. Assuming the state fidelity is described by F = (F BSB ) M , the fidelity of a single sideband operation is given by F BSB = 99.23(7)%.
Apart from the state preparation and measurement, the full circuit is implemented using only single-ion sideband operations. In these operations, the phonon-ion coupling scales with the Lamb-Dicke parameter η, which heavily depends on the trap geometry, or more precisely, on the overlap of the incident laser beam and the trap axis. In our setup we measure η = 0.038, implying a requirement of up to 2mW of peak laser power to implement the sideband operations, which in turn induce AC stark shifts ∆ by coupling to the carrier transition with a strength on the order of ∆ ≈ 5 kHz. These shifts are actively compensated by a second, far detuned laser beam [75]. Since this requires substantial power in the compensation beam and with the total available laser power being finite, the achievable Rabi-frequency is also limited. We found the optimum to be at a sideband coupling strength of ηΩ = 2π · 8 kHz, such that a π rotation |0, n = 0 ←→ |1, n = 1 is performed in ≈ 125 µs.
In contrast to other entangling schemes like the Mølmer-Sørensen gate [50], the sideband operations executed in our circuits actively entangle the ions with the phonon mode. Uncontrollable interactions with the environment cause the qubits to depolarize, mainly due to heating and motional dephasing as a consequence of, predominantly, electric field noise in the trap [76]. In our setup, we measure the heating rate of Γ H = 27 (2) phonons per second in an 8 ion crystal. Furthermore, we measure the motional coherence time τ COM of the COM mode via Ramsey spectroscopy. For a single qubit we obtain τ COM = 101.9(1)ms, which is comparable to the coherence time of the laser given by T 2 = 107(15)ms. However, the motional coherence time will decrease with the number of ions in the trap. On an 8 ion string we measure τ COM = 21.9(2)ms, which is by a factor of 5 lower than T 2 . As such, we identify heating and motional dephasing as our main decoherence mechanisms. Consequently, in our setup, it is paramount, that each sideband sequence is finished well within the characteristic times τ H = 1 /ΓH and τ COM to ensure the faithful implementation of the target state.

V. VARIATIONAL OPTIMIZATION
The objective of our variational quantum eigensolver is to prepare the many-body ground state of the target Hamiltonian H eSSH of Eq. (5) via closed-loop optimization. For a given set of couplings t − ∈ [−1, 1] and δ > 0 the protocol proposes sets of trial parameters θ and evaluates the energy functional via data obtained directly on the ion trap quantum simulator. A simplified sketch highlighting the full VQE is shown in Fig. 4(a). For each parameter set θ we implement the circuit Fig. 1(b) and measure collectively in the X, Y and Z basis. The measurement outcome is then fed back to the classical computer, which evaluates Eq. (9) and makes an improved guess for new input parameters, thus iteratively minimizing E(θ).

A A A A A A A
The classical optimization algorithm that minimizes the energy functional over the parameters θ is based on the pattern search algorithm [77,78]. This local search algorithm moves around in parameter space by polling nearby points to a candidate solution θ c . The polling points are organized according to a stencil, centered at the candidate solution and comprised of orthogonal vectors in each of the possible search directions. Based on the experimentally measured cost function values, sampled at each of the polling points, the algorithm decides to move the stencil to a new candidate solution θ c+1 . If none of the polling points provided an improvement, the size of the stencil is decreased. Contrary, upon a successful energy lowering move, the stencil size is increased. The stencil is rotated such that the first polling vector is oriented along the direction of the last successful move.
Since the cost function landscape is only sampled through noisy projective measurements, some additions to the algorithm have been made. Firstly, a gaussian process model [79] is fitted to the data to provide a better estimate of the cost function in the neighbourhood of the current candidate point. This is to be compared with standard gradient based algorithms, that fit a linear model to the locally obtained function values and move the candidate solution accordingly. Here instead we fit a global model, taking into account all previous measurement outcomes. A second modification for dealing with noisy cost function values is the option for the algorithm to request additional samples at already polled points, in cases where the error bars on the energy estimates are deemed too large to be able to make a good decision on where the stencil should be moved. In this refinement step, elements of optimal computational budget allocation are employed [80].
We now briefly discuss details related to the implementation in the experiment. Each set of trial parameters θ is defined by 14 individual angles θ k in units of π, which are transpiled into a sequence of sideband operations adhering the proposed pattern in Fig. 1(b) [46]. However, the experimental setup imposes a lower limit on the θ k due to electronic components. Any angle yielding a laser pulse duration below 10µs cannot be faithfully implemented and is thus automatically dropped by the transpiler. With a target Rabi frequency of Ω ≈ 2π·8kHz this pulse length translates to an angle ≈ 0.03π. Thus, the actual circuit might differ from Fig. 1(c). Nevertheless, simulations show that such small angles have negligible impact on the implemented quantum state.
An example of the closed-loop optimization presenting real data for the parameters t − = 0 and δ = 0 -the critical point in the thermodynamic limit -is shown in Fig. 4(b). Starting with an initial guess θ = {0} the algorithm evaluates Eq. (9) for each optimization step. The solid line indicates the minimum for E(θ) in the experimental data achieved thus far. For comparison we show the value from numerical simulations of the circuit for each of the proposed parameter sets. We show our 'best guess' reaching E(θ opt. ) ≈ −8 compared to the exact value of E min. = −9.52.

VI. RESULTS
In this section we present the experimental verification of different phases in the eSSH model by means of manybody topological invariants as phase detectors. For each data set we first run the variational optimization algorithm and, upon convergence, obtain an optimal parameter set θ opt. . In our experiment, the quantum state of the bulk region -i.e., the four ion state ρ 3,4,5,6 in the center of the string, as shown in Fig. 2(d) -is measured via full state tomography and reconstructed using maximum-likelihood methods [81,82]. We calculate the MBTIsZ R andZ T according to Eqs. (6) and (7) from the reconstructed density matrices and calculated the measurement errors via bootstrapping. We also compare the experimental results with (i) the target values from the exact finite-size ground states, (ii) values from a ground state close to the thermodynamic limit, and (iii) values from the states obtained by numerically simulating the circuits with the experimentally obtained optimal parameter sets θ opt. .
First we discuss the case δ = 0, where Eq. (5) is equiv-  alent to the plain SSH model. In the thermodynamic limit, the transition between the topological and trivial phases occurs with increasing t − at t − = 0 and is indicated by the abrupt change of the MBTIs values from −1 to 1, as given by solid black lines in Fig. 5. However, in the finite-size ground states, the MBTIs, shown by dashed black lines, change smoothly. Moreover, the center of the slope is shifted from t − = 0 to the negative values due to finite-size effects, and specifically, that the number of odd-even terms in the eSSH Hamiltonian, Eq (5), dominates the number of even-odd terms by one.
The data obtained in our experiment shown by orange square markers is in good agreement with the described behavior and clearly demonstrates the transition of the MBTIs from the negative values in the topological phase to the positive in the trivial phases. We observe the most pronounced deviation of the experimental values from the target values and from the numerical simulation in the topological phase, especially with respect to the timereversal symmetry in Fig. 5(b). This observation is attributed to decoherence processes in the COM mode and residual Stark shifts, resulting in depolarization of the prepared states along the quantization axis of the qubits, which both MBTIs are sensitive to. A more detailed simulation analysis of this qualitative argument is included in Appendix A. Considering the example states shown in Fig. 2(d), it can be shown that both the MBTIs are quite fragile, under dephasing, in the topological phase, while they are relatively robust in the trivial phase. This effect is due to the algebraic properties of the MBTI itself, given that we consider a bulk of four sites. In fact, in the topological phase, losing coherences can causeZ R to increase all the way up to zero. Conversely, in the trivial phase,Z R can decrease only to 1/ √ 2. Such asymmetry is reflected in the quality of the MBTIs in either phase.
Having explored the integrable case, we now consider the case δ = 4, which exhibits the symmetry-broken phase (an Ising Antiferromagnet) between the topological and trivial phases due to the sharp anisotropy favoring interactions in the z-direction. This phase is indicated by the abrupt change of the MBTIs to 0 value, as shown in Fig. 6. However, the exact finite-size ground state indicates almost no presence of the symmetry-broken phase. Interestingly, the symmetry-broken phase can be clearly observed from the experimental data. We attribute this behavior to the fact that at a finite size, the exact model resolves the ground state degeneracy, and, in the finitesize unique ground state, the parity symmetry is not spontaneously broken. In contrast, due to imperfections and noise, experimental VQE polarizes into a physical ground state with spontaneous symmetry breaking, see middle panel in Fig. 2(d). This state is less entangled, has smaller finite-size effects, and is thus closer to the values of the thermodynamic limit. For more details, see Appendix A. Similar to the case δ = 0, we also observe deviation of the experimental values from the target values in the topological phase.

VII. CONCLUSIONS AND OUTLOOK
In the present work, we implemented a variational quantum eigensolver in an ion trap quantum device, capable of targeting tensor network states. We demonstrated that our technique is efficient at preparing entangled ground states of gapped Hamiltonians, including symmetry-protected topological phases. Our strategy encodes the tensor network variational ansatz in a quantum circuit which explicitly includes the COM vi-brational mode as an entanglement mediator. The native interactions of the ions with the COM mode were used in our circuit as variational resource operations, and the target pure state was approximated in ions via a variational quantum eigensolver. We carried out our experiments in traps with 8 ions, and successfully prepared each gapped phase of the non-integrable interacting extension of the SSH model at zero temperature.
We view our work as a first step to realize a scalable tensor-network simulator on an ion trap platform. The quality of our experiment has been improved by spectroscopic decoupling to suppress cross-talk between the ions -one of the major imperfections in ion traps. We observed that one of the limitations on the system size N in our circuit is finite coherence time τ COM of the COM mode. Since the coherence time of the COM mode scales roughly as 1/N and the depth of the TN circuit is proportional to N , one can estimate the size limit to N ∝ √ τ COM . Scalable simulation can be achieved with several sideband sequences separated by intermediate recooling of the phonon mode. In such a scheme, the ions must be variationally disentangled from the phonon mode at the end of each sequence. Additionally, implementing ion traps exhibiting local phonon modes will also circumvent this restriction [83,84]. Our TN circuit employed translation invariance in the bulk, thus we expect that the complexity of the optimization problem will not increase with the system size but only with the entanglement growth in the target state. Another challenge for large-scale simulation of many-body phases is their verification. In ion traps, this problem can be tackled by adapting randomized measurement techniques [53,[85][86][87]. A crucial step towards a useful quantum simulation is the extension of the TN-VQE towards 2D lattice systems [43,[88][89][90] and tensor network geometries for 2D, such as Projected Entangled-Pair States (PEPS) [91,92]. Precisely, recent work [90] discusses quantum circuits for sequential generation of plaquette-PEPS -a subclass of PEPS, that is believed to include a large class of 2D phases, including the topological ones [88]. In Appendix C, we speculate on a potential implementation of 2D TN states in scalable ion trap architectures that allow ion-crystal reconfiguration. In particular, we discuss a detailed schedule of operations for a 5x5 lattice in a well-tested microstructured ion trap [83]. Alternatively, we are confident that scalable realization of the TN circuit for 2D systems in simpler trap architectures can be achieved by implementing in-sequence projective measurements and reset of individual ions [93]. Indeed, assessing the optimal implementation scheme for a given trap architecture will require further research.
The demonstrated experimental capabilities open opportunities for ion trap implementation of a variety of protocols that make use of 1D variational tensor-network circuits. These include a protocol to study infinite 1D systems [47], imaginary and real time evolution [94], and quantum machine learning [95]. Finally, the resonant in-teractions with one or several phonon modes potentially can be used to construct circuits beyond the TN ansätze, e.g., to address problems in quantum chemistry [96,97]. The design of the appropriate variational circuit might be obtained in a closed-loop optimization on a quantum device itself using recent hybrid algorithms such as adaptive algorithm [98] or reinforcement learning [99]. In this section, we examine and compare different error models to explain the deviation between the experimental data and the predicted values forZ R/T in Fig. 2, Fig. 5 and Fig. 6 of the main text. We numerically simulate the variational circuit for δ = 4 with the experimentally obtained optimal parameter sets θ opt. , considering the full circuit as a sequence of operations on the blue sideband -errors are modelled to occur after each individual sideband operation. We quantify the performance of the considered error models via the residual sum of squares (RSS)

VIII. ACKNOWLEDGEMENTS
with i labelling each point obtained forZ R/T . From the following analysis, we conclude that, in our task, the initial temperature of the COM mode can significantly affect only the symmetry-broken phase. Because of the finite temperature, VQE can prepare a low-energy symmetry-broken state instead of the symmetric exact ground state. However, this imperfection is not relevant for investigating condensed matter phases. We identified the dominant sources of errors as the COM mode heating, fluctuations of the tip voltages of the trap electrodes, and imperfect compensation of the Stark shift, which effectively result in depolarization of the prepared states along the quantization axis of the qubits. Particularly, the coherence time of the COM mode should scales roughly as 1/N with N the number of ions. Since the We investigate the effect of finite temperature of the COM mode in the initial state caused by imperfect cooling. In Fig. 7 we show the results of the numerical simulation of the VQE where we consider the mean phonon numbersn = 0.00 andn = 0.05 in the thermal state of the COM mode. Here we do not include heating and dephasing, as these mechanisms are considered individually in the next section. For both values ofn we observe no shift ofZ R in the topological phase, which, in contrast, is present in our experimental data. However, in the symmetry-broken phase, unlike forn = 0, the case for n = 0.05 demonstrates results close to the experimentally obtained values. Forn = 0.05 we obtain a residual sum of squares RSS = 1.11 (96) forZ R and RSS = 1.92 (48) for Z T -when averaged, we get RSS = 1.52 (22)). In contrast,n = 0.00 yields larger deviations, namely RSS = 2.41(76) forZ R and RSS = 3.23 (46) for Z T (averaged RSS = 2.82 (61)).
We attribute the significant difference in the MBTIs for the casesn = 0.00 andn = 0.05 in the symmetrybroken phase to the ground state degeneracy. At finite size, the exact model resolves this degeneracy. In the unique finite-size ground state the parity symmetry is not spontaneously broken, as shown at the top left panel of Fig. 8, which gives the reduced states of 4 middle qubits in the ground state of the eSSH model at δ = 4 and t − = 0.333. In contrast, due to the finite temperature of the COM mode, the VQE algorithm polarizes into a physical ground state with spontaneous symmetry breaking as indicated by the top right panel in Fig. 8. This state is less entangled, has smaller finite-size effects and is thus closer to the valueZ R = 0 in the thermodynamic limit.
The fragility of the global energy minimum corresponding to the exact ground state of the eSSH model at δ = 4, t − = 0.333 is shown in Fig. (7). We simulated the variational circuit for a range of parameters θ(α) = (1 − α)θ opt. |n =0 + αθ opt. |n =0.05 , which interpo-lates between the optimal parameters θ opt. |n =0 found by VQE for the cases withn = 0 and the optimal parameters θ opt. |n =0.05 found forn = 0.05. For the respective states we calculate the relative energy error ∆ rel = [E targ − E(θ)]/E targ with respect to the ground state energy and compute the corresponding values of Z R . At α = 0 and α = 1 ∆ rel exhibits two local minima. Forn = 0 the global minimum corresponds to the symmetric state at α = 0, while atn = 0.05 we observe a substantial energy shift. In contrast, the symmetry-broken state at α = 1 demonstrates almost no energy shift when n increases, thus its energy becomes the global minimum forn = 0.05. This robustness is a result of the smaller amount of entanglement in the reduced bulk state compared to the exact ground state -this is highlighted by the top panels in Fig. 8).
Breaking of the symmetry can also be achieved in the target eSSH model by introducing sufficiently large pinned staggered magnetic field B on the outermost ions For B = −3 the exact ground state is in good agreement with the experimental data as shown in Fig. 7 by the red dashed line -this is also evident from the density matrices in Fig. 8 when directly comparing the top right and bottom panels. With increasing system size, the gap between the ground state and the first excited states decreases, in turn also decreasing the required pinning field B until it vanishes in the thermodynamic limit. Our analysis shows, that for fixed circuit parameters θ moderate temperature of the initial state of the COM mode does not affect the MBTIs significantly. Instead it can cause the VQE to obtain different optimal parameters θ opt , preparing a physical low-energy state with broken symmetry instead of the exact symmetric ground state in the symmetry-broken phase. This imperfection, however, is not important when investigating phases of quantum matter, since -in the thermodynamic limitthe same symmetry breaking occurs spontaneously regardless.

b. Heating of motional modes
Residual noise, most importantly due to electric field fluctuations in the ion trap, disturbs the motional state of the ion string, causing the string to heat up. In our setup, these perturbations occur at a rate Γ H = 27(1) phonons per second. This effect is modelled by the channel inducing an error with probability p = tΓ H on the state ρ depending on the time t 1 /ΓH. We simulate heating in the circuits with the initial states having different temperatures, namelyn = 0.00 andn = 0.05. To establish a notion of time, we assume each blue sideband operation implementing |S, n = 0 → |D, n = 1 to require a laser pulse of length t π = 125 µs. In Fig. 10(a/b) we observe a significant effect on the MBTIsZ R/T in the topological phase t − → −1, which is less pronounced in the trivial phase t − → +1. However, we see deviations inZ T , which tend to increase as we approach t − = 1. These can be attributed to different lengths T required to execute the pulse sequences as shown in Fig. 11. As T grows when approaching the edge cases t − → ±1, these sequences are more affected by heating effects. Considering a perfectly cooled initial state with a mean phonon numbern = 0, the residual sum of squares are given by RSS = 0.84 (21) forZ R and RSS = 0.45 (11) forZ T -when averaged, we obtain RSS = 0.65 (16). Increasing the initial temperature ton = 0.05 yields RSS = 0.71 (18) forZ R and RSS = 0.35 (10) for Z T (average RSS = 0.65 (16)). However, as discussed in the previous subsection, the temperature of the initial state does not affect the MBTIs significantly for fixed circuit parameters.

c. Depolariziation & weighted Pauli errors
Decoherence processes on the motional modes translate to errors on the addressed qubits. Such errors can be modelled by a depolarizing channel, where each spatial axis is equally affected by a product of single qubit Pauli errors σ (x,y,z) i , with probability p /3. In reality, it is unlikely that all errors are distributed equally -thus, we must consider probabilities for each axis p x , p y and p z , such that for a system of N qubits For simplicity we assume p x = p y = p xy -this assumption was shown to be valid in simulations, where exchanging p x for p y and vice versa yields identical results. We evaluate the model for bothZ R andZ T for all combinations p xy,z = {0, 0.001, 0.005, 0.01, 0.02, 0.03, 0.04} and estimate the 'goodness of fit' via the RSS to deduce the most likely source of error in our setup. Note, that this analysis also considers the case without any errors, i.e. p xy = p z = 0.
The results of our comparison are shown in Fig. 12(a/b) -the best fit is achieved, when p xy = 0 and p z = 0.03 for both MBTIs. We obtain a RSS = 0.17 (7) forZ R and RSS = 0.098 (35) inZ T (RSS = 0.13 (5) when averaged) and conclude, that the circuit is predominantly impacted by dephasing, resulting in Z-type errors. From a purely experimental point of view, this argument is expected for two reasons: first, even in the presence of laser intensity and laser phase noise we are able to control the rotation angles θ and φ of single qubit gates to an extend, such that we can achieve average gate fidelities of beyond 99.99%. Second, the VQE algorithm is supposed to take systematic overand underrotations into account when converging to an optimal set of circuit parameters θ opt. . The observed Z-errors are a consequence of various sources: we expect contributions due to heating and dephasing of the motional states as well as residual detuning from the first blue sideband, which is a consequence of imperfect Stark shift compensation.
We also consider a finite temperature of the initial state with mean phonon number ofn = 0.05 and evaluate the models performance. The results are shown in Fig. 13; we find a best fit for the probabilities p xy = 0 and p z = 0.03 in agreement with the data obtained for the ideally cooled initial state. The residual sum of squares is RSS = 0.15(7) forZ R and RSS = 0.081(32) for Figure 11. Lengths T of the individual laser pulse sequences required to implement the desired target state. T increases as t− approaches ±1, giving the system more time to be affected by heating. Figure 12. 'Goodness of fit' estimation via residual sum of squares for the weighted Pauli error model. We compare the model with the data forZR andZT individually in (a). As we expect X and Y errors to occur with the same probability, we assume px = py = pxy. The best fit is achieved for pxy = 0 and pz = 0.03 with RSS = 0.17 (7) inZR and RSS = 0.098 (35) inZT shown in (b), respectively.
Z T -when averaged, we find an overall RSS = 0.12 (5). This analysis yields a better agreement than the previous analysis forn = 0.00, however, we found the RSS to overlap, within their margin of errors.

d. Comparison of error models
To conclude this section we compare the individual models via their respective RSS in Tab I: heating of the phonon mode with a rate Γ H = 27(1) phonons per second and a model considering dephasing of the qubits with p z = 0.03, neglecting X and Y errors. We consider both cases of a perfectly cooled initial staten = 0.00 and a finite temperature staten = 0.05. As it is evident from both models, the initial temperature in the considered range has only a small effect. Forn = 0.05 we find an average RSS = 0.53 (14) assuming errors are only related to heating of the phonon mode, compared to RSS = 0.12 (5) for the weighted Pauli model. Dephasing of the qubits is attributed to various mechanisms: first, heating yields to the occupation of different motional states, each evolving with a different frequency depending on the number of phonons. Second, fluctuations of the tip voltages of the trap electrodes induce phases on the COM mode, in turn inducing a phase shift on the qubits upon each sideband operation. A third mechanism is imperfect compensation of the Stark shift, resulting in a small detuning from the first blue sideband of the COM mode. All of these sources manifest as Z errors. As increasing the system size N extends the length of the circuit while simultaneously decreasing the coherence time of the COM mode, these errors will become evermore relevant when adding more qubits. In contrast, systematic X and Y errors do not impact the measurement outcomes, as the VQE tunes the rotation angles to implement the desired operations.  We test our blue sideband sequences on the smallest instance of the SSH model, namely on a system of four qubits. The circuit follows the scheme shown in Fig. 1(a). However, the final 'box' θ D is executed right after θ C on qubits 3 and 4. We obtain the optimal parameter set θ opt. from theory; the 14 angles (in units of π) required to implement the circuit are given in Tab Table II. Rotation angles for the individual boxes of the circuit shown in Fig. 1(a) for the 4 qubit test bed system; all angles are given in units of π.
Appendix C: Proposed implementation of a 2D TN-VQE in a novel ion trap In this section we speculate on a potential preparation of plaquette-PEPS [90] with two-dimensional (2D) TN-VQE, using a novel quantum technology of linear ion trap as quantum hardware [83]. We demonstrate a minimal scheme, involving 25 ions, which can exploit bulktranslation invariance in the target state; the TN diagram is shown in Fig. 14(a). Similar to the MPS diagram in Fig. 1(c), the 'plaquettes' encode a sequence of unitary operations on local sets of ions. The architecture of choice is shown in Fig. 14(b) and has already been tested under experimental conditions [83]. It incorporates four 'storage' branches, which are used to load and to store ions. Via 'shuttling' zones the ions can be moved to the central 'quantum region', where they reside in a common potential and thus share a (local) phonon mode. Similar to the scheme in Fig. 1(a/b) we use single-ion sideband operations to implement a variational unitary. However, we need to consider controlled reconfigurations of ion string(s) in the QR to realize the different plaquettes; these reconfigurations are implemented by modulation of individual trap electrodes. We employ three operations, namely (1) 'shuttling' operations, where ions are moved in between different trap zones, (2) 'splits', which divide an ion string into substrings and (3) 'merges', where substrings are joined into one, larger string.
We will now show how the TN in Fig. 14(a) can be mapped to the chosen trap architecture. Considering the lower left corner of the TN diagram, we sketch the ion string reconfiguration operations required to realize the first three plaquettes A, B and C in detail in Fig. 14(c); optical single qubit sideband operations are not explicitly shown. The ions are represented by colored circles, where empty grey circuits are considered 'fresh' resources, i.e., ions that have not yet participated in the circuit. Solid grey circles indicate ions that have already been addressed, but must be kept in one of the memory regions to be used in later plaquettes. Finally, solid black circles show those ions, that have already 'dropped out' of the circuit; they are thus parked in the lower left branch of the trap. To better identify the movement of individual ions we assign tuples of integers (n, m), which correspond to the column and row of the particle in the TN diagram.
1. We implement the first plaquette A by shuttling four ions from the resource register to the QR region. Note, that the initial arrangement of the qubits is chosen as such, that the number of required reconfigurations required in later plaquettes are minimized.
2. Continuing with B, we can move qubit (1, 1) to the lower left branch -this ion will not participate in subsequent plaquettes. Ion (2, 1) is split from the substring (1, 2) & (2, 2) and shuttled to the left memory. Two more ions (1, 3) and (2, 3) join the QR and merge with the aforementioned substring.
In total, we need to implement 15 reconfiguration operations to implement the first three plaquettes; the individual operations are listed in Tab. III below. We extend our investigation to the full circuit, ending up with a total of 152 operations. However, the feasibility of such Step Shuttling Split Merge Total  1  1  1  0  2  2  3  4  1  7  3  4  3  2  7   Table III. Required number of shuttling, split and merge operations to realize each of the first three plaquettes of the TN diagram for the given trap configuration, see Fig. 14.
demanding circuits remains an open question. Recently a fault-tolerant parity readout scheme has been realized, requiring more than 40 split-and-merge operations and 110 shuttlings including state preparation and readout [100]. However, the trap geometry used in this experiment does not feature branches as memory regions; an implementation in the trap architecture described in this section might greatly simplify the circuit. Tensor-Network-VQE in a possible architecture of a 2D ion trap. (a) TN representation of a 2D circuit in a system of 5 × 5 ions (circles). Each plaquette couples four qubits to a local phonon mode via sideband operations similarly to Fig. 1. Plaquettes labeled by identical letters can share variational parameters to enforce the approximate translational invariance. (b) A suitable linear trap architecture [83] -the quantum region in the center is connected to four branches, which serve as loading and memory zones. Ions can be shuttled in between regions by precise control of local trap potentials. (c) Possible implementation of the 2D circuit via sideband and shuttling operations -we show the first three steps in the circuit with all participating ions shuttled to their respective positions. Here, empty grey circles indicate ions that have thus far not participated in the circuit, while those in solid grey have already been addressed and must be kept for later use in different plaquettes. Black circles highlight ions that have already dropped out of the circuit and remain from this point on untouched.