Quantum Simulation of Pairing Hamiltonians with Nearest-Neighbor Interacting Qubits
aa r X i v : . [ qu a n t - ph ] N ov Quantum Simulation of Pairing Hamiltonians with Nearest-Neighbor Interacting Qubits
Zhixin Wang, Xiu Gu,
1, 2
Lian-Ao Wu, and Yu-xi Liu
1, 2, 4, ∗ Department of Microelectronics and Nanoelectronics, Tsinghua University, Beijing 100084, China Institute of Microelectronics, Tsinghua University, Beijing 100084, China Ikerbasque, Basque Foundation for Science, 48011 Bilbao and Department of Theoretical Physics and History of Science,The Basque Country University (UPV/EHU), PO Box 644, 48080 Bilbao, Spain Tsinghua National Laboratory for Information Science and Technology (TNList), Beijing 100084, China (Dated: September 14, 2018)Although a universal quantum computer is still far from reach, the tremendous advances in controllablequantum devices, in particular with solid-state systems, make it possible to physically implement “quantumsimulators”. Quantum simulators are physical setups able to simulate other quantum systems efficiently that areintractable on classical computers. Based on solid-state qubit systems with various types of nearest-neighborinteractions, we propose a complete set of algorithms for simulating pairing Hamiltonians. Fidelity of thetarget states corresponding to each algorithm is numerically studied. We also compare algorithms designedfor different types of experimentally available Hamiltonians and analyze their complexity. Furthermore, wedesign a measurement scheme to extract energy spectra from the simulators. Our simulation algorithms mightbe feasible with state-of-the-art technology in solid-state quantum devices.
PACS numbers: 03.67.Lx, 03.67.Mn, 74.20.Fg
I. INTRODUCTION
Classical computers fail to efficiently simulate quantumsystems with complex many-body interactions due to the ex-ponential growth of variables for characterizing these sys-tems [1, 2]. For instance, N parameters are required forthe complete description of a quantum system composed of N entangled spin-1/2 particles. In the 1980s, quantum sim-ulation was proposed to solve such an exponential explosionproblem using a controllable quantum system [3]. In 1996,it was shown that a quantum computer only containing few-particle interactions can be used to efficiently simulate many-body quantum Hamiltonians [4]. Recently quantum simu-lation has attracted extensive attention in condensed matterphysics [5, 6], high energy physics [7] and quantum chem-istry [8] due to the rapid progress on coherent control of quan-tum systems aimed at quantum information processing [9–13]. Quantum simulators using trapped ions [14–17], coldatoms [18, 19] and photons [20] have already been experi-mentally demonstrated to some extent.Quantum simulators are classified into analog and digi-tal ones [1] . An analog quantum simulator is a control-lable quantum system mimicking the behaviors of the targetquantum system whose evolution can be effectively mappedonto the simulator [6], while a digital quantum simulator nor-mally imitates the time evolution operator of the target systemthrough the implementation of a series of elementary quan-tum gates [16]. Practically, these two approaches are oftencombined. A simulation task can be performed through thefree evolution of quantum simulators combined with externallogic gates at given time instants. Quantum simulators con-taining tens of qubits are practically useful to carry out classi-cally intractable tasks while the number of qubits required for ∗ Electronic address: [email protected] practical quantum computing is much larger [1]. Therefore,in comparison with universal quantum computing, quantumsimulation is more feasible in the near future, and the simu-lation algorithms for certain task using existed Hamiltoniansare strongly desired.Pairing Hamiltonians, for example, BCS Hamiltonian inconventional superconductors, feature long-range many-bodyinteractions which are generally intractable on classical com-puters. Nevertheless, large-scale numerical calculations basedon pairing Hamiltonians are of great importance, for instancein mesoscopic condensed matter, ultrasmall metallic grainsand heavy nuclei [21]. To tackle this problem, a polynomial-time quantum algorithm based on a nuclear magnetic reso-nance (NMR) quantum computer was proposed [21], and hasbeen demonstrated experimentally [22, 23]. However, liquidNMR has several constrains that make NMR quantum com-puter not scalable [24]. Therefore, the large-scale implemen-tation of the NMR-based quantum algorithm is unlikely withthe state-of-the-art technology.Rapid progress of superconducting quantum circuits hasbeen witnessed in the past decades [24–28]. This enablesthem to be one of the most promising candidates towardspractical quantum information processing. Based on nearest-neighbor coupled superconducting qubits, single-qubit andtwo-qubit gates with fidelity at the threshold of fault tolerantquantum computing has been realized recently [29]. Varioustheoretical and experimental explorations on quantum simu-lation have been carried out using this approach [30–32]. Theunique flexibility on design and fabrication of superconduct-ing circuits enables wide tunability in extensive Hamiltonianparameter ranges and the techniques for scaling up are com-patible with those for modern integrated circuits. Both ofthe above aspects are significant advantages for superconduct-ing quantum circuits to serve as practical quantum simulators.Moreover, the research on other solid-state qubit devices, e.g.quantum dot in semiconductors [33] and defect systems [34],has also made significant progress in past years. Therefore, itwould be of great interest to update the algorithm in NMR sys-tem [21] to simulate the paring Hamiltonian using solid-statequbit systems.Here we propose a complete set of simulation algorithmsand a measurement scheme to simulate many-body pairingmodels using various Hamiltonians existed in solid-state qubitsystems. The algorithms are suitable for a wide range of quan-tum systems, especially superconducting quantum circuitsand semiconducting qubit systems. Sec. II is a brief descrip-tion of our simulation task and available theoretical Hamilto-nians for quantum simulation. General simulation algorithmsbased on qubits with various types of nearest-neighbor inter-actions are presented in Sec. III. The fidelity of the algorithmsand its variation with various parameters are numerically stud-ied in Sec. IV. The complexity of the algorithms as well as itsparameter dependence is analyzed in Sec. V. We have alsodesigned an effective measurement scheme based on the en-tanglement of a single qubit in the simulator and an ancillaryqubit in Sec. VI, through which crucial information in the en-ergy spectrum could be extracted. We finally summarize ourresults in Sec. VII.
II. PRELIMINARIES
For the completeness of this paper, first we briefly reviewpairing Hamiltonians and outline our simulation tasks.
A. Pairing Hamiltonians and qubit representation
The general BCS pairing Hamiltonian is extensively usedin condensed matter physics and nuclear physics, and has theform: H BCS = N X m =1 ǫ m n Fm + n F − m ) + N X l =1 N X m =1 V ml c † m c †− m c − l c l , (1)where c †± m and c ± l are fermionic creation and annihilationoperators, and n F ± m = c †± m c ± m are fermionic number op-erators. As has been analyzed (see for example Ref. [21]),the BCS pairing Hamiltonian made by fermionic pair oper-ators can be mapped onto qubit operators σ xm , σ ym and σ zm through the transformation { σ xm , σ ym , σ zm } = { c † m c †− m + c − m c m , ic − m c m − ic † m c †− m , n Fm + n F − m − } . With the map-ping, we can rewrite the Hamiltonian in Eq. (1) as, H p = N X m =1 ε m σ zm + X m Solid-state qubits can be coupled through various types ofinteractions. For instance, superconducting qubits can be cou- TABLE I: Summary for various interaction Hamiltonians in solid-state systems with single-qubit Hamiltonians H = P Nl =1 ( ω l σ zl / .Note that J xl = J yl = J l for the Heisenberg model in the table.Interaction types Interaction HamiltoniansLongitudinal Ising model H Ising,L = H + P N − l =1 J l σ zl σ zl +1 Transverse Ising model H Ising,T = H + P N − l =1 J l σ xl σ xl +1 XY model H XY = H + P i = x,y P N − l =1 J l σ il σ il +1 Heisenberg model H H = H + P i = x,y,z P N − l =1 J il σ il σ il +1 pled to their nearest neighbors through capacitances [35], in-ductances [36] or Josephson junctions [37]. Different inter-action models [38] resulted from different coupling schemescan be classified into four categories of commonly used inter-action Hamiltonians: longitudinal Ising types [39–41], trans-verse Ising types [42–44], XY types [45–47] and Heisenbergtypes [48–52]. These four different types of nearest-neighborHamiltonians can be unified as H = H + H I (3)with H denoting the single-qubit Hamiltonian H = N X l =1 ω l σ zl (4)and H I denoting the interaction Hamiltonian H I = N − X l =1 ( J xl σ xl σ xl +1 + J yl σ yl σ yl +1 + J zl σ zl σ zl +1 ) . (5)Here σ xl , σ yl , σ zl are Pauli matrices in the basis of σ zl , and l denotes the l th qubit. Parameters J l ( l = 1 , . . . , N − denotethe coupling strength between the l th and ( l + 1) th qubits.Each of these four types of Hamiltonians is a special caseof Eq. (3) with parameters being properly chosen. The Hamil-tonian (3) can be thus reduced to (i) longitudinal Ising Hamil-tonian for parameters J xl = J yl = 0 and J zl = J l ; (ii) thetransverse Ising Hamiltonian for parameters J yl = J zl = 0 and J xl = J l ; (iii) the XY Hamiltonian for parameters J zl = 0 and J xl = J yl = J l ; (iv) the Heisenberg Hamiltonian for pa-rameters J xl = J xl = J yl = J l . For clarity, we summarizethem in Table I.Assuming that qubit systems with nearest-neighbor cou-plings in Eq. (3) are available experimentally, we can henceuse them to simulate dynamical behaviors of pairing Hamil-tonians in Eq. (2) with the help of single-qubit operations.These operations will be done by applying external fields F = P Nl =1 ( f xl σ xl + f yl σ yl + f zl σ zl ) to individual qubits. It isclear that the pairing models in Eq. (2) do not share the sameform of Hamiltonians with in Eq. (3). Therefore, it is neces-sary to design algorithms to simulate these pairing Hamilto-nians using the four types of interaction Hamiltonians men-tioned above. C. Quantum simulations Our goal is to simulate the BCS-type pairing Hamiltonianin Eq. (2) which can be rewritten as H p = H p + H pI , (6)where H p = P Nm =1 ε m σ zm / is the single-qubit Hamilto-nian and H pI = P Nl>m V ml ( σ xm σ xl + σ ym σ yl ) / is the interac-tion Hamiltonian.We now come to simulate parameters ε m ( m = 1 . . . N ) and V ml ( m, l = 1 . . . N, m < l ) in Eq. (6) with each ofthe Hamiltonians in Table I by using the decoupling and re-coupling techniques in Ref. [53]. The single-qubit terms ε m σ zm / m = 1 , . . . , N ) and two-qubit interaction Hamil-tonian terms V ml ( σ xm σ xl + σ ym σ yl ) / l, m = 1 , . . . , N, l In this section, we will give detailed discussion on how tosimulate pairing Hamiltonian using Ising-type, or XY-type, orHeisenberg-type interaction Hamiltonians as listed in Table I. A. Algorithm for Simulators with Longitudinal IsingHamiltonian We first study an algorithm to simulation pairing Hamilto-nians using longitudinal Ising Hamiltonian. The simulationalgorithm needs two steps [21]. The first step is to simulate H p and the second one is to mimic H pI . We then combine H p and H pI to obtain the complete pairing Hamiltonian. Thedetailed description for these steps is given as below. 1. Simulation Algorithm for Individual Terms H p We now use the simulators with longitudinal Ising-type in-teraction H Ising , L in Table I to simulate H p . The operator U Ising,L ( τ ) = e − iτH Ising,L (7) denotes the time evolution operator of the quantum simula-tor. Let us consider the time evolution operator U zl ( τ ) =exp( − iτ ω l σ zl / of the l th qubit. We show that U zl ( τ ) canbe given as U zl ( τ ) = O j ′′ = l e i π σ xj ′′ U Ising,L (cid:16) τ (cid:17) O j ′ e i π σ xj ′ U Ising,L (cid:16) τ (cid:17) O j ′′ = l e − i π σ xj ′′ U Ising,L (cid:16) τ (cid:17) O j ′ e − i π σ xj ′ U Ising,L (cid:16) τ (cid:17) (8)using U Ising , L ( τ / . Here j ′ and j ′′ denote the j ′ th and the j ′′ th qubits in the simulator. j ′ is an even (odd) number and j ′′ is an odd (even) number if l is an odd (even) number, where j ′′ = l . Figure 1 shows the corresponding quantum circuits. .. …… .. … .. … … .. … … … FIG. 1: The quantum circuit to realize U zl ( τ ) from U Ising,L . The redbold line stands for the l th qubit in the system, whose individualevolution term is going to be extracted through this step. Black dotson both sides stand for the periodic extension of the logic gates to therest of the system in the given direction. Here, the period is two. X ± stand for external gates e ± i π σ x . For a given U zl ( τ ) , we can extract the time evolution H Ising,L of each single qubit and then simulate H p through thecircuit e − itH p = N Nm =1 U zm ( τ ) with τ = ε m t/ω m . There-fore, individual U zl ( τ ) is the building block when simulating H p . In the following subsections we will show that U zl ( τ ) can be simulated on simulators with other types of interac-tions. 2. Simulation Algorithm for Interactions H pI and H p Before simulating H pI of the pairing Hamiltonian in Eq. (6)we first consider a time evolution operator U zzl,m ( τ ) =exp( − iτ J l σ zl σ zm ) . If the simulator has longitudinal Isingnearest neighbor couplings, U zzl,l +1 ( τ ) can be directly obtainedfrom U Ising,L ( τ ) through the circuit as in Fig. 2, which can beexpressed as U zzl,l +1 ( τ ) = O j ′ 1. Simulation Algorithm for H p Based on Trotter’s formula, we can decompose the timeevolution operator with Heisenberg type interaction Hamilto-nian into U H ( τ ) = e − iτH H ≈ exp N X l =1 ω zl σ z ! exp − iτ N X l =1 H zzl ! × exp − iτ N − X l =1 H xyl ! + O ( J τ ) , (13)with H zzl = J l σ zl σ zl +1 .When the quantum simulators possess Heisenberg type in-teraction, we can first simulate longitudinal Ising Hamiltonianusing Heisenberg Hamiltonian and then H p in terms of lon-gitudinal Ising Hamiltonian. The latter step has already beensolved. Assume that error up to O ( J τ ) is tolerable, wecould obtain the longitudinal Ising Hamiltonian through U Ising,L ( τ ) ≈ O j e i π σ z j U H (cid:16) τ (cid:17) O j e − i π σ z j U H (cid:16) τ (cid:17) , (14)Here j denotes that the qubits with even indices are rotated ± π/ around the z direction, the qubits with odd indices re-main unchanged. We can also rotate all qubits with odd in-dices ± π/ around z axis to obtain the same result. Figure 3shows the corresponding quantum circuit. 2. Simulation Algorithm for Interaction H pI We now design the algorithm to simulate the long range XYinteraction terms in the pairing Hamiltonian in Eq. (6) usinga quantum simulator with Heisenberg interaction. Let us firstconsider two time evolution operators U zz ( τ ) = exp( − iτ N − X l =1 J l σ zl σ zl +1 ) , (15) … .. … .. …… FIG. 3: Quantum circuit to simulate U Ising,L ( τ ) from U H ( τ ) . Z ± stand for external gates e ± i π σ z . Black dots on both sides standfor the periodic extension of the logic gates starting from any sin-gle qubit. Here, the period is two. and U zzl,m ( τ ) = exp( − iτ J l σ zl σ zm ) . (16)As shown in Fig. 4(a), U zz ( τ ) can be simulated through U H ( τ ) with error up to order O ( J τ ) . Here U H ( τ ) is thetime evolution operator of simulators with Heisenberg inter-actions, U zz ( τ ) ≈ O j e i π σ x j O j e i π σ y j +1 U H (cid:16) τ (cid:17)O j e − i π σ x j O j e − i π σ y j +1 U H (cid:16) τ (cid:17) . (17)Figure 4(b) shows that U zzl,l +1 ( τ ) can be simulated through U zz ( τ ) , that is, U zzl,l +1 ( τ ) = O j ′ 1. Simulation Algorithm for H p If simulators with the XY interaction are available and thetransition frequencies of all qubits in the simulators are iden-tical or almost identical, the time evolution operator is separa-ble, U XY ( τ ) = e − iτH XY = exp( − iτ H ) exp " − iτ N − X l =1 H xyl . (22)Here H and H xyl are given in Eq. (4) and Eq. (11) respec-tively. Figure 6(a) shows when the Heisenberg Hamiltonian inFig. 3 is replaced by the XY type Hamiltonian, we can obtain U z ( τ ) = exp( − iτ P Nl =1 ω l σ zl / , which can be expressed as U z ( τ ) = O j e i π σ z j U XY (cid:16) τ (cid:17) O j e − i π σ z j U XY (cid:16) τ (cid:17) . (23)The quantum circuit for U zl ( τ ) is shown in Fig. 6(b) and canbe expressed as U zl ( τ ) = O j = l e i π σ xj U z (cid:16) τ (cid:17) O j = l e − i π σ x j U z (cid:16) τ (cid:17) . (24)The single-qubit Hamiltonian H p can be simulated in thesame way from U zl ( τ ) as those for longitudinal Ising andHeisenberg interaction. … … …… (a) … …… … (b) FIG. 6: The quantum circuit for achieving (a) U z ( τ ) from U XY ( τ ) and (b) U zl ( τ ) from U z ( τ ) . The red bold line stands for the l th qubitin the system. The period of extension represented by black dots istwo in (a), while is one in (b). 2. Simulation Algorithm for Interaction H pI If we use simulators with XY interactions, then U xx + yy ( τ ) can be simulated in the following way U xx + yy ( τ ) = O j e i π σ xj U XY (cid:16) τ (cid:17) O j e − i π σ xj U XY (cid:16) τ (cid:17) , (25)which is shown in Fig. 7(a). … … …… (a) … … … … (b) FIG. 7: The quantum circuits for simulating (a) U xx + yy ( τ ) from U XY ( τ ) and (b) U xxl,l +1 ( τ ) directly from U xx + yy ( τ ) . The red boldlines stand for the l th and the ( l + 1) th qubits in the system. Blackdashed lines on both sides stand for the periodic extension of thelogic gates. In (a) the period is one while in (b) the period is two. Figure 7(b) shows that we can acquire U xxl,l +1 ( τ ) from U xx + yy ( τ ) in the following way U xxl,l +1 ( τ ) = O j ′′ e i π σ xj ′′ O j ′ 1. Simulation Algorithm for Individual H p In contrast to the above cases, U Ising,T ( τ ) = e − iτH Ising,T forthe Hamiltonian with transverse Ising interactions is not ex-actly separable. Even so, as shown in Fig. 8(a), U z ( τ ) can beextracted by the same quantum circuit as it was in Fig. 6(a)and Fig. 6(b), if error up to order O ( ωJτ ) is tolerable. Here ω is the typical value of qubit frequency. However, we willshow in Sec. IV that this approximation due to Trotter’s for-mula leads to the fluctuation in fidelity. 2. Simulation Algorithm for Interaction H pI Let us first consider a time evolution operator U z + xx ( τ ) = exp " − iτ N X l =1 ω l σ zl + N − X l =1 J l σ xl σ xl +1 ! , (29) … … …… (a) … … …… (b) FIG. 8: The quantum circuit for simulating (a) U z ( τ ) from U Ising,T and (b) U z + yy ( τ ; ω , . . . , ω N ) from U z + xx ( τ ; ω , . . . , ω N ) . Blackdots on both sides stand for the extension of the logic gates withperiod one. √ Z ± stand for external gates e ± i π σ z . which is exactly the same as U Ising,T ( τ ) with tunable parame-ters ω l ( l = 1 , . . . , N ) . Similarly, we also consider U z + yy ( τ ) = exp " − iτ N X l =1 ω l σ zl + N − X l =1 J l σ yl σ yl +1 ! . (30)The operation U z + yy ( τ ) can be acquired from U z + xx ( τ ) by U z + yy ( τ ) = N O j =1 e i π σ zj U z + xx ( τ ) N O j =1 e − i π σ zj , (31)which is graphically shown in Fig. 8(b). If all qubit resonancefrequencies in U z + xx ( τ ) are set to be ω and those in U z + yy ( τ ) to be − ω , we have U xx + yy ( τ ) ≈ U z + xx ( τ ) U z + yy ( τ ) (32)up to order O ( J τ ) . Here, the negative qubit resonancefrequencies mean that inverse external fields are used to flipthe sign of P Nl =1 ω l σ zl / . Thus, the simulation using thetransverse Ising Hamiltonian is converted to the the cases ofHeisenberg or XY Hamiltonian as studied before. IV. NUMERICAL STUDY OF FIDELITY We now come to study numerically the error introduced byTrotter’s formula, which are typically in the order O ( J t ) or O ( ωJt ) . Our focus will be the fidelity of the simulatedHamiltonian using the simulation algorithms.In Sec. III, we assume for generality that all of the al-gorithms are applicable to simulators with constant cou-pling strengths and frequencies. However, if simulators withtunable parameters are available, we can not only simplifythe simulation process but also reduce digital errors signifi-cantly. For instance, e − iτH Ising,T = exp( − iτ P Nl =1 ω l σ zl / P N − l =1 J l ( σ xl σ xl +1 + σ yl σ yl +1 + σ zl σ zl +1 ) does not hold ex- | ω t | / 2 π F i d e lit y (a) | ω t | / 2 π (b) FIG. 9: Digital fidelity of (a) e − iτH Ising,T ≈ exp( − iτ P Nl =1 ω l σ zl / P N − l =1 J l ( σ xl σ xl +1 + σ yl σ yl +1 + σ zl σ zl +1 ) and (b) process shownin Fig. 8(a) for a 4-qubit quantum simulator. ε l = 2 × Hz, ω l =5 × Hz ( l = 1 , . . . , N ) , V l = − × Hz and J l = 3 × Hz ( l = 1 , . . . , N − . The fluctuations on the curves originate fromthe trigonometric terms in Eq. (33) with period π/ω in (a) and π/ω in (b) respectively. actly in general. Detailed calculation indicates exp − iτ N X l =1 ω σ zl ! exp − iτ N − X l =1 Jσ xl σ xl +1 ! = exp ( − iτ N X l =1 ω σ zl + N − X l =1 Jσ xl σ xl +1 ! + iτ N − X l =1 J (cid:0) σ xl σ xl +1 − σ yl σ yl +1 (cid:1) − i ω N − X l =1 J (cid:0) σ xl σ xl +1 − σ yl σ yl +1 (cid:1) sin(2 ωτ )+ i ω N − X l =1 J ( σ yl σ xl +1 + σ xl σ yl +1 ) [cos(2 ωτ ) − ) (33)where iτ P N − l =1 J ( σ xl σ xl +1 − σ yl σ yl +1 ) / is the error that willaccumulate when simulation processes. These terms contain-ing sin(2 ωτ ) and cos(2 ωτ ) are the origins of the fluctuationsof period π/ω on the fidelity curve, as shown in Fig. 9(a).Numerical simulation shows that the fidelity of the circuit inFig. 8(a) is plotted in Fig. 9(b). The period of fluctuation be-comes π/ω because U Ising,T ( τ / instead of U Ising,T ( τ ) is in-volved. However, if the interaction in H Ising,T can be turnedoff when H p is simulated, those fluctuations can be avoidedas shown in Fig. 10(a).The fidelity of a simulation algorithm increases when tun-able parameters of the simulator increase, as exemplified inFig. 10. The effect of tunable parameters is significant exceptfor simulator with longitudinal Ising Hamiltonian, in whichall terms commute with each other.The use of Trotter’s formula, or the short time approxima-tion improves the simulation fidelity significantly, as illus-trated in Fig. 11. For simulators with longitudinal or trans-verse Ising interactions, tuneable XY or Heisenberg inter-actions, we find that the more time steps the total simula-tion time is divided into, the higher simulation fidelity is ob-tained, as in Figs. 11 (a), (b), (c), and (e). We believe that |Vt| F i d e lit y (a) |Vt| (b) |Vt| F i d e lit y (c) |Vt| (d) FIG. 10: Increase of simulation fidelity contributed by tunable pa-rameters of 4-qubit simulators containing (a)(b) transverse Isingtype, (c) XY type and (d) Heisenberg type nearest-neighbor inter-action. In all four subfigures, ε l = 2 × Hz, ω l = 5 × Hz ( l = 1 , . . . , N ) , V l = − × Hz. In (a)(b), J l = 3 × Hz ( l =1 , . . . , N − and the whole simulation task is divided into M = 20 intervals, while in (c)(d) J l = 3 × Hz ( l = 1 , . . . , N − andthe whole simulation task is divided into M = 10 intervals. Fur-thermore, the subprocess in (c) and (d) for simulators with constantparameters to simulate e ± iτ ( σ xl σ xl +1 + σ yl σ yl +1 ) π/ ( l = 1 , . . . , N − is divided into 200 intervals in order to reduce the error introducedby short-time approximation. In (a), red solid curve stands for theaverage fidelity when all the parameters of the simulator are con-stant. Magenta solid curve stand for the average fidelity when only ω l ( l = 1 , . . . , N ) are tunable. Green dot curves stand for the highfrequency fluctuation of simulation fidelity originated from Eq. (33)associated with the above two cases. Blue solid curve gives the fi-delity when all the parameters of the simulator can be turned on andoff during the simulation process. The detailed shape of the fidelityfrom a simulator with constant parameters over a short interval isshown in (b). In both (c) and (d), red solid curve stands for the av-erage fidelity when all the parameters of the simulator are constantwhile blue solid curve gives the fidelity when all the parameters ofthe simulator can be turned on and off. the error introduced by short time approximation can be re-duced in this way. Nevertheless, for simulators with fixed-parameter XY or Heisenberg interaction Hamiltonians, we no-tice that the overall fidelity is sensitive to that of the sub-circuitused to simulate U xx + yyl,l +1 [ π/ (4 J l )] = e ± iτ ( σ xl σ xl +1 + σ yl σ yl +1 ) π/ ( l = 1 , . . . , N − (see Figs. 11 (d) and (f)). On the otherhand, on the contrary to the error due to short-time approxi-mation, the error from simulating U xx + yyl,l +1 [ π/ (4 J l )] increaseswith the running times of the sub-circuit. Figures 11(d) and(f) show that when the simulation time is relatively short andthe error from the Trotter’s formula is small, the overall sim-ulation fidelity is mainly determined by that of simulating U xx + yyl,l +1 [ π/ (4 J l )] . Therefore, we conclude that in this casethe more time steps the total time for completing the wholetask is divided into, the more time U xx + yyl,l +1 [ π/ (4 J l )] simula- |Vt| F i d e lit y M = 3 M = 5 M = 10 M = 20 (a) |Vt| M = 3 M = 5 M = 10 M = 20 (b) |Vt| F i d e lit y M = 1 M = 2 M = 3 M = 5 M = 10 (c) |Vt| M = 1 M = 2 M = 3 M = 5 M = 10 (d) |Vt| F i d e lit y M = 1 M = 2 M = 3 M = 5 M = 10 (e) |Vt| M = 1 M = 2 M = 3 M = 5 M = 10 (f) FIG. 11: Effect of changing the number of intervals M that the wholesimulation process is divided into on the fidelity for 4-qubit simula-tors containing (a) longitudinal Ising type, (b) transverse Ising type,(c)(d) XY type and (e)(f) Heisenberg type nearest-neighbor interac-tion. In (c) and (e), all parameters of the simulators can be turned offand on during the process, while in (d) and (f) they are constant. Inall six subfigures, ε l = 2 × Hz, ω l = 5 × Hz ( l = 1 , . . . , N ) , V l = − × Hz, while in (a)(b) J l = 3 × Hz and in(c)(d)(e)(f) J l = 1 × Hz. The subprocess in (d) and (f) for sim-ulators with constant parameters to simulate e ± i ( σ xl σ xl +1 + σ yl σ yl +1 ) π/ ( l = 1 , . . . , N − is divided into G = 200 intervals in order toreduce the error introduced by short-time approximation. tion process is run, the larger the error is. When the simu-lation time becomes longer, error from theTrotter’s formulabecomes dominant. Therefore, as shown in Figs. 11(d) and(f), all curves are almost flat at the short simulation time anddrop when the time increases. Moreover, the fidelity of simu-lation with smaller number M of intervals is higher at shortersimulation time but lower at longer simulation time.Figure 12 suggests that if other conditions are the same,the fidelity of simulation algorithms decreases when the qubitnumber in quantum simulators increases. In practice, this er-ror can be compensated by reducing the error due to the short-time approximation. |Vt| F i d e lit y N = 4 N = 5 N = 6 N = 7 N = 8 FIG. 12: Fidelity of the algorithm run by quantum simulator withtransverse Ising type nearest-neighbor interactions containing differ-ent number of qubits. Every term in the Hamiltonian of the simu-lator can be turned on and off during the simulation process. ε l =2 × Hz, ω l = 5 × Hz ( l = 1 , . . . , N ) , V l = − × Hz, J l = 3 × Hz ( l = 1 , . . . , N − . The whole simulation task isdivided into M = 20 intervals. V. COMPLEXITY ANALYSIS Let us now analyze complexity of the algorithms, which isdetermined by the total number of external single-qubit logicgate required for the simulation process. We will show thatthe algorithms are polynomial.The number of single-qubit gates in any of quantum cir-cuits in the previous sections increases linearly with the num-ber of qubits N in the simulator. One can extract U zl ( τ ) =exp( − iτ ω l σ zl / in the complexity O ( N ) . The simulation ofindividual e − itH p = N Nm =1 U zm ( ε m t/ω m ) requires O ( N ) external gates. We find that U zzl,l +1 ( τ ) = exp( − iτ J l σ zl σ zl +1 ) or U xx + yyl,l +1 ( τ ) = exp[ − iJ l ( σ xl σ xl +1 + σ yl σ yl +1 )] can also besimulated within O ( N ) . As shown in [21], it can be shown e i π X l,l +1 U zzm,l ( τ ) e − i π X l,l +1 = U zzm,l +1 ( τ ) , (34)for m < l , where X l,m = ( σ xl σ xm + σ yl σ ym ) / . Thus itneeds O ( N ) to simulate an arbitrary long-range interaction U zzm,l ( τ ) . Since N ( N + 1) / terms are included in the pair-ing Hamiltonian Eq. (6), simulation of all these interactionsrequires O ( N ) logic gates, or the complexity of the wholealgorithm is O ( N ) .In Fig. 13 we calculate total numbers of external gates forfour types of simulators with all constant parameters when N ≤ . It shows that longitudinal Ising simulator has thelowest complexity. On the other hand if the parameters of thesimulators are tunable, the complexity is significantly reducedwhen the scale of the simulator increases, as shown in Fig. 14.Furthermore, although higher simulation fidelity can be ob-tained through the short-time approximation by dividing the0 N N u m b e r o f G a t e s Longitudinal IsingTransverse IsingXYHeisenberg FIG. 13: Total number of external gates in the simulation pro-cess for four types of simulators with constant parameters. Here thecomplexity is shown for a single running process, in which both thesimulation task and subprocess to simulate e ± i ( σ xl σ xl +1 + σ yl σ yl +1 ) π/ ( l = 1 , . . . , N − are not further divided into any subintervals( M = 1 , G = 1 ). N N u m b e r o f G a t e s Constant XYTunable XYConstant HeisenbergTunable Heisenberg FIG. 14: Effects of tunable parameters on the complexity for XYtype and Heisenberg type simulators. Solid lines stand for the casewhen all parameters in the simulator are constant, while dashed linesgive the result when they are all tunable. Here the total number of ex-ternal gates is shown for a single running process, in which both thesimulation task and subprocess to simulate e ± i ( σ xl σ xl +1 + σ yl σ yl +1 ) π/ ( l = 1 , . . . , N − are not further divided into any subintervals( M = 1 , G = 1 ). simulation process into M > intervals, this is at the expenseof higher complexity, since the number of logic gates is M times of the gate number without the short-time approxima-tion. Moreover, although Trotter’s formula can be applied tothe subprocess of simulating exp[ ± iπ ( σ xl σ xl +1 + σ yl σ yl +1 ) / l = 1 , . . . , N − , the complexity grows linearly with the G N u m b e r o f G a t e s Constant XYTunable XYConstant HeisenbergTunable Heisenberg FIG. 15: Increasing of complexity due to more subintervals theprocesses for simulating e ± i ( σ xl σ xl +1 + σ yl σ yl +1 ) π/ ( l = 1 , . . . , N − are divided into. XY and Heisenberg simulators are taken as theexamples. Solid lines stand for the case when all parameters in thesimulator are constant, while dashed lines give the result when theyare all tunable. Here the whole simulation task is not further dividedinto any subintervals ( M = 1 ). number G of intervals that this subprocess is divided into,which is shown in Fig. 15. We can also find that if simula-tors with tunable parameters are available, the effect of G onthe total number of external gates will be reduced. In particu-lar, when all parameters of an XY type simulator are tunable,the growing of G has no effect on the simulation complexity. VI. MEASUREMENT SCHEME Now we come to the measurement or readout approach asthe last step of our simulation algorithm. The measurementcircuit is shown in Fig. 16 directly after the simulation forparing Hamiltonian. We use an ancillary qubit, denoted by ared color and entangled with the simulator. We measure theancilla at the end of the simulation process. We use one qubit,marked by a red color , to directly interact with the ancillaryqubit in the whole simulator. |±i = ( | ↑ i ± | ↓ i ) / √ areeigenstates of operator σ x for the ancillary qubit.The readout processes as follows. First, the ancillary qubitis prepared to the state | + i , while the simulator, includingqubits to N involved in the algorithms, is prepared to a state | ψ i . The whole system is initially at | Ψ i = | + i ⊗ | ψ i . Wethen apply a controlled-NOT gate denoted by CNOT → to qubit and qubit so that the system state becomes | Ψ i = CNOT → | Ψ i , where the ancilla serves as thecontrol qubit and the qubit serves as the target qubit. Werun the complete simulation algorithms for a time interval t ,represented by U p ( t ) = e − itH p such that the system becomes | Ψ i = [ I ⊗ U p ( t )] | Ψ i . A new controlled-NOT gate is ap-plied to the qubit and the qubit again, where the ancilla1 … … … … or 10 (ancilla) or FIG. 16: Measurement scheme for the simulation algorithms in or-der to extract information such as energy spectrum from the simula-tor. | + i and | ψ i stand for the initial state of the ancillary qubit andthe quantum simulator respectively. Index 0 and 1 denote the acillaand the single qubit in the simulator directly coupled with the ancillarespectively. | Ψ i , | Ψ i , | Ψ i and | Ψ i are the states of the wholesystem during the simulation process. U p ( t ) is the time evolutionoperator including the complete simulation algorithms, which is ex-plicitly equivalent to the pairing Hamiltonian. Two CNOT gates inwhich qubit 0 (ancilla) serves as the control qubit and qubit 1 servesas the target qubit are graphically shown. Measurement is done to theancilla on | + i or |−i basis, giving probability P +0 ( t ) or P − ( t ) asthe result respectively, where t corresponds to the evolution time forBCS system in U p ( t ) . again serves as the control qubit and the qubit serves asthe target qubit. The state of the whole system ends up with | Ψ i = CNOT → | Ψ i . Finally, the ancilla is measured on {| + i , |−i } basis. The probabilities of obtaining the states |±i are P ± ( t ) , respectively, which vary with time t . We canuse either P +0 ( t ) and P − ( t ) to extract spectrum informationof the simulator. For example, P +0 ( t ) is calculated as P +0 ( t ) = 12 + 14 [ h ψ | σ x ( t ) σ x (0) | ψ i + c.c. ] , (35)where σ x ( t ) = U † p ( t ) σ x U p ( t ) is in the Heisenberg picture.We could now expand the initial state with a complete set ofquantum numbers n , i and β i , | ψ i = P n,i,β i B n,i,β i | n, i, β i i .Here n is the total number of spin-up qubits in our simulator. i denotes the energy level E n,i with a given n and β i denotesthe degeneracy for a given energy level E n,i . Physically, n is also the number of Cooper pairs in the simulated pairingHamiltonian. States with different quantum number n are mu-tually orthogonal, h m, j, β j | n, i, β i i = δ nm h m, j, β j | n, i, β i i .In order to simplify the calculation, we can specially prepare qubit on the spin-up state, such that σ z | ψ i = | ψ i and h ψ | σ x ( t ) σ x (0) | ψ i = X n,i,j ˜ C n,i,j e it ( E n,j − E n − ,i ) (36)where ˜ C n,i,j = P β i ,β j B ∗ n,j,β j B n,i,β i h n, j, β j | n, i, β i i . Tostudy the measurement result in frequency domain, we cantake the Fourier transform ˜ ρ +0 ( ω ) of P +0 ( t ) , ˜ ρ +0 ( ω ) = πδ ( ω )+ π X n,i,j ˜ C n,i,j δ ( ω + E n,j − E n − ,i )+ π X n,i,j ˜ C ∗ n,i,j δ ( ω − E n,j + E n − ,i ) . (37)The spectrum of P +0 ( t ) is symmetric on amplitude and an-tisymmetric on phase because P +0 ( t ) is real. In amplitudespectrum, there should be sharp peaks, ideally δ -shape peaksat frequencies ω ± n,i,j = ± ( E n,j − E n − ,i ) . Energy spectraof the pairing Hamiltonian could be extracted from these fre-quencies. For instance, the energy gap between the groundstate and the first excitation state with n -Cooper pairs can beobtained by n = E n, − E n, = ( E n, − E n − ,i ) − ( E n, − E n − ,i )= ω + n,i, − ω + n,i, . (38)Generally, low-lying energy spectra of a pairing Hamilto-nian are of great interest. Using the adiabatic method devel-oped in Ref. [21], we can prepare initial states that involvea few eigenstates of the simulated pairing Hamiltonian to im-prove the efficiency of our measurement approach. VII. CONCLUSION In summary, we study algorithms for simulating the pairingHamiltonians based on various nearest-neighbor interactions,e.g., longitudinal Ising Hamiltonian, transverse Ising Hamil-tonian, XY Hamiltonian and the Heisenberg Hamiltonian,which are available in the solid-state quantum devices. Withcurrent experimental advances on gate fidelity and coherencetime, our proposal might be feasible in superconducting qubitcircuits. Since the four types of nearest-neighbor interactionsare shared by various quantum systems [48, 49, 52], our algo-rithms are adaptable to various solid-state systems. [1] I. Buluta and F. Nori, Science , 108 (2009).[2] I. M. Georgescu, S. Ashhab, and Franco Nori, Rev. Mod.Phys. , 153 (2014).[3] R. Feynman, Int. J. Theor. Phys. , 467 (1982).[4] S. Lloyd, Science , 1073 (1996).[5] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk and J. J. Bollinger, Nature ,489 (2012).[6] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss and M.Greiner, Nature , 307 (2011).[7] R. Gerritsma, G. Kirchmair, F. Z¨ahringer, E. Solano, R. Blattand C. F. Roos, Nature , 68 (2010). [8] B. P. Lanyon, J. D. Whitfield, G. G. Gillett, M. E. Goggin,M. P. Almeida, I. Kassal, J. D. Biamonte, M. Mohseni, B. J.Powell, M. Barbieri, A. Aspuru-Guzik and A. G. White, NatureChem. , 106 (2010).[9] C. Monroe and J. Kim, Science , 1164 (2013).[10] D. D. Awschalom, L. C. Bassett, A. S. Dzurak, E. L. Hu, Jason.Petta, Science , 1174 (2013).[11] R. Blatt and D. Wineland, Nature , 1008 (2010).[12] I. Bloch, Nature , 1016 (2010).[13] R. Hanson and D. D. Awschalom, Nature , 1043 (2010).[14] R. Blatt and C. F. Roos, Nature Phys. , 277 (2012).[15] J. T. Barreiro, M. M¨uller, P. Schindler, D. Nigg, T. Monz, M.Chwalla, M. Hennrich, C. F. Roos, P. Zoller and R. Blatt, Na-ture , 486 (2011).[16] B. P. Lanyon, C. Hempel, D. Nigg, M. M¨uller, R. Gerritsma, F.Z¨ahringer, P. Schindler, J. T. Barreiro, M. Rambach, G. Kirch-mair, M. Hennrich, P. Zoller, R. Blatt, C. F. Roos, Science ,57 (2011).[17] K. Kim, M.-S. Chang, S. Korenblit, R. Islam, E. E. Edwards,J. K. Freericks, G.-D. Lin, L.-M. Duan and C. Monroe, Na-ture , 590 (2010).[18] I. Bloch, J. Dalibard and S. Nascimb`ene, Nature Phys. , 267(2012).[19] J. Struck, C. ¨Olschl¨ager, R. Le Targat, P. Soltan-Panahi, A.Eckardt, M. Lewenstein, P. Windpassinger, K. Sengstock, Sci-ence , 996 (2011).[20] A. Aspuru-Guzik and P. Walther, Nature Phys. , 285 (2012).[21] L.-A. Wu, M. S. Byrd and D. A. Lidar, Phys. Rev. Lett. ,057904 (2002).[22] K. R. Brown, R. J. Clark, and I. L. Chuang, Phys. Rev. Lett. ,050504 (2006)[23] X. Yang, A. M. Wang, F. Xu and J. Du, Chem. Phys. Lett. ,20 (2006).[24] T. D. Ladd, F. Jelezko, R. Laflamm, Y. Nakamura, C. Monroeand J. L. OBrien, Nature , 45 (2010).[25] M. H. Devoret and R. J. Schoelkopf, Science , 1169 (2013).[26] J. Clarke and F. K. Wilhelm, Nature , 1031 (2008).[27] J. Q. You and F. Nori, Phys. Today (11), 42 (2005).[28] J. Q. You and F. Nori, Nature , 589 (2011).[29] R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey,T. C.White, J. Mutus, A. G. Fowler, B. Campbell, Y. Chen,Z. Chen, B. Chiaro, A. Dunsworth, C. Neill, P. O’Malley, P.Roushan, A. Vainsencher, J. Wenner, A. N. Korotkov, A. N.Cleland and J. M. Martinis, Nature , 500 (2014).[30] A. A. Houck, H. E. T¨ureci and J. Koch, Nature Phys. , 292(2012).[31] D. Underwood, W. E. Shanks, J. Koch and A. A. Houck, Phys.Rev. A , 023837 (2012). [32] P. Roushan, C. Neill, Yu Chen, M. Kolodrubetz, C. Quintana, N.Leung, M. Fang, R. Barends, B. Campbell, Z. Chen, B. Chiaro,A. Dunsworth, E. Jeffrey, J. Kelly, A. Megrant, J. Mutus, P. J.J. O’Malley, D. Sank, A. Vainsencher, J. Wenner, T. White, A.Polkovnikov, A. N. Cleland and J. M. Martinis, Nature , 241(2014).[33] M. A. Eriksson, S. N. Coppersmith, and M. G. Lagally, MRSBulletin , 794 (2013).[34] L. Gordon, J. R. Weber, J. B. Varley, A. Janotti, D. D.Awschalom, and C. G. Van de Walle, MRS Bulletin , 802(2013).[35] T. Yamamoto, Yu. Pashkin, O. Astafiev, Y. Nakamura and J. S.Tsai, Nature , 941 (2003).[36] J. B. Majer, F. G. Paauw, A. C. J. ter Haar, C. J. P. M. Harmans,and J. E. Mooij, Phys. Rev. Lett. , 090501 (2005).[37] J. Siewert, R. Fazio, G.M. Palma, and E. Sciacca, J. Low Temp.Phys. , 795 (2000).[38] N. Schuch and J. Siewert, Phys. Rev. A , 032301 (2003).[39] T. P. Orlando, J. E. Mooij, L. Tian, C. H. van der Wal, L. S. Lev-itov, S. Lloyd, and J.J. Mazo, Phys. Rev. B , 15398 (1999).[40] Y. Makhlin, G. Sch¨on, and A. Shnirman, Nature , 305(1999).[41] J. Q. You, J. S. Tsai, and F. Nori, Phys. Rev. Lett. , 197902(2002).[42] Yu. A. Pashkin, T. Yamamoto, O. Astafiev, Y. Nakamura, D. V.Averin and J. S. Tsai, Nature , 823 (2003).[43] R. McDermott, R. W. Simmonds, Matthias Steffen, K. B.Cooper, K. Cicak, K. D. Osborn, Seongshik Oh, D. P. Pappas,John M. Martinis, Science , 1299 (2005).[44] P. Bertet, C. J. P. M. Harmans, and J. E. Mooij, Phys. Rev. B ,064512 (2006).[45] A. Imamo ¯ glu, D. D. Awschalom, G. Burkard, D. P. DiVincenzo,D. Loss, M. Sherwin, and A. Small, Phys. Rev. Lett. , 4204(1999).[46] D. Mozyrsky, V. Privman, and M. Glasser, Phys. Rev. Lett. ,5112 (2001).[47] L. Quiroga and N. F. Johnson, Phys. Rev. Lett. , 2270 (1999).[48] D. Loss and D. DiVincenzo, Phys. Rev. A , 120 (1998).[49] B. Kane, Nature , 133 (1998).[50] G. Burkard, D. Loss, and D. P. DiVincenzo, Phys. Rev. B ,2070 (1999).[51] X. Hu and S. Das Sarma, Phys. Rev. A , 062301 (2000).[52] R. Vrijen, E. Yablonovitch, K. Wang, H. W. Jiang, A. Balandin,V. Roychowdhury, T. Mor and D. DiVincenzo, Phys. Rev. A ,012306 (2000).[53] D. W. Leung, I. L. Chuang, F. Yamaguchi, and Y. Yamamoto,Phys. Rev. A61