Tensor-Network-Based Distributed Quantum Dynamics on Independent Quantum Computers Anurag Dwivedi,1, 2 Melissa C. Revelle,3 Daniel S. Lobser,3 Brian K. McFarland,3 Edward C. Tortorici,3 Christopher G. Yale,3 Susan M. Clark,3 Philip Richerme,4, 2 and Srinivasan S. Iyengar2, 1, ∗ 1
arXiv:2606.11579v1 [quant-ph] 10 Jun 2026
2
Department of Chemistry, Indiana University, Bloomington, Indiana 47405, USA Indiana University Quantum Science and Engineering Center, Bloomington, Indiana 47405, USA 3 Sandia National Laboratories, Albuquerque, New Mexico 87123, USA 4 Department of Physics, Indiana University, Bloomington, Indiana 47405, USA (Dated: June 11, 2026)
We present an approach based on tensor networks for distributed quantum computing simulation of chemical wavepacket dynamics in a continuous variable representation. The central idea is that the tensor-network representation of the multidimensional time-evolution operator naturally induces an elevated Hilbert space in which the global dynamics decomposes into a collection of independent lower-dimensional propagations. This transformation converts an entangled quantum evolution into a set of parallel computational tasks that can be executed asynchronously across heterogeneous quantum and classical computing architectures. The resulting formalism establishes a direct connection between tensor-network decompositions, uniformly controlled quantum circuits, and asynchronous distributed quantum computing. The approach is developed with a goal towards hybrid quantum/classical implementation, and is appropriate for a general heterogeneous mixture of quantum hardware systems; it is demonstrated here on ion-trap quantum computers. The experimental realization of the asynchronously distributed quantum processes that arise from the tensor-network decomposition are carried out on the Sandia National Laboratories’ trapped-ion quantum computer, where the circuits are compiled using native partial-entangling XX(θ) gates, reducing the expected two-qubit gate infidelity by more than 30% relative to conventional fully entangling decompositions. We demonstrate the methodology by quantum computing the vibrational spectral properties of a small protonated water cluster that shows critical quantum nuclear behavior. Such water cluster systems have been found to be challenging for classical computation and for experimental action spectroscopy and here, for the first time, we provide results for vibrational spectroscopy that are in −1 agreement with the respective classical results to within 4cm , thus allowing for the potential for spectroscopic accuracy from quantum computations. To compute molecular vibrations efficiently, we also introduce here a modified version of the phase estimation algorithm to directly obtain energy differences (instead of absolute energies) that are most relevant as spectroscopic observables. More broadly, we note that while the methods developed here have been demonstrated for chemical dynamics and vibrational spectroscopy, the general prescription of introducing an elevated space where the dynamical components are decoupled and treated as a parallel tasks, with individual components being represented as continuous variables, is also appropriate for other continuous systems such as entangled optical beams and correlated fluid elements. These extensions provide a pathway toward scalable distributed quantum simulation on future heterogeneous quantum-computing architectures.
I.
INTRODUCTION
Simulating the quantum dynamics of interacting many-body systems remains a grand challenge[1–3] in multiple areas of the physical sciences. This difficulty is particularly pronounced in molecular systems, with coupled nuclear degrees of freedom, such as in hydrogenbond networks in water, as well as in hydrogen, hydride, and proton transfer reaction problems in biological and atmospheric systems, where an accurate treatment of high-dimensional wavepacket dynamics on correlated potential energy surfaces is needed. Despite significant advances in classical computing algorithms, such as the multiconfigurational time-dependent Hartree (MCTDH) approximation[4–10] the exponential scaling of both com-
∗
putation and storage, with system size[10–12], continues to limit predictive simulations of such complex systems. At the same time, quantum computing appears to offer a fundamentally different route to address such challenges by encoding the wavefunction directly on a quantum register and evolving these under unitary dynamics[13–15]. In principle, algorithms based on quantum phase estimation (QPE)[13] enable access to spectral features through time evolution and Fourier transforms. However, practical implementations of such approaches on near-term quantum devices are severely restricted by circuit depth[16, 17], entangling gate counts[18, 19], and limited qubit connectivity[20]. These limitations are particularly acute for quantum dynamical simulations, where repeated application of time-evolution operators may be needed. In this publication, we demonstrate a tensor-networkbased framework for quantum wavepacket dynamics and show that this formalism yields an asynchronous and
2
[1]
ϕα (x1 )
[1]
α
β
[2]
ϕα (x2 )
[2]
Uβ
Uβ
Ion-trap-1
Ion-trap-2
FIG. 1: The essence of the distributed algorithm presented here which appears from a tensor network formalism. There is one set of parallel streams of computational tasks for each value of the pair of variables (α, β). Thus, there are [α × β × number of nuclear dimensions] parallel computing tasks, each executed here as a separate ion-trap computation, possible due to elevation of the entanglement variables. This yields an enormous computational advantage as the number of dimensions increase. Future developments will involve a heterogeneous set of classical and hardware systems all runing these parallel set of tasks.
parallelizable quantum computing algorithm, which we implement on Sandia National Labs’ QSCOUT ion-trap quantum computing system[21, 22]. This approach enables a decomposition of the time-evolution operator into a structured set of lower-dimensional operations suitable for distributed quantum execution. By expressing both the wavefunction and propagator in a tensor network form, we show that the global unitary evolution can be recast as a block-diagonal operator in an elevated, higherdimensional space, whose block-components correspond to conditionally applied product-state time-evolution operators. This representation establishes a direct mapping between tensor network contractions, uniformly controlled quantum circuits[23, 24], and distributed/hybrid quantum computing architectures where entanglement indices are used to create an elevated space and are also used as control quantum registers. As a result, the overall dynamics is exactly decomposed into a family of independent, reduced-dimensional quantum tasks that may be executed in parallel[16, 17, 21, 25, 26] across multiple quantum processors with potential for use as a hybrid quantum/classical algorithm. Figure 1 captures the essence of this idea. In addition to the above, we have also utilized here a quantum resource-optimized variant of the phase estimation algorithm tailored for quantum dynamical simulations on NISQ environments. Here, the time-evolution stage is implemented on quantum hardware, while the Fourier transform is performed classically. This hybrid approach significantly reduces circuit complexity while preserving access to spectroscopic observables, including energy differences and vibrational features derived from time-correlation functions. Together with the tensornetwork decomposition, this enables a practical pathway
for performing quantum dynamical simulations on nearterm devices with limited coherence and connectivity. We demonstrate this framework through the simulation of vibrational dynamics in a protonated water-wire system, implemented on a distributed set of trappedion quantum processors. The water-cluster systems have been considered to be a major challenge to both experiment[27, 28] and theory[29–34]. The study of protonated water clusters have deep fundamental as well as applied implications due to prevalence in biological ion-channels[35], inside enzyme active sites[36–38], in polymer electrolyte fuel cells[39] and in the earth’s atmosphere[40–54]. The study of such systems is complicated by, multi-dimensional quantum nuclear effects arising from hydrogen bonded networks as evidenced by the Grothhuss mechanism of proton transfer[55]. By decomposing the multidimensional Hamiltonian for such systems into a collection of effective one-dimensional components as formally allowed by the aforementioned tensor network formalism, we realize parallel quantum simulations whose combined results reproduce the key spectral features of the system. Indeed, the agreement between quantum hardware results and exact classical simulations is of the order of ≈ 4cm−1 . This is at the spectroscopic accuracy level of agreement between the quantum simulation and classical simulation and validates the approach and suggests its potential for scaling to more complex systems in future. This paper is organized as follows: In Sections II and III, we present the theoretical framework underlying the tensor-network-based parallel quantum dynamics approach and the phase estimation algorithm developed for resource-optimized quantum vibrational spectroscopy. In particular, Section II introduces the tensornetwork decomposition of the multidimensional propagator, the resulting block-diagonal structure in the elevated representation, and the distributed quantum execution strategy for parallel implementation on quantum hardware. Section III describes the modified phase estimation protocol for extracting vibrational spectra with reduced circuit depth requirements. The experimental implementation of the trapped-ion quantum simulations is described in Section IV. In Section V, we demonstrate the framework through simulations of proton-transfer vibrational dynamics in a protonated water-wire system implemented on a trapped-ion quantum processor at Sandia National Laboratories. In particular, we present the reduced-dimensional potential energy surfaces, the quantum wavepacket dynamics obtained from trapped-ion quantum hardware, vibrational spectra extracted from the quantum dynamics, and a detailed comparison between the quantum-computed spectra and exact multidimensional calculations, including a quantitative error analysis. Finally, conclusions and future directions are discussed in Section VI.
3 II. TENSOR NETWORK BASED PARALLEL QUANTUM DYNAMICS SIMULATIONS ON QUANTUM COMPUTERS
ϕ
[1]
α1
ϕ
′
ᾱ
j=1
ᾱ
where, the dependence on the discretized version of the continuous position basis, x̄ ≡ (x1 , x2 , . . . , xN ), is encoded in tensor elements parameterized by xj in the
U
β̄
j=1
β̄
where the indices β̄ ≡ (β0 , β1 , . . . , βN ), bounded by η̄ ≡ (η0 = 1, η1 , . . . , ηN −1 , ηN = 1), encode the entanglement structure of the operator across different degrees
αN -2
···
βN -2
ϕ
[1]
U
[N -1]
αN -1
ϕ
′
x2 β1
x1
′
xN -1 β2
[2]
U
x2
[N ]
xN
[N -1]
βN -1
U
xN -1
[N ]
xN
FIG. 2: Time evolution of a wavepacket in a tensornetwork representation as given by Eq. (3). The initial state is parameterized by entanglement variables ᾱ, while the time-evolution propagator is characterized by β̄. |β1 ⟩
|β1 ⟩ |χ⟩
[1]
[2]
[1]
|χ⟩
U1,β1 ⊗ Uβ1 ,1
[2]
U1,1 ⊗ U1,1
[1]
[2]
U1,2 ⊗ U2,1
(b)
(a)
FIG. 3: (a) Uniformly controlled gate-like structure with single control qubit encoding the bond-dimension index β1 . (b) An illustration of Figure (a) with maximum bond dimension η1 = 2 (i.e., upper bound to the quantity β1 ). Also see Eq. (5). Critical to note that the [1] [2] [1] [2] individual blocks U1,1 ⊗ U1,1 and U1,2 ⊗ U1,2 are direct product operators that act on product states of dimensions “[1]” and “[2]”. of freedom. Each tensor U [j] represents the reduceddimensional effective propagator that evolves the corresponding lower-dimensional function ϕ[j] in Eq. (1) and ′ Uβ̄x̄,x̄ represent the respective product state operators. These operators propagate the MPS state as per N Z ′ ′ X Y [j]x x [j]x dx′j U j j ϕαj−1j ,αj χt (x̄) = βj−1 ,βj
j=1
ᾱ,β̄
[j]x
superscript of, ϕαj−1j ,αj . These tensors ϕ[j] , can thus be interpreted as reduced-dimensional functions associated with each coordinate xj . The superscripts in square x̄ brackets [·] label the tensor cores and the terms ϕᾱ refer to the associated product states. The summation indices ᾱ ≡ (α0 , α1 , . . . , αN ) are commonly referred to as the bond dimensions (or equivalently, entanglement dimensions or Schmidt ranks). The time-evolution operator Û is expressed in the matrix product operator (MPO) form [58, 59] as η̄ η̄ N ′ X Y X ′ [j]x x x̄′ Û x̄ = Uβj−1j ,βjj = Uβ̄x̄,x̄ , (2)
···
′
x1
The central challenge in quantum dynamical simulations lies in the efficient representation and implementation of the time-evolution operator acting on a high-dimensional Hilbert space. For systems with multiple coupled degrees of freedom, such as wavepackets for molecular systems, the dimensionality of the Hilbert space grows exponentially with system size[1, 2, 10], rendering direct implementations of the time-evolution operator to be prohibitively expensive on classical hardware architectures. To address this challenge, we begin with a tensornetwork-based decomposition of quantum dynamics that exposes a structured representation of the propagator amenable to distributed quantum execution. The key idea is to exploit the entanglement structure of both the wavefunction and the propagator to create a global evolution strategy that decomposes exactly into a collection of lower-dimensional operations and may then be executed in parallel on multiple hardware systems. We begin by representing the initial wavepacket, χ0 , as a matrix product state (MPS) [56]. (Although we use MPS states here, this is not a restriction and other graph based tensor-network schemes[57] can also be used within our formalism.) In this representation, a wavepacket with N nuclear degrees of freedom is expressed as a multiconfigurational expansion [6, 8–10] that may be graphically represented as in the top row (orange) in Figure 2: N X x̄ X Y [j]x ϕᾱ , (1) χ0 (x̄) = ϕαj−1j ,αj =
α2
[2]
=
XZ
′
′
x̄ dx̄′ Uβ̄x̄,x̄ ϕᾱ .
(3)
ᾱ,β̄
This expression corresponds to a contraction of the MPO with the MPS, resulting in a new state whose bond structure is governed jointly by ᾱ and β̄. These equations are illustrated in Figure 2. A.
Block-diagonal structure of the tensor-network propagator in an elevated representation
The propagator in Eq. (2) can be recast as a blockdiagonal matrix acting on an enlarged Hilbert space. We illustrate this idea here for the case N = 2 and β1 = 2, where the time-evolution operator, [1]
[2]
[1]
[2]
U1,1 U1,1 + U1,2 U2,1
(4)
4
ϕ[1] (x1 , 0)
[1]
ϕ[2] (x2 , 0)
[2]
[N ]
U1,1 ⊗ U1,1 ⊗ · · · ⊗ U1,1
0 .. . 0
β1
[1]
U1,β1
QC-3
QC-4
U1,(β1 =1)
[1]
U1,(β1 =2)
[1]
U(β1 =1),1
[2]
U(β1 =2),1
|q0 ⟩ |q1 ⟩ |q2 ⟩
|q0 ⟩ |q1 ⟩ |q2 ⟩
|q0 ⟩ |q1 ⟩ |q2 ⟩
|q0 ⟩ |q1 ⟩ |q2 ⟩
(β1 =1),1
[1] [2]
[1]
···
0
..
.. .
.
···
[1] U1,η1 ⊗
[2] Uη1 ,η2 ⊗
[N ]
· · · ⊗ UηN −1 ,1
product operators where each operator acts on lower dimensional function as noted in Eqs (8) and (3).
β1 =1,1
[2]
φβ1 ,1 (x2 , t)
FIG. 4: The circuit in Fig. 3 can be distributed, in principle, across four quantum computers, where each constituent sub-circuit can be executed simultaneously. In this paper, this is all done sequentially on one quantum computer using three qubits. The aggregation (the double vertical arrow at the bottom of the figure) is also done classically here through measurement, and the final state inherits the entanglement from the unitary above. As described in Section II C, such entanglement arises from the potential energy surface.
may be rewritten in an enlarged Hilbert space as, " # [1] [2] U1,1 ⊗ U1,1 0 Û = [1] [2] 0 U1,2 ⊗ U2,1
(5)
Here, the size of the resultant vector space has been increase by a factor corresponding to the extent of entanglement, defined here by the upper bound to the quantity β1 , say η1 . That is, the size of the unitary in Eq. (5) is now D2 ∗ η1 for D basis discretizations per dimension. Thus, iwhen Eq. (5) acts on a product initial state h ϕ[1]x1 ϕ[2]x2 , we obtain
0
U U(β 1 =2),1
β1 =1,1
β1
0
···
[2][2]
U U(β 1 =1),1
(β1 =1),1
φ1,β1 (x1 , t)
[2]
[2][2]
UU1,(β1 =2)
.. .
···
[N ] ⊗ U1,1
FIG. 5: Block diagonal form of the overall unitary in arbitrary dimensions. The dimensionality of the probQ lem is expanded by a factor of [ η̄] (see Section II B) which results in a block diagonal representation with each block acting on a product state. As noted ear ′ ′ QN [j]xj xj x̄,x̄ lier, the individual blocks, Uβ̄ ≡ U j=1 βj−1 ,βj are
QC-2
[1] [2]
[2] U2,1 ⊗
[2]
Uβ1 ,1
QC-1
UU1,(β1 =1)
0 [1] U1,2 ⊗
The block-diagonal operator in Eq.(5) can be implemented directly as a quantum circuit using uniformly controlled gates (UCG)[23] with a single set of control qubits encoding the bond index (β1 ), as shown in Fig. 3. However, an important feature of Eq.(5) is that each diagonal block is itself a tensor product of lower-rank unitary operators. See Eq. (6). This tensor-product structure exposes an additional level of parallelism beyond that provided by the UCG representation. Specifically, the constituent operators within each diagonal block can be executed as independent computational streams, as illustrated in Fig. 4. While the UCG formulation provides a compact circuit-level representation, it still requires the implementation of the corresponding controlled operations. In contrast, the block-diagonal decomposition naturally lends itself to a distributed execution environment where the constituent operators within each diagonal block may be evolved independently and in a completely asynchronous manner on separate hardware resources. Bipartite systems with greater entanglement are discussed in Appendix A. The principal advantage of this formulation is therefore its compatibility with large-scale parallel execution. Multiple quantum processors, or quantum processors coupled to classical high-performance computing (HPC) resources, can be employed to evaluate the independent streams associated with the different tensor-product sectors. This makes the approach particularly well suited for hybrid quantum–HPC architectures.
B.
Entanglement Entropy as a computational
[1]x [2]x [1] resource: Formal Scaling of Distributed Quantum [2] 1 2 Execution from Tensor-Network Decompositions ϕβ1 =1 ϕβ1 =1 0 U1,1 ⊗ U1,1 ϕ[1]x1 ϕ[2]x2 = [1]x1 [2]x2 [1] [2] ϕβ1 =2 ϕβ1 =2 ϕ[1]x1 ϕ[2]x2 0 U1,2 ⊗ U2,1 The main features of the algorithm presented above are [1] [2] depicted in Figures 5 6 and 7 for arbitrary dimensions. As U1,1 ϕ[1]x1 ⊗ U1,1 ϕ[2]x2 before, these figures are obtained from the block-diagonal (6) = [2] [1]x1 [1] [2]x2 nature of the propagator in an elevated Hilbert space U1,1 ϕ ⊗ U1,2 ϕ that contains both the system degrees of freedom and the entanglement degrees of freedom. In this section, we where the elevated vector space is clearly marked based present a formal analysis of the reduction in complexity on the additional index β1 .
5 ...
1 .. .
k
(N − 1) log2 η controls
.. .
.. .
.. .
(N − 1) log2 η − 1
...
(N − 1) log2 η
...
[j] j Uβj-1 βj
N
|χ⟩
(a) [1]
[2]
[N ]
[1]
U1,1 ⊗ U1,1 ⊗ · · · ⊗ U1,1
t
[2]
[N ]
...
U1,2 ⊗ U2,1 ⊗ · · · ⊗ U1,1
[1]
[N ]
[2]
U1,η1 ⊗ Uη1 ,η2 ⊗ · · · ⊗ UηN −1 ,1
(b)
FIG. 6: (a) Uniformly controlled gate (UCG) representation of the block-diagonal operator shown in Fig. 5.
[1] ϕ
U
[1]
···
[2] ϕ
β 1 U
β 2
[2]
[N -1] ϕ
β N -2
···
···
U
[N ] ϕ
β N -1
[N -1]
U
[N ]
QC-(1,1,1)
QC-(2,1,1)
QC-(2,1,2)
QC-(N-1,1,1)
QC-(N-1,1,2)
QC-(N,1,1)
[1] U 1,(β1=1)
[2] U (β1=1),(β2=1)
[2] U (β1=1),(β2=2)
[N −1] U (β1=1),(β2=1)
[N −1] U (β1=1),(β2=2)
U
q2
q2
q2
q2
q2
q2
q1
[1]
q1
U1,1 [1]
U1,1
[2]
q1
U1,1 [2]
U1,1
[2]
q1
U1,2 [2]
U1,2
[N −1]
q1
U1,1
[N −1]
U1,1
[N −1]
U1,2
[N −1]
U1,2
[N ] (βN −1=1),1
q1
[N ]
U1,1 [N ]
U1,1
q0
q0
q0
q0
q0
q0
QC-(1,1,2)
QC-(2,2,1)
QC-(2,2,2)
QC-(N-1,2,1)
QC-(N-1,2,2)
QC-(N,2,1)
[1] U 1,(β1=2)
[2] U (β1=2),(β2=1)
[2] U (β1=2),(β2=2)
[N −1] U (β1=2),(β2=1)
[N −1] U (β1=2),(β2=2)
U
q2
q2
q2
q2
q2
q2
q1
[1]
q1
U1,2 [1]
U1,2
q0
[2]
q1
U2,1 [2]
U2,1
q0
φ
[1]
[2]
q1
U2,2 [2]
U2,2
q0
β 1 φ
[2]
[N −1]
q1
U2,1
[N −1]
U2,1
q0
β 2
···
[N −1]
U2,2
[N −1]
U2,2
q0
β N -2 φ
[N -1]
[N ] (βN −1=2),1
q1
[N ]
U2,1 [N ]
U2,1
q0
β N -1 φ
[N ]
FIG. 7: The circuit in Fig. 6 can, in principle, be distributed across several quantum computers, and this figure shows the full-scale parallel generalization of our algorithm. The ion-traps are indexed here using the dimension (first index), and the two entanglement indices following the dimension. So for example, “QC-(2,1,1)” implies the specific quantum hardware system (trapped ions in our case) that computes the second dimension for entanglement P 2 variables β1 = β2 = 1. As is clear, Q the number of parallel streams (dotted boxes) is equal to i ηi . By contrast the number of blocks in Fig. 5 is [ η̄] as is the number of constraint operations in Figure 6.
and the extent of parallelism afforded by this algorithm. In general, for a system with N dimensions and D basis discretizations per dimension, the formal computational scaling of quantum propagation involves matrix-vector operations that are potentially O(DN ). The algorithm presented here provides an efficient quantum computing strategy for this problem that automatically reduces to
an asynchronous and distributed set of quantum processes that can be executed on a hybrid set of quantum and classical hardware platforms. For an N -dimensional system, the block-diagonal structure for the evolution operator is illustrated in Fig. 5. As before, this structure maps onto a uniformly controlled gate (UCG) configuration where the entangle-
6 ment indices β̄ are encoded in control registers. The associated quantum circuit is provided in Figure 6.Q The numberQ of required control qubits may scale as {ln [ η̄]}, where [ η̄ ≡ η1 η2 · · · ], that is product of the individual maximum bond dimensions, or the maximumm entanglement volume within this representation.The quantity η̄ is the upper bound in Eq. (2). Q For the case where all bonds are cut, the quantity {ln [ η̄]} represents the maximum possible entanglement entropy of the associated tensor network state. Thus, the number of control qubits may scale as the maximum possible entanglement entropy of the system. The individual operations in each Q block may be viewed as a collection of [ i ηi ] independent unitary blocks, where each block corresponds to a ′ ′ x̄ single operation (Uβ̄x̄,x̄ ϕᾱ ) acting on a product-state sector of dimension DN . The key observation, however, is that each block is itself composed of a direct product of N lower-rank unitary operators, one associated with each physical dimension. As shown in Eqs. (3) and (6) (see also Eq. (12) below), these operators act independently on product states and therefore admit a highly decomposed computational structure. However, when the individual one-dimensional streams are propagated in parallel, one can parse these our separately and the number of such individual streams is bounded by 2η1 + 2ηN +
N −1 X i=2
ηi ηi+1 ≈
X 2 ηi
(7)
i
for large N , as shown in Fig. 7. This is because while Qthe block-diagonal form of Fig. 5 appears to contain N i ηi independent unitary operations that could be executed in parallel, many of these operations are repeated across different blocks and therefore need not be recomputed when treated in an asynchronous manner. For example, [1] the operator U1,1 appears in every block for which the remaining bond indices vary while β1 = 1. Consequently, this operator need only be evaluated once and reused wherever it appears. Accounting for such redundancies, the number of distinct computational tasks is reduced to P a maximum of i ηi2 operations, each corresponding to an O(D1 ) operation of the form Z ′ ′ [j]x x [j]x dx′j Uβj−1j ,βjj ϕαj−1j ,αj (8) which also follows directly from Eq. (3). A similar reduction may also be possible in the UCG implementation, but will require additional control gates which in current implementations may be impractical. This parallel and hybrid aspect is emphasized in Figure 7. The family of operators may be used to spawn a set of asynchronous tasks that can be performed on separate ion-trap quantum computers. In the current implementation these O(D1 ) operations are performed on a single ion-trap, sequentially, but this is not the only possible implementation. In principle, these separate
tasks can be performed, as noted, in a completely asynchronous way, on multiple quantum computers, some of which may involve ion-traps and some may involve other quantum architectures. The labeling scheme of quantum hardware in Figure 7 suggests this generality. For example, the label “QC-(2,2,1)” implies the specific quantum hardware system that computes the action of the circuit corresponding to the second dimension for entanglement variables β1 = 2 and β2 = 1. In summary, a uniformly controlled gate implementation of Figure 6 may require the number Q of ancilla to be equal to up to a maximum of ln [ η̄], which is the maximum possible entanglement entropy that can be supported by the tensor network. By contrast, an asynchronousPand distributed implementation needs a maximum of i ηi2 number of parallel streams, which may in general scale as the area of the entanglement hyperspace defined by edges {ηi } to capture the quantum dynamics, as can be seen in Figure 7. As hardware systems evolve, future versions of this algorithm can include both aspects by constructing parallel tasks that may contain a small set of uniformly controlled gate implementations. We wish to add one caveat to this argument. Although the general trend is clear, one may ask how the method scales for systems that truly contain exponential entanglement. In such cases {ηi } may grow exponentially with system size, and while the current approach still partitions the full problem into manageable chunks, the number of such independent processes could grow exponentially. For such problems, a hybrid approach, such as the one mentioned at the end of the previous paragraph, could be useful to consider. However, for most practical problems[60] we expect the proposed methodology to produce efficient algorithms. Thus, the decomposition introduced here provides a pathway for mitigating two of the primary limitations of near-term quantum devices: (a) circuit depth and (b) limited qubit connectivity. By transforming a global, potentially highly entangled unitary into a family of ′ ′ x̄ product-state operations (Uβ̄x̄,x̄ ϕᾱ ) that may be thought to be conditioned by the entanglement index β̄, the approach reduces both the entangling gate count and the effective circuit depth of each sub-task. The resulting structure is inherently parallel, making it well-suited for emerging distributed quantum computing architectures and hybrid Quantum-HPC in which multiple quantum (and classical) processors are coordinated to perform large computations. Figure 6 shows the extent parallelism possible within this formalism.
C.
Obtaining Uβ̄
from electronic structure
As specified, the formalism here pertains to chemical dynamics, where the overall Hamiltonian may be written as a sum of the nuclear kinetic energy and potential energy operators represented here on a grid representation,
7
′
Ĥ(x̄, x̄ ) =
X
K(xi , x′i ) + δ(x̄ − x̄′ ) V̂ (x̄),
(9)
i
where V̂ (x̄) is the (local) potential energy operator, diagonal in the coordinate representation and obtained from electronic structure calculations. The varianble x̄ now represents a nuclear geometry where this electronic structure is computed. In contrast to the kinetic term, the potential energy is naturally not separable across dimensions, which introduces additional complexity. The time evolution operator Û = e−iĤ∆t/ℏ is approximated using a first-order Trotter decomposition [61, 62]: ( ) Y ′ e−iĤ∆t/ℏ = e−iV̂ (x̄)∆t/ℏ e−iK̂(xi ,xi )∆t/ℏ + O(∆t2 ) i
(10) and a key complication arises from the non-separable nature of the potential energy operator in the coordinate representation. To address this, we express the potential propagator also as an MPS, η̄ N Y X [j]x Vβj−1j ,βj , (11) e−iV̂ (x̄)∆t/ℏ = β̄
j=1
and using this expression, Eq. (2) may now be rewritten for a product state initial wavepacket as N Z ′ X Y ′ ′ [j]x ′ −i K̂(x ,x )∆t/ℏ [j]x i i dxj e Vβj−1j ,βj ϕ j χt (x̄) = β̄
=
XZ
j=1 ′
′
dx̄′ Uβ̄x̄,x̄ ϕx̄ .
(12)
β̄
It is thus clear that all the entanglement in the propagated state arises from the potential. Again, as noted above, this paper deals with an MPS representation of wavepackets and potential propagators, but also allows more general graph-based tensor-network treatments as discussed in Ref. 57. (Ref. 57 provides general expressions for tensor network treatment of potential propagators when graph-based molecular fragmentation techniques are used to obtain the potential energy surface in high dimensions[63] with post-Hartree-Fock accuracy.) However, the factorization of the potential propagator above brings about an extremely subtle aspect that ′ [j]xj needs to be further discussed. The individual Vβj−1 ,βj propagators may be rewritten as ′ ′ [j]x [j]x Vβj−1j ,βj = exp −ıVβj−1j,βj t/ℏ , wherethe effective reduced dimensional entangled poten′ [j]xj tials Vβj−1 ,βj may now be complex thus making the in-
due to the entanglement enforced by the potential. For n ′o x̄,x̄ example, if the Uβ̄ turn out to be completely unitary, there is no effective population redistribution between the various modes after each propagation step, thus rendering all dimensions to be orthogonal. Thus non-unitarity of the individual propagated branches of the tensor network is a critical feature in this algorithm.
D.
Key experimental advances that make the action of Uβ̄ efficient on ion-traps
In addition to the theoretical aspects presented here, it is critical to mention an experimental advance that makes the associated calculations possible. Experimental details are presented in Section IV. As might be clear, now the key elements of this algorithm to a fam reduces ′ [j]xj xj [j]xj ily of parallel operations given by Uβj−1 ,βj ϕαj−1 ,αj . In Ref. [64] we present a general Quantum Shannon Decomposition (QSD) algorithm to perform these individual operations. In this QSD circuit though, arbitrary decompositions of three-qubit unitary propagators obtained during the implementaion of operations [j]x x
′
[j]x
Uβj−1j ,βjj ϕαj−1j ,αj in this paper, require 24 fully-entangling CNOT gates or XX(π/2) gates. In Section IV and in Appendix D, we show that we extend the QSD decomposition to allow XX gates with arbitrary angle θ and reconstruct the desired unitary at each timestep. The result is a circuit for each dimension with 6 fully-entangling XX(π/2) gates, and 18 partial-entangling XX(θ) gates, with 0 < θ < π/2. Partial-angle gates are implemented by keeping the gate time constant (to satisfy phase-space closure constraints) and scaling the overall laser amplitude to generate the desired entanglement[65]. This gate count is lower than that obtained from the standard QSD algorithm[64]. In the following section, we combine this tensornetwork decomposition with a resource-optimized variant of quantum phase estimation to construct a practical algorithm for extracting dynamical and spectroscopic information from quantum simulations.
III. MODIFIED PHASE ESTIMATION ALGORITHM FOR RESOURCE OPTIMIZED QUANTUM VIBRATIONAL SPECTROSCOPY
′
dividual operators, Uβ̄x̄,x̄ potentially non-unitary. However, this non-unitarity is critical in that it only arises
In quantum phase estimation, illustrated in Fig. 8, a unitary time-evolution operator is used to construct the
8 Prop
Init |a1 ⟩
where |e⟩a represents the computational basis states associated with the discretized energy values Ee . of |Final⟩ therefore yields the states n Measurement o
Final
H
|a2 ⟩
H
|a3 ⟩
H
|q⟩ ≡ |χ0 ⟩
2
QF T
U
U
2
U
4
FIG. 8: Illustration of the phase estimation algorithm for quantum wavepacket dynamics. U = exp{−ıH∆t/ℏ}. |q⟩ refers to a family of qubits operated upon by powers of U .
final quantum state as a
q
1
2X −1 2X −1
1
2X −1 2X −1
|Final⟩ = √ a q |e⟩a |xk ⟩ 2 2 e=0 k=0 a 2X −1 1 √ a eı(e∆E)(m∆t)/ℏ χ(xk ; m∆t) × 2 m=0 a
q
|e⟩a |xk ⟩ ξ(xk , Ee ), =√ a q 2 2 e=0 k=0
(13)
where ξ(xk , Ee ) denotes the time-to-energy Fourier transform of χ(xk , tm ). Here, the multidimensional quantum wavepacket is represented on a discrete coordinate grid mapped onto the computational basis of the q-qubit register according to xk −→ |k⟩q ,
(14)
with n n oo |k⟩q ∈ |00 · · · 00⟩q , |00 · · · 01⟩q , · · · , |2q − 1⟩q , (15) where xk labels the discrete coordinate grid points of the wavepacket representation. Appendix B provides a detailed dervation of Eq. (13). More details on the general map between the continuous representation |x⟩ and discrete qubit representation for the nuclear dynamics problem treated here in discussed in Section V B. Here, the index j is reserved exclusively for labeling the physical dimensions of the multidimensional system. Similarly, the ancilla register, after Fourier transform, encodes a discrete energy grid, {|e⟩a ∈ {|00 · · · 00⟩a , |00 · · · 01⟩a , · · · , |2a − 1⟩a }} , (16)
|k⟩q ; |e⟩a with probability |ξ(xk , Ee )| , corresponding to a specific coordinate grid point xk and energy value Ee . The resulting probability distribution forms the Fourier power spectrum associated with the vibrational eigenstates of the system. This standard phase-estimation procedure requires a large number of controlled unitary operations, which remain challenging to implement efficiently on current quantum hardware. In the following, we therefore introduce an alternative approach that substantially reduces the number of controlled operations while naturally integrating with the tensor-network formalism developed above. In this publication, we provide two effective ways to help implement the above algorithm, using tensor networks on current quantum computers. (a) We first split the algorithm depicted in Figure 8 into two parts. The first part does the quantum propagation to obtain multiple time samples of the quantum state and arrives at “Prop” (for propagated) in Figure 8. The second part is the Fourier transform. In our current implementation, we perform the Fourier transform on a classical computer whereas the first part, quantum propagation, is done on a quantum computer. (b) We introduce an alternative algorithm that yields eigen-energy differences, as opposed to absolute energies as these are experimentally observable. This second aspect is a key distinction between the phase estimation algorithm and the approach used here. As noted, the phase estimation algorithm provides an approach to compute the eigenstates of a Hamiltonian from time-propagation followed by Fourier transforms. But, in most chemical applications, including vibrational spectroscopy, energy differences are experimentally observed and have chemical meaning, but absolute energies do not have any physical meaning. For example, in molecular spectroscopy, it is energy differences between eigenvalues that are measured and not absolute eigenenergies. Compare this aspect with Eq. (13) which results in spectral intensities at absolute energy values given by Ee . Hence we begin our approach here towards energy differences to arrive at a reduced quantum resource phase estimation using the Fourier transform of the density-density autocorrelation function and arrive at the following final expression: Z Z +∞ 2 ıωt P(ω) = dx dt e ρ(x, x; t) (17) −∞ 2
Z =
dx
X
δ(ω − (Ei − Ej ))ci (0)c∗j (0)ϕi (x)ϕj (x) ,
i,j
(18) As can be seen from above, this expression yields spectral energy differences directly. Furthermore, using ten-
9 sor networks, we use an approximate form of the expression above given by, P(ω) ≈
η̄ Y N Z X
Z +∞
ıωt
dt e
dxj
2 [j]x ϕαj−1j ,αj
2
(19)
−∞
ᾱ j=1
Let’s define, Pα[j]j−1 ,αj (ω) =
Z +∞
Z dxj
[j]x
dt eıωt ϕαj−1j ,αj
2
2
(20)
−∞
Then, Eq. (19) can be written as:
P(ω) ≈
η̄ Y N X
Pα[j]j−1 ,αj (ω)
(21)
ᾱ j=1
As we will see, this expression provides a ≈ 4cm−1 -level agreement between the quantum computing results and classical computing results.
IV.
EXPERIMENTAL IMPLEMENTATION
Experiments are performed using 171 Yb+ ions confined in a surface electrode trap (Sandia Peregrine [66]). Qubit states are encoded using the |F = 0, mF = 0⟩ ≡ |0⟩ and |F = 1, mF = 0⟩ ≡ |1⟩ hyperfine levels of the 2 S1/2 manifold [67]. Doppler cooling and state initialization into |0⟩ are performed using near-resonant laser beams at 369 nm, and resolved sideband cooling prepares the ions near their ground motional states. State-dependent fluorescence from the ions at 369 nm is captured using a high numerical aperture objective (NA=0.6) and imaged onto a multimode fiber array to provide site-resolved readout of the ion qubit state [68]. State preparation and measurement (SPAM) errors are estimated to be 0.7%. Quantum gate operations are engineered using pairs of Raman beams at 355 nm, with a propagation direction aligned to one of the radial modes of the ion crystal. The Raman beams have two path options, one of which is shaped elliptically and focused to an 8 µm×160 µm spot size at the center of the trap, such that it globally illuminates all ions. The other passes through a multichannel acousto-optic modulator [69] such that each ion may be individually addressed with independent laser amplitudes, frequencies, and phases. The Raman beams can be operated in either a co-propagating configuration, such that both Raman tones pass along the same path, or a counter-propagating configuration that is sensitive to ion motion. Single qubit gates are driven by tuning the beatnote frequency between Raman beams onto resonance with the hyperfine transition freq of 12.642819 GHz. Depending on the chosen phase and laser pulse duration, this executes rotations on the Bloch sphere of form Rx (θ) or Ry (θ). Rotations around the z−axis of the Bloch
sphere are implemented virtually. Rotations around the x− and y−axes are performed in the co-propagating Raman configuration with 99.5(3)% typical fidelity for π/2 rotations. Due to the large number of π/2 gates in the QSD of our unitary propagator, we dynamically correct these gates for small amplitude and frequency errors in the laser pulse [70–73] to improve the overall circuit coherence. Two qubit gates are implemented between pairs of ions by individually addressing them with counterpropagating Raman beams to drive Mølmer-Sørensen interactions [74]. These gates take the Ising-type form XX(θ), which are equivalent to CNOT gates (up to single-qubit rotations) when θ is set equal to π/2 [75]. The amplitude of the laser addressing beam is modulated with a Gaussian envelope to minimize phase-space nonclosure errors and to reduce high-frequency drive components which may off-resonantly couple to nearby transitions [76]. In addition, wrapper gates are added around the Mølmer-Sørensen gates to maintain phase correlation with the co-propagating single qubit rotations [65]. Typical gate times for full entanglement are 200 µs with two-qubit gate fidelities of 98.0(3)%. A key aspect of our experimental implementation is the use of partial-angle two-qubit entangling gates, which can be performed at higher fidelity than fully-entangling gates [65]. In the standard QSD circuit, arbitrary decompositions of three-qubit unitary propagators require 24 fully-entangling CNOT gates or XX(π/2) gates. Here, we extend the decomposition to allow XX gates with arbitrary angle θ and reconstruct the desired unitary at each timestep. The result is a circuit with 6 fullyentangling XX(π/2) gates, and 18 partial-entangling XX(θ) gates, with 0 < θ < π/2. Partial-angle gates are implemented by keeping the gate time constant (to satisfy phase-space closure constraints) and scaling the overall laser amplitude generate the desired entanglement. Careful calibrations have been performed between π/32 ≤ θ ≤ π/2 to minimize effects of gain non-linearity and variable light shifts. To improve the overall fidelity, single- and two-qubit gates with rotation angles θ < π/32 and θ > 31π/32 are approximated as θ = 0 and θ = π/2, respectively. We estimate that our use of partial-entangling gates has led to a > 30% infidelity reduction when implementing QSD circuits on our trapped ion hardware. Also see discussion in Appendix D.
V.
QUANTUM WAVEPACKET DYNAMICS OF PROTONATED WATER CLUSTERS ON ION-TRAP QUANTUM COMPUTERS
Water clusters provide a uniquely controlled environment for understanding how hydrogen bonding, manybody interactions, and quantum effects give rise to the unusual properties of water [29, 32, 77–89]. Additionally, protonated water wires and water clusters are an important class of molecular systems found in many con-
10 strained environments such as biological membranes and enzyme active sites [36–38], ions channels [35, 90], carbon nanotubes [91–93], and fuel cells [39]. Water wires are also present in the photosynthetic reaction center of Rhodobactersphaeroides where they are responsible for proton transfer to a secondary quinone group [37]. Furthermore, the lightweight hydrogen nucleus makes quantum nuclear effects important in such cases [94–97]; additionally the multidimensional quantum nuclear effects in such systems is also known to be critical[27, 98]. Consequently, they have been studied extensively using high-resolution infrared (IR), terahertz vibration–rotation–tunneling (VRT), and rotational spectroscopy across cluster sizes ranging from the water dimer to larger aggregates[27, 28, 98–100]. These experiments reveal highly fluxional structures, dense manifolds of tunneling states, and vibrational spectra that are strongly influenced by cooperative hydrogen bonding[27, 28, 98– 100]. Even small water clusters exhibit signficantly complex quantum effects [89, 101] and vibrational mode couplings[33, 102, 103], and their spectra are extremely sensitive to both the underlying potential energy surface and the treatment of multi-dimensional quantum nuclear effects[28–31, 98, 100, 103–109]. Despite decades of experimental and theoretical effort, significant discrepancies between measured spectra and theoretical predictions have persisted for many protonated as well as hydroxide-rich water clusters[28, 110, 111], leading to long-standing controversies regarding structural assignments and spectral interpretation[99, 112]. A notable example relates to the spectroscopy of the hydrated excess proton, where competing models historically favored either the Eigen-like H9 O+ 4 structure or the Zundel-like H5 O+ 2 motif [94–97, 112, 113]. Subsequent studies revealed that the hydrated proton actually samples a highly fluxional environment in which these structures represent limiting configurations connected by large-amplitude proton transfer and hydrogenbond rearrangements[33, 114–118]. More broadly, the difficulty arises because water clusters are intrinsically multidimensional quantum nuclear systems in which collective intermolecular vibrations strongly influence the spectra[99, 100]. Accurately describing these effects requires quantum nuclear dynamics on high-dimensional potential energy surfaces that also capture electron correlation, many-body polarization, and dispersion[87, 119, 120]. The simultaneous need for high-level electronic structure and multidimensional quantum dynamics has made quantitative agreement between theory and experiment extremely challenging, and water cluster spectroscopy therefore remains an important benchmark for developing predictive methods quantum nuclear dynamics[99, 112, 114]. The system considered in this work is a protonated water wire, H7 O3 + [106]. This system acts a bridge between the key Zundel and Eigen cations and we show here that these kinds of systems can now be studied on current quantum computers, thus taking a critical step towards
FIG. 9: Molecular structure of the H7 O3 + molecule (left) and the corresponding two-dimensional potential energy surface (right) plotted along the shared hydrogen coordinates x1 and x2 . Each coordinate represents the proton shared between the water molecules.
TABLE I: Geometric parameters for the H7 O3 + system. The atom numbering follows that presented in the left panel of Fig. 9. Donor-acceptor distances
Value
Hydrogen bond OO distance (O1 –O2 ) Hydrogen bond OO distance (O1 –O3 )
2.715 Å 2.715 Å
Hydrogen bond angles ◦
Hydrogen bond OHO angle (O1 –H1 –O2 ) Hydrogen bond OHO angle (O1 –H2 –O3 )
175.27 ◦ 177.69
Water molecule bond angles ◦
Water bond angle (H4 –O2 –H5 ) Water bond angle (H6 –O3 –H7 )
107.11 ◦ 106.89
future predictive vibrational spectrocopy calculations on quantum computers, beyond the harmonic approximation.
TABLE II: Characteristics of the grid over which the two-dimensional potential surface is created. Parameter No. of grid points along the x1 -dimension No. of grid points along the x2 -dimension Grid size along the x1 -dimension Grid size along the x2 -dimension Level of theory
Value 8 8 0.70 Å 0.70 Å B3LYP/6-311++G(d,p)
11 TABLE III: Initial wavepacket characteristics. The initial state is defined on a discrete grid using the Kronecker delta δi,j , where δi,j = 1 if i = j and 0 otherwise. Wavepacket Formula χ0 (x1 , x2 ) = δx1 ,0.25 √12
δx2 ,−0.15 + δx2 ,−0.05
TABLE IV: Simulation parameters for the water wire system H7 O3 + . Parameter
Average Energy
1
40.99 kcal/mol
H7 O3
alytic banded Toeplitz distributed approximating functional (DAF) [127–130] representation for the grid representation of the kinetic energy operator:
+
1 fs 150 fs 5 1000
( 2 ) x − x′ −ℏ2 K(x, x ) =K( x − x ) = √ exp − 4mσ 3 2π 2σ 2 n MDAF /2 X x − x′ −1 1 √ H . (22) 4 n! 2n+2 2σ n=0
Molecular geometry and reduced dimensional potential energy surfaces
The analytic banded-Toeplitz representation of the DAF approximation for the kinetic energy operator is one where nthe matrix Kij ≡ K(|i − j|), The quan elements, o
Time step ∆t Total time T Qubit counts D Number of shots
A.
Bond Dimension
′
′
′
To investigate the structural, dynamical, and vibrational properties of the H7 O3 + water-wire system (Fig. 9), we perform quantum wavepacket dynamics simulations on Born–Oppenheimer potential energy surfaces constructed at the density functional theory level using the B3LYP/6-311++G(d,p) method. The optimized geometrical parameters associated with these surfaces are summarized in Table I, while the details of the multidimensional grid used to represent the wavepacket are provided in Table II. The characteristics of the initial wavepackets employed in the simulations are listed in Table III and simulation parameters are given in Table IV. For the H7 O3 + water-wire system, the hydrogenbonded network formed by three water molecules supports two coupled proton-transfer coordinates. The shared proton stretch coordinates are the only degrees of freedom treated quantum mechanically in this study and future studies will work towards scaling up to the full system with 33 nuclear degrees of freedom. The treatment of the proton stretches here gives rise to an effective twodimensional PES, as shown in Fig. 9. The coordinates x1 and x2 describe proton motion along the respective donor–acceptor oxygen axes associated with each hydrogen bond. The interaction between these coordinates introduces coupling between the hydrogen-bonding degrees of freedom, resulting in a multidimensional energy landscape. The PES is discretized on a direct-product grid points along each coordinate, symmetrically distributed about the respective grid centers along the donor–acceptor axes. The kinetic energy operator may be approximated in a number of ways. One approach is to recognize that this operator is diagonal in the momentum representation and hence fast Fourier transforms are commonly employed[121–126]. In this paper, we employ an an-
√ tities H2n+2 x−x in Eq. (22) are the even order 2σ Hermite polynomials that only depend on the spread separating the grid basis vectors, |x⟩ and x′ , and MDAF and σ are parameters that together determine the accuracy and efficiency of the resultant approximate kinetic energy operator. The system is subsequently investigated through quantum wavepacket propagation implemented using quantum computing algorithms. The time evolution of the wavepacket is obtained using a quantum algorithm implemented on the QSCOUT platform at Sandia National Laboratories, following the approach described in Section III. (Mapping the grid representation to the qubit representation is discussed in Section V B.) This framework enables the simulation of multidimensional quantum nuclear dynamics and the extraction of vibrational information directly from the propagated quantum states. All the quantum computing results here are compared with the exact treatment of time-evolution operator represented in the eigenstates basis of the nuclear Hamiltonian, with a detailed error analysis provided in Section V D. All calculations are performed using a Pythonbased classical and quantum algorithm software developed within the Iyengar group.
B. Mapping the continuous grid representation onto the discrete computational basis from qubits
To enable simulation on a gate-based quantum computer, the discretized nuclear-coordinate basis is mapped onto the computational basis of a qubit register. Let D
G = {|gi ⟩}i=1
(23)
12 denote the set of D grid basis points obtained from the discretization of the continuous one-dimensional nuclearcoordinate space for the evolution of each unitary in ′ [j]xj xj Uβj−1 ,βj . Each grid point is assigned a unique integer
TABLE V: First 10 eigenenergies of the H7 O3 + molecule, reported in kcal/mol. Energy (kcal/mol)
1 2 3 4 5 6 7 8 9 10
8.3884 15.3090 15.9444 21.0859 21.4919 23.4887 27.4780 27.7994 28.9209 30.2794
label “i”, which is subsequently encoded in binary form using n = ⌈log2 D⌉ qubits. The resulting correspondence is |gi ⟩ ←→ i ←→ |i⟩q ,
(24)
where |i⟩q = |bn−1 bn−2 · · · b0 ⟩, is a Hamming space vector, and bk ∈ 0, 1 are the binary components of the integer (i). The computational basis states belong to the Hamming space Bn = {|00 · · · 00⟩ , |00 · · · 01⟩ , . . . , |11 · · · 11⟩} .
a
and the discussion here establishes a one-to-one correspondence between the discretized nuclear configurations and the computational basis states of the qubit register. Using this mapping, an arbitrary nuclear wavepacket, [j]x say ϕαj−1j ,αj in Eqs. (3) or (8), that is used in the quantum propagation formalism discussed here, can be represented as ϕ[j] αj−1 ,αj
E
=
D X i=1
ci |gi ⟩ =
D X
ci |i⟩q ,
a
Eigenstate
−1
1kcal/mol = 349.76 cm = 10.49 THz. Energy differences here are generally greater than the resolution imposed by choice of total propagation time of 150fs apart from the difference between eigenstates 2 and 3.
(25)
i=1
and {|gi ⟩} are now the discritized grid D representation E for dimension xj over which the state
xj ϕ[j] αj−1 ,αj
≡
(a)
(b)
[j]x
This allows the initial ϕαj−1j ,αj is represented. wavepacket to be prepared directly on the qubit register and subsequently evolved under the encoded unitary dynamics. The molecular nuclear dynamics can therefore be simulated entirely within the computational basis of the quantum processor. C.
Quantum Wavepacket Dynamics on a Trapped-Ion Quantum Computer
This section presents the implementation of multidimensional quantum wavepacket dynamics on distributed trapped-ion quantum hardware. Each one-dimensional wavepacket was propagated on the quantum computer for a total of 150 fs using a timestep of 1 fs, as summarized in Table IV. These parameters were selected based on the following practical considerations. Firstly, for hydrogen bond dynamics, time-steps of the order of 1fs are commonly used given that oscillation periods are of the order of ≈10fs in such systems. This choice provides good resolution for the higher frequency hydrogen bond oscillations. The second choice of total propagation time of 150fs corresponds to a frequency resolution of approximately 6.67 THz (222.5 cm−1 or 0.64kcal/mol) in the Fourier- spectrum. This resolution window, as we will see, is sufficient to resolve
FIG. 10: (a) Potential energy surface (PES) on 8 × 8 grid, with the initial discrete wavepacket overlaid on a uniform grid. The wavepacket is defined using a Kronecker delta distribution, where filled (black) markers denote grid points with unit amplitude and open (white) markers indicate zero amplitude. The grid spans the interval [−0.35, 0.35] along both x1 and x2 , discretized into 8 points per dimension (∆x = 0.10). The wavepacket is localized at x1 = 0.25 and forms an equal superposition at x2 = −0.15 and x2 = −0.05, consistent with ψ(x1 , x2 , 0) = δx1 ,0.25 √12 (δx2 ,−0.15 + δx2 ,−0.05 ). (b) Projection of the initial wavepacket onto the eigenstates of the system, shown as |⟨ϕn |ψ(0)⟩|2 versus eigenstate index n. The five largest contributions are explicitly labeled.
all energy differnces in the current vibrational problem (see Table V). Other problems may require longer propagation times. The 1 fs sampling interval yields a Nyquist frequency of 1000 THz, which very much exceeds the largest energy gap among all the eigenenergies of the lower-dimensional subsystems. Consequently, the chosen temporal discretization captures the full range of physically relevant vibrational frequencies without aliasing while providing sufficient spectral resolution to identify
13 individual vibrational modes. For the quantum dynamics presented here, the initial wavepacket, χ0 (x1 , x2 ), illustrated in Fig. 10(a) and summarized in Table III, is chosen to be a product state, ψ(x1 , x2 , 0) = δx1 ,0.25 √12 (δx2 ,−0.15 +δx2 ,−0.05 ). The wavepacket is spatially localized in the proton-transfer coordinate x1 and prepared as an equal superposition of two neighboring grid points along x2 . This initialization places the proton predominantly within one channel of the double-well landscape, enabling us to directly monitor population transfer and coherence between the two channels during the subsequent time evolution. The spectral composition of the initial state is shown in Fig. 10(b), which presents the overlap probabilities |⟨ϕn |ψ(0)⟩|2 with the exact vibrational eigenstates of the system. The localized nature of the wavepacket results in a distribution of amplitudes over multiple eigenstates rather than a single vibrational level, with the five largest contributions explicitly indicated in the Fig. 10(b). The resulting non-equilibrium state generates the characteristic oscillatory dynamics observed in the proton-transfer process and provides a suitable initial state for a robust assessment of the accuracy of the distributed quantum simulation framework. In contrast, the potential energy surface exhibits a modest degree of inter-dimensional correlation, characterized by the presence of only two terms in the summation in Eq. (11). Within this framework, each product potential propagator term in Eq. (11) generates two distinct one-dimensional propagators per dimension. For the present two-dimensional proton-transfer system, the resulting dynamics can therefore be distributed across four independent quantum processors and propagated simultaneously. More generally, the achievable level of parallelization is determined jointly by the number of nuclear dimensions and the degree of correlation encoded within the tensor-network representation of the multidimensional potential energy surface. To assess the accuracy of the quantum simulations, benchmark dynamics are first computed using numerically exact propagation, where the time-evolution operator is evaluated in the eigenbasis of the full nuclear Hamiltonian. These classical calculations provide reference trajectories against which the experimentally obtained quantum dynamics are compared. The quantum propagators are compiled into quantum circuits as described in Section IV and executed on the QSCOUT trapped-ion quantum computing platform at Sandia National Laboratories. As shown in Fig. 11, the populations associated with the time-evolved wavepackets obtained from the trappedion quantum hardware (blue markers) are compared directly with the numerically exact classical simulations (solid lines) for each grid basis state of the multidimensional system. The time-dependent dynamics are resolved into the effective one-dimensional subsystems generated through the tensor-network decomposition of the multidimensional propagator as described in Section II C.
In particular, Figs. 11(a) and 11(b) present the time traces associated with the two entanglement components of dimension x1 and similarly, Figs. 11(c) and 11(d) display the dynamics associated with the two entanglement components of dimension x2 . (See Fig. 9.) Although significant deviations in the oscillation amplitudes are observed between the quantum-hardware results and the exact classical propagation, as we will see in the next section, the overall waveform structure and oscillation frequencies remain in strong agreement throughout the evolution. The preservation of the characteristic dynamical frequencies indicates that the tensor-network decomposition together with the trapped-ion implementation captures the essential multidimensional quantum coherence and coupled vibrational motion of the H7 O3 + system. This agreement becomes particularly evident in the Fourier-domain analysis presented later in this work, where the vibrational spectra extracted from the quantum simulations reproduce the dominant spectral features of the exact multidimensional dynamics.
D.
Vibrational spectral results from wavepacket dynamics on quantum computers
The time evolution of the shared-proton wavepacket provides direct access to the vibrational structure of the multidimensional system. In particular, vibrational frequencies are extracted from the Fourier transform of the density–density time-correlation function defined in Eq. (18). The resulting spectra exhibit distinct peaks whose positions correspond to energy differences between the underlying vibrational eigenstates. Within the tensor-network framework employed in this work, the spectral analysis is performed independently for each entanglement copy, that is for each effective one-dimensional subsystem generated from the decomposition of the multidimensional propagator. Specifically, the Fourier transforms are computed using the Eq. (20), thereby allowing the vibrational response associated with each tensor-network leg to be analyzed separately. A detailed derivation of the corresponding spectral expression is provided in Appendix C. The frequency spectra obtained from the timedependent populations shown in Fig. 11 are presented in Fig. 12. In particular, Figs. 12(a) and 12(b) display the Fourier spectra associated with the two entanglement components of dimension x1 , while Figs. 12(c) and 12(d) display the Fourier spectra associated with the two entanglement components of dimension x2 . Thus, panels (a)–(d) in Fig. 12 represent the Fourier transforms of the corresponding time-domain dynamics shown in panels (a)–(d) of Fig. 11. The spectra extracted from the trapped-ion quantum simulations are shown as magenta curves and are obtained directly from the Fourier transform of the experimentally measured density dynamics. For comparison, the solid vertical black lines denote the vibrational tran-
14
(a)Quantum (blue) and classical (yellow) time-traces for |x1 ⟩ (b)Quantum (blue) and classical (yellow) time-traces of |x1 ⟩
(c)Quantum (blue) and classical (yellow) time-traces of |x2 ⟩
(d)Quantum (blue) and classical (yellow) time-traces of |x2 ⟩
FIG. 11: Time-dependent populations of the propagated wavepacket projected onto individual grid basis states, comparing numerically exact classical simulations (solid lines) with trapped-ion quantum hardware results (blue markers). The dynamics are resolved into effective one-dimensional subsystems generated through the tensornetwork decomposition of the multidimensional propagator. Panels (a) and (b) display the dynamics corresponding to the two entanglement components associated with the vibrational coordinate x1 along the proton-transfer direction connecting the two oxygen atoms. Panels (c) and (d) show the dynamics corresponding to the two entanglement components associated with the orthogonal vibrational coordinate x2 , as illustrated in Fig. 9. Although significant deviations in the oscillation amplitudes are observed between the quantum-hardware results and the exact classical propagation, the overall waveform structure remain in strong agreement throughout the evolution.
sition frequencies obtained from exact diagonalization of the corresponding Hamiltonians. The observed agreement between the quantum-computed spectra and the exact vibrational energy differences demonstrates that the distributed trapped-ion implementation accurately captures the essential vibrational structure and coherent proton-transfer dynamics of the multidimensional H7 O3 +
system. Additionally, a quantitative comparison of the vibrational transition frequencies is presented in Fig. 13. Panel (a) shows the transition energies obtained from exact diagonalization of the full two-dimensional Hamiltonian along the horizontal axis, plotted against the transition energies extracted from the trapped-ion simulations
15
(a)
(b)
(c)
(d)
FIG. 12: Frequency spectra obtained from the Fourier transform of the time-dependent wavepacket populations shown in Fig. 11. Panels (a) and (b) correspond to the two entanglement components associated with the protontransfer coordinate x1 connecting the two oxygen nuclei, while panels (c) and (d) correspond to the entanglement components associated with orthogonal vibrational coordinate x2 , shown in Fig. 9. The magenta curves represent spectra extracted from trapped-ion quantum simulations using Eq. (20), and the solid vertical black lines denote vibrational transition frequencies obtained from exact diagonalization. The close agreement between the quantumcomputed spectra and exact results demonstrates that the distributed trapped-ion implementation accurately captures the vibrational structure and multidimensional proton-transfer dynamics of the H7 O3 + system.
of the effective entangled one-dimensional propagations along the vertical axis. The dashed diagonal line represents the guide for perfect agreement, and very clearly from the plots, the agreement is at an extremely high level. Panel (b) provides a magnified view of the lowfrequency region, where deviations between the exact multidimensional spectrum and the reduced-dimensional quantum-computed spectra can be resolved more clearly. The comparison demonstrates that the dominant vibrational transition energies obtained from the trappedion simulations remain in extremely close agreement with the exact multidimensional eigenenergy differences. In particular, the mean absolute deviation between the quantum-computed and exact transition energies, is 4.17 cm−1 , which is within spectroscopic accuracy. These results indicate that the distributed quantum implementation of the tensor-network decomposition preserves the essential spectral structure of the full proton-transfer Hamiltonian.
VI.
DISCUSSION AND OUTLOOK
We introduce a tensor-network based distributed quantum algorithm and implement the algorithm on ion-trap quantum computers. The method is demonstrated by computing vibrational spectra, from wavepacket dynamics, constructed on accurate potentials obtained from electronic structure. The approach also utilizes a new
resource-optimized modified phase estimation algorithm, where we directly find energy differences, rather than absolute energies and use these to compute vibrational frequencies. At the end, the agreement between the quantum computed frequencies and classically computed frequencies are of the order of 4 cm−1 , that is within spectroscopoc accuracy. Thus, we believe, this paper constitutes an important step in the area of quantum computation of chemical dynamics. The present work also establishes a unified connection between tensor network representations, structured quantum circuit decompositions, and distributed quantum computing architectures and implements the same on real quantum computing devices. This perspective provides a pathway to leverage the entanglement structure as a computational resource for hardware parallelism and to extend the reach of quantum simulations beyond the limits imposed by current device capabilities. The present work is closely connected to several recent advances in quantum simulations of molecular dynamics and spectroscopy. In particular, related progress has been made in the quantum simulation of vibronic spectra and nonadiabatic molecular processes [131–134], the study of model systems involving conical intersections and geometric phase effects [135–138], and the development of quantum algorithms for reduced-dimensional reactive scattering and molecular dynamics simulations [139, 140]. The tensor-network-based distributed quantum dynamics framework introduced here is also ex-
16 [1] [2] U1,1 ⊗ U1,1 0 0 0 [1] [2] 0 U1,2 ⊗ U2,1 0 0 [1] [2] 0 0 U1,3 ⊗ U3,1 0 [1] [2] 0 0 0 U1,4 ⊗ U4,1 (a)
|β⟩ (1D)
(a)∆E (Quantum) vs. (2D) ∆E (Classical)
(b)Zoomed low-frequency region
[1]
|χ⟩
[2]
U1,β1 ⊗ Uβ1 ,1
(b)
FIG. 13: Comparison between the vibrational transition frequencies computed using trapped ion quantum hardware and extracted from the effective one-dimensional Hamiltonians and those obtained from exact diagonalization of the full two-dimensional Hamiltonian. Figure (a) Transition energies extracted from the effective onedimensional Hamiltonians simulated on trapped-ion quantum hardware, ∆E (1D) , plotted against the corresponding transition energies of the full two-dimensional Hamiltonian, ∆E (2D) , obtained through exact diagonalization. Figure (b) presents a magnified view of the low-frequency region, highlighting that the tensornetwork-based quantum simulation accurately reproduces the dominant vibrational transition energies of the full H7 O3 + Hamiltonian.
pected to benefit substantially from ongoing developments in quantum electronic-structure algorithms, including fermion-to-qubit mappings, variational quantum eigensolvers, adaptive ansatz constructions, correlationdriven methods, and fragmentation-based approaches for molecular systems [16, 102, 141–167]. These developments, together with continuing progress in experimental quantum hardware implementations [149, 152, 156, 168– 171], provide an important foundation for scalable quantum simulations of multidimensional molecular dynamics and vibrational spectroscopy.
VII.
|β1 ⟩
|χ⟩
[1]
[2]
U1,1 ⊗ U1,1
[1]
[2]
U1,2 ⊗ U2,1
[1]
[2]
U1,3 ⊗ U3,1
[1]
[2]
U1,4 ⊗ U4,1
(c)
FIG. 14: (a) Block diagonal form of the overall unitary in the elevated Hilbert space including entanglement, (b) associated uniformly controlled gate-like representation of Figure (a) with two control qubits encoding β1 , and (c) explicit form of the circuit corresponding to Figure (b). As in Figure 3, the individual blocks [1] [2] U1,β1 ⊗ Uβ1 ,1 are direct product operators that act on product states of dimensions “[1]” and “[2]”.
operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. SAND202622099O
ACKNOWLEDGMENT
This research was supported by the National Science Foundation under Grant CHE-2102610 awarded to SSI, including support provided through a Special Creativity Extension. Author AD acknowledges support from the Siedle Materials Fellowship Foundation and Lynne L. Merritt Fellowship Foundation. These fellowships were awarded to AD by the department of chemistry, Indiana University. Author PR is supported by the Gordon and Betty Moore Foundation, grant DOI 10.37807/GBMF12963. The QSCOUT open-access testbed is funded by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research Quantum Testbed Program. Sandia National Laboratories is a multimission laboratory managed and
Appendix A: Bipartite systems with more entanglement
For N = 2 and and with maximum value of the entanglement index, β1 = 4, the operator consists of four blocks, requiring two control qubits to select the corresponding channel. The associated matrix and circuit representation are shown in Fig. 14, where each computational basis state of the control register encodes a specific value of β1 . The parallel streams that arise from such a situation are given in Figure 15. In both cases discussed above, the number of ancilla needed for a uniformly controlled gate implementation is ln η1 and the number of parallel streams is 2η1 .
17
[1]
[2]
ϕ (x1 , 0)
ϕ (x2 , 0)
β1
[1]
U1,β
[2]
Uβ ,1
1
QC-1 [1]
U1,(β =1) 1
QC-2 [1]
U1,(β =2) 1
1
QC-3 [1]
U1,(β =3) 1
QC-4
QC-5
[1]
[2]
U1,(β =4)
[1]
[1]
[1] Uβ=2 U1,(β 1 =2)
[1]
[1] Uβ=1 U1,(β 1 =3)
U(β =2),1
1
1
[1]
[2]
[2]
[2] U(βU1β=1 =1),1
[2] U(βU1β=2 =2),1
β1
[1]
[2]
U(β =3),1 1
QC-8 [2]
U(β =4),1 1
|q0 ⟩|q1 ⟩|q2 ⟩ |q0 ⟩|q1 ⟩|q2 ⟩ |q0 ⟩|q1 ⟩|q2 ⟩ |q0 ⟩|q1 ⟩|q2 ⟩
[1] Uβ=2 U1,(β 1 =4)
φ1,β1 (x1 , t)
QC-7
[2]
U(β =1),1
1
|q0 ⟩|q1 ⟩|q2 ⟩ |q0 ⟩|q1 ⟩|q2 ⟩ |q0 ⟩|q1 ⟩|q2 ⟩ |q0 ⟩|q1 ⟩|q2 ⟩ [1] Uβ=1 U1,(β 1 =1)
QC-6
[2]
[2] U(β U1β=1 =3),1
[2]
[2] U(β U1β=2 =4),1
[2]
φβ1 ,1 (x2 , t)
FIG. 15: The circuit in Fig. 14 can, in principle, be distributed across eight quantum computers, with each constituent sub-circuit being executed in parallel as explained in Figure 4. But, as noted in Figure 4, for simplicity, in this paper, these operations are performed sequentially on a single quantum device.
Appendix B: A Green’s function description of the Phase estimation algorithm: relation to vibrational spectroscopy from wavepacket dynamics
In quantum phase estimation, illustrated in Figure 8, a specific unitary is used to construct the time-series that is stored as “Prop” on the quantum computer, and can be written as, a
2 −1
h i 1 X |Prop⟩ = √ a |j⟩a U j |χ0 ⟩q . 2 j=0
(B1)
Here the kets, {|j⟩a } represent the states of the ancilla, that is
|q⟩ are mapped onto a multi-dimensional grid basis representation to describe the quantum nuclear wavepacket. More details on the general map between the continuous representation |x⟩ and discrete qubit computational basis representation for the nuclear dynamics problem treated here are presented in discussed in Section V B. In general, for a discrete grid index representation, |k⟩q → xk , with n o |k⟩q | ∀q ∈ {00 · · · 00, 00 · · · 01, · · · , 2q − 1} (B3) and xk being a specific multi-dimensional grid point, Eq. (B1) may be explicitly written as a
a
{|j⟩a | ∀j ∈ {00 · · · 00, 00 · · · 01, · · · , 2 − 1}} (B2) h i and for each state on the ancilla, the state U j |χ0 ⟩q is stored in the qubits thus representing an entangled ancilla-qubit state at the |Prop⟩ stage. The state |χ0 ⟩q is an initial wavepacket state that is represented on the computational basis of q (bottom set of qubits in Figure 8) and here, the computational basis corresponding to
q
1
2X −1 2X −1
1
2X −1 2X −1
h i |Prop⟩ = √ a q |j⟩a |k⟩q ⟨k|q U j |χ0 ⟩ 2 2 j=0 k=0 a
q
h i ≡√ a q |j⟩a |xk ⟩ ⟨xk | U j |χ0 ⟩ (B4) 2 2 j=0 k=0 where in the last equation we have explicitly noted the grid representation and time representation at “Prop”.
18 where
Thus, if U = exp{−ıH∆t/ℏ}, then a
1
q
2X −1 2X −1
i h −ıH[j∆t] |χ0 ⟩ |Prop⟩ = √ a q |j⟩a |xk ⟩ ⟨xk | e ℏ 2 2 j=0 k=0 a
1
a
1 √ a eı(e∆E)(j∆t)/ℏ χ(xk ; j∆t) 2 j=0 Z ≡ dt exp{ıEt/ℏ}χ(x, t) (B9)
ξ(x, E) =
q
2X −1 2X −1
=√ a q |j⟩a |xk ⟩ χ(xk ; j∆t) 2 2 j=0 k=0
2X −1
(B5)
Finally, the Fourier transform simply rotates this vector from the time representation to the conjugate energy representation and hence the operation being performed at the QFT stage is given by,
The measurement of |Final⟩ gives you the state {|k⟩ ; |e⟩}, 2 with probability, |ξ(xk , Ee )| at a specific value of xk (system qubit computational basis value) and Ee (ana a q 2X −1 2X −1 2X −1 cilla computational basis value). In passing, we also note 1 |e⟩a ⟨e|a |j⟩a |xk ⟩ χ(xk ; j∆t) that ξ(x, E) can also be written as |Final⟩ = √ a a q 2 2 2 e=0 j=0 k=0 a
1
q
2X −1 2X −1
|e⟩a |xk ⟩ =√ a q 2 2 e=0 k=0 a 2X −1 1 √ a eı(e∆E)(j∆t)/ℏ χ(xk ; j∆t) 2 j=0
Z ξ(x, E) =
dt exp{ıEt/ℏ} exp{−ıHt/ℏ}χ(x, 0)
=δ(E − H)χ(x, 0)
(B10)
(B6) where δ(E − H) is a Green’s function, also known as the spectral density operator, and is the difference between the advanced and retarded Green’s functions δ(E −H) = G + (E) − G − (E)[172, 173].
where the last line above is of course the Fourier transform, and in rotating the ancilla from time to energy, we have also introduced an energy representation: {|e⟩a | ∀e ∈ {00 · · · 00, 00 · · · 01, · · · , 2a − 1}} ,
Appendix C: Qubit resource optimized phase estimation using tensor networks
(B7)
hence the state of the system prior to the measurement can also be written as a
Using the tensor-network representation introduced above, the density–density correlation spectrum in Eq. (18) can be expressed in terms of the tensor-network wavefunction components as
q
2X −1 2X −1
1 |Final⟩ = √ a q |e⟩a |xk ⟩ ξ(xk , Ee ) 2 2 e=0 k=0
Z +∞
Z P(ω) =
dx̄
(B8)
2
dt eiωt ρ(x̄, x̄; t)
−∞ 2
Z =
dx̄
X X Z +∞ ᾱ
β̄
iωt
dt e
−∞
The expression above contains both diagonal and off-
N Y [j]x ,t [j]x ,t ∗ ϕαj−1j ,αj ϕβj−1j βj .
(C1)
j=1
diagonal tensor-network channels. Separating these contributions yields
19
X Z +∞
Z P(ω) =
dx̄
−∞
ᾱ
+
dt eiωt
N Y
2
N Y [j]x ,t ∗ [j]x ,t . ϕαj−1j ,αj ϕβj−1j βj
iωt
dt e
−∞
P (ω) =
network legs. The corresponding reduced spectral function is therefore written as
η̄ Z +∞ X
Z dx̄
iωt
dt e
−∞
ᾱ
N Y
ϕ[j]xj ,0 .
N Y
2 2 [j]x ,t ϕαj−1j ,αj
.
(C3)
j=1
gator,
As shown in Eq. (C3), restricting the analysis to the diagonal tensor-network channels modifies the relative spectral intensities while preserving the locations of the vibrational transition frequencies. To illustrate the resulting spectral structure, we consider the special case of an initially separable wavepacket corresponding to the product-state limit of Eq. (1), χ0 (x̄) =
(C2)
j=1
In the present work, we specifically evaluate the diagonal contribution associated with the individual tensor-
′
2
j=1
X X Z +∞ ᾱ β̸̄=ᾱ
[j]x ,t
ϕαj−1j ,αj
[j]x ,t
ϕαj−1j ,αj =
Z
′
′
[j]x x
dx′j Uαj−1j ,αj j ϕ[j]xj ,0 ,
(C5)
The initial wavepacket associated with the jth tensornetwork leg is expanded in the eigenbasis of the corresponding effective Hamiltonian,
(C4)
X
ϕ[j]xj ,0 =
j=1
[j]
[j]
clj,α
j−1 ,αj
lj,α
φlj,α
j−1 ,αj
(xj ).
(C6)
j−1 ,αj
Each tensor-network leg then evolves independently under its corresponding effective one-dimensional propa-
Z
[j]x x
The eigenstates satisfy
[j]
′
−iEl
dx′j Uαj−1j ,αj j φlj,α [j]
j−1
(x′j ) = e ,α
[j]x x
[j]
φlj,α
j−1 ,αj
(xj ),
From Eqs. (C5) to (C7), the propagated tensornetwork leg therefore becomes
is a key feature in
[j]
[j]x ,t
X
ϕαj−1j ,αj =
[j]
clj,α
j−1
lj,α
(C7)
this algorithm.)
′
Uαj−1j ,αj j
t/ℏ
j
where, for clarity, we retain only the unitary component of the tensor-network propagator. (As noted in Section II C the non-unitarity of
j,αj−1 ,αj
−iEl
[j]
φlj,α ,α j
j−1
t/ℏ
j,αj−1 ,αj
(xj ) e ,α
.
(C8)
j
j−1 ,αj
The corresponding density associated with the jth tensor-network leg is therefore 2 [j]x ,t ϕαj−1j ,αj =
X lj,α
j−1 ,αj
X mj,α
[j] clj,α
j−1
[j] φlj,α j−1 ,αj
c[j]∗ mj,α ,α j
j−1
(xj ) φ[j]∗ mj,α ,α j
j−1 ,αj
[j] −i El
(xj )e
j,αj−1 ,αj
[j] j,αj−1 ,αj
−Em
t/ℏ
.
j−1 ,αj
(C9)
20 Hence, the time-dependent density associated with an individual tensor-network leg can be expressed as a coherent superposition of oscillatory contributions arising from
[j]
2
[j]x ,t
ϕαj−1j ,αj
energy differences between the eigenstates of the corresponding effective entangled lower-dimensional Hamiltonians,
X
=
lj,α
j−1 ,αj
[j]
where the coefficients Clj ,mj,α
Clj ,mj,α
mj,α
j−1
j
(xj ) e ,α
j,αj−1 ,αj
t/ℏ
,
(C10)
j
j−1 ,αj
the corresponding eigenstate pair,
(xj ) contain the spa-
j−1 ,αj
−i∆El ,m
[j]
X
tial overlap and population amplitudes associated with
[j]
[j]
Clj ,mj,α
j−1 ,αj
(xj ) = clj,α
j−1 ,αj
[j]
[j]
∆Elj ,mj,α
j−1 ,αj
[j]
The quantity ∆Elj ,mj,α
j−1 ,αj
= Elj,α
j−1 ,αj
j−1 ,αj
[j] − Em j,α
[j]x ,t
ϕαj−1j ,αj
=
j=1
where
j−1
j−1 ,αj
,mj,α
j−1 ,αj
N Y
{lj,α
P
{lj,α
X
}
,α ,mj,α j
j−1
,α } j
dt e −∞
N Y iωt
j−1 ,αj
j=1
j=1
{lj,α
,mj,α } j−1 ,αj j−1 ,αj
(xj ),
(C11) (C12)
(xj ) exp−
N it X
ℏ j=1
[j]
∆Elj ,mj,α
j−1 ,αj
,
N Y
Equation (C14) shows that the multidimensional Fourier spectrum contains peaks corresponding both to sums of transition energies across multiple tensornetwork legs and to individual transition energies associated with a particular subsystem. In particular, when
(C13)
evolution contains oscillatory components whose frequencies are determined by sums of transition energies originating from the effective one-dimensional subsystems. The Fourier transform of the multidimensional density then becomes
N X 1 [j] [j] ∆Elj ,mj,α ,α . Clj ,mj,α ,α (xj ) δ ω − j−1 j j−1 j ℏ j=1 j=1
X
j−1 ,αj
.
[j]
Clj ,mj,α
denotes the summation
2 [j]x ,t ϕαj−1j ,αj =
(xj ) φ[j]∗ mj,α
over all eigenstate pairs for every tensor-network leg. Equation (C13) shows that the multidimensional time
Z +∞
j−1 ,αj
network leg and αj−1 , αj entanglement index. For the full tensor-network wavefunction, the multidimensional density is obtained by combining the contributions from all tensor-network legs,
therefore represents the
2
φlj,α
j−1 ,αj
transition energy between the eigenstates lj and mj of the effective Hamiltonian associated with the jth tensor-
N Y
[j]
c[j]∗ mj,α
(C14)
lj,αj−1 ,αj = mj,αj−1 ,αj for all dimensions except a specific leg j1 , the spectral contribution reduces to [j ]
ω=
∆Elj 1 mj ,αj −1 ,αj 1
1
ℏ
1
1
,
21 yielding peaks corresponding directly to the eigenenergy differences of the effective one-dimensional Hamiltonian associated with the j1 th tensor-network leg. Consequently, the multidimensional Fourier spectrum naturally
P ′ (ω) =
Z dx̄
η̄ X
X
N Y
ᾱ {lj,α
j−1
,α ,mj,α j
j−1
,α } j
[j]
Clj ,mj,α
j=1
The expression above explicitly connects the spectral peaks observed in the tensor-network quantum dynamics to the transition energies of the effective one-dimensional Hamiltonians. As demonstrated in the Results section, these transition frequencies closely reproduce the vibrational energy differences obtained from exact diagonalization of the full multidimensional Hamiltonian.
Appendix D: Decomposing the family of
′ [j]x x Uβj−1j ,βjj
into native ion-trap gates
This appendix provides the full decomposition used to implement the operations
[j]x x
′
[j]x
Uβj−1j ,βjj ϕαj−1j ,αj
discussed
in Section II.
M Ry
M Rz A1
B1
M Rz A2
contains both collective multidimensional excitations and the individual vibrational transition frequencies of the effective one-dimensional subsystems. Substituting Eq. (C14) into Eq. (C3) yields
B2
FIG. 16: Three-qubit QSD decomposition grouped into alternating two-qubit and multicontrolled three-qubit subroutines.
[1] Richard P. Feynman, “Simulating physics with computers,” International Journal of Theoretical Physics, International Journal of Theoretical Physics 21, 467–488 (1982). [2] Richard Phillips Feynman, JG Hey, and Robin W Allen, Feynman Lectures on Computation (Addison-Wesley Longman Publishing Co., Inc.75 Arlington Street, Suite 300 Boston, MAUnited States, 1998). [3] R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals (McGraw-Hill Book Company, New York, 1965). [4] H-D Meyer, U Manthe, and L S Cederbaum,
j−1 ,αj
(xj ) δ ω −
N 1X
ℏ j=1
2 [j]
∆Elj ,mj,α
j−1 ,αj
.
(C15)
The QSD approach described in Ref. 174 requires 24 fully-entangling two-qubit gates to implement an arbitrary three-qubit unitary. For quantum hardware platforms that admit partial-angle entangling gates, such as trapped ions, the circuit depth and complexity may be further reduced. We begin by writing the three-qubit QSD in terms of 7 subroutines, as shown in Fig. 16. In Fig. 16, the four subroutines A1,2 and B1,2 operate only on two qubits and may be fully decomposed using the KAK method[175]. This yields a minimal quantum circuit with 15 single-qubit gates and 3 fully entangling two-qubit operations. In our implementation, the full two-qubit entanglement gates are replaced by partialangle XX(θ) gates, with the surrounding single-qubit gates adjusted accordingly. The resulting sequences for these two-qubit unitary subroutines is shown in Fig. 17. The full QSD additionally requires three multicontrolled operations acting on all three qubits. In Fig. 16, subroutines M Rz and M Ry indicate Rz and Ry rotations (respectively), with rotation angles controlled by the four possible states of the top two qubits. In Fig. 18, we provide the full decomposition of these subroutines using partial-angle entangling gates. Compared to a typical multicontrolled gate on three qubits, which requires four fully-entangling operations, our decomposition executes the same unitary using two fully-entangling and two partially-entangling operations. In total, six adjustable angles are used to specify the desired M Rz and M Ry subroutines.
“The multi-configurational time-dependent hartree approach,” Chem. Phys. Lett. 165, 73–78 (1990). [5] U. Manthe, H. D. Meyer, and L. S. Cederbaum, “Wavepacket dynamics within the multiconfiguration hartree framework: General aspects and application to NOCl,” J. Chem. Phys. 97, 3199–3213. [6] M H Beck, A Jäckle, G A Worth, and H-D Meyer, “The multiconfiguration time-dependent hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets,” Phys. Rep. 324, 1–105 (2000). [7] Hans-Dieter Meyer, Fabien Gatti, and Graham A Worth, Multidimensional Quantum Dynamics:
22 Z(θ1 )
Y (θ3 )
Z(θ5 )
Y ( π2 )
X(θ8 ) XX(θ7 )
Z(θ2 )
Y (θ4 )
Z(θ6 )
Z(− π2 )
Z(θ1 )
Y (θ3 )
Z(θ5 )
Z(− π2 )
X(θ8 )
Z( π2 )
Y ( π2 )
X(θ8 )
Z( π2 )
Y ( π2 )
XX(θ7 ) Z(θ2 )
Y (θ4 )
Z(θ6 )
X(θ10 )
Y (− π2 )
Z(− π2 )
XX(θ9 ) X(θ10 )
Y (− π2 )
X(θ8 )
Z( π2 )
X(θ10 )
Y (− π2 )
X(θ10 )
Y (− π2 )
X(θ12 )
X(θ12 )
XX(θ9 ) Y ( π2 )
X(θ12 )
Z(θ13 )
Y (θ15 )
Z(θ17 )
Z(θ14 )
Y (θ16 )
Z(θ18 )
Z(θ13 )
Y (θ15 )
Z(θ17 )
Z(θ14 )
Y (θ16 )
Z(θ18 )
XX(θ11 )
XX(θ11 ) Z(− π2 )
X(θ12 )
Z( π2 )
FIG. 17: Full gate sequence for subroutines A1,2 (top) and B1,2 (bottom) used in the partial-angle QSD decomposition. Each circuit has 18 adjustable angles θi which are used to specify the desired two-qubit unitary. Z(θ1 )
Y ( π2 )
X(θ3 )
Y (− π2 )
X(− π2 )
Y ( π2 )
XX12 (θ2 ) Y ( π2 )
XX13 ( π2 )
X(θ3 )
Y ( π2 )
Y (θ1 )
X(θ3 )
Z( π2 )
X(− π2 )
XX12 (θ2 )
Y ( π2 )
Y (− π2 )
X(θ5 )
Y (− π2 )
X(− π2 )
Z(θ6 ) XX13 ( π2 )
Y (− π2 )
X(π)
Z(− π2 ) Y ( π2 )
X(θ5 ) XX12 (θ4 )
Z(− π2 )
X(θ5 )
Z( π2 )
X(θ5 )
Y (− π2 )
X(− π2 )
Y (θ6 )
XX12 (θ4 ) X(θ3 )
XX13 ( π2 ) X(π)
XX13 ( π2 ) Y (− π2 )
FIG. 18: Gate decomposition for multicontrolled three-qubit subroutines M Rz (top) and M Ry (bottom). Each circuit contains 6 adjustable angles θi used to specify the target multicontrolled unitary.
MCTDH Theory and Applications (John Wiley & Sons, 2009). [8] Haobin Wang and Michael Thoss, “Multilayer formulation of the multiconfiguration time-dependent hartree theory,” J. Chem. Phys. 119, 1289–1299 (2003). [9] Robert Binder and Irene Burghardt, “First-principles quantum simulations of exciton diffusion on a minimal oligothiophene chain at finite temperature,” Faraday Discuss. 221, 406–427 (2020). [10] Frank Otto, “Multi-layer potfit: An accurate potential representation for efficient high-dimensional quantum dynamics,” J. Chem. Phys. 140, 014106 (2014). [11] Nicole DeGregorio and Srinivasan S Iyengar, “Adaptive dimensional decoupling for compression of quantum nuclear wave functions and efficient potential energy surface representations through tensor network decomposition,” J. Chem. Theory Comput. 15, 2780–2796 (2019). [12] Samuel M Greene and Victor S Batista, “Tensor-train split-operator fourier transform (tt-soft) method: Multidimensional nonadiabatic quantum dynamics,” Journal of chemical theory and computation 13, 4034–4042 (2017). [13] Michael A Nielsen and Isaac L Chuang, Quantum computation and quantum information (Cambridge University Press, Cambridge, New York, NY, 2000). [14] John Preskill, “Quantum computing and the entanglement frontier,” arXiv:1203.5813 [quant-ph] (2012). [15] John Preskill, “Quantum computing 40 years later,” (2021), arXiv:2106.10522 [quant-ph].
[16] Juncheng (Harry) Zhang and Srinivasan S. Iyengar, “Graph-|Q⟩ ⟨C|: A graph-based quantum-classical algorithm for efficient electronic structure on hybrid quantum/classical hardware systems: Improved quantum circuit depth performance,” J. Chem. Theory Comput. 18, 2885 (2022). [17] Srinivasan S. Iyengar, Juncheng Harry Zhang, Debadrita Saha, and Timothy C. Ricard, “Graph-|Q⟩ ⟨C|: A quantum algorithm with reduced quantum circuit depth for electronic structure,” J. Phys. Chem. A 127, 9334 (2023). [18] Muhammad AbuGhanem and Hichem Eleuch, “Twoqubit entangling gates for superconducting quantum computers,” Results in Physics 56, 107236 (2024). [19] Yuxuan Yan, Zhenyu Du, Junjie Chen, and Xiongfeng Ma, “Limitations of noisy quantum devices in computing and entangling power,” npj Quantum Information 11, 188 (2025). [20] Pei Yuan, Jonathan Allcock, and Shengyu Zhang, “Does qubit connectivity impact quantum circuit complexity?” IEEE Transactions on computer-aided design of integrated circuits and systems 43, 520–533 (2023). [21] Anurag Dwivedi, A. J. Rasmussen, Philip Richerme, and Srinivasan S. Iyengar, “Quantum nuclear dynamics on a distributed set of ion-trap quantum computing systems,” J. Am. Chem. Soc. 146, 29355–29363 (2024). [22] Anurag Dwivedi, Miguel Angel Lopez-Ruiz, and Srinivasan S. Iyengar, “Resource optimization for quantum dynamics with tensor networks: Quantum and
23 classical algorithms,” The Journal of Physical Chemistry A 128, 6774–6797 (2024), pMID: 39101545, https://doi.org/10.1021/acs.jpca.4c03407. [23] Ville Bergholm, Juha J Vartiainen, Mikko Möttönen, and Martti M Salomaa, “Quantum circuits with uniformly controlled one-qubit gates,” Physical Review A—Atomic, Molecular, and Optical Physics 71, 052330 (2005). [24] José Alex de Carvalho, Carlos Batista, Tiago de Veras, Israel Araujo, and Adenilton José da Silva, “Quantum multiplexer simplification for state preparation,” ACM Transactions on Quantum Computing 6, 1–12 (2025). [25] John Preskill, “Quantum Computing in the NISQ era and beyond,” Quantum 2, 79 (2018). [26] Nai-Hui Chia, Kai-Min Chung, and Ching-Yi Lai, “On the need for large quantum depth,” in Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, STOC 2020 (Association for Computing Machinery, New York, NY, USA, 2020) p. 902–915. [27] J.-W. Shin, N. I. Hammer, E. G. Diken, M. A. Johnson, R. S. Walters, T. D. Jaeger, M. A. Duncan, R. A. Christie, and K. D. Jordan, “Infrared signature of struc+ tures associated with the H (H2 O)n (n = 6 to 27) clusters,” Science 304, 1137 (2004). [28] Conrad T. Wolke, Joseph A. Fournier, Laura C. Dzugan, Matias R. Fagiani, Tuguldur T. Odbadrakh, Harald Knorke, Kenneth D. Jordan, Anne B. McCoy, Knut R. Asmis, and Mark A. Johnson, “Spectroscopic snapshots of the proton-transfer mechanism in water,” Science 354, 1131 (2016). [29] Oriol Vendrell, Fabien Gatti, and Hans-Dieter Meyer, “Dynamics and infrared spectroscopy of the protonated water dimer,” Ang. Chem. Intl. Ed. 46, 6918 (2007). [30] O. Vendrell, F. Gatti, and H.-D. Meyer, “Full dimensional (15d) quantum-dynamical simulation of the protonated water-dimer ii: Infrared spectrum and vibrational dynamics dynamics and infrared spectroscopy of the protonated water dimer,” J. Chem. Phys. 127, 184303 (2007). [31] O. Vendrell, F. Gatti, and H.-D. Meyer, “Strong isotope effects in the infrared spectrum of the zundel cation,” Angew. Chem. Intl. Ed. 48, 352 (2009). [32] S. M. Dietrick and S. S. Iyengar, “Constructing periodic phase space orbits from ab initio molecular dynamics trajectories to analyze vibrational spectra: Case study + of the zundel (H5 O2 ) cation,” J. Chem. Theory and Comput. 8, 4876 (2012). [33] S. S. Iyengar, M. K. Petersen, T. J. F. Day, C. J. Burnham, V. E. Teige, and G. A. Voth, “The properties of ion-water clusters. i. the protonated 21-water cluster,” J. Chem. Phys. 123, 084309 (2005). [34] S. S. Iyengar, “Further analysis of the dynamically averaged vibrational spectrum for the “magic” protonated 21-water cluster,” J. Chem. Phys. 126, 216101 (2007). [35] Toby W Allen, Turgut Baştuğ, Serdar Kuyucak, and Shin-Ho Chung, “Gramicidin A channel as a test ground for molecular dynamics force fields,” Biophys. J. 84, 2159–2168 (2003). [36] J. F. Nagle and H. J. Morowitz, “Molecular mechanisms for proton transport in membranes,” Proc.Natl.Acad.Sci. 75, 298– (1978). [37] Laura Baciou and Hartmut Michel, “Interruption of the water chain in the reaction center from Rhodobacter sphaeroides reduces the rates of the proton uptake and
of the second electron transfer to QB,” Biochemistry 34, 7967–7972 (1995). [38] Haibo Guo and Amanda S Barnard, “Proton transfer in the hydrogen-bonded chains of lepidocrocite: a computational study,” Phys. Chem. Chem. Phys. 13, 17864 (2011). [39] Yun-Sheng Ye, John Rick, and Bing-Joe Hwang, “Water soluble polymers as proton exchange membranes for fuel cells,” Polymers 4, 913–963 (2012). [40] S Solomon, R. R. Garcia, F. S. Rowland, and D. J. Wuebbles, “On the depletion of antarctic ozone,” Nature 321, 755 (1986). [41] B. J. Finlayson-Pitts and J. N. Pitts, Jr., “Chemistry of the upper and lower atmosphere: Theory, experiments, and applications,” (Academic, San Diego, 2000). [42] R. P. Turco, O. B. Toon, and P. Hamil, “Heterogeneous physicochemistry of the polar ozone hole,” Journal Of Geophysical Research 94, 16493 (1989). [43] S. Solomon, “Stratospheric ozone depletion: A review of concepts and history,” Reviews of Geophysics 37, 275 (1999). [44] Dwayne E. Heard and Michael J. Pilling, “Measurement of oh and ho2 in the troposphere,” Chem. Revs. 103, 5163 (2003). [45] J. Troe, “Detailed modeling of the temperature and pressure dependence of the reaction h+o2 (+ m)→ ho2 (+ m),” Proceedings of the combustion institute 28, 1463 (2000). [46] J. A. Logan, M. J. Prather, S. C. Wofsy, and M. B. McElroy, “Tropospheric chemistry: A global perspective,” J. Geophys. Res. 86, 7210 (1981). [47] M. J. McEwan and L. F. Phillips, “Chemistry in the upper atmosphere,” Accounts of Chemical Research 3, 9 (1970). [48] A. W. Castleman and R. G. Keesee, “Ionic clusters,” Chem. Rev. 86, 589 (1986). [49] M. J. McEwan and L. F. Phillips, Chemistry of the Atmosphere (Eward Arnold:London, 1975). [50] R. P. Wayne, Chemistry of the Atmosphere (Clarendon Press: Oxford, 1994). [51] J. S. Aloisio, S.; Francisco, “Radical water complexes in earth’s atmosphere,” Accounts of Chemical Research 33, 825 (2000). [52] D. F. Harris, D. A. Lukoyanov, S. Shaw, P. Compton, M. Tokmina-Lukaszewska, B. Bothner, N. Kelleher, D. R. Dean, B. M. Hoffman, and L. C. Seefeldt, “The Mechanism of N2 Reduction Catalyzed by FeNitrogenase Involves Reductive Elimination of H2 ,” Biochemistry 57, 701–710 (2018). [53] R. N. F. Thorneley and D. J. Lowe, “Molybdenum enzymes,” (Wiley-Interscience: New York, 1985) Chap. Kinetics and Mechanism of the Nitrogenase Enzymatic System, pp. 221–284. [54] D. V. Yandulov and R. R. Schrock, “Reduction of dinitrogen to ammonia at a well-protected reaction site in a molybdenum triamidoamine complex,” J. Am. Chem. Soc. 124, 6252 (2002). [55] N. Agmon, “The Grotthuss Mechanism,” Chem. Phys. Lett. 244, 456 (1995). [56] Román Orús, “A practical introduction to tensor networks: Matrix product states and projected entangled pair states,” Ann. Physics 349, 117 – 158 (2014). [57] Anup Kumar, Nicole DeGregorio, Timothy Ricard, and Srinivasan S. Iyengar, “Graph-theoretic molecular frag-
24 mentation for potential surfaces leads naturally to a tensor network form and allows accurate and efficient quantum nuclear dynamics,” J. Chem. Theory Comput. 18, 7243 (2022). [58] Michael L Wall and Lincoln D Carr, “Out-ofequilibrium dynamics with matrix product states,” New Journal of Physics 14, 125015 (2012). [59] B Pirvu, V Murg, J I Cirac, and F Verstraete, “Matrix product operator representations,” New Journal of Physics 12, 025012 (2010). [60] M B Hastings, “An area law for one-dimensional quantum systems,” Journal of Statistical Mechanics: Theory and Experiment 2007, P08024 (2007). [61] M. F. Trotter, “On the product of semi-groups of operators,” Proc. Am. Math. Soc. 10, 545 (1959). [62] E. Nelson, “Feynman integrals and the schrödinger equation,” J. Math. Phys. 5, 332 (1964). [63] Xiao Zhu and Srinivasan S. Iyengar, “Large language model-type architecture for high-dimensional molecular potential energy surfaces,” Phys. Rev. X 16, 011012 (2026). [64] Debadrita Saha, Philip Richerme, and Srinivasan S Iyengar, “Quantum circuit and mapping algorithms for wavepacket dynamics: case study of anharmonic hydrogen bonds in protonated and hydroxide water clusters,” Journal of Chemical Theory and Computation 21, 3814–3831 (2025). [65] Christopher G Yale, Ashlyn D Burch, Matthew NH Chow, Brandon P Ruzic, Daniel S Lobser, Brian K McFarland, Melissa C Revelle, and Susan M Clark, “Realization and calibration of continuously parameterized two-qubit gates on a trapped-ion quantum processor,” IEEE Transactions on Quantum Engineering (2025). [66] Melissa C. Revelle, “Phoenix and peregrine ion traps,” arXiv preprint arXiv:2009.02398 (2020). [67] S. Olmschenk, K. C. Younge, D. L. Moehring, D. N. Matsukevich, P. Maunz, and C. Monroe, “Manipulation and detection of a trapped yb[sup +] hyperfine qubit,” Phys. Rev. A 76, 052314 (2007). [68] Susan M. Clark, Daniel Lobser, Melissa C. Revelle, Christopher G. Yale, David Bossert, Ashlyn D. Burch, Matthew N. Chow, Craig W. Hogle, Megan Ivory, Jessica Pehr, Bradley Salzbrenner, Daniel Stick, William Sweatt, Joshua M. Wilson, Edward Winrow, and Peter Maunz, “Engineering the quantum scientific computing open user testbed,” IEEE Transactions on Quantum Engineering 2, 1–32 (2021). [69] L3 Harris, “Acousto optic solutions,” https://www.l3harris.com/all-capabilities/acoustooptic-solutions (2020). [70] Evangelos Piliouras, Hisham Amer, Susan M. Clark, Melissa C. Revelle, Edward C. Tortorici, Matthew N. H. Chow, Brandon Ruzic, Daniel S. Lobser, Brian K. McFarland, Christopher G. Yale, Edwin Barnes, and Sophia E. Economou, “High-performance gates on trapped ion qubits using counterpropagating pulseshaped laser beams,” (2026), manuscript in preparation. [71] Evangelos Piliouras, Dennis Lucarelli, and Edwin Barnes, “An automated geometric space curve approach for designing dynamically corrected gates,” npj Quantum Information 12, 46 (2026). [72] Hisham Amer, Evangelos Piliouras, Edwin Barnes, and Sophia E. Economou, “Implementing and benchmark-
ing dynamically corrected gates on superconducting devices using space curve quantum control,” arXiv preprint arXiv:2504.09767 (2025). [73] Edwin Barnes, Fernando A. Calderon-Vargas, Wenzheng Dong, Bikun Li, Junkai Zeng, and Fei Zhuang, “Dynamically corrected gates from geometric space curves,” Quantum Science and Technology 7, 023001 (2022). [74] Klaus Mølmer and Anders Sørensen, “Multiparticle entanglement of hot trapped ions,” Phys. Rev. Lett. 82, 1835 (1999). [75] Dmitri Maslov, “Basic circuit compilation techniques for an ion-trap quantum machine,” New Journal of Physics 19, 023035 (2017). [76] Brandon P Ruzic, Matthew NH Chow, Ashlyn D Burch, Daniel S Lobser, Melissa C Revelle, Joshua M Wilson, Christopher G Yale, and Susan M Clark, “Leveraging motional-mode balancing and simply parametrized waveforms to perform frequency-robust entangling gates,” Physical Review Applied 22, 014007 (2024). [77] R. Pomes and B. Roux, “Structure and dynamics of + a proton wire: A theoretical study of H translocation along the single-file water chain in the gramicidin a channel,” Biophys J. 71, 19 (1996). [78] M. J. McEwan and L. F. Phillips, Chemistry of the Atmosphere (Eward Arnold: London, 1975). [79] R. P. Wayne, Chemistry of the Atmosphere (Clarendon Press: Oxford, 1994). [80] Regis Pomes and Benoı̂t Roux, “Theoretical study of H+ translocation along a model proton wire,” J. Phys. Chem. 100, 2519–2527 (1996). [81] H. Decornez, K. Drukker, and S. Hammes-Schiffer, “Solvation and hydrogen-bonding effects on proton wires,” J. Phys. Chem. A 103, 2891 (1999). [82] M. L. Brewer, U. W. Schmitt, and G. A. Voth, “The formation and dynamics of proton wires in channel environments,” Biophys J. 80, 1691 (2001). [83] MM Teeter, “Water structure of a hydrophobic protein at atomic resolution: Pentagon rings of water molecules in crystals of crambin,” Proceedings of the National Academy of Sciences 81, 6014–6018 (1984). [84] Stephen Neidle, Helen M. Berman, and H. S. Shieh, “Highly structured water network in crystals of a deoxydinucleoside-drug complex,” Nature 288, 129 (1980). [85] Leigh Ann Lipscomb, Mary Elizabeth Peek, Fang Xiao Zhou, Jay A Bertrand, Donald VanDerveer, and Loren Dean Williams, “Water ring structure at DNA interfaces: hydration and dynamics of DNA-anthracycline complexes,” Biochemistry 33, 3649–3659 (1994). [86] Chingkuang Tu, Roger S Rowlett, Brian C Tripp, James G Ferry, and David N Silverman, “Chemical rescue of proton transfer in catalysis by carbonic anhydrases in the β-and γ-class,” Biochemistry 41, 15429– 15435 (2002). [87] Andreas Heßelmann, “Correlation effects and manybody interactions in water clusters,” Beilstein journal of organic chemistry 14, 979–991 (2018). [88] Vasiliy S Znamenskiy and Michael E Green, “Quantum calculations on hydrogen bonds in certain water clusters show cooperative effects,” Journal of chemical theory and computation 3, 103–114 (2007). [89] Frank H Stillinger, “Water revisited,” science 209, 451–
25 457 (1980). [90] Alberte Pullman and Catherine Etchebest, “The gramicidin A channel: the energy profile for single and double occupancy in a head-to-head β6. 33, 3-helical dimer backbone,” FEBS Lett. 163, 199–202 (1983). [91] David J. Mann and Mathew D. Halls, “Water alignment and proton conduction inside carbon nanotubes,” Phys. Rev. Lett. 90, 195503 (2003). [92] Woochul Song, Himanshu Joshi, Ratul Chowdhury, Joseph S. Najem, Yue-Xiao Shen, Chao Lang, Codey B. Henderson, Yu-Ming Tu, Megan Farell, Megan E. Pitz, and et al., “Artificial water channels enable fast and selective water permeation through water-wire networks,” Nat. Nanotechnol. 15, 73–79 (2020). [93] G. Hummer, J. C. Rasaiah, and J. P. Noworyta, “Water conduction through the hydrophobic channel of a carbon nanotube,” Nature 414, 188–190 (2001). [94] M. E. Tuckerman, P. J. Ungar, T. Vonrosenvinge, and M. L. Klein, “Ab initio molecular dynamics simulations,” J. Phys. Chem. 100, 12878 (1996). [95] M. Tuckerman, K. Laasonen, M. Sprik, and M. Parrinello, J. Phys. Chem. 99, 5749 (1995). [96] Udo W Schmitt and Gregory A Voth, “The computer simulation of proton transport in water,” J. Chem. Phys. 111, 9361–9381 (1999). [97] Udo W Schmitt and Gregory A Voth, “Multistate empirical valence bond model for proton transport in water,” J. Phys. Chem. B 102, 5547–5551 (1998). [98] J. M. Headrick, E. G. Diken, R. S. Walters, N. I. Hammer, R. A. Christie, J. Cui, E. M. Myshakin, M. A. Duncan, M. A. Johnson, and KD Jordan, “Spectral signatures of hydrated proton vibrations in water clusters,” Science 308, 1765 (2005). [99] Frank N Keutsch and Richard J Saykally, “Water clusters: Untangling the mysteries of the liquid, one molecule at a time,” Proceedings of the National Academy of Sciences 98, 10533–10540 (2001). [100] Kn Liu, JD Cruzan, and RJ Saykally, “Water clusters,” Science 271, 929–933 (1996). [101] Andrey M Tokmachev, “Networks of hydrogen bond networks in water clusters,” The Journal of Physical Chemistry A 128, 2763–2771 (2024). [102] Srinivasan S. Iyengar, Debadrita Saha, Anurag Dwivedi, Miguel A. Lopez-Ruiz, Anup Kumar, Juncheng (Harry) Zhang, Timothy C. Ricard, Philip Richerme, and Amr Sabry, “Quantum algorithms for the study of electronic structure and molecular dynamics: Novel computational protocols,” in Comprehensive Computational Chemistry (Elsevier, 2023). [103] Timothy C Ricard, Anup Kumar, and Srinivasan S Iyengar, “Embedded, graph-theoretically defined manybody approximations for wavefunction-in-dft and dft-indft: Applications to gas- and condensed-phase ab initio molecular dynamics, and potential surfaces for quantum nuclear effects,” Int. J. Quantum Chem. 120, e26244 (2020). [104] Anup Kumar, Nicole DeGregorio, and Srinivasan S. Iyengar, “Graph-theory-based molecular fragmentation for efficient and accurate potential surface calculations in multiple dimensions,” J. Chem. Theory Comput. 17, 6671–6690 (2021). [105] Xiao Zhu and Srinivasan S. Iyengar, “Graph theoretic molecular fragmentation for multidimensional potential energy surfaces yield an adaptive and general transfer
machine learning protocol,” J. Chem. Theory Comput. 18, 5125–5144 (2022). [106] Joseph A. Fournier, Conrad T. Wolke, Mark A. Johnson, Tuguldur T. Odbadrakh, Kenneth D. Jordan, Shawn M. Kathmann, and Sotiris S. Xantheas, “Snapshots of proton accommodation at a microscopic water surface: Understanding the vibrational spectral signatures of the charge defect in cryogenically cooled + H(H2 O)n =2–28 clusters,” J. Phys. Chem. A 119, 37 (2015). [107] Eric G. Diken, Jeffrey M. Headrick, Joseph R. Roscioli, Joseph C. Bopp, Mark A. Johnson, and A. B. McCoy, “Fundamental excitations of the shared proton in the − + H3 O2 and H5 O2 complexes,” J. Phys. Chem. A 109, 1487 (2005). [108] Knut R. Asmis, Nicholas L. Pivonka, Gabriele Santambrogio, Mathias Brümmer, Cristina Kaposta, Daniel M. Neumark, and Ludger Wöste, “Gas-phase infrared spectrum of the protonated water dimer,” Science 299, 1375 (2003). [109] Joseph A. Fournier, Conrad T. Wolke, Christopher J. Johnson, Mark A. Johnson, Nadja Heine, Sandy Gewinner, Wieland Schöllkopf, Tim K. Esser, Matias R. Fagiani, Harald Knorke, and Knut R. Asmis, “Site-specific vibrational spectral signatures of water molecules in the + + magic h3 o (h2 o)20 and cs (h2 o)20 clusters,” Proceedings of the National Academy of Sciences 111, 18132– 18137 (2014). [110] D. Asthagiri, L. R. Pratt, J. D. Kress, and M. A. − Gomez, “Hydration and mobility of ho (aq),” Proc. Natl. Acad. Sci. 101, 7229 (2004). [111] X. Li, V. E. Teige, and S. S. Iyengar, “Can the fourcoordinated, penta-valent oxygen in hydroxide water clusters be detected through experimental vibrational spectroscopy?” J. Phys. Chem. A 111, 4815 (2007). [112] Paul B Calio, Chenghan Li, and Gregory A Voth, “Resolving the structural debate for the hydrated excess proton in water,” Journal of the American Chemical Society 143, 18672–18683 (2021). [113] G. Zundel, Hydration and Intermolecular Interaction: Infrared Investigations with Polyelectrolyte Membranes (Academic Press, 2012). [114] Dominik Marx, “Proton transfer 200 years after von grotthuss: Insights from ab initio simulations,” ChemPhysChem 7, 1848–1870 (2006). [115] D. Marx, M. E. Tuckerman, and M. Parrinello, “Solvated excess protons in water: Quantum effects on the hydration structure,” J. Phys,: Condens. Matter 12, A153 (2000). [116] M. E. Tuckerman, D. Marx, and M. Parrinello, “The nature and transport mechanism of hydrated hydroxide ions in aqueous solution,” Nature 417, 925 (2002). [117] M. E. Tuckerman, D. Marx, and M. Parrinello, “The nature and transport mechanism of hydrated hydroxide ions in aqueous solution,” Nature 417, 925– (2002). [118] D. Marx, M. E. Tuckerman, J. Hutter, and M. Parrinello, “The nature of the hydrated excess proton in water,” Nature 397, 601– (1999). [119] Nathan I Hammer, Eric G Diken, Joseph R Roscioli, Mark A Johnson, Evgeniy M Myshakin, Kenneth D Jordan, Anne B McCoy, Xinchuan Huang, Joel M Bowman, and Stuart Carter, “The vibrational predissoci+ ation spectra of the H5 O2 .RGn (RG= Ar, Ne) clus-
26 ters: Correlation of the solvent perturbations in the free OH and shared proton transitions of the Zundel ion,” J. Chem. Phys. 122, 244301 (2005). [120] David Eisenberg and Walter Kauzmann, The structure and properties of water (Oxford University Press, 2005). [121] M. D. Feit and J. A. Fleck, “Solution of the schrödinger equation by a spectral method ii: Vibrational energy levels of triatomic molecules,” J. Chem. Phys. 78, 301 (1982). [122] M. D. Feit and J. A. Fleck, J. Chem. Phys. 79, 301 (1983). [123] M. D. Feit and J. A. Fleck, J. Chem. Phys. 80, 2578 (1984). [124] P. DeVries, “Fast fourier transform techniques in numerical simulations of intense pulse-molecule interactions,” in Atomic and Molecular Processes with Short Intense Laser Pulses, NATO ASI Series B; Physics, Vol. 171, edited by A. D. Bandrauk (Plenum Press, New York, 1988) p. 481. [125] D Kosloff and R Kosloff, “A fourier method solution for the time dependent schrödinger equation as a tool in molecular dynamics,” Journal of Computational Physics 52, 35–53 (1983). [126] D. Kosloff and R. Kosloff, J. Chem. Phys. 79, 1823 (1983). [127] Donald J. Kouri, Youhong Huang, and David K. Hoffman, “Iterated real-time path integral evaluation using a distributed approximating functional propagator and average-case complexity integration,” Phys. Rev. Lett. 75, 49 (1995). [128] D. K. Hoffman, N. Nayar, O. A. Sharafeddin, and D. J. Kouri, “Analytic banded approximation for the discretized free propagator,” J. Phys. Chem. 95, 8299 (1991). [129] Debadrita Saha, Srinivasan S Iyengar, Philip Richerme, Jeremy M Smith, and Amr Sabry, “Mapping quantum chemical dynamics problems to spin-lattice simulators,” J. Chem. Theory Comput. 17, 6713–6732 (2021). [130] S. S. Iyengar, “Ab Initio dynamics with wave-packets and density matrices,” Theo. Chem. Accts. 116, 326 (2006). [131] Chris Sparrow, Enrique Martı́n-López, Nicola Maraviglia, Alex Neville, Christopher Harrold, Jacques Carolan, Yogesh N Joglekar, Toshikazu Hashimoto, Nobuyuki Matsuda, Jeremy L O’Brien, David P Tew, and Anthony Laing, “Simulating the vibrational quantum dynamics of molecules using photonics,” Nature 557, 660–667 (2018). [132] Christopher S. Wang, Jacob C. Curtis, Brian J. Lester, Yaxing Zhang, Yvonne Y. Gao, Jessica Freeze, Victor S. Batista, Patrick H. Vaccaro, Isaac L. Chuang, Luigi Frunzio, Liang Jiang, S. M. Girvin, and Robert J. Schoelkopf, “Efficient multiphoton sampling of molecular vibronic spectra on a superconducting bosonic processor,” Phys. Rev. X 10, 021060 (2020). [133] Joonsuk Huh, Gian Giacomo Guerreschi, Borja Peropadre, Jarrod R McClean, and Alán Aspuru-Guzik, “Boson sampling for molecular vibronic spectra,” Nature Photonics 9, 615–620 (2015). [134] Yuchen Wang, LeeAnn M Sager-Smith, and David A Mazziotti, “Quantum simulation of bosons with the contracted quantum eigensolver,” New Journal of Physics 25, 103005 (2023). [135] Christopher S. Wang, Nicholas E. Frattini, Benjamin J.
Chapman, Shruti Puri, S. M. Girvin, Michel H. Devoret, and Robert J. Schoelkopf, “Observation of wave-packet branching through an engineered conical intersection,” Phys. Rev. X 13, 011008 (2023). [136] Christophe H Valahu, Vanessa C Olaya-Agudelo, Ryan J MacDonell, Tomas Navickas, Arjun D Rao, Maverick J Millican, Juan B Pérez-Sánchez, Joel YuenZhou, Michael J Biercuk, Cornelius Hempel, Ting R Tan, and Kassal Ivan, “Direct observation of geometricphase interference in dynamics around a conical intersection,” Nature Chemistry 15, 1503–1508 (2023). [137] Jacob Whitlow, Zhubing Jia, Ye Wang, Chao Fang, Jungsang Kim, and Kenneth R Brown, “Quantum simulation of conical intersections using trapped ions,” Nature Chemistry 15, 1509–1514 (2023). [138] Yuchen Wang and David A Mazziotti, “Quantum simulation of conical intersections,” Physical Chemistry Chemical Physics 26, 11491–11497 (2024). [139] Xiaodong Xing, Alejandro Gomez Cadavid, Artur F Izmaylov, and Timur V Tscherbul, “A hybrid quantumclassical algorithm for multichannel quantum scattering of atoms and molecules,” The Journal of Physical Chemistry Letters 14, 6224–6233 (2023). [140] Sumit Suresh Kale and Sabre Kais, “Simulation of chemical reactions on a quantum computer,” The Journal of Physical Chemistry Letters 15, 5633–5642 (2024). [141] P. Jordan and E. Wigner, “Über das paulische Äquivalenzverbot,” Z. Phys. 47, 631–651 (1928). [142] G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, “Quantum algorithms for fermionic simulations,” Phys. Rev. A 64, 022319 (2001). [143] Sergey B. Bravyi and Alexei Yu. Kitaev, “Fermionic quantum computation,” Ann. Physics 298, 210–226 (2002). [144] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, “Simulated quantum computation of molecular energies,” Science 309, 1704 (2005). [145] William M. Kirby and Peter J. Love, “Variational quantum eigensolvers for sparse hamiltonians,” Phys. Rev. Lett. 127, 110503 (2021). [146] Alba Cervera-Lierta, Jakob S. Kottmann, and Alán Aspuru-Guzik, “Meta-variational quantum eigensolver: Learning energy profiles of parameterized hamiltonians for quantum simulation,” PRX Quantum 2, 020329 (2021). [147] P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, “Scalable quantum simulation of molecular energies,” Phys. Rev. X 6, 031007 (2016). [148] Scott E. Smart and David A. Mazziotti, “Quantum solver of contracted eigenvalue equations for scalable molecular simulations on quantum computing devices,” Phys. Rev. Lett. 126, 070504 (2021). [149] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M Chow, and Jay M Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum
27 magnets,” Nature 549, 242 (2017). [150] Rongxin Xia and Sabre Kais, “Quantum machine learning for electronic structure calculations,” Nat. Commun. 9, 4195 (2018). [151] Dylan J Gorman, Boerge Hemmerling, Eli Megidish, Soenke A Moeller, Philipp Schindler, Mohan Sarovar, and Hartmut Haeffner, “Engineering vibrationally assisted energy transfer in a trapped-ion quantum simulator,” Phys. Rev. X 8, 011038 (2018). [152] Yunseong Nam, Jwo-Sy Chen, Neal C Pisenti, Kenneth Wright, Conor Delaney, Dmitri Maslov, Kenneth R Brown, Stewart Allen, Jason M Amini, Joel Apisdorf, Kristin M. Beck, , Aleksey Blinov, Vandiver Chaplin, Mika Chmielewski, Coleman Collins, Shantanu Debnath, Kai M. Hudek, Andrew M. Ducore, Matthew Keesan, Sarah M. Kreikemeier, Jonathan Mizrahi, Phil Solomon, Mike Williams, Jaime David Wong-Campos, David Moehring, Christopher Monroe, and Jungsangothers Kim, “Ground-state energy estimation of the water molecule on a trapped-ion quantum computer,” npj Quantum Inf. 6, 33 (2020). [153] Bi-Xue Wang, Ming-Jie Tao, Qing Ai, Tao Xin, Neill Lambert, Dong Ruan, Yuan-Chung Cheng, Franco Nori, Fu-Guo Deng, and Gui-Lu Long, “Efficient quantum simulation of photosynthetic light harvesting,” npj Quantum Inf. 4, 1–6 (2018). [154] A. W. Chin, E. Mangaud, O. Atabek, and M. DesouterLecomte, “Coherent quantum dynamics launched by incoherent relaxation in a quantum circuit simulator of a light-harvesting complex,” Phys. Rev. A 97, 063823 (2018). [155] Anton Potočnik, Arno Bargerbos, Florian A. Y. N. Schröder, Saeed A. Khan, Michele C. Collodo, Simone Gasparinetti, Yves Salathé, Celestino Creatore, Christopher Eichler, Hakan E. Türeci, Alex W. Chin, and Andreas Wallraff, “Studying light-harvesting models with superconducting circuits,” Nat. Commun. 9, 1–7 (2018). [156] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Alán Aspuru-Guzik, and Jeremy L O’brien, “A variational eigenvalue solver on a photonic quantum processor,” Nat. Commun. 5, 4213 (2014). [157] Harper R Grimsley, Sophia E Economou, Edwin Barnes, and Nicholas J Mayhall, “An adaptive variational algorithm for exact molecular simulations on a quantum computer,” Nat. Commun. 10, 1–9 (2019). [158] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, Sergio Boixo, Michael Broughton, Bob B. Buckley, David A. Buell, Brian Burkett, Nicholas Bushnell, Yu Chen, Zijun Chen, Benjamin Chiaro, Roberto Collins, William Courtney, Sean Demura, Andrew Dunsworth, Edward Farhi, Austin Fowler, Brooks Foxen, Craig Gidney, Marissa Giustina, Rob Graff, Steve Habegger, Matthew P. Harrigan, Alan Ho, Sabrina Hong, Trent Huang, William J. Huggins, Lev Ioffe, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Cody Jones, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Seon Kim, Paul V. Klimov, Alexander Korotkov, Fedor Kostritsa, David Landhuis, Pavel Laptev, Mike Lindmark, Erik Lucero, Orion Martin, John M. Martinis, Jarrod R. McClean, Matt McEwen, Anthony Megrant, Xiao Mi, Masoud Mohseni, Wojciech Mruczkiewicz, Josh Mutus, Ofer
Naaman, Matthew Neeley, Charles Neill, Hartmut Neven, Murphy Yuezhen Niu, Thomas E. O’Brien, Eric Ostby, Andre Petukhov, Harald Putterman, Chris Quintana, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger, Vadim Smelyanskiy, Doug Strain, Kevin J. Sung, Marco Szalay, Tyler Y. Takeshita, Amit Vainsencher, Theodore White, Nathan Wiebe, Z. Jamie Yao, Ping Yeh, and Adam Zalcman, “Hartree-fock on a superconducting qubit quantum computer,” Science 369, 1084–1089 (2020), https://science.sciencemag.org/content/369/6507/1084.full.pdf. [159] Robert M. Parrish, Edward G. Hohenstein, Peter L. McMahon, and Todd J. Martinez, “Quantum Computation of Electronic Transitions Using a Variational Quantum Eigensolver,” Phys. Rev. Lett. 122, 230401 (2019). [160] Nikolay V. Tkachenko, James Sud, Yu Zhang, Sergei Tretiak, Petr M. Anisimov, Andrew T. Arrasmith, Patrick J. Coles, Lukasz Cincio, and Pavel A. Dub, “Correlation-informed permutation of qubits for reducing ansatz depth in the variational quantum eigensolver,” PRX Quantum 2, 020337 (2021). [161] Alba Cervera-Lierta, Jakob S. Kottmann, and Alán Aspuru-Guzik, “Meta-variational quantum eigensolver: Learning energy profiles of parameterized hamiltonians for quantum simulation,” PRX Quantum 2, 020329 (2021). [162] William J Huggins, Jarrod R McClean, Nicholas C Rubin, Zhang Jiang, Nathan Wiebe, K Birgitta Whaley, and Ryan Babbush, “Efficient and noise resilient measurements for quantum chemistry on near-term quantum computers,” npj Quantum Inf. 7, 1–9 (2021). [163] Jarrod R McClean, Nicholas C Rubin, Kevin J Sung, Ian D Kivlichan, Xavier Bonet-Monroig, Yudong Cao, Chengyu Dai, E Schuyler Fried, Craig Gidney, and Brendan Gimby, “Openfermion: the electronic structure package for quantum computers,” Quantum Sci. Technol. 5, 034014 (2020). [164] Mario Motta, Tanvi P Gujarati, Julia E Rice, Ashutosh Kumar, Conner Masteran, Joseph A Latone, Eunseok Lee, Edward F Valeev, and Tyler Y Takeshita, “Quantum simulation of electronic structure with a transcorrelated hamiltonian: improved accuracy with a smaller footprint on the quantum computer,” Phys. Chem. Chem. Phys. 22, 24270–24281 (2020). [165] Robert A. Lang, Ilya G. Ryabinkin, and Artur F. Izmaylov, “Unitary transformation of the electronic hamiltonian with an exact quadratic truncation of the baker-campbell-hausdorff expansion,” J. Chem. Theory Comput. 17, 66–78 (2021). [166] Ilya G. Ryabinkin, Tzu-Ching Yen, Scott N. Genin, and Artur F. Izmaylov, “Qubit coupled cluster method: A systematic approach to quantum chemistry on a quantum computer,” J. Chem. Theory Comput. 14, 6317– 6326 (2018). [167] Artur F. Izmaylov, Tzu-Ching Yen, Robert A. Lang, and Vladyslav Verteletskyi, “Unitary partitioning approach to the measurement problem in the variational quantum eigensolver method,” J. Chem. Theory Comput. 16, 190–195 (2020). [168] Benjamin P Lanyon, James D Whitfield, Geoff G Gillett, Michael E Goggin, Marcelo P Almeida, Ivan Kassal, Jacob D Biamonte, Masoud Mohseni, Ben J Powell, Marco Barbieri, et al., “Towards quantum chem-
28 istry on a quantum computer,” Nature chemistry 2, 106–111 (2010). [169] Dawei Lu, Nanyang Xu, Ruixue Xu, Hongwei Chen, Jiangbin Gong, Xinhua Peng, and Jiangfeng Du, “Simulation of chemical isomerization reaction dynamics on a nmr quantum simulator,” Physical Review Letters 107, 020501 (2011). [170] Cornelius Hempel, Christine Maier, Jonathan Romero, Jarrod McClean, Thomas Monz, Heng Shen, Petar Jurcevic, Ben P. Lanyon, Peter Love, Ryan Babbush, Alán Aspuru-Guzik, Rainer Blatt, and Christian F. Roos, “Quantum chemistry calculations on a trappedion quantum simulator,” Phys. Rev. X 8, 031022 (2018). [171] Google AI Quantum, Collaborators*†, Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Sergio Boixo, Michael Broughton, Bob B Buckley, et al., “Hartree-fock on a su-
perconducting qubit quantum computer,” Science 369, 1084–1089 (2020). [172] R. G. Newton, Scattering Theory of Waves and Particles (Springer-Verlag, New York, 1982). [173] S. S. Iyengar, D. J. Kouri, and D. K. Hoffman, “Particular and homogeneous solutions of time-independent wavepacket schrodinger equations: Calculations using a subset of eigenstates of undamped or damped hamiltonians.” Theor. Chem. Acc. 104, 471 (2000). [174] Debadrita Saha, Philip Richerme, and Srinivasan S. Iyengar, “Quantum circuit and mapping algorithms for wavepacket dynamics: Case study of anharmonic hydrogen bonds in protonated and hydroxide water clusters,” Journal of Chemical Theory and Computation 21, 3814–3831 (2025). [175] Robert R Tucci, “An introduction to cartan’s kak decomposition for qc programmers,” arXiv preprint quantph/0507171 (2005).