Efficiently simulating the quantum dynamics of physical and molecular systems represents an important near-term application of quantum computers [1,2]. Several quantum algorithms have been proposed for simulating quantum Hamiltonian dynamics, such as Lie-Trotter-Suzuki approximations [3–5] (known as Trotterizations), quantum random walks [6,7], multi-product formula [8], truncated Taylor series [9,10], qubitization [11], and quantum signal processing [12]. Among them, Trotterization can be regarded as a basic algorithm in quantum dynamics simulations. The accuracy of Trotterization, however, appears to be orders of magnitude looser than its prediction error bounded by numerical simulations [13,14]. Many algorithms have been proposed to improve the Trotterizations, e.g., product order randomization [15], time-dependent Hamiltonian simulation [16], and truncated Taylor series [9,10]. The order of Lie-Trotter-Suzuki approximation determines the error of Trotterization. For example, the first-order formula requires an O(t2/ϵ)-step Trotterization with an error of ϵ. Expectably, utilizing a higher-order formula can result in a substantial improvement of simulation precision at the expense of more quantum operations and higher circuit depth [17,18]. Another efficient multi-product algorithm was proposed by Childs and Wiebe [8], which can reach the same simulation precision as the high-order Trotterization but with a lower circuit depth, representing an efficient and precise algorithm to simulate quantum dynamics. Implementing this multi-product algorithm, however, requires a digital quantum simulator with high controllability and high precision at the hardware level.

The quantum representation for implementing the multi-product algorithm is a linear combination of unitaries (LCU) [8]. LCU can be achieved using auxiliary multi-level systems or sampling [19,20], and it supports both unitary and non-unitary operators essential for various quantum computing tasks like Hamiltonian simulation [8,9], Harrow-Hassidim-Lloyd algorithm [21], passive quantum error correction [22], and simulation of the Yang-Baxter equation [23]. The realization of LCU circuits in quantum devices requires the implementation of a sequence of multi-qubit controlled-unitary operations [24]. Realizing such multi-qubit controlled-unitary operations, in which the unitary can be arbitrarily controllable, is generally a challenging task, though multi-qubit gates have been reported in different quantum systems, e.g., superconducting qubits [25], trapped ions [26], and photons [27]. In photonic systems, controlled-unitary operations between qubits or qudits have been demonstrated [28–32]. Silicon-photonics quantum technologies [33], which can integrate photon-pair sources [34], multiphoton multidimensional entangled devices [35], large-scale programmable quantum circuits [36], and efficient single-photon detectors [37], could provide a versatile platform for LCU-based Hamiltonian dynamics simulation.

Here, we demonstrate an LCU-based quantum simulator on a programmable four-qubit silicon-photonic chip, and we further implement a modified multi-product algorithm for quantum simulation of Hamiltonian dynamics. The algorithm allows a higher accuracy of dynamics simulation compared to conventional Trotterizations as well as the original multi-product [8], and it also enables a near-deterministic success probability by combining with the oblivious amplitude amplification (OAA). To benchmark the capabilities of the LCU-based quantum simulator, we simulate a general Rabi-type Hamiltonian of a coupled nuclear spin and electron spin system. To the best of our knowledge, this work reports the first LCU-based quantum simulator, on which we show the first implementation of the multi-product algorithm for Hamiltonian dynamics simulation.

We first overview the multi-product algorithm [8] and then discuss how to improve it with the assistance of the OAA [10]. When considering the second-order Trotter approximation of a unitary operator of exp(−i∑j=1mHjt) given exp(−iHjt), we have e−i∑j=1mHjt=S1l(t/l)+E3t3/l2+E5t5/l4+⋯,(1)where S1(t)=∏j=1me−iHjt/2∏k=m1e−iHkt/2 is the second-order Trotter product, E3 and E5 are high-order error terms that are functions of Hj, t is the evolution time of the quantum system, and l is the number of iterations. Simply increasing l, one can reduce the high-order errors but cannot get rid of them. The key idea of the multi-product algorithm is to properly combine multiple low-order Trotterizations linearly so that the low-order error terms can be cancelled. For example, by properly choosing the coefficients of cn=n2/(n2−m2) and cm=m2/(m2−n2), and linearly combining two second-order Trotterizations with n and m iterations, respectively, one can directly cancel the third-order error term E3 and thus improve the approximation to the fifth-order error: e−i∑j=1mHjt=cnS1n(t/n)+cmS1m(t/m)−E5t5/n2m2+O(t7).(2)That being said, the error of quantum dynamics simulation grows in O(t5) by a linear combination of two second-order Trotterizations. Experimentally, this could be realized using an LCU circuit of the two Trotterizations. In general, one can scale up the procedure and reach an error of O(t2d+1) by the LCU of d second-order Trotterizations with a proper choice of their coefficients [38]: e−i∑j=1mHjt=∑q=1dcqS1L(q)(t/L(q))+(−1)d+1E2d+1t2d+1∏q=1d1L2(q)+O(t2d+3),(3)where cq is chosen as ∏p≠q,p=1d(L2(q)/(L2(q)−L2(p))), the cq satisfies ∑q=1dcq=1, and L(q) is the number of iterations for the q-th second-order formula [without loss of generality, we assume L(q) a monotonically increasing function].

Figure 1 illustrates the scheme for the multi-product algorithm and LCU circuit where Aj represents S1L(j) (j=1,2,⋯,d) in Eq. (3). The LCU circuit can be realized with a sequence of quantum controlled-unitary operations [24]. Implementing the multi-product Trotterization could significantly improve the simulation precision. While obtaining the required linear combination B1=∑j=1dcjAj, the LCU circuits will also return the results of other linear combination terms (B2,…,Bd) that should be in spam. These unwanted combinations are caused by the entanglement between the ancillary and data register. It will result in a low success probability of the LCU circuits. The success probability drops quickly as the number of unitaries grows. This is because ∑q=1dcq=1, and cq and cq+1 are different in sign, making the success probability 1/(∑q=1d|cq|)2 less than the unit. To overcome this, one can design the function L(q) to meet ∑q=1d|cq|≈1 to improve the probability. However, there is a trade-off between high success probability and low algorithm error, since if one wants cq to satisfy ∑q=1d|cq|≈1 and ∑q=1dcq=1 simultaneously, the only way (without loss of generality) is to choose 1≈|c1|≫|c2|≫|c3|⋯. Considering cq=∏p≠q,p=1d(L2(q)/(L2(q)−L2(p))), the original multi-product algorithm chooses L(q)={qifq≤d−1eγdifq=d(4)to be the number of iterations, where L(d−1)=eγd can be a much larger iteration number than others. This means that the Trotterization corresponding to coefficient c1 contributes most of the output while the other Trotterizations contribute little. More importantly, with such a selection, the iteration numbers for the L(1),⋯,L(d−1) orders are very small, leading to larger errors in the corresponding Trotterizations. When these are combined linearly, the errors propagate to the final result, adversely affecting the overall accuracy. Therefore, aiming for a high success probability would reduce the precision of the multi-product algorithm and result in deeper circuits. Additionally, while this selection method can satisfy the condition 1≈|c1|≫|c2|≫|c3|⋯, it will simultaneously lead to an exponential decrease in the coefficients c2,c3,⋯, requiring high controllability and high precision at the hardware level. Since the errors of Trotterizations are devised to interfere with each other destructively, their contributions to the output are preferably similar (at least intuitively).

We modified the multi-product algorithm by combining it with the OAA, as shown in Fig. 1. The OAA can amplify the success probability of LCU to near 100% in a similar way to the Grover’s algorithm. The quantum circuit of OAA is plotted in Fig. 2(b). The OAA is enabled by applying a UOAA=(−WRW†R)NW circuit on the |00⟩|ψ⟩ input state, where W operator represents the LCU circuit, R is the amplitude flip operator defined as (I−2|00⟩⟨00|)⊗I, and N refers to the number of amplification iterations. When applying a single W operator on the state, we have W|00⟩|ψ⟩=sinθ|00⟩U|ψ⟩+cosθ|Φ⊥⟩, where the operator U is the target operator of LCU, sinθ indicates the complex amplitude of the expected state, and (|00⟩⟨00|⊗I)|Φ⊥⟩=0. If we project this state with a projector |00⟩⟨00|⊗I, it will turn into a typical LCU algorithm with a success probability of sin2θ. When applying the OAA operator on the state, given U is unitary, the OAA circuit will amplify θ and return the output state as sin((2N+1)θ)|00⟩U|ψ⟩+cos((2N+1)θ)|Φ⊥⟩. This process is very similar to the Grover’s algorithm. By doing this, the OAA circuit can thus enhance the success probability from P to sin2((2N+1)arcsin(P)). The OAA circuit requires larger circuit resources compared with the LCU. In fact, it requires O((2N+1)χ) resources, where O(χ) is the resource for LCU.

One notices that, if the LCU returns a unitary operator with a success probability of 25%, the OAA can amplify the probability to 100% with N=1; if the LCU returns a non-unitary and its success probability is not exactly 25%, the OAA can also allow an efficient amplification [9]. To certify the observation, we take an example of an iteration number selection with L(q)=2q. In this case, the ∑q=1d|cq| quickly converges to 1.969 as d increases (e.g., ∑q=17|cq|≈1.969), which results in P≈25.79%. Having OAA allows a fast amplification of the success probability to be 99.93%. Note that any L(q)∝2q yields the same cq; therefore the success probability of L(q)=a·2q (a is a positive integer that can be chosen arbitrarily) can also be amplified to nearly 100%. So, given by Eq. (3), we choose the iteration number of the modified multi-product algorithm as L(q)=a·2q,(5)and amplify the output of LCU with OAA (N=1). This choice of L(q) corresponds to LCU with success probability ≈25%; OAA can thus amplify it to about 100%. With such a selection of iteration numbers, the simulation accuracy of lower-order Trotterizations has been improved, further enhancing the overall simulation precision. Taking the non-unitary effect of LCU into account, we can both achieve an error of ‖|ψ0⟩−|ψ′⟩‖≤‖E2d+1‖3t2d+1a2d2d2+d,(6)which decreases much faster than the original multi-product algorithm as d grows (|ψ0⟩ and |ψ′⟩ are the exact state and the output state of algorithms, respectively, d is the number of combined unitary gates, and ‖·‖ denotes the L2 norm for vector and spectral norm for matrix), and a nearly deterministic success probability.

Figure 3 plots the simulation errors for different algorithms. We choose L(q)=2·2q,q∈{1,2,3,4} for the modified multi-product algorithm, L(q)=1,2,3,96 for the original multi-product, and l=96 for Trotterization to make them work with the same circuit depths. It indicates that the modified multi-product can significantly improve the precision of dynamics simulation when t<40. Such improvement becomes more dominant for longer-time dynamics when a larger number of Trotterizations and a deeper circuit depth (larger q) are available. Technologically, the integrated quantum photonic platform is capable of enabling more Trotterizations and deeper circuits.

Figures 2(c) and 2(d) show the scheme of a programmable four-qubit silicon-photonic quantum simulator. One photograph of a fabricated simulator is shown in the inset. We implement the on-chip mapping of quantum dit (qudit) states [34] into multiple quantum bit (qubit) states. This allows us to implement high-fidelity and arbitrary entangling operations between the qubits, which are key to the implementation of multi-product algorithms. The LCU circuit is enabled by translating multidimensional multiphoton Greenberger-Horne-Zeilinger (GHZ) entanglement [35] to a sequence of controlled-unitaries [31]. Figure 2 shows a four-qubit version of LCU circuits.

We now discuss how to implement the LCU circuit in the quantum simulator (when N=0). The LCU is enabled by a sequence of programmable controlled–controlled–unitary operation (CCU) gates; see Fig. 2(a). The photonic quantum simulator in Fig. 2(d) integrates several key functional modules, including arbitrary local single-qubit preparation, quantum CCU, and local projective measurement. The quantum simulator monolithically integrates 451 optical components, and it was fabricated on the silicon-on-insulator platform using complementary metal-oxide-semiconductor (CMOS) processes.

We map qudit states into multiple qubit states. For example, in our photonic chip, we have the following mapping: {|0⟩qudit↔|00⟩qubit,|1⟩qudit↔|01⟩qubit,|2⟩qudit↔|10⟩qubit,|3⟩qudit↔|11⟩qubit.(7)

In this work, the notation of |0⟩, |1⟩ refers to the qubit representation.

As shown in Fig. 2(d), a pair of single photons is created in integrated spontaneous four-wave mixing (SFWM) sources. By coherently pumping four sources and wavelength demultiplexing, it returns a four-dimensional Bell state [34]. Using the mapping in Eq. (7), we rewrite it as a four-qubit entangled state: (|0000⟩+|0101⟩+|1010⟩+|1111⟩)/2,(8)where {|0⟩,|1⟩} represent the logical basis of each qubit. The basic idea is to coherently translate the entanglement in the source module to the CCU entangling gates [28–30].

First, we apply the space expansion procedure (|ij⟩→|ij,ψ⟩, |ij⟩ becomes the layer qubit, i,j=0,1) to create a six-qubit entangled state: 12(|00⟩|00,ψ⟩+|01⟩|01,ψ⟩+|10⟩|10,ψ⟩+|11⟩|11,ψ⟩),(9)where the first two qubits represent the ancillary-register states, and the two-qubit state |ψ⟩ represents the data register and it can be arbitrarily initialized by module 3. The |ij,ψ⟩ denotes the state |ψ⟩ encoded in the (i,j)-th layer qubits (will be erased eventually) that are entangled with the ancillary-register qubits.

Second, we implement local operations {A1,A2,A3,A4} on the |ψ⟩, where A represents two-qubit unitary realized by module 2. We obtain a state of 12(|00⟩A1|00,ψ⟩+|01⟩A2|01,ψ⟩+|10⟩A3|10,ψ⟩+|11⟩A4|11,ψ⟩).(10)

Third, we implement projective measurement M (module 4) on the qubits of the ancillary register and post-select its |00⟩⟨00| outcomes. This process projects the data register into an output of m1A1|00,ψ⟩+m2A2|01,ψ⟩+m3A3|10,ψ⟩+m4A4|11,ψ⟩,(11)where {m1,m2,m3,m4} correspond to the first column of the C operator in Fig. 2(a).

Finally, the coherent compression modules (modules 5, 6) implement the C′ operator, as well as the projector of |00⟩p⟨00| acting on the layer qubits |ij⟩ in |ij,ψ⟩. Then the output state is given as $$\begin{gathered} m1|00⟩p⟨00|C′A1|00,ψ⟩+m2|00⟩p⟨00|C′A2|01,ψ⟩ \\ +m3|00⟩p⟨00|C′A3|10,ψ⟩+m4|00⟩p⟨00|C′A4|11,ψ⟩, \end{gathered}$$(12)where the unitaries {C′,A1,A2,A3,A4} and complex numbers {m1,m2,m3,m4} can be arbitrarily controlled and reprogrammed. The above post-selection only takes account of the elements in the first row of the matrix representation of gate C′. Similar procedures can be applied to other projective post-selections.

Assuming the first row of C′ is {m1′,m2′,m3′,m4′} (arbitrary complex numbers), we obtain the output state of data-register as |Ψ⟩data=(m1m1′A1+m2m2′A2+m3m3′A3+m4m4′A4)|00,ψ⟩=(m1m1′A1+m2m2′A2+m3m3′A3+m4m4′A4)|ψ⟩,(13)where we have removed the layer index for clarity. Though mi and mi′ can be chosen arbitrarily, the combined coefficients mimi′=ci/∑i=14|ci| can be optimized to improve the success probability of the algorithm. When mi=mi′, the success probability of LCU becomes the largest. For this reason, we choose C′ to be CT [see Fig. 2(a)], so that the elements in the first column of C and the first row of CT are the same.

When integrating LCU and OAA techniques, the implementation can also be achieved using the configuration illustrated in Fig. 2. To achieve this, it suffices to set N and iterate the quantum circuit enclosed within the parentheses. The OAA was implemented on chip in the following way: it consists of alternating LCU circuits and |00⟩ ancillary qubits amplitude-flip gates [Fig. 2(b)]. We first fed the initial state into the simulator that had been programmed as LCU circuits and then performed quantum tomography on the output state. After that, we performed a |00⟩ ancillary qubits amplitude-flip operation on the most probable pure state in the measured state and re-fed it into the simulator that had already been programmed for the next LCU circuit. By repeating this offline process, we were able to implement the OAA circuit experimentally. Though the device comprises a single layer of the repetitive structure, it is sufficient to validate the algorithm.

We benchmark the modified multi-product algorithms in our LCU quantum simulator to simulate the dynamics of a Hamiltonian describing the Rabi-type interactions of an electron spin and a nuclear spin in the rotating frame of external electromagnetic field frequency. The Hamiltonian H is written as H=Ω2σxe+δ2σze+|1⟩e⟨1|⊗(E1|0⟩n⟨0|+E2|1⟩n⟨1|),(14)where σx and σz are Pauli matrices, subscript and superscript of n and e represent the nuclear and electron spins, respectively, |0⟩n and |1⟩n represent nucleus spin up and down states with hyperfine coupling strength of E1 and E2, and Ω and δ are the amplitude of external electromagnetic field and fine energy splitting, respectively. The Hamiltonian is defined by the following parameters in the experiment: E1=0.3,E2=0.7,Ω=0.2,δ=0.5. The parameters are chosen so that distinct deviations between the time evolution of different eigenstates can be observed in our quantum simulator, but they do not refer to actual quantum systems. This way, we can efficiently test the implementation of multi-product algorithms in experiments. We use the LCU quantum simulator to emulate the unitary U(t)=exp(−iHt). The operators {A1,A2,A3,A4} are chosen as {S12(t/2), S14(t/4), S18(t/8), S116(t/16)}, respectively. The S1(t) is a second-order Trotter product of e−iH1t/2e−iH2te−iH1t/2, where H1 is the first two terms in Eq. (14); H2 is the last term.

Figure 4(a) shows experimental probability distributions for the |11⟩ state measured in the LCU simulator, emulating the dynamic evolution of a |01⟩ initial state governed by the Hamiltonian. We compare the performance of the modified multi-product algorithm, which combines S12(t/2), S14(t/4), S18(t/8), and S116(t/16), with the standard Trotterization [S116(t/16)] and the original multi-product algorithm [combining S1(t), S12(t/2), S13(t/3), and S116(t/16)]. This comparison is made without OAA to validate the enhanced precision achieved through the new iteration numbers L(q), which is independent of whether OAA is used or not. It shows that, for t≤30, the modified multi-product algorithm can approximate the exact solution more precisely than the other two algorithms. Their discrepancies to the exact result are plotted in Fig. 4(b). For t>30, both original and modified multi-product algorithms show deviation from the exact solution due to high-order error and the non-unitary effect.

After confirming the enhanced accuracy of the modified multi-product algorithm, we further assess the improvement of its integration with OAA on success probability. Figures 5(a)–5(h) report experimental dynamics for the four logical states measured in the quantum simulator, when the initial state is chosen as 0.3|00⟩+0.7|01⟩, i.e., 0.3|0⟩e|0⟩n+0.7|0⟩e|1⟩n. Both modified multi-product algorithms before OAA [Figs. 5(a)–5(d)] and after OAA [Figs. 5(e)–5(h)] are implemented. Experimental results well agree with their theoretical predictions, and means of classical statistical fidelities of Fc=0.996±0.005 and Fc=0.996±0.004 are obtained, respectively, for the modified multi-product algorithm with and without OAA. The Fc is defined as (∑ipiqi)2, where pi and qi are measured distribution and theoretical distribution, respectively. Moreover, we measured the success probability of algorithms, as shown in Fig. 5(i). Experimental results show that when t≤30, the success probability (Pr) is about 26.33% in the modified multi-product algorithm before the OAA step. We then implement OAA and obtain the amplified probability. Note in this case, OAA requires only one step of amplification, i.e., N=1. The OAA circuit in Fig. 2(b) can be simplified to WRW†RW, and compiled in the simulator. By implementing the OAA, we amplify the success probability to 99.79% at t≤30, approaching a deterministic probability. All measured probabilities are consistent with theoretical predictions. When t>30, due to the non-unitary LCU, the high-level success probability cannot be retained.

Figure 6 shows a generalized scheme for the implementation of LCU circuits: ∑j=1dcj∏i=1n⊗Aij|ψ⟩i,(15)where |ψ⟩=∏i=1n⊗|ψ⟩i is the initial state and Aij is the unitary operator. The LCU circuit in Fig. 6(a) can be prepared by using the d-dimensional (n+1)-photon Greenberger-Horne-Zeilinger (GHZ) state, which is defined as |GHZ⟩n+1,d=∑i=0d−1|i⟩⊗(n+1).(16)

The scheme employs a sequence of processes, referred to as “space expansion-local operation-coherent compression,” as shown in Fig. 6(b). Utilizing advanced silicon photonics, a 10-photon LCU-quantum simulator is feasible [31]. Moreover, iteratively repeating the colored segment in Fig. 6(b) by M iterations, like the structure within the brackets in Fig. 2(c), and applying swapping operations between the replicated structures to exchange the states |ψ⟩i generated by P^ and |ϕ⟩i, one can construct the extended circuit as Fig. 6(c). This circuit features an alternating cascade of M layers of generalized controlled operations. If the target qudit remains unchanged and only single-qudit gates are inserted between the controlled gates, we can set all the O^ operations within this layer to be the same.

We estimate the required physical resources, that is, the number of on-chip phase-shifters and classical controls when scaling up the LCU-quantum simulator. This scheme requires several nd multi-qubit state preparation units, nd local unitary operation units, n layer information erasure units (nd multi-qubit measurement units), and n+1 multi-qubit measurement units. Each unitary operation unit needs Nunit=d2−d phase shifters, and each multi-qubit state preparation (measurement) unit requires Npre=Nmes=2(d−1) phase shifters. Given the state |GHZ⟩n+1,d, we can estimate the number of the classical controls for the M-layer LCU device as N=ndNpre+nMdNunit+(nd+n+1)Nmes=(d−1)(nMd2+4nd+2n+2).(17)

Overall we have proposed and demonstrated a modified multi-product algorithm. By combining the multi-product and OAA algorithms, it can improve the simulation precision and relax the success probability constraints. The algorithm is certificated experimentally in a small-scale silicon-photonic quantum simulator. Our experiments, even at a small scale, unveil error propagation’s influence on quantum computation, shaping strategies to optimize quantum device controls for large implementations. Moreover, the linearly combined operations are not necessarily unitary (mostly, they are non-unitary), implying that LCU can simulate non-Hermitian systems’ dynamics. The quantum algorithms and quantum hardware could be scaled up to provide an efficient solution to simulate the dynamics of physical and molecular systems.