Variational quantum state preparation via quantum data buses

We propose a variational quantum algorithm to prepare ground states of 1D lattice quantum Hamiltonians specifically tailored for programmable quantum devices where interactions among qubits are mediated by Quantum Data Buses (QDB). For trapped ions with the axial Center-Of-Mass (COM) vibrational mode as single QDB, our scheme uses resonant sideband optical pulses as resource operations, which are potentially faster than off-resonant couplings and thus less prone to decoherence. The disentangling of the QDB from the qubits by the end of the state preparation comes as a byproduct of the variational optimization. We numerically simulate the ground state preparation for the Su-Schrieffer-Heeger model in ions and show that our strategy is scalable while being tolerant to finite temperatures of the COM mode.


I. INTRODUCTION
Realizing many-body quantum states on quantum devices offers an experimental pathway for studying the equilibrium properties of interacting lattice models [1][2][3], quench dynamics [4][5][6][7], or it can be viewed as a quantum resource, e.g., for sensing [8][9][10][11]. Adiabatic state preparation has been experimentally demonstrated, e.g., on trapped ions [12] and atoms in optical tweezers [13]. However, this strategy is currently limited by the finite coherence times of quantum devices, incompatible with the long annealing times required by the adiabatic criterion [14]. An alternative, promising pathway towards quantum ground state preparation is given by feedbackloop quantum algorithms, and especially by the Variational Quantum Eigensolver [15][16][17][18][19][20][21] (VQE). In VQE, parametrized − generally non-universal − quantum resource operations, available on a given quantum device, are arranged in a variational sequence, or circuit ansatz, to generate entangled trial states. The trial states approximate a target ground state as a result of a closedloop optimization of the variational parameters, where an energy cost function is minimized. For trapped ion qubits, VQE was recently demonstrated using as entangling resource either the Mølmer-Sørensen gate [22], or a programmable analog quantum simulator of the longrange XY spin model [23]. Both these entangling resources generate effective interactions between ion qubits by coupling them to vibrational modes (phonons). Implementing the optical couplings off-resonantly, one effectively achieves robustness of these resource operations to the temperatures of the phonon modes. However, this strategy yields relatively slow rates [24,25], limiting performance of VQE, as the quantum processing suffers from decoherence.
Here we propose a strategy to improve the efficiency of VQE on programmable quantum devices where interaction between qubits is mediated by auxiliary degrees of freedoms: Quantum Data Buses (QDBs). This strategy employs interactions between qubits and QDBs, naturally available in these devices, as variational resource operations in the quantum algorithm. We illustrate our approach on trapped ions, where resonant coupling of ion qubits to the axial Center-Of-Mass (COM) phonon mode [26] − acting as a single QDB − yields faster rates with respect to the off-resonant strategies, and is thus less prone to decoherence. Unlike logical gates designed with the resonant coupling, such as the C-phase gate [27], our variational approach is tolerant to finite temperatures of the COM mode, which do not limit the state prepara-tion fidelity. Moreover, our strategy does not require to disentangle the QDB after each resource operation, but only at the end of the state preparation, ultimately yielding faster processing. The task of disentangling the QDB from the qubits is carried out by the optimization algorithm itself: As the output qubit state approaches the unique ground state of a target non-degenerate Hamiltonian, the QDB becomes disentangled without ever being measured.
Via numerical simulations, we investigate the scaling and robustness properties of our strategy for 1D lattice models. As a benchmark goal, we aim to variationally prepare the ground state, |ψ targ , of the Su-Schrieffer-Heeger (SSH) Hamiltonian [28][29][30] on ion qubits in a linear trap. Using only blue-detuned sideband optical pulses [31] as resource operations, we design the variational circuit ansatz shown in Fig. 1, which can efficiently realize Matrix Product States (MPS). This 'QDB-MPS circuit' can incorporate various symmetries of the target model for enhanced performance, including approximate translational invariance in the bulk. In this paper, we will show that our approach is scalable, as we can approximate 1D ground states at saturating precision in the system size N ions , without increasing the number of variational parameters N p [see Fig. 2(a)]. Additionally, we compare the results from the QDB-MPS circuit with other VQE strategies, which employ closed-system (coherent) ion qubit resources, e.g., as in Ref. [23]. We demonstrate that our strategy can realize high-fidelity ground states even for QDB initialized at finite temperatures. Finally, we show that the QDB-MPS circuit can be up-scaled beyond the single-trap limit by being implemented in modular ion traps [32,33]. Overall, we consider our strategy a viable, efficient route towards quantum state preparation in ion traps as well as in other quantum platforms that rely on QDBs, such as atom qubits coupled via photons in waveguides [34,35] or cavities [36][37][38][39], or superconducting qubits coupled by microwave resonators [40][41][42][43].
The paper is organized as follows. In section II, we describe the ion trap resource operations which we employ in the QDB-MPS circuit. We also review the resources used by other VQE strategies in ion traps, which we consider for comparison. We briefly introduce the target SSH model in section III, and we discuss in detail the VQE strategies. We exhibit the results of our numerical simulation in section IV. Finally, in section V, we summarize our conclusions and argue future perspectives for our approach.

II. QUANTUM RESOURCES IN TRAPPED IONS
We consider a setup of N ions ions in a linear Paul trap as a relevant quantum platform. The system qubit, or spin 1/2, at site j is encoded by two internal electronic levels of ion j [31], labeled |↓ j and |↑ j . Each ion, driven by laser fields in the transition between levels |↓ j ↔ |↑ j , experiences a recoil associated with absorption or emission of photons. This mechanism couples the internal ion levels with vibrational (phonon) modes, either axial or radial, depending on the optical setup [31], and allows phonons to be used as bosonic QDBs for the ion qubits. In ion traps, as well as in other quantum hardware using QDBs [34-42, 44, 45], various strategies have been developed to construct entangling quantum operations between qubits, while ensuring that the QDBs are disentangled from the qubits by the end of each operation. In what follows, we review some of these strategies.

Resonant
Sideband Resources − One possible approach to exploit QDBs in quantum hardware is to transfer information from one qubit to another by populating QDBs with real excitations [27,41,42,44,45]. In ion traps, such transfer of the information can be realized by coupling an individual ion to the COM mode using sideband optical pulses [27] implemented in the Lamb-Dicke regime [31] (where the Lamb-Dicke parameter η is small, η 1). These pulses have frequencies detuned from the atomic transition exactly by the frequency of the COM mode ν. For the case considered further in this paper, it is sufficient to consider a resonant sideband with positive (blue) detuning, described by the Anti-Jaynes-Cummings where Ω is the Rabi frequency, and achievable rates are of the order ηΩ ∝ 25 kHz [46,47].
Here a is the destruction operator for the COM mode, and σ − j = (σ + j ) † = |↓ ↑| j is the lowering operator for ion j. In Eq. (1), we used a gauge transformation a → e iϕ a to highlight the anti-unitary symmetry (σ ± → σ ∓ , a → a † , t → −t) in the interaction Hamiltonian.
In Ref. [27], it was proposed to use resonant sideband pulses to realize the C-phase gate in ion traps. However, the main hindrance in using gates built out of resonant interactions with a bosonic QDB, as Eq. (1), is the requirement for the QDB to be perfectly initialized in a pure state, typically the vacuum |0 a , where index a refers to the QDB mode. This requirement stems from the incommensurate rates √ l + 1ηΩ of the transitions |↓ j |l a ↔ |↑ j |l + 1 a , depending on the number of excitations l stored in the QDB, as highlighted in Fig. 1(a). As a result, the gate fidelity of resonant sideband gates in ion traps decreases roughly linearly with temperature T of the initial thermal state ρ 0 (T ) = (1 − e −1/T ) exp(−a † a/T ) (in units of k B = = 1) of the COM mode.
Off-Resonant Sideband Resources − To address this hindrance, alternative strategies to realize entangling operations were developed, in which couplings of the qubits to QDBs are implemented off-resonantly, and the QDBs are only virtually populated [38,[48][49][50][51]. A prototypical example in ion traps, the Mølmer-Sørensen (MS) gate [48], can be implemented by coupling a subset S of ions to the COM mode by a bichromatic laser beam with two frequencies, detuned from the qubit frequency by ± (ν − δ) with δ ν. The resulting Hamiltonian is an all-to-all interaction among qubits in S, with σ x j = σ + j + σ − j . The condition δ ηΩ [52] ensures that the qubits interact with each other via second-order processes, populating the COM mode only virtually and keeping it disentangled from the qubits. This approach is insensitive to the temperature of the COM mode. However, its rates are intrinsically lower than the rates of the resonant sideband resources (1), (ηΩ) 2 /δ ηΩ. Analog Spin-model Resource − Recent experiment [23] showed that VQE in ion traps can also employ a programmable analog quantum simulator as entangling resource operation. Such analog quantum simulator is realized by using a bichromatic laser beam that offresonantly couples all the ions to all the transverse phonon modes, ultimately implementing the XY Hamiltonian [49,53] where σ {x,y,z} j are the Pauli operators acting on qubit j and B is a uniform magnetic field. Thanks to the off-resonant couplings, the interaction (3) is also insensitive to the temperature of the phonon modes. However, typical rates J 0 ∼ 0.05 − 0.5 kHz [24,25] are, again, much slower than the rates of resonant sideband interactions (1).

III. VARIATIONAL QUANTUM EIGENSOLVERS
VQE enables the experimental realization of ground states of arbitrary Hamiltonians, providing the best approximation for the given resource operations. Once the optimization is converged, the optimal quantum state can then be re-prepared on demand. For example, this state can be used as a resource to initialize the system for digital quantum simulation [54], which was recently demonstrated on ion traps in the context of lattice gauge theories [6].
Target model: Su-Schrieffer-Heeger − As an example, in this work, we consider VQE for the ground state preparation of the SSH model on ion qubits. The SSH is an exactly solvable Hamiltonian [28][29][30], and a free-fermion prototype for symmetry-protected topological order in 1D, which reads in open boundary conditions. H SSH is real in the canonical basis, exhibits z-magnetization symmetry [H SSH , j σ z j ] = 0, and has zero-magnetization ground states with real amplitudes. We added a small staggered magnetic fieldB = 0.1 at the boundaries, to lift the four-fold degeneracy arising in the topological phase (t < 0). With this prescription, H SSH is also CP symmetric (invariant under global spin-flip plus spatial reflection), translationally invariant in the bulk with a period of two sites, and has a critical point at t = 0.
VQE with Closed System Analog resources (CSA) − In Ref. [23], VQE was performed on trapped ions using Closed-System programmable Analog resources as in Eq. (3). We report the variational circuit ansatz considered in that work in Fig 2(b). This circuit alternates global resource operations, generated by H XY , and single-qubit rotations exp(−iθσ z j ). It was shown that this strategy is potentially scalable and capable of approximating ground states of 1D lattice models, including lattice gauge theories, with high fidelity. Here, we consider the CSA circuit as a comparison benchmark for the QDB-MPS circuit in preparing of the ground state of H SSH . We briefly review some details of the CSA strategy to highlight the differences with respect to the QDB-MPS.
The variational CSA circuit built with closed-system resource interactions H k ) |ψ in which are ideally pure. The system qubits are initialized in a Néel product state |ψ in = |↓ 1 |↑ 2 |↓ 3 ... |↑ Nions which can be prepared with high fidelity. For each given set of trial parameters θ, multiple instances of |ψ out (θ) are prepared, and projectively measured in various local bases to reconstruct the energy cost function H SSH θ = ψ out (θ)| H SSH |ψ out (θ) from one-and two-point correlation functions σ a i , σ a i σ a j . The cost function H SSH θ , function of θ, is evaluated at finite precision and minimized with a search algorithm capable of taking into account statistical errors. Upon convergence, H SSH θopt reaches its global minimum, and the optimal set of parameters θ opt can be used to re-prepare the (approximate) non-degenerate ground state |ψ out (θ opt ) on demand.
Convergence of the search algorithm can be enhanced by incorporating symmetries of the target model in the trial states. For the CSA circuit, the zero zmagnetization of the initial Néel state |ψ in is protected both by the entangler H XY and rotations exp(−iθ j σ z j ). CP symmetry is protected by H XY , and can be enforced on single-qubit transformation 'layers' by using correlated rotations angles θ j = −θ Nions+1−j . can be approximated as well, by repeating odd-site and even-site rotation angles, θ j = θ j+2 , away from the system edges [23].
While VQE provides the best approximate state preparation for a given set of resources in the presence of actual imperfections and noise [10], it is still affected by decoherence. By considering a different VQE circuit ansatz for trapped ions, which uses the high-rate resonant sideband resources of Eq. (1), we aim to tackle decoherence by speeding-up the preparation of trial states.
VQE with Quantum Data Buses (QDB-MPS)− We design a variational circuit ansatz U (θ) = k exp(−iθ k H j(k) ) built on ion-COM mode interaction pulses (1), according to Fig. 1(b). Since now the COM mode becomes entangled to the ion qubits, the state of the qubits is generally mixed, both during the processing and in the output trial state ρ out (θ) = Parameters {θ k } can be controlled by the duration of the sideband laser pulses. Similar to the closed-system case, the cost function H SSH θ = Tr[ρ out (θ)H SSH ] is reconstructed from correlation functions evaluated on the trial qubit state, while the output state of the COM mode is discarded. As H SSH θ is minimized to approach the ground energy, ρ out (θ) approximates the unique, pure ground state |ψ targ ψ targ |, and the output state of the COM mode becomes automatically disentangled.
In our design, resource operations are arranged such that the circuit ultimately generates a variational Matrix Product State (MPS). In fact, MPS efficiently approximate ground states of 1D quantum lattice Hamiltonians, including the SSH model [55,56]. In contrast to CSA resources, which are inherently global, blue-sideband operations (1) offer single-site addressability. This property allows us to progressively tailor the MPS from one edge to the other, as shown by the Tensor Network diagram [55,56] in Fig. 1(c), in analogy to quantum circuits developed in the context of quantum machine learning [57]. Previous strategies for realizing MPS on quantum devices were proposed, both for deterministic state preparation [58,59], as well as for variational state preparation [60] which employed long-range interactions. In contrast, our variational circuit requires only local ionphonon interactions, and, therefore, can be up-scaled to modular ion trap architectures (see next section). Additionally, we introduce an approximate bulk-translational invariance to the ansatz, to reduce its variational complexity: In our circuit, interactions in the bulk can be virtually grouped in 'boxes' (drawn in Fig. 1(b) as grey rectangles at the edges, and orange and green for the bulk). Each box is defined by a set of variational parameters, which are repeated in the bulk boxes of the same color. At the same time, the maximum achievable bond dimension D ∼ 2 l−1 of the generated MPS can be controlled by increasing the size of the bulk boxes l (see Appendix A).
As in the CSA case, we again incorporate symmetries of the target model in the trial states. First of all, since H j is fully imaginary, the operation U (θ) is real as well as |ψ in and ρ 0 (T ), and, thus, the generated states ρ out (θ) have real amplitudes. Moreover, operations (1) protect an 'extended' magnetization symmetry 1/2 j σ z j − a † a, forcing the output qubits state to have well-defined z-magnetization upon qubits and QDB becoming variationally disentangled (see Appendix B for details). We have analytically verified that, in the theoretical limit of infinite-depth circuits, the resources given by Eq. (1) are sufficient to generate ground states of any real, magnetization conserving Hamiltonian (proof in Appendix C). Controllability for an arbitrary lattice model (e.g., without symmetries) can be achieved by adding to the resources a resonant carrier operation σ z j for at least one ion.
We stress that resources of the form (1) are typical in various experimental platforms with bosonic QDBs, such as superconducting qubits connected via a resonator [41,42,45]. As such, the numerical analysis presented below can be generalized to these platforms.
CSD-MPS Circuit Ansatz − We also compare our QDB-based approach with a variational MPS circuit ansatz realized with Closed System Digital resource operations (CSD-MPS circuit), as presented in Fig 2(c). This circuit is built with single-qubit σ z i rotations and site-filtered MS gates (2) applied to a local set of ions. Such resource operations protect z-magnetization in the trial state.

IV. RESULTS
Performance− We compare the performance of the variational circuits by numerically simulating VQE for preparation of the ground state |ψ targ of Hamiltonian (4) in the gapped phase with t = 0.5. Trial states are obtained by evolving the initial quantum state in terms of density matrixes. For the QDB-MPS circuit, we consider the COM mode to be prepared in a thermal state ρ 0 (T ) with realistic average number of excitations n 0 (T ) = Tr[a † a ρ 0 (T )] = 1/(e 1/T − 1) = 0.01 [61,62]. Other sources of imperfections are not taken into account, and, therefore, in contrast to the QDB-MPS circuit, simulations of the CSA and CSD-MPS circuits are free of imperfections. We consider the theoretical limit of unlimited measurement budget [10,23], i.e., the cost function is the exact expectation value H SSH θ . To optimize parameters θ, we use a gradient-based optimization algorithm (see Appendix D), which can be adapted for noisy cost functions [63][64][65]. We emphasize that we consider a low number of free variational parameters, N p ∼ 10 − 25, which enables global search methods, such as the dividing rectangles algorithm [66,67], already adopted in VQE [23]. For the analog resource entangler H XY , realistic values α = 1.34 and B = 20 [23] are used. Further details on the circuit ansätze and their simulation are given in Appendix E.
To assess the accuracy of the ground state preparation, we employ various figures of merit. Besides the fidelity and the excitation energy (defined later on), we also consider a correlation-based error function calculated using (parallel) two-point correlators Here C ij is evaluated on the optimized state ρ out (θ opt ) and C DMRG ij is the exact correlator of |ψ targ obtained by DMRG method [68]. We consider this error function in Fig. 2(a), to show that the proposed QDB-MPS circuit is able to prepare states with F err saturating in N ions at fixed number of variational parameters N p . In particular, the insets show that these states have correlations C ij accurately capturing C DMRG ij . By contrast, F err arising from the CSA circuit grows rapidly for this specific task, and the correlators C ij decay slower than in |ψ targ . We argue that this follows from the necessity to variationally remove the strong (power-law) correlations established by the long-range interactions in H XY (3) to recover the exponential decay of correlations in the ground state. Moreover, the SSH ground state exhibits considerable (short-range) entanglement dimerization, which in turn seems to favor quasi-local resource circuits over global resource circuits. While the results for CSD-MPS and QDB-MPS circuits given in Fig 2(a) are compatible, the QDB-MPS is realized with higher-rate operations, and, therefore, is expected to be less prone to decoherence.
Scalability − We study the scalability of the circuits in terms of the excitation energy in units of the energy gap, which bounds infidelity of the prepared states 1 − F ≤ ε = ( H SSH θopt − E 0 )/∆ (Appendix F), and is thus a figure of accuracy. Here, F = ψ targ |ρ out (θ opt ) |ψ targ and ∆ = E 1 − E 0 is the energy gap obtained by DMRG, which converges to a finite non-zero value in the thermodynamical limit N ions → ∞ since the target model is non-critical. Figure 3(a) shows that, for the QDB-MPS and CSD-MPS circuits, ε grows linearly with N ions , in contrast to the CSA circuit. This is compatible with the previous observation that the two-point correlation functions carry, on average, an error saturating in N ions , and the energy functional contains a linear amount of such correlators.
Scaling of the accuracy ε with N p , for the various methods and at fixed N ions , is plotted in Fig. 3(b). Again, we observe that MPS circuits offer asymptotically better accuracy (or better parametric efficiency) than the CSA global resources. While increasing the number of parameters N p , we observe the CSD-MPS circuit reaches slightly better accuracies than the QDB-MPS circuit; we attribute this effect to the finite temperature of the COM mode we employed for the simulations. This argument is supported by the gray curves in Fig. 3(b), showing relative excitation energies ε q obtained with the same QDB-MPS circuit and the same parameter values as the black plot but with different initial pure states of the COM mode |q a for q ∈ {0, 1}. For n 0 = 0.01, we have ε = ε 0 p 0 + ε 1 p 1 + O(n 2 0 ), with p q = q| a ρ 0 (T )|q a ∼ n q 0 . The plot shows that the two quantum trajectories are both variationally optimized in such a way that their error contributions are commensurate, ε 0 p 0 ∼ ε 1 p 1 .
Tolerance to temperature − Here, we further investigate tolerance of the QDB-MPS approach to finite temperatures T of the COM mode. As T increases, more amplitudes p q are populated in ρ 0 (T ) = p q (T ) |q a q|,  and, therefore, more equations U (θ opt ) |ψ in |q a ∼ |ψ targ |q a need to be simultaneously satisfied, requiring higher N p to achieve the same accuracies. In Fig. 4, we consider the fidelity of the optimized output state of N ions = 6 qubits, as a function of the average number of excitations n 0 (T ) and number of variational parameters N p . Panel (a) shows that, for fixed N p , F decreases linearly with n 0 (similarly to the C-phase gate [27] built with sidebands), but with a slope that flattens for increasing N p . Therefore, to achieve a fidelity threshold, N p must be increased with n 0 . Remarkably, for the case under consideration, it seems that the N p required to obtain a static fidelity threshold scales logarithmically with n 0 . This behavior is illustrated in panel (b), where the density plot shows the infidelity of the optimized state as a function of n 0 and N p . The red curve highlights a fixed fidelity threshold 1 − F ≤ 0.002. This property sets our QDB-based VQE apart from the approaches using sideband interactions to construct logical gates [27], in which imperfect cooling imposes a hard bound on the attainable accuracies. Scalability on a modular ion trap architectures − The proposed QDB-MPS circuit 1(a) relies on interactions between individual ions and the COM mode. However, the number of ions which can be coupled to the COM mode in a single trap is limited by ∼ 10 − 100 [33]. To further scale the ion trap quantum processors, modular ion trap architectures were proposed [32,33]. Following these proposals, we describe ground state preparation with QDB-MPS circuit across multiple traps, as shown in Fig. 5. We consider coherent shuttling of few ions [32,69] (or their internal states [33]) as a quantum communication channel among the traps, and we assume it to be errorfree. We consider each trap containing a small number of ions, such that each ion is addressable via blue sideband pulses, realizing interactions as Eq. (1) with the COM mode of the corresponding trap. For simplicity, we assume that the initial thermal states, ρ 0 , of the COM modes are identical in each trap, as well as and their effective coupling to the ions.
The variational circuit ansatz is constructed sequentially along the network of traps, i.e., by applying operations in trap k + 1 after all operations in trap k are completed. Since each trap has an independent COM mode as a bosonic QDB, we introduce a variational 'interface' operation U I (θ 5 ), constructed with resource interactions (1), which is applied to a set of l − 1 ions (with l the size of the bulk boxes) as the last operation for trap k. This operation can be optimized in order to disentangle the COM mode of trap k, effectively implementing U I ρ qubits,k U † I = ρ qubits ⊗ ρ 0 , with ρ qubits,k the common state of ions and the COM mode. After transferring the set of l − 1 ions to trap k + 1, we apply the unitary U † I to the same ions plus the COM mode of trap k + 1 prepared in state ρ 0 , effectively reconstructing the state ρ qubits,k+1 = U † I (ρ qubits ⊗ ρ 0 )U I = ρ qubits,k . Ultimately, this operation removed the 'impurity' introduced by the interface between the two traps.
In Fig. 6, we compare the excitation energies, from optimizing the QDB-MPS circuit, with a single trap and with a modular trap network. Both circuits have the same number of variational operations in the boxes in the edges and in the boxes of the bulk, and the interface operation U I contains 2 variational sideband operations. The error emerging from the optimization of the ion-trap network deviates only slightly from the error in the single trap case, and we recover linear scaling of the error with the system size. This result is an even stronger indicator of scalability of the QDB-MPS approach for VQE, as it can treat hardware interfaces as impurities, and efficiently compensate for them with small variations of the circuit ansatz.

V. OUTLOOK
We proposed an approach for VQE which uses qubit-QDB interactions, available in currently-existing programmable quantum devices, as variational resource operations. In our strategy, the QDB is disentangled from the optimized states only at the end of the state preparation, as a byproduct of the optimization. For trapped ions, we discussed how to realize our strategy with resonant blue-sideband pulses, potentially yielding faster state-preparation rates compared to existing strategies, and thus helping to mitigate decoherence effects. The variational ansatz circuit we designed is tailored to efficiently prepare MPS, leading to a potential improvement in scalability over current VQE strategies.
We numerically simulated the ion trap implementation of VQE based on the QDB-MPS ansatz circuit, and showed strong evidence of scalability of our approach, even when it is realized in modular ion traps. While the ions-COM mode sideband interactions operate in an extended Hilbert space, we demonstrated that these resource operations can incorporate symmetries in the variational circuit, similar to the closed-system VQE approaches. Although here we demonstrated robustness of the method to the initial temperature of the COM mode, further investigations are required to study other realistic imperfections, including (i) qubits decoherence, (ii) fluctuations of the control parameters, and (iii) shot noise resulting from a finite budget of measurements to reconstruct the cost function.
The developed QDB-based VQE can be readily used in other currently available experimental platforms, such as quantum dots [70] or atoms coupled by (chiral) waveguides [71]. Moreover, we expect that the scaling and parametric efficiency of our results, found for the MPS circuit ansatz, extend to other efficiently-contractible tensor network states [57,60], such as MPS 2 [72,73], thus allowing us to explore higher dimensions. Finally, we envision the possibility of including controllable dissipative elements to the QDB resources [74] in order to construct optimized cooling algorithms beyond stabilizer codes [75,76].

APPENDIX Appendix A: Maximum bond dimension
In the section, we study the maximum bond dimension, D M , of MPS generated with the QDB-MPS circuit using a bosonic QDB and built with resource operations (1) described by the Anti-Jaynes-Cummings Hamiltonian. Here we distinguish two sets of MPS accessible by the circuit. The first set contains all pure states which can be generated by the QDB-MPS circuit. The states of this set can be generated probabilistically [58] after measuring (and therefore disentangling) the QDB at the end. It defines the search space for VQE and, therefore, is desired to be limited. The second set -a subset of the first set -contains deterministically generated translationally invariant MPS, which are disentangled from the QDB by the end of the circuit without measuring the EM. This set, in opposite to the first one, is required to be big enough to include a good approximation of the target state.

Set of probabilistically generated MPS
We start with the set of MPS generated probabilistically. We characterize this set of MPS by considering a possible range of their bond dimensions. Bond dimension D M of MPS of size N is the maximum of bond dimensions D n between two blocks of qubits [1, n] and [n + 1, N ]. D n corresponds to entanglement entropy between the blocks [56], such that the Schmidt decomposition of MPS |Φ in the basis of two blocks {|φ The entanglement entropy of the blocks is S [1,n] = −Tr ρ [1,n] log(ρ [1,n] ) ≤ logD n , with ρ [1,n] the state of qubits [1, n].
Using singular value decomposition, one can show that, for an arbitrary state of size N , it is valid that D n ≤ 2 min(n,N +1−n) , and, therefore, D M ≤ 2 N/2 . By using a not universal set of resource operations (B1) we limit D M and, therefore, restrict the set of achievable states. As presented in Appendix B, interactions (1) preserve an extended symmetry Z = 1/2 N j=1 σ z j − a † a. In the following, we explain that this symmetry imposes D n n/2 · 2 l−1 , where l is the size of the bulk boxes in the QDB-MPS circuit, and limits the maximum bond dimension as D M 2/3 · N 2 l−1 . Therefore, in comparison with the general case, the bond dimension of the MPS generated by the QDB-MPS circuit is exponentially suppressed. Figure 7, demonstrates entanglement distribution in the QDB-MPS circuit between the qubits by a subsystem, highlighted by orange color, which includes the QDB and l − 1 qubits shared by neighbor virtual boxes. Such a subsystem can be considered as a virtual ancilla that sequentially interacts with the qubits distributing entanglement among them [58]. The effective dimension of this ancilla (dimension of the populated Hilbert space) bounds bond dimension D n .
Let us consider the QDB initially prepared in state |p 0 a = |0 a , and the qubits prepared in state |ψ in = |↓ 1 |↑ 2 |↓ 3 ... |↑ N . Then, because of the preserved symmetry, the maximum possible population of the QDB after interaction with n ions can not be higher than | n/2 a . Multiplying the growing dimension of the QDB by the dimension of l − 1 qubits in the ancilla, we obtain that D n n/2 ·2 l−1 . From the other side, as explained in [58], D n has to decrease as n approaches to N to guarantee that the ancilla (with QDB measured in |p N a ) is decoupled by the end of the state preparation. Therefore, there exists some n M , for which D n M achieves its maximum. To obtain n M , we study how D n can decrease while the ancilla interacts with the qubits after qubit n M . Intuitively, this can be done other way around by considering the growth of bond dimensionD N −n if MPS would be generated from the other side of the qubits chain, starting from n = N . Since the measured output state |p N a can contain number of excitations p N ≥ 0, the number of excitations in the QDB can not only increase but also decrease as the QDB interacts with the qubits, populating Hilbert space with the twice bigger dimension than if the QDB would be initially prepared in |0 a . Therefore,D N −n ≈ 2Dñ ñ=N −n (N − n)2 l−1 . Solving D n M =D N −n M we obtain that, after n M 2/3 · N , D n has to decrease with n. Thus, we find that D M = D n M 2/3 · N 2 l−1 .

Set of deterministically generated translationally invariant MPS
In order to generate MPS with translationally invariant D n in the bulk, the virtual ancilla has to own dimension which repeats itself with a target period along the state preparation. This includes the ancilla dimension after realization of the last box in the bulk, after which the edge box is aimed to disentangle the QDB from l − 1 ancilla qubits as required by deterministic preparation of a pure state. This can be done only if the effective dimension of the ancilla does not exceed dimension 2 l−1 of the qubits. Hence, the limit of the bond dimension of the translationally invariant MPS deterministically generated in the considered QDB-MPS circuit is ∼ 2 l . Here, we defined the limit not exactly, since the actual bond dimension can vary within one period of the bulk.

Appendix B: Symmetries
In this section, we discuss the symmetries protected by the operations withΩ the Rabi frequency, a is the destruction operator for the bosonic QDB, and σ − j = (σ + j ) † = |↓ ↑| j the lowering operator for qubit j, with {|↓ j , |↑ j } the internal logical states of the qubit.
Magnetization − A number symmetry, forming an Abelian U (1) group of symmetric transformations W (ϕ) = exp(−iϕZ), is generated by the integer operator where N is the number of qubits, M = 1/2 j σ z j , and σ z j = |↑ ↑| j − |↓ ↓| j . It is straightforward to check that Z is a symmetry since [H j , Z] = 0 for every site j ∈ {1..N }. Z thus defines a quantum number, which is conserved during the controlled variational dynamics. Consider the initial state of the system of qubits plus the QDB as |Ψ in = |ψ in ⊗ |s a , with |ψ in = |↓ 1 |↑ 2 |↓ 3 ... |↑ N the state of the qubits and |s a the Fock state of the QDB, a † a|s a = s|s a , s ∈ N. This state has well-defined quantum number Z |Ψ in = s |Ψ in . Due to the symmetry of Eq. (B2), the evolution of this state will preserve its quantum number, and the state will always be of the form |Ψ out = q,m,t δ m,q−s c m,t |m, t m ⊗ |q a , where the qubit states are labeled by the magnetization m (eigenvalues of M ) and degeneracy t m . At the end of the variational quantum state preparation, after optimization of the process, the output state will be disentangled between phonon and qubits, and reads |Ψ opt = ( t c m,t |m, t m ) ⊗ |m + s a . The final qubit state will thus have a well-defined magnetization m ∈ N. This requirement identifies the class of target Hamiltonians H targetable with the chosen resources: which must contain the magnetization symmetry: [H, M z ] = 0.
Conjugation − A second symmetry is nested within our controls, and it becomes evident after gauging the resource Hamiltonians as in Eq. (1). In this format, the resource Hamiltonians are purely imaginary H j = −H j when expressed in the canonical basis, and thus antisymmetric H j = −H t j and traceless Tr[H j ] = 0. It follows that the unitary operators U j (t) = e −iHj t/ of time evolution under H j must be real, as they are imaginary exponentials of imaginary matrices, and have determinant det(U j ) = e −itTr[Hj ] = 1. This means that, with the chosen resources, we can only perform transformations in the special orthogonal group SO(2 d ), a subgroup of the unitaries U(2 d ), which, in turn, limits the set of states we can prepare to real states (vectors with real coefficients in the canonical basis). In conclusion, the class of models targetable with our variational scheme is restricted to models expressed by real, or symmetric, Hamiltonians, which always possess real eigenbases. We stress that the complex conjugation symmetry can be relaxed by adding a detuning δ j = ν ± to one or more of the control lasers, which can be helpful if we aim to quantum simulate a non-real Hamiltonian H T .

Appendix C: Controllability
We now discuss the controllability problem and show that, within the constraints set in the symmetry section, the selected controls are theoretically able to prepare any quantum state, in the limit of infinite available resources. First of all, we argue that starting from the state |ψ in ⊗ |q a it is possible to shift the system qubits state to any magnetization sector. This is shown recursively, by highlighting that the state |↓ j |q a under the action of H j for a time equal to t = π/(Ω √ q + 1) is completely mapped into the state |↑ j |q + 1 a , ultimately increasing the magnetization by +1. Secondly, we show that our controls H j provide strong controllability within each magnetization sector m (of dimension d m ) of real states, meaning that the Lie algebra generated by {H j } contains the Lie algebra of imaginary, hermitian matrices onto the m sector, generators of SO(d m ), as we will show later in this section. From this observation, it follows that any two real states |ψ 1 |q a and |ψ 2 |q a can be dynamically connected, i.e., |ψ 2 = e −iH t/ |ψ 1 for some t, simply because the effective Hamiltonian H = i|ψ 1 ψ |−i|ψ ψ 1 |, with |ψ ∝ |ψ 2 −|ψ 1 ψ 1 |ψ 2 , can be generated with our resources.
We now prove that, by means of linear combination and commutation, we can generate the whole algebra of imaginary Hermitian matrices for the sector m by starting from the bare resources H j = i(σ + j a † − σ − j a).
We start by considering the commutator i[H j , H j ] = i(σ − j σ + j − σ + j σ − j ), for sites j = j , which is an endomorphism of magnetization sector m. Moreover, it is easy to show that i[H j , i[H j , H j ]] = σ z j H j , for j = j . By recursion and locality, it follows that we can generate any string operator of the form for any binary string q k . We run one additional commu- where all sites are different. By recursion, we conclude that we can generate operatorsH which are tensor product of one operator at a pair of sites j = j , any amount of σ z , and any amount of (σ − k σ + k + σ + k σ − k ) (each one acting on separate sites). We remark that these operators form a basis for the imaginary Hermitian endomorphisms of the symmetry sector m. To see this, consider the canonical basis i(|q 1 q 2 q 3 . . . r 1 r 2 r 3 . . . | − H.c.), where q and r are any two-bit strings with total magnetization m. This operator can be easily decomposed in operators of the formH: for bits j that do not flip (q j = r j = ±1) just insert a (1 ± σ z j ) in the tensor product string, while for pairs of bits that flip-flop, insert a (σ − . This concludes the proof.

Appendix D: Parameters optimization
In the section, we describe the method used for the optimization of variational parameters θ in the numerical simulations of VQE with the circuits given in the main text. Since the purpose of the paper is to study the achievable accuracy of the proposed approaches, we employ the exact value of cost function H SSH θ in the optimization. The expectation value H SSH θ is obtained by simulating the circuits in terms of density matrixes without considering shot noise caused by the finite number of measurements of the trial ions states. To optimize parameters θ we use a gradient-based optimization algorithm, which relies on numerically accessible ap- To be more precise, we use the basin-hopping method [77,78], which is a two-phase method that combines a global stepping algorithm with local minimization at each step. As a local minimization algorithm, we use the quasi-newton method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS) [79][80][81][82].
To decrease noise in the data presented in the result section, we ensured that H SSH θ is monotonically decreases with decrease of temperature T of the COM mode (or equivalently number of excitations n 0 ), with decrease of state size N ions , and with growth of the number of variational parameters N p . We begin by independently optimizing parameters for all data points (defined by one or two values of n 0 , N ions , and N p ), each of which starts from several sets of initial parameters θ with values ∼ 0.001 − 0.1. Then, the optimized parameters of the best result in one data point is used as the initial parameters set for the neighbor data point, and the obtained result is compared with the existing ones. In the case when the neighbor data point requires less initial parameters, the extra parameters are removed from the inserted parameters set; otherwise, parameters with 0 can be inserted. To obtain the presented data, we iteratively re-optimize the parameters from one side of data points range to another and, afterward, in the reverse direction until convergence.

Appendix E: Gates sequence in the circuit
In the section, we consider in detail the sequences of variational operations, Fig. (8), used to construct circuits considered in our work. The number of variational operations in the edge boxes of the QDB-MPS and CSD-MPS circuits is fixed to 2 and 3 correspondingly, as shown in Fig. (8). In these circuits, the number of variational parameters is increased by progressively inserting operations in the first bulk box (orange) or second (green), alternately. In the QDB-MPS circuit, we increase this number by 1, arranging the operations such that the last operation in one box and the first in the next box do not act to the same qubit. In the CSD-MPS circuit, we add a single layer, consisting of a local MS operation and σ z i rotations at each ion. To increase the number of variational operations in the CSA circuit, we add a single layer consisting of H XY operations and σ z i per ion. The trial states in the CSA circuit are obtained as |ψ out (θ) = i exp{−iθ i H i } the resource operations. Afterward, the required correlation functions can be obtained from |ψ out (θ) . Since this method requires numerical operation with the full quantum state, the number of the qubits in the CSA circuit in our study is limited by 12. On the contrary, one can see from Fig. 8 that, in order to obtain correlation functions of the states generated by the QDB-MPS and CSD-MPS, the full state is not required to be kept, but only a reduced density matrix of at maximum l qubits plus the QDB, with l the size of the bulk boxes. The qubits which do not more participate in the evolution until the end of the circuit can be measured in the required bases, and the next qubit can be added to the remained qubits state by using the Kronecker product operation.

Appendix F: Energy bounds on fidelity and purity
In this section, we show that the final energy H = Tr[Hρ(θ)], in respect to a target Hamiltonian H, measured on the optimized variational quantum state ρ(θ) imposes a bound on its fidelity F = ψ targ |ρ(θ)|ψ targ , with |ψ targ the exact ground state. Establishing this bound requires a good estimate of the two lowest energy levels E 0 and E 1 , which can be estimated on the quantum device via a subspace expansion technique for mixed states, illustrated in detail in the Supplementary Information of Ref. [83]. This analysis provides selfverification of the variational quantum state preparation algorithm [23]. Alternatively, for 1D lattice Hamiltonians, E 0 and E 1 can be obtained using numerical method DMRG.
Let us consider the following inequality stating that the state fidelity F is bound from below by the ratio In conclusion, as long as H ≤ E 1 , we have that approaching Tr[ρ 2 ] = 1 (certified pure state) in the limit H → E 0 .