Probabilistic eigensolver with a trapped-ion quantum processor

PHYSICAL REVIEW A 101, 052333 (2020) Probabilistic eigensolver with a trapped-ion quantum processor Jing-Ning Zhang ,1,2,*Iñigo Arrazola ,3Jorge Casanova ,3,4Lucas Lamata ,3,5Kihwan Kim ,2 and Enrique Solano 3,4,6,7,† 1Beijing Academy of Quantum Information Sciences, Beijing 100193, China 2Center for Quantum Information, Institute for Interdisciplinary Information Sciences, Tsinghua University, Beijing 100084, People’s Republic of China 3Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apdo. 644, 48080 Bilbao, Spain 4IKERBASQUE, Basque Foundation for Science, Maria Diaz de Haro 3, 48013 Bilbao, Spain 5Departamento de Física Atómica, Molecular y Nuclear, Universidad de Sevilla, 41080 Sevilla, Spain 6International Center of Quantum Artificial Intelligence for Science and Technology (QuArtist) and Department of Physics, Shanghai University, 200444 Shanghai, China 7IQM, Munich, Germany (Received 4 April 2020; accepted 7 April 2020; published 18 May 2020) Preparing the eigenstate, especially the ground state, of a complex Hamiltonian is of great importance in quantum simulations. Many proposals have been introduced and experimentally realized, among which are quantum variational eigensolver and heat-bath algorithmic cooling, with the former hindered by local minima and the latter lacking of complex system Hamiltonians. Here we introduce a dissipative quantum-classical hybrid scheme, the probabilistic eigensolver. The scheme repeatedly uses an ancilla qubit to acquire information on the system, based on which it postselectively lowers the average energy of the system. The optimal reduction is achieved through classical optimization with a single variational parameter. We describe the implementation of the probabilistic eigensolver with trapped-ion systems and demonstrate the performance by numerically simulating the ground-state preparation of several paradigmatic models, including the Rabi and the Hubbard models. We believe the scheme would enrich the functionalities of universal quantum simulators and be useful as a module for various quantum-computation tasks. DOI: 10.1103/PhysRevA.101.052333 I. INTRODUCTION Simulating complex quantum systems is known to be an intractable task if accomplished with classical computing resources. This is because the dimension of the Hilbert space needed to describe a quantum system increases exponentially with its number of constituents. Quantum simulation [1], i.e., the use of a fully controlled quantum system that simulates the dynamics of another one, was originally proposed to overcome this problem, and has led to plenty of developments in different quantum platforms [2]. On the other hand, quantum computing is a more general approach that allows us not only to simulate quantum dynamics, but also to solve systems of linear equations [3], linear differential equations [4–6], or the eigenvalue problem of complex Hamiltonians [7,8], with an exponential speedup. Regarding the latter, it is important to remark that the knowledge of the energy spectrum of a system is crucial in different fields such as quantum chemistry and condensed-matter physics. The phase estimation algorithm [7] is a prominent method for finding eigenvalues and preparing eigenstates of target Hamiltonians; however, other approaches have also been designed. Among them, we can mention the use of adiabatic evolution [9,10] and variational quantum eigensolvers [11–13]. *[email protected] †enr[email protected] The adiabatic method starts from a simple Hamiltonian and slowly evolves to a complex Hamiltonian of interest. During the evolution the adiabaticity condition has to be fulfilled. This can be challenging, especially when critical points corresponding to quantum phase transitions are present. Here, the energy gap between two or more eigenstates becomes exponentially small as the size of the system increases. On the other hand, variational eigensolvers are designed to prepare the ground state (or excited state [11,12]) of a complex Hamiltonian, which interestingly combine a quantum processor and a classical optimization algorithm. The variational eigensolvers have been implemented experimentally in state-of-the- art quantum platforms such as photonics [11], superconducting circuits [14], or trapped ions [15,16], and their accuracy in eigenstates preparation, as well as in the computing of the associated eigenvalues, highly depends on the flexibility of the ansatz [14,15,17]. In this respect, having a more complex and flexible input state, or ansatz, implies that more variational parameters have to be optimized simultaneously, which could be hinderedby the presence of local minima. In addition, the preparation of the ansatz may involve nontrivial entanglement operations among the available quantum registers. Thus, it is important to explore alternative algorithms to be applied in near-future quantum processors [18]. State preparation can also be achieved by dissipatively pumping entropy into the environment. For example, heatbath algorithmic cooling [19,20],whichiswidelyusedin nuclear magnetic resonance quantum computing, transfers 2469-9926/2020/101(5)/052333(11) 052333-1 ©2020 American Physical Society JING-NING ZHANG et al. PHYSICAL REVIEW A 101, 052333 (2020) entropy from the system to ancillary spins that rapidly thermalize through the environment. This method is usually used for preparing the spin system in a fully polarized product state. Dissipative methods to prepare entangled states also exist [21], with proof-of-principle experiments performed in trapped ions [22,23] and superconducting circuits [24]. Although powerful and robust, these state-oriented proposals assume complete prior information of the target state and thus cannot be used to find solutions encoded in some complex Hamiltonians. A closely related proposal, called dissipative quantum computation [25], designs a dissipative process whose unique steady state is the ground state of the target Hamiltonian; however, it only works for frustration-free systems. In this article, we propose a hybrid classical-quantum eigensolver that uses projective measurements on a single ancilla qubit to probabilistically prepare the eigenstates of general multipartite Hamiltonians. Unlike the quantum variational solver or adiabatic state preparation, the only requirement in our method is that the initial state has a non-negligible overlap with the desired target state. Previously, there were theoretical attempts on probabilistic preparation of the ground state of a single qubit [26]orann-level quantum system [27–29], with the former experimentally demonstrated with quantum optics [26]. Although the mechanism of our probabilistic eigensolver relies on arguments similar to those of Refs. [26–29], here we consider the implementation of this method in trapped-ion setups, which are capable of realizing various Hamiltonians consisting of multiple qubits and harmonic motional modes. In addition, our proposed eigensolver is featured by the introduction of a classical optimization module to reduce the coherent circuit depth. Moreover, we extend the method to prepare excited states with minor modifications. In Sec. II, we introduce the general method for preparing an arbitrary target eigenstate, given that the overlap between the target and the initial state is nonzero. Afterwards, we present how the method can be combined with a classical optimization algorithm, resulting in a significant reduction on the number of steps required for the preparation, while maintaining acceptable success probability. In Sec. III,we present the toolkit that trapped-ion quantum simulators offer to implement our method. Finally, in Sec. IV, we present some models in which our probabilistic eigenstate preparation could be applied by using a trapped-ion quantum processor and demonstrate its performance with numerical simulations. II. GENERAL FRAMEWORK Consider a target quantum system described by the Hamiltonian ˆ Hsuch that ˆ H|j=Ej|j.(1) Here, |j(Ej)isthe jth eigenstate (eigenenergy) and the energy spectrum {Ej}is bounded from below and sorted in ascending order. To prepare one of the eigenstates of ˆ H,we consider a quantum system isomorphic to the target system and an additional ancillary qubit. The evolution of the whole system is described by a unitary operator acting on the product Hilbert space HS⊗HA, with HS(HA) being the Hilbert space of the target system (the ancilla qubit). For preparing the jth eigenstate |j, our method requires knowledge of all eigenenergies below Ej. As, in general, this is not the case, we have to start with the ground-state preparation. In this manner, this first step can be considered as a cooling protocol. To this end, we apply the following unitary operator on the whole system: ˆ Wγ(τ)=exp −i(ˆ HS+γ)ˆσx Aτ.(2) Here, ˆσx Ais the x Pauli matrix acting on the ancillary qubit and the parameters γand τare real numbers to be determined. ˆ Wγ(τ) can be viewed as an evolution operator of the whole system with τbeing the effective evolution time. We use γto shift the energy spectrum of ˆ HSsuch that E0+γ⩾0 (i.e., the shifted spectrum is always positive), where E0is the groundstate energy of the target system. Our protocol is an iterative method where each repetition involves, first, initializing the ancillary qubit to a reference state |0Asuch that ˆσz A|0A=−|0A, and second, applying ˆ Wγ(τ). In the kth iteration we obtain the state |k, which is |k=ˆ Wγ(τ)|ψk−1|0A =ˆ Cγ(τ)|ψk−1|0A−iˆ Sγ(τ)|ψk−1|1A.(3) Here, |ψk−1is the normalized state vector in HSthat results from the k−1 iteration, and the nonunitary operators are defined as ˆ Cγ(τ)=cos[( ˆ H+γ)τ] and ˆ Sγ(τ)= sin[( ˆ H+γ)τ]. Finally, we perform a projective measurement on the ancillary qubit and make postselection on the measurement results to select certain the output state in HS. The two possible normalized output states, conditioned on the measurement results, are |ψk0= ˆ Cγ(τ)/(P0 k)1/2|ψk−1and |ψk1=ˆ Sγ(τ)/(P1 k)1/2|ψk−1, and their probabilities are P0 k=j|j|ψk−1|2cos2[(Ej+γ)τ] and P1 k=j|j|ψk−1|2sin2[(Ej+γ)τ], respectively. We would like to point out that, although we use a pure input state |ψk−1to illustrate the process, our method can be straightforwardly applied to mixed input states without any modification. Now, if we compare the average energy of the two possible output states with the average energy of the input state, we find that, for positive and small τ, the following relation holds (see Appendix A): ˆ H(0) k⩽ˆ H(0) k−1<ˆ H(1) k.(4) Here ·(i) k≡ψ(i) k|ˆ H|ψ(i) kand the superscript i=0,1 denotes the outcome of the projective measurement on the ancillary qubit. In other words, the average energy of the system is lowered compared with the one of the input state, i.e., the postselected state from the previous round, if we successfully project the ancillary qubit in |0A. Here, it is noteworthy to comment that the probability of projecting on the reference ancilla state |0A(i.e., on the state that leads to an effective average energy reduction) is close to one for a sufficiently small τ(note that P0 k=1ifτ=0). An intuitive way to understand this energy reduction for small τ is provided when studying the transformation of the state. If |ψk−1=jck−1 j|jis the system input state before the kth iteration, the output state after the application of ˆ Wγ(τ) and 052333-2 PROBABILISTIC EIGENSOLVER WITH A TRAPPED-ION … PHYSICAL REVIEW A 101, 052333 (2020) postselection of |0Awill be |ψk=|ψk0≡jck j|j, where the probability distribution of the state has changed as ck j 2=1 P0 k cos2[(Ej+γ)τ]ck−1 j 2.(5) As cos2(·) is a monotonically decreasing function in the vicinity of τ=0, the probability amplitudes of those eigenstates with lower eigenenergies Ejwill be enhanced [note that a larger value of Ejimplies a lower value for the corresponding cos2(·) function] compared with those with higher eigenenergies, leading to the decrease of the average energy—an effective cooling effect. After the kth iteration, we have to measure the average energy of the resulting state ¯ Ek≡ˆ H(0) k, and compare it with the average energy of the previous stage ¯ Ek−1.If|¯ Ek−1− ¯ Ek|< holds (with being a small positive value), we conclude that that the ground-state within a certain precision has been reached and stop the protocol. If the state does not converge at the kth iteration, saying |¯ Ek−1−¯ Ek|>,the output state |ψk0will be sent through another iteration to further reduce its energy. An important feature of our protocol is that, if the initial state |ψ0has no overlap with the desired eigenstate, it is not possible to prepare it. It can be seen from Eq. (5) that, if the (k−1)st input state has no overlap with a given eigenstate |j (i.e., ck−1 j=0), the resulting kth state will also share the same condition ck j=0. Although this may seem a disadvantage, it can actually be used to prepare excited eigenstates. To this end, we start from an arbitrary initial state |φ0and apply the unitary operation ˆ Us=exp −iπ 2Esˆ Hˆσx A(6) and measure the ancilla qubit. The action of the unitary operator ˆ Uson a general state |φis to conditionally produce a state |ψ(0)where the contribution of the basis state |sgoes to zero (s|ψ(0)=0) as long as the ancilla qubit is projected on |0A. As a result, the output state after postselection has no overlap with the eigenstate |s. For example, to prepare the first-excited state of ˆ Hstarting from an arbitrary initial state |φ0, one starts by applying the unitary ˆ U0and projecting the ancilla on the |0Astate. If the projection succeeds, one obtains a state |ψ0that has zero overlap with the ground state of ˆ H, which serves as input state for the cooling stage. Then, one applies the previously described protocol (see Fig. 1), until the state converges to the lowest energy eigenstate that has a nonzero overlap with the initial state, which in this case is the first-excited state. To prepare an arbitrary eigenstate |jstarting from |φ0, one first removes any overlap with lower-energy eigenstates by sequentially applying ˆ Ujand the following projective measurement with j=0,1,..., j−1. If all of the postselection operations succeed, the resulting state is sent to the cooling stage, which would probabilistically converge to the desired eigenstate |j. Note that once a postselection operation fails, the whole procedure should be started all over again. We consider now two different variants for the cooling stage, namely, the fixed-step and variational approaches. 1.0 0.0 01234 N-1N j 1.0 0.0 01234 N-1N |c0 j|2 j 1.0 0.0 01234 N-1N j ˆ U0 ˆ Wγ(τk) Restart |ψ0 |φ0 |w0 j|2 A A Restart 0A 1A 0A1A ×k |ψk |ck j|2 FIG. 1. Scheme for preparing the first-excited state of an N- dimensional Hamiltonian. First, we apply ˆ U0and measure the ancilla qubit. Then we use postselection to remove the ground-state contribution from the initial state |φ0. After that, we apply the cooling unitary ˆ Wγ(τk) and subsequent ancilla measurement ktimes, until we reach convergence. The fixed-step approach treats τas a small fixed parameter. The operational procedure is as follows: After each implementation of ˆ Wγ(τ), we perform a projective measurement on the ancilla qubit. We continue with the next iteration if the outcome is |0A, otherwise we restart the whole process. At the kth stage of the process, we estimate the average energy of the output state ¯ Ek. The process stops if |¯ Ek−1−¯ Ek|⩽, where the precision is a small positive value. Otherwise, we continue the protocol by moving to the (k+1)st stage. In the variational approach, the value of the parameter τ is optimized in each stage such that it minimizes the average energy. This hybrid optimization is done by feeding the average energy obtained from the quantum platform to a classical optimization algorithm, which provides new trial value for τ to the quantum platform based on the given information (see Fig. 2). More specifically, the classical algorithm uses Brent’s method to find a minimum in the interval 0 <τ k<1. After several optimization steps (trial steps), the classical algorithm will converge to τk, the optimal value of τat the kth stage, and the quantum platform will estimate the average energy ¯ Ek. As before, the process stops if |¯ Ek−1−¯ Ek|⩽.Atthe end, we obtain a set of optimized values {τ1,τ 2,...,τ k}. Compared to the fixed-step approach, the variational approach allows the probabilistic preparation of the desired eigenstate with a substantial reduction on the circuit depth, i.e., with a qualitatively fewer gates and postselections. Both the fixed-step and the variational approaches are probabilistic because we need to restart the whole procedure depending on the results of the projective measurements on the ancilla qubit and provide an estimation of the eigenenergy of the target eigenstate within the preset precision . However, the variational protocol additionally provides a recipe for shallow-depth probabilistic preparation of the target state, as described by the set of optimized parameters {τ1,...,τ k}.The overall success probability Psuc of both protocols is ultimately 052333-3 JING-NING ZHANG et al. PHYSICAL REVIEW A 101, 052333 (2020) FIG. 2. Scheme for obtaining the optimal value of τfor the kth stage in the variational approach to the probabilistic eigensolver. First, prepare the state at the (k−1)st stage by applying the effective cooling k−1 times with optimized parameters {τ1,τ 2,...,τ k−1}. Then, search for the optimal value τkwith the help of a classical optimization algorithm. Specifically, the quantum simulator estimates the average energy and feeds it into the classical algorithm, which determines whether the optimal value is reached or alternatively the next trial value for τ. bounded by the overlap between the initial state and the desired eigenstate, i.e., Psuc ⩽|ψ0|j|2. On the other hand, to maximize Psuc for the preparation of eigenstate |j, one should choose the value of γto be the closest possible to −Ej. However, it is possible that we do not know the exact value of Ej, in which case one can always choose γto be minus the lower bound of the norm of the Hamiltonian, which can be analytically approximated in general. Finally, we comment that the unitary operator ˆ Wγ(τ) required in our method can be decomposed as ˆ Wγ(τ)≡ ˆ W(τ)RA x(2γτ), with ˆ W(τ)=exp(−ˆ Hˆσx Aτ) and RA α(θ)≡ exp(−iθ 2ˆσα A), α=x,y,z, being the single qubit rotation acting on the ancilla qubit. Now, the ˆ W(τ) operator that entangles the system and the ancilla qubit, can be trotterized into a set of basic operations [30]. More specifically, suppose that the system Hamiltonian is composed of Mcomponents ˆ Hmthat do not commute with each other, namely ˆ H=M m=1ˆ Hm. Then, the unitary evolution in Eq. (2) can be approximately written by using the second-order Trotter-Suzuki expansion, ˆ W(τ)=ˆ W1τ 2r··· ˆ WM−1τ 2rˆ WMτ r ׈ WM−1τ 2r··· ˆ W1τ 2rr +Oτ3 r2,(7) where ris the number of Trotter steps and ˆ Wm(τ)= exp(−iˆ Hmˆσx Aτ) with m=1,2,...,M. Note that, by increasing the number of Trotter steps, arbitrary precision can be achieved in principle. III. TOOLKIT IN TRAPPED-ION PLATFORMS Trapped ions have been demonstrated to be suitable for digital [31] and analog [32,33] quantum simulations, having access not only to the manipulation of the ion-qubits, but also to the coherent control of their collective vibrational modes [34]. The digital-analog approaches which combine both resources could enhance the computational power of the quantum platform for solving hard problems [35], for example, nontrivial many-body models or even problems in quantum field theories [36]. We envision four types of primitive system-ancilla couplings in the toolkit of trapped-ion systems, including (1) ˆσz iˆσx A,(2)ˆa† mˆamˆσx A,(3)ˆσx i1ˆσx i2ˆσx A, and (4) (ˆa† m+ˆam)ˆσx iˆσx A. Using these interactions and the measurement on the ancilla, the protocol described in the previous section could be implemented in trapped-ion systems for a wide variety of models. The free Hamiltonian ˆ H0of an array of Nions trapped in a linear Paul trap is written as follows: ˆ H0= N  n=1 ¯hω0 2ˆσz n+ N  m=1 ¯hωmˆa† mˆam+1 2,(8) where ˆσx,y,z nare the Pauli matrices acting on the nth ion with ω0being the frequency splitting of the two involved internal levels, and ˆam(ˆa† m) is the annihilation (creation) operator of the mth collective motional mode with the mode frequency ωm. For a better readability, from now on we omit the subscript n(m)ofˆσx,y,z n(ˆam) when there is no ambiguity. In the following, we provide a brief explanation of how to engineer the four types of system-ancilla couplings in a trapped-ion quantum simulator: First, spin-spin couplings can be implemented by the Mølmer-Sørensen (MS) type of interaction, which takes the form of σx iσx jand can be engineered between any pair of ions embedded in a linear chain of ions [37]. With the appropriate single qubit rotations on the first qubit, we can extend the coupling to the form ˆσnσx A, which would correspond to a simulated local spin-1/2 operator ˆ HS∝ˆσn, where n= (nx,ny,nz)Tis a unit vector. Second, the ˆa†ˆaˆσx Ainteraction term involving the mth mode could be realized with a detuned red-sideband coupling, for which the Hamiltonian is written, in the rotating frame with respect ˆ H0and after the rotating wave approximation, as ˆ HI=¯hA 2[(ˆσ+eit+ˆσ−e−it)] +i¯hηA 2(ˆaˆσ+e−iδt−ˆa†ˆσ−eiδt),(9) where Ais the ancilla Rabi frequency, η≡ηA m= δk√¯h/(2Mωm) is the Lamb-Dicke (LD) parameter associated with the mth mode and the ancilla qubit, with δkbeing the 052333-4 PROBABILISTIC EIGENSOLVER WITH A TRAPPED-ION … PHYSICAL REVIEW A 101, 052333 (2020) (a) (b) FIG. 3. Quantum circuits for the trapped-ion toolkit. (a) Scheme for implementing the unitary exp (−iφ 2ˆσx i1ˆσx i2ˆσx A)by using single-qubit rotations and global Mølmer-Sørensen gates MS =exp(−iπ 4i<jˆσx iˆσx j), where the summation runs over all possible pairs of qubits. Here, the single-qubit rotations are defined as Rα(φ)=exp(−iφ 2ˆσα) with α=x,y,orzand we use Rα(π 2)≡Rα for simplicity. (b) Scheme for implementing the unitary operator exp [−iφ 2(ˆa+ˆa†)ˆσx iˆσy A] using MS gates, single-qubit gates, and the analog block ˆ UR(φ) involving the mth motional mode. transferred wave vector of the laser and Mbeing the mass of the ion, and and δare the detunings to the carrier and the first red-sideband transitions, respectively. In the dispersive regime, i.e., Aand ηAδ, the effective Hamiltonian can be written as follows [38]: ˆ Heff =¯h2 A 4ˆσz+¯hη22 A 4δˆa†ˆa+1 2ˆσz,(10) where the AC-Stark shift term can be absorbed into the free Hamiltonian, and the term that depends on the bosonic operators can be transformed into ˆa†ˆaˆσx Aby appropriate single-qubit rotations on the ancilla. The third term ˆσz 1ˆσz 2ˆσz 3can be implemented by combining the controlled-not (CNOT) gates [8]ortheMSgates[39] with single-qubit zrotations. Moreover, global MS gates [40] can be used to reduce the number of CNOTs or two-qubit MS gates, increasing the process fidelity. Figure 3(a) shows a recipe to construct this nonlocal spin operation with global MS gates. Finally, similarly to the previous quantum circuit, we show that the dynamics governed by (ˆa†+ˆa)ˆσx iˆσx Acan be implemented by combining CNOT gates (or MS gates [41]) with analog blocks that involve the dipolar coupling between the mth motional mode with the ancilla qubit. In Fig. 3(b) we show the complete circuit to achieve the desired interaction. The central part of the circuit is a unitary operator produced by a simulated dipolar coupling that can be implemented by simultaneously driving the red-sideband and blue-sideband transitions of the ancilla qubit [42]. The associated evolution operator ˆ UR(φ) in the rotating frame defined by the free Hamiltonian ˆ H0is written as follows: ˆ UR(φ)=exp −iφ 2(ˆa+ˆa†)ˆσx A,(11) where φis proportional to the Rabi frequency A,theLamb- Dicke parameter ηA m, and the interaction time. Another qubit will be involved by a pair of CNOT gates and all the other single-qubit rotations are used for changing the basis of the Pauli matrices. With the described toolkit available on trapped-ion quantum platforms, the probabilistic eigensolver can be extended to a wide range of models involving both qubits and bosonic modes. The latter include spin-spin interaction models like the Ising, XY, or Heisenberg models [43], quantum Rabi models [44–47], the Dicke models [48], or second-quantization Hamiltonians like the Holstein model [49]. Also, more general methods such as those for computing n-time correlation functions [50] or to simulate dissipative processes [51] could benefit from the presented eigensolver as the latter can be used to prepare arbitrary eigenstates. In the next section, we numerically investigate the performance of the probabilistic eigensolver with several examples, including a single harmonic oscillator, the quantum Rabi model, and the Hubbard model. IV. EXAMPLES A. Harmonic oscillator One of the simplest and most important models of quantum mechanics is the quantum harmonic oscillator, which gives a mathematical description of various physical phenomena, including mechanical oscillators in harmonic potentials and electromagnetic fields [52]. The Hamiltonian of a single harmonic oscillator is written as follows: ˆ H=¯hωh.o.ˆa†ˆa,(12) where ˆa(ˆa†) is the annihilation (creation) operator of the mode with the mode angular frequency ωh.o.. The Hamiltonian in Eq. (12) provides an equal-spacing spectrum with an infinite number of levels. One could think on a simple example of the protocol presented in Sec. II, where the application of a unitary of the form ˆ W(τ)= exp [−iτˆa†ˆaˆσx A] (see Sec. III) in a single ion, combined with measurements on the ancilla, could lead to a probabilistic preparation of the motional ground state. Note that the ground-state energy for a harmonic oscillator is positive, thus γis set to zero. Starting from, for example, a thermal state, each time the ancilla is measured to be in the state |0A after applying ˆ W(τ), the motional state is effectively “cooled down” (i.e., the average energy is lowered). In Fig. 4, we show the results for the probabilistic groundstate preparation of a harmonic oscillator with ωh.o.=1. The initial state is naturally chosen as a thermal equilibrium state characterized by the thermal average number ¯nth =0.5. Both the fixed-step cooling with τ=0.3 [Fig. 4(a)] and the variational cooling [Fig. 4(b)] show a monotonically decreasing average energy ˆ Hwith respect to the number of iterations k. The overall success probability Psuc(k)uptothekth stage is defined as Psuc(k)= k  k=1 P0 k,(13) where P0 kis the probability of projecting the ancilla to |0A at the kth iteration. It is shown in Fig. 4that both in the fixed-step and variational cases, Psuc(k) decreases towards a saturate value around 0.65. In this example, the fixed-step 052333-5 JING-NING ZHANG et al. PHYSICAL REVIEW A 101, 052333 (2020) FIG. 4. Probabilistic ground-state preparation of a single harmonic mode. (a) Average energy ˆ H(black dots) and the overall success probability Psuc(k) (blue dots) as functions of the number of iterations kfor the fixed-step protocol. The normalized step size τis fixed to be 0.3. (b) Average energy ˆ H(red dots over the black curve) and the overall success probability Psuc(k) (red dots over the blue curve) as functions of the number of iterations kfor the optimized variational protocol. Black dots represent results for all the different trial steps that are needed until reaching convergence for the optimum τk, which are {0.8487,0.5044,0.9919,0.9919,0.9910,0.7430,0.4194,0.9881}. and the variational protocols reach convergence in 18 and 8 steps, respectively, where a precision of =10−3is pursued. Figure 4(b) shows that the total number of trial steps (black dots) needed to obtain the optimized variational protocol are around 80 (≈10 trials per step). However, once you have the set of optimized values {τk}, the variational approach requires fewer steps than the fixed-step approach, which reduces the circuit depth and alleviates the requirement for long-term quantum coherence. B. Quantum Rabi model The quantum Rabi model [53] describes the dipolar interaction between a magnetic dipole and an oscillating magnetic field. Recently, its quantum simulation has been realized experimentally with highly tunable parameters in many physical platforms, including superconducting circuits [54,55] and trapped ions [56]. The Hamiltonian of the quantum Rabi model can be separated into two noncommuting parts, ˆ H=ˆ H1+ˆ H2, with ˆ H1=¯hω0 2ˆσz+¯hωˆa†ˆa, (14) ˆ H2=¯hg(ˆa+ˆa†)ˆσx, where ˆσxand ˆσzare the Pauli matrices acting on the two-level system with the energy splitting ω0,ˆa(ˆa†) is the annihilation (creation) operator of the harmonic mode with frequency ω, and gis the coupling strength of the dipolar interaction. We implement ˆ W(τ) through the second-order Trotter- Suzuki expansion, ˆ Wγ(τ)=ˆ W1δτ 2ˆ W2(δτ)ˆ W1δτ 2r ,(15) where δτ =τ/ris the size of the Trotter step and ˆ Wj(τ)= exp(−iˆ Hjˆσx Aτ) with j=1 and 2. ˆ W1could be implemented by combining the first and second interactions described in Sec. III, while ˆ W2requires the implementation of the fourth interaction. Overall, the simulation requires two qubits and a motional mode. Note that, depending on the step size τ,more than one Trotter step may be needed. Figure 5shows the results for the ground-state preparation of the quantum Rabi model in the deep strong-coupling regime [57]. We choose the state |↓|0as the initial state, for it can be easily prepared by optical pumping and sideband cooling. Comparing both cases we clearly see how the variational approach offers a faster route to the ground-state preparation, with an overall success probability around 60%, very similar to the fixed-step approach. C. Hubbard model The Hubbard model [58] is one of the most important models in solid-state physics. In spite of the simplicity in form, it describes the quantum phase transition between superconducting and Mott-insulator phases. The analytical solution has not been available for arbitrary dimensions and the exact numerical treatment is believed to be hard for classical computers because of the sign problem. To illustrate the performance of our protocol in a system with only qubits, we consider the one-dimensional Hubbard model with Lsites and open boundary conditions, whose Hamiltonian is written as follows: ˆ H=−t L−1  i=1 σ (ˆc† i,σ ˆci+1,σ +H.c.) +U L  i=1 ˆc† i,↑ˆci,↑ˆc† i,↓ˆci,↓, (16) where tand Uare the nearest-neighbor hopping strength and the on-site interaction strength, respectively. This model can be mapped to a nearest-neighbor spin chain of 2Nspins through the Jordan-Wigner transformation [59], which establishes a mapping between a fermionic operator and a set of spin operators. We choose the one-dimensional Hubbard model as an example, while the extension to two- and threedimensional (2D and 3D) models can be done in trapped ions by using nonlocal interactions [36]. The operator in Eq. (2) can be constructed by using the three-body spin operators described in Sec. III and the Suzuki-Trotter expansion. For 052333-6 PROBABILISTIC EIGENSOLVER WITH A TRAPPED-ION … PHYSICAL REVIEW A 101, 052333 (2020) FIG. 5. Probabilistic ground-state preparation of the quantum Rabi model with ω0=1.2, ω=0.8, and g=1.0. The energy is in units of gand the number of stages is denoted by k. The number of Trotter steps is chosen as r=3. (a) Fixed-step probabilistic state preparation. The step size is chosen as τfix =0.3. (b) Variational probabilistic state preparation. The total number of trial steps is denoted by n(black dots). The trial steps for different stages are discriminated by different background colors. Red dots represent the energy and overall probability for different stages using the optimized τvalues {0.4762,0.9839,0.9032,0.5575}. the simulations, we choose t=1 and U=2 if not explicitly specified otherwise. Figure 6shows the numerical simulations for the groundstate preparation of the two-site and three-site Hubbard models using the variational approach for probabilistic state preparation. As initial states, we choose |↑↓↑↓and |↓↑↓↓↑↓for the two-site and three-site models, respectively. In addition, the ˆ Wγ(τ) operation is implemented by using the symmetric Trotter expansion with r=3 Trotter steps. The results show that, with less than 10 steps, one can obtain the ground state for both models. The overall success probability is in this case around 15% which is still an acceptable value. V. CONCLUSIONS We have presented a probabilistic eigensolver which is capable of preparing arbitrary eigenstates of Hamiltonians that are implementable with a trapped-ion quantum processor. The method, applicable to a digital or digital-analog quantum simulator, requires control operations and measurements on an extra ancilla qubit. Moreover, we provide a recipe to enhance the performance of the probabilistic eigenstate prepa- FIG. 6. Probabilistic ground-state preparation of (a) two-site and (b) three-site Hubbard models. The black dots indicate the average energy of the trial states after successfully projecting the ancilla qubit to the reference state |0A, while the blue dots represent the overall success probability. Red dots represent the average energy and overall success probability after the kth stage, which is done by applying the ˆ Wγ(τ) unitary with the optimized τk,usingr=3Trotter steps. The optimized taus are {0.3692,0.3959,0.3684,0.3959}and {0.2694,0.3661,0.3247,0.3973,0.2726,0.3959}for the two-site and three-site models, respectively. ration by means of a hybrid classical-quantum optimization algorithm. We describe a basic toolbox natural in trapped-ion quantum platforms, which can be used as building blocks to implement the method for complex Hamiltonian models. Finally, we numerically simulate the method for some interesting examples that could be implemented in state-of-the-art trapped-ion setups. ACKNOWLEDGMENTS This work was supported by the National Key Research and Development Program of China under Grants No. 2016YFA0301900 and No. 2016YFA0301901 and the National Natural Science Foundation of China Grants No. 11374178, No. 11574002, No. 11504197, and No. 11974200. We also acknowledge financial support from Spanish Government via PGC2018-095113-B-I00(MCIU/AEI/FEDER, UE), Basque Government via IT986-16, as well as from QMiCS (820505) and OpenSuperQ(820363) of the EU 052333-7 JING-NING ZHANG et al. PHYSICAL REVIEW A 101, 052333 (2020) Flagship on Quantum Technologies, and the EU FET Open Grant Quromorphic. I.A. acknowledges support from Basque Government Ph.D. Grant No. PRE-2015-1-0394. J.C. acknowledges support from the UPV/EHU grant EHUrOPE. APPENDIX A: PROBABILISTIC COOLING EFFECT For an arbitrary input state |ψ=jcj|j, the state for the system plus the ancilla qubit after applying ˆ Wγ(τ)is denoted as |≡ ˆ Wγ(τ)|ψ|0A. The values of the average energy with respect to the input state and conditioned output states are ˆ H=ψ|ˆ HS|ψ= j Ej|cj|2,(A1) ˆ H(0) =|ˆ HS⊗|0A0A|| =1 P0 j Ej|cj|2cos2[(Ej+γ)τ], ˆ H(1) =|ˆ HS⊗|1A1A|| =1 P1 j Ej|cj|2sin2[(Ej+γ)τ], where P0=j|cj|2cos2[(Ej+γ)τ] and P1=j|cj|2 sin2[(Ej+γ)τ]. For small τ, i.e., (Ej+γ)τ1, we expand ˆ H(0) to the second order in τand calculate the difference between ˆ H(0) and ˆ Has follows: ˆ HS−ˆ HS(0) = j (|cj|2−|cj|2)Ej,(A2) where the new probability distribution is related with the old one as cj cj=1 P0 cos [(Ej+γ)τ] ≈1 P01−1 2(Ej+γ)2τ2+O(τ4).(A3) Because E0+γ⩾0, the factor 1 −1 2(Ej+γ)2τ2decreases for higher excited states, meaning that the higher the energy of the state, the lower the probability amplitude at the output state. Similarly, we can calculate the difference between ˆ HS(1) and ˆ HSto the leading order and obtain the following expression: ˆ HS(1)−ˆ HS= j (|c j|2−|cj|2)Ej,(A4) where now the new probability distribution is related with the old one as c j cj=1 P1 sin [(Ej+γ)τ] ≈1 P1 (Ej+γ)τ+O(τ3).(A5) Notice that, in this case, the lower the energy of the state, the lower the ratio between the new and the old probability amplitude. APPENDIX B: SCALING ANALYSIS In this Appendix, we study the scaling behavior of the probabilistic eigensolver in the preparation of the ground state of a generic Hamiltonian. Specifically, we focus on the dependence of the number of cooling steps on the infidelity of the prepared ground state. To simplify the scenario, we analytically obtain the scaling behavior of the number of cooling steps in a cooling process with constant cooling strength, which serves as an upper bound for the number of cooling steps in a classically assisted cooling process with optimized cooling strength for each cooling step. A generic input state is |ψin=jcj|j. If there are k steps and we get kzeros, we get the following unnormalized state: |ψout= j cjcoskφj|j,(B1) with φj=Ejτ+γ. Let us define the condition that the state is convergent to the ground state as M j=1Pjcos2kφj P0cos2kφ0 ⩽, (B2) with Pj≡|cj|2. In the parameter regime that 0 ⩽φ0<φ 1<···<φ M⩽ π 2, thus cos φ1⩾cos φjfor j>2 and it is evident that the following inequality is valid: M j=1Pjcos2kφj P0cos2kφ0 ⩽(1−P0) P0cos φ1 cos φ02k ,(B3) where we make use of the fact that, in the parameter regime, 1 ⩾cos φ1>cos φ2>···>cos φM⩾0. Comparing the right-hand sides of the above two equations, we obtain the scaling for the number of steps kas follows: k⩾1 2ln P0 1−P0ln cos φ1 cos φ0 .(B4) The denominator can be further processed as ln cos φ1 cos φ0=ln 1−tan φ0τ −1 22τ2+O(3τ3) =− tan φ0τ +1 22τ2+O(3τ3),(B5) where =E1−E0is the energy gap between the ground and the first-excited state. Finally, we obtain k⩾−1 2ln P0 1−P0tan φ0τ +1 22τ2−1 .(B6) APPENDIX C: NOISE ANALYSIS All realistic quantum computations are subjected to noise, leading to decoherence. Here we study how dephasing or depolarizing errors can harness the presented cooling protocol. First, we consider the case in which each qubit is subjected to a dephasing channel only while performing multiqubit gates. This simplified scenario is reasonable because dephasing affects more to multiqubit operations, as these are slower than single-qubit gates. Let us consider the quantum circuit in Fig. 3(a) as our ideal process. The quantum map 052333-8 PROBABILISTIC EIGENSOLVER WITH A TRAPPED-ION … PHYSICAL REVIEW A 101, 052333 (2020) (a) (b) FIG. 7. Noisy analysis. (a) The average fidelity between the ideal and noisy processes for the quantum circuit in Fig. 3(a).The horizontal axis is the single-qubit dephasing rate per multiqubit gate. (b) The excited-state probability P1,kafter ksuccessful steps of cooling, when the probabilistic eigensolver is subjected to a depolarizing environment characterized by the error rate δ=0.1. The input state is chosen as the completely mixed state, and the cooling strength is characterized by φ1=π/8, π/4, and π/2. corresponding to the ideal process is Eid =UR† z◦UMS†◦UR† y(φ)◦UMS ◦URz,(C1) where ◦means the concatenation of quantum maps, and the unitary channels are UO(ρ)=OρO†with Obeing single- or multiqubit operations. With spatially localized dephasing noise, the quantum map becomes E(p)=UR† z◦E⊗3 deph ◦UMS†◦UR† y(φ) ◦E⊗3 deph ◦UMS ◦URz,(C2) where the single-qubit dephasing channel is characterized by Edeph(ρ)=(1 −p)ρ+pˆσzρˆσz, with pbeing the single-qubit dephasing rate. To show the effect of the dephasing channel, we calculate the average fidelity ¯ F(Eid,E(p))between the ideal and the noisy process. As shown in Fig. 7(a), the average fidelity decreases as the single-qubit dephasing rate increases such that it is over 98% when the dephasing rate p=0.01. We also consider the effect of a depolarizing error acting on both the system and the ancillary qubits. For simplicity, we consider ground state preparation in a two-level system. The depolarizing error acts with probability δper unit of time. In this case, γand τare carefully chosen to make φ0=0 and φ1⩽π 2. We consider the completely mixed state ˆρin = 1 2(|00|+|11|) as the initial state, although the protocol should work for any initial state with a nonzero overlap with the desired state. After ksteps, the unnormalized output state (a) (b) FIG. 8. Cooling of a two-level system under dephasing. The system Hamiltonian is the one in Eq. (C5)withbx=bz=1. (a) Average energy after a one cooling step versus the effective evolution time τ. The initial state is |0S0Aand the average is taken over the postselected final state. The applied cooling unitary Wγ(τ)is depicted in the inset. Solid markers are results obtained in the ideal situation while hollow ones are obtained with dephasing noise. The error rate per MS gate is p=0.01. (b) Average energy (solid red) and ground-state fidelity (dashed blue) versus the dephasing strength p. becomes Ek de(ˆρin )=˜ P0,k|00|+˜ P1,k|11|,(C3) where the unnormalized probabilities are ˜ P0,k=(1−δ)P0,k−1+δ 4, (C4) ˜ P1,k=(1−δ)P1,k−1cos2φ1+δ 4, with initial conditions P0,0=P1,0=1 2. The probabilities for the ground and the excited states can then be derived as Pi,k= ˜ Pi,k/(˜ P0,k+˜ P1,k) with i=0or1. Figure 7(b) shows that the effective cooling operation changes the stationary state of the system. Instead of retrieving the completely mixed state, which is the stationary state without effective cooling, we end up with a new stationary state, where the excited-state population has been suppressed due to the effective cooling protocol. Finally, we investigate the effect of dephasing on the cooling of a two-level system. The system Hamiltonian of a spin-1/2 system under a static magnetic field is ˆ Hs=−bz 2ˆσz−bx 2ˆσx,(C5) where bxand bzare magnetic field components along the x and zaxes. We first initialize the system and ancilla qubits to the state |0S0A. Then, we implement ˆ Wγ(τ)forsome fixed τand γ. Next, we measure the ancilla qubit. The 052333-9