Digital Quantum Simulation of Wavepacket Correlations in a Chemical Reaction

TL;DR

This paper demonstrates how digital quantum circuits can accurately simulate wavepacket correlations in a chemical reaction, showing promise for quantum chemistry applications.

Contribution

A new hybrid quantum-classical method for computing wavepacket correlation functions using digital quantum simulation and a reduced-depth MFE protocol.

Findings

Quantum-estimated correlation functions match classical simulations across the full time domain.

The MFE protocol achieves accurate results with reduced circuit depth.

The method is validated on the H + H2 exchange reaction with high accuracy.

Abstract

We present hybrid quantum–classical algorithms to compute time-dependent Møller wavepacket correlation functions via digital quantum simulation. Reactant and product channel wavepackets are encoded as qubit states, evolved under a discretized molecular Hamiltonian, and their correlation is reconstructed using both a modified Hadamard test and a multi-fidelity estimation (MFE) protocol. The method is applied to the collinear H + H2 exchange reaction on a London–Eyring–Polanyi–Sato potential energy surface. Quantum-estimated correlation functions show quantitative agreement with high-resolution classical wavepacket simulations across the full time domain, reproducing both short-time scattering peaks and long-time oscillatory dynamics. The ancilla-free MFE protocol achieves matching results with reduced circuit depth. These results provide a proof of principle that digital quantum circuits can be used to accurately calculate the wavepacket correlation functions for a benchmark chemical reaction system.

Figures

Full text

1. Introduction

Understanding how quantum states evolve in time lies at the foundation of modern chemical physics, condensed matter theory, and quantum information science. Many dynamical processes of physical interest—including energy transfer and wavepacket propagation—can be characterized through time-dependent correlation functions, which quantify the overlap between quantum states as they evolve under a given Hamiltonian. These correlation functions provide direct insights into coherence, interference, and the temporal structure of molecular motion, making them essential tools for analyzing quantum dynamical systems.

Classically computing such correlations often requires propagating high-dimensional wavefunctions on strongly coupled potential energy surfaces, a task that becomes increasingly demanding as the system size grows [1]. Quantum computers offer a natural platform for simulating unitary dynamics, motivating the development of digital quantum algorithms capable of estimating dynamical overlaps with scalability beyond classical limits [2,3,4]. Recent progress in quantum hardware and circuit design has made it possible to explore such algorithms in realistic molecular contexts, even within the constraints of near-term devices [5,6,7,8,9].

In this work, we present a hybrid quantum–classical framework for computing time-dependent correlation functions using digital quantum simulation. Our method implements the time-dependent Møller operator construction of reactant and product wavepackets on quantum circuits, allowing the correlation between two wavepackets to be expressed as a measurable quantum expectation value. We employ two complementary correlation estimation strategies: the modified Hadamard test, which uses an ancilla qubit to extract complex-valued overlaps, and the multi-fidelity estimation (MFE) protocol, which performs overlap estimation without controlled operations [10]. These approaches provide a flexible set of tools with reasonable circuit depth while retaining full access to the real and imaginary components of the correlation function.

To demonstrate the method, we consider the collinear H + H_2_ reaction system, modeled through wavepacket propagation on a London–Eyring–Polanyi–Sato (LEPS) potential energy surface. Vibrational motion is represented using Morse oscillator eigenfunctions, while translational motion is encoded in Gaussian wavepackets. The full molecular Hamiltonian is discretized on a two-dimensional grid in bond coordinates and mapped to qubits. Time evolution is implemented through unitary gates, and the resulting quantum-estimated correlation functions are benchmarked against high-resolution classical simulations of standard quantum calculations.

The results show close agreement between the quantum circuit and classical simulations of standard quantum calculations of correlation functions across the full time domain, capturing the expected long-time behavior and phase evolution of the wavepackets. This demonstrates that digital quantum simulation can reliably reproduce dynamical features encoded in time-dependent overlaps [11].

Overall, the present study establishes correlation function computation as a practical and informative intermediate step toward the quantum simulation of molecular dynamics. The combination of ancilla-based and ancilla-free estimation schemes provides circuit flexibility compatible with near-term quantum hardware. These results highlight the potential of digital quantum simulation to serve as a scalable framework for studying quantum dynamics and lay the groundwork for future extensions to higher-dimensional systems and more complex molecular processes.

Previous quantum simulation studies have addressed molecular ground-state energies, vibrational spectra, and short-time dynamics; however, the direct computation of time-dependent correlation functions using Møller operators has remained largely unexplored within a fully digital quantum framework [10,12]. The present work fills this gap by introducing a circuit-level formulation of channel-resolved wavepacket correlations and benchmarking two complementary quantum circuit estimators against classical computer simulations of standard quantum mechanical correlation calculations.

From an information-theoretic viewpoint, the Møller correlation $[eqn]$ is a Loschmidt echo-type overlap that quantifies how much dynamical information about the reaction is retained in the evolving wavepacket. Its magnitude $[eqn]$ is equivalent to the quantum fidelity between reactant and product preparations [13].

To summarize, in this work, we (i) formulate time-dependent Møller wavepacket correlation functions as quantum-measurable overlaps; (ii) conceptually experiment with a modified Hadamard test and a multi-fidelity estimation (MFE) protocol as prototype implementations, primarily to confirm the consistency between the measured data and the predicted theoretical bounds for the H + H_2_ reaction; and (iii) perform conceptual data confirmation by comparing the resulting channel-resolved correlations with high-resolution standard quantum calculations on classical computer wavepacket propagation and correlation. Our aim is proof of principle: we show that chemically meaningful overlaps can be computed on digital quantum circuits and are compatible with near-term hardware constraints.

2. Model: Collinear Hydrogen Exchange Reaction

$[eqn]$ serves as a fundamental prototype for studying reactive molecular dynamics [14,15,16,17,18]. Owing to its minimal number of degrees of freedom and well-established potential energy surfaces, the model provides an ideal setting to test new numerical and quantum simulation frameworks. Despite its simplicity, the system exhibits key dynamical features found in more complex reactions, including barrier crossing, transient energy redistribution, and vibrational state selectivity.

2.1. Bond and Jacobi Coordinates

In a collinear triatomic system, shown in Figure 1, two coordinate systems are commonly used to represent nuclear motion.

The first is Jacobi coordinates, which naturally describe individual channels. For the H + H_2_ system, the reactant arrangement channel is described by $[eqn]$ , where $[eqn]$ is the bond length of the diatomic H_2_, and $[eqn]$ is the separation between the incoming H atom and the diatomic center of mass. Similarly, the product arrangement channel is described by $[eqn]$ , where $[eqn]$ is the bond length of the newly formed H_2_, and $[eqn]$ is the distance from the scattered H atom to the product diatomic center of mass.

Channel definition: In this work, channel refers to a choice of vibrational quantum numbers associated with the asymptotic reactant and product arrangements. We label the reactant vibrational quantum number $[eqn]$ and the product vibrational quantum number $[eqn]$ . Throughout this study, the initial reactant is prepared in its vibrational ground state $[eqn]$ , and we consider two product channels corresponding to $[eqn]$ (vibrationally elastic scattering) and $[eqn]$ (transition to the first vibrationally excited product state).

The second is Bond coordinates, defined simply as the interatomic distances

$$[eqn]$$

This representation provides a single $[eqn]$ grid on which both reactant and product wavepackets can be expressed simultaneously.

The two descriptions are related by a straightforward linear transformation. Specifically, the bond coordinates $[eqn]$ can be obtained from either the reactant Jacobi coordinates $[eqn]$ or the product Jacobi coordinates $[eqn]$ as

$$[eqn]$$

In practice, the initial wavepackets are prepared in Jacobi coordinates, reflecting their channel character, and then transformed into the common bond coordinate grid $[eqn]$ . All subsequent propagation and correlation function calculations are carried out in this bond coordinate representation, which unifies the treatment of reactant and product channels.

2.2. Structure of the Wavepacket

For each channel, the wavepacket is constructed as a product of three factors:

$$[eqn]$$

Here,

1. $[eqn]$ describes the internal vibration of the diatomic molecule (H_2_);2.The Gaussian envelope $[eqn]$ ensures the localization of the wavepacket in the translational coordinate R, centered at $[eqn]$ ;3.The carrier plane wave $[eqn]$ gives the packet a mean momentum $[eqn]$ , directed either inward (reactant) or outward (product) relative to the interaction region.

The prefactor is chosen to guarantee square-integrability and to match the conventional normalization of Gaussian wavepackets:

$$[eqn]$$

Vibrational Component: Morse Oscillator

The vibrational eigenstates of H_2_ are modeled using the Morse oscillator, which provides a realistic one-dimensional description of bond stretching. The eigenfunction for vibrational quantum number n is

$$[eqn]$$

where $[eqn]$ are generalized Laguerre polynomials. The dimensionless coordinate z and parameter $[eqn]$ are defined as

$$[eqn]$$

Here,

$[eqn]$ is the dissociation energy of the bond;a is an inverse range parameter controlling the curvature of the potential at equilibrium; $[eqn]$ is the equilibrium bond length; $[eqn]$ is the reduced mass of the diatomic fragment (for H-H, $[eqn]$ in atomic units).

The normalization constant $[eqn]$ is fixed by the orthonormality of Morse eigenstates:

$$[eqn]$$

2.3. Interaction Potential

To model the nuclear motion in the collinear H + H_2_ reaction, we employ a London–Eyring–Polanyi–Sato (LEPS)-type potential energy surface. In this construction, each pairwise H-H interaction is described by an effective diatomic potential, and the three pair potentials are then combined into a three-center form that interpolates smoothly between reactant and product valleys while generating a single barrier along the exchange path. The LEPS ansatz captures the essential barrier–well topology of the H_3_ system with relatively few parameters and has long served as a benchmark surface for time-dependent studies of hydrogen exchange. In the following, we summarize the two-body building blocks and their London-type three-body combination, which are used to define the interaction potential in our simulations.

2.3.1. Two-Body Forms: Morse and Anti-Morse Potentials

Each pairwise interaction is represented by two effective potential forms:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the standard Morse potential describing the covalent bond well, and $[eqn]$ is an anti-Morse counterpart that emphasizes the short-range repulsive character while retaining exponential decay. Following the London construction, these are combined into symmetric and antisymmetric forms, interpreted as approximate Coulomb and exchange contributions:

$$[eqn]$$ $$[eqn]$$

2.3.2. Three-Center London Combination

For the collinear triatomic system, we evaluate Q and $[eqn]$ on all three bonds. The London-type three-body combination is

$$[eqn]$$

where $[eqn]$ sets the sign of the exchange-like term. In this work, we take $[eqn]$ , which lowers the energy when one bond is near its covalent minimum while the others are elongated, producing a barrier–well profile consistent with the H + H_2_ exchange reaction.

The resulting surface exhibits a single saddle corresponding to the collinear exchange barrier, with reactant and product valleys along the asymptotic channels. The qualitative topology of $[eqn]$ agrees with benchmark LEPS representations used in classical and semiclassical scattering studies [19,20,21]. Historically, this construction corresponds to the London–Eyring–Polanyi–Sato (LEPS) potential energy surface, which was evaluated empirically for H + H_2_ by Cashion and Herschbach [22] and remains a benchmark model for testing reaction models [14,15].

For our chemical reaction model, we use (in a.u.)

$$[eqn]$$

which reproduce a qualitatively reasonable covalent well for each pair in Equation (8) and generate a barrier profile consistent with collinear exchange under Equation (12).

The spatial structure of the reactant and product wavepackets in bond coordinates $[eqn]$ is shown in Figure 2, and the corresponding LEPS potential energy surface governing the interaction dynamics is shown in Figure 3.

2.4. Discrete Kinetic Energy Operator in Bond Coordinates

The total Hamiltonian on the $[eqn]$ bond coordinate grid requires, in addition to the London potential, a representation of the nuclear kinetic energy operator. For a collinear triatomic system with equal masses (hydrogen atoms), the kinetic part of the Hamiltonian in bond coordinates can be expressed as

$$[eqn]$$

where $[eqn]$ is the mass of a hydrogen atom in atomic units. The first two terms correspond to the standard Laplacian in X and Y, while the mixed derivative $[eqn]$ arises due to the non-orthogonal nature of the bond coordinate system.

3. Methods: Møller Operator Formulation of Wavepacket Propagation

In the time-dependent formulation of wave propagation theory, the central objects are the Møller operators, which connect asymptotic free states to interacting scattering states. Let $[eqn]$ denote the channel Hamiltonian that governs the non-interacting asymptotic dynamics, and let $[eqn]$ be the full Hamiltonian including the interaction potential V. The Møller operators are defined as

$$[eqn]$$

Physically, $[eqn]$ maps a free state prepared in the distant past to the corresponding incoming reactant state, while $[eqn]$ maps a free state prepared in the distant future to the corresponding outgoing product state. Thus, the ‘+’ sign encodes incoming boundary conditions and the ‘−’ sign encodes outgoing boundary conditions.

Acting with $[eqn]$ on an asymptotic channel basis state $[eqn]$ , characterized by a relative momentum $[eqn]$ and an internal state label $[eqn]$ , produces the interacting scattering state

$$[eqn]$$

These states retain the asymptotic character of the channel $[eqn]$ far from the interaction region but evolve according to the full Hamiltonian H near the reaction zone.

In the present work, the term channel refers to a choice of vibrational quantum numbers for the reactant and product H-H bond. We label the reactant vibrational level $[eqn]$ and the product vibrational level $[eqn]$ . Throughout, we consider an initial reactant prepared in its vibrational ground state, $[eqn]$ , and we focus on two product channels: $[eqn]$ (vibrationally elastic scattering) and $[eqn]$ (transition to the first vibrationally excited product state) [23]. The corresponding channel-resolved time-dependent correlation function is

$$[eqn]$$

For the choices above, $[eqn]$ measures the time-dependent probability amplitude for the system to remain in the vibrational ground state, while $[eqn]$ measures the amplitude for an inelastic transition into the first excited vibrational channel. The peaks and long-time oscillations of these quantities, shown later in the Section 6, encode the dynamical signatures of barrier crossing and vibrational recurrences in the H + H_2_ exchange reaction [16,24,25].

In other words, $[eqn]$ is the overlap between a reactant-like wavepacket $[eqn]$ propagated forward in time by t and a product-like wavepacket $[eqn]$ propagated backward in time by t, with both wavepackets evaluated at a common reference time. Varying t traces the dynamical progression of the scattering process, with the transient reaction event reflected in the behavior of $[eqn]$ near $[eqn]$ . This correlation function contains the same physical information as the traditional S-matrix elements, but it is naturally expressed in terms of time evolution operators and overlaps, which are directly implementable on a quantum computer as unitary circuits. It is precisely this quantity that we encode and measure with the quantum algorithms described in the following sections.

4. Standard Quantum Calculations

To obtain the results reported here, we used a uniform bond coordinate grid with $[eqn]$ points spanning $[eqn]$ , giving a Hilbert-space dimension $[eqn]$ . The time-dependent correlations were evaluated on a uniform time grid $[eqn]$ with step size $[eqn]$ , which resolves both the short-time scattering peak and long-time vibrational recurrences [15,25,26,27,28].

These standard quantum calculations on classical computer simulations, shown in Figure 4, provide a reference for assessing the quantum algorithms. In particular, they show that physically relevant information about the reaction and product overlap and crossing, partial transmission into the product channel, and long-time vibrational recurrences is encoded in the time-resolved overlaps $[eqn]$ . The key question that we address in the following sections is whether digital quantum circuits can reproduce these features with sufficient accuracy to serve as reliable surrogates for standard quantum calculations on classical computer wavepacket propagation.

5. Quantum Algorithms and Implementation

5.1. Hamiltonian Construction in the Grid Basis

The total Hamiltonian for the collinear H + H_2_ reaction is written as

$$[eqn]$$

where X and Y denote the bond coordinates, $[eqn]$ and $[eqn]$ are the kinetic energy operators, and $[eqn]$ is the interaction potential. The asymptotic Hamiltonian is given by

$$[eqn]$$

corresponding to free propagation in the absence of interaction. The kinetic energy operator in bond coordinates contains both diagonal and mixed derivative terms,

$$[eqn]$$

where $[eqn]$ and $[eqn]$ denote the second derivatives with respect to X and Y, respectively, and $[eqn]$ is the mixed kinetic term arising from the coordinate transformation from Jacobi to bond coordinates. This mixed term encodes the kinematic coupling between the two bond lengths and is essential in accurately describing the collinear three-body dynamics.

Both Hamiltonians are represented on a uniform two-dimensional grid in $[eqn]$ . The kinetic energy operators are discretized using second-order finite difference approximations, while the interaction Hamiltonian is obtained by evaluating the potential energy surface $[eqn]$ pointwise on the grid. This representation captures both the translational motion along the reaction coordinate and the vibrational dynamics of the diatomic fragment.

The resulting Hamiltonian matrices are then mapped to qubit operators by expressing them as linear combinations of Pauli strings, which are subsequently used to construct the time evolution operator $[eqn]$ in the quantum circuit.

5.2. Grid-Based Representation of Channel Wavepackets

The reactant and product channel states, denoted by $[eqn]$ and $[eqn]$ , incorporate both the translational motion along the reaction coordinate and the vibrational dynamics of the diatomic fragment. In our implementation, these states are represented using a first-quantized, grid-based encoding on a two-dimensional bond coordinate grid $[eqn]$ , following the standard real-space discretization approach used in quantum simulations of molecular scattering [4].

Each channel wavepacket is initially defined in its natural Jacobi coordinates $[eqn]$ as a separable product

$$[eqn]$$

where $[eqn]$ labels the reactant and product channels, respectively. The vibrational component $[eqn]$ is taken to be an eigenstate of the Morse potential, while the translational component $[eqn]$ is chosen as a Gaussian wavepacket with central momentum $[eqn]$ and momentum width $[eqn]$ ,

$$[eqn]$$

To place both channels on a common computational grid, the Jacobi coordinates are expressed in terms of the bond coordinates $[eqn]$ appropriate for the collinear H + H_2_ geometry. For the reactant channel,

$$[eqn]$$

while, for the product channel,

$$[eqn]$$

Using these channel-dependent coordinate maps, the continuous wavepackets are evaluated directly on the $[eqn]$ grid,

$$[eqn]$$

The grid is discretized uniformly as $[eqn]$ and $[eqn]$ , with spacings $[eqn]$ and $[eqn]$ . After discretization, each channel wavepacket is normalized according to

$$[eqn]$$

5.3. Qubit Encoding and State Preparation

For the grid size used here ( $[eqn]$ ), the system dimension $[eqn]$ maps to $[eqn]$ qubits. The modified Hadamard test therefore uses $[eqn]$ qubits (including the ancilla), whereas the MFE protocol uses $[eqn]$ qubits without an ancilla.

Each bond coordinate is represented on a discrete grid of size $[eqn]$ in the $[eqn]$ bond coordinate plane. The joint Hilbert space of these grid points is mapped to a register of $[eqn]$ qubits.

In this proof-of-principle work, we load these states using a dense state preparation unitary obtained from the classical amplitude vectors, i.e., we assume access to a circuit that maps $[eqn]$ to $[eqn]$ exactly. The Gaussian envelope parameters, such as the central momentum $[eqn]$ and width $[eqn]$ , are chosen so that the initial wavepacket is localized near a chosen bond length $[eqn]$ and carries sufficient kinetic energy to traverse the barrier on the LEPS potential.

This approach of loading states using a dense state preparation unitary is restricted to small Hilbert spaces and is used here exclusively for validation and benchmarking. We do not claim this method to be scalable; rather, it allows us to isolate and test the algorithmic structure of the correlation estimation protocols. Scalable implementations would instead employ structured or approximate state preparation techniques.

5.4. Modified Hadamard Test

To evaluate the time-dependent correlation function,

$$[eqn]$$

we implement a modified Hadamard test circuit augmented with controlled state preparation unitaries for both the reactant and product wavepackets.

Circuit Structure

Figure 5 shows the layout. The ancilla qubit is first prepared in a superposition state via a Hadamard gate. Conditioned on the ancilla, we then apply the following:

• 1.Controlled preparation of the reactant state $[eqn]$ ;
• 2.Controlled short-time Møller operator blocks $[eqn]$ and $[eqn]$ built from the non-interacting Hamiltonian $[eqn]$ ;
• 3.A controlled real-time propagator $[eqn]$ ;
• 4.Controlled preparation of the product state $[eqn]$ , implemented by flipping the ancilla into $[eqn]$ , applying the preparation;
• 5.Additional $[eqn]$ blocks acting only when the ancilla is in $[eqn]$ , ensuring the proper phase ordering in Equation (27).

Finally, another Hadamard gate on the ancilla and a projective measurement in the computational basis yield the expectation value

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the measured probabilities of finding the ancilla in states $[eqn]$ and $[eqn]$ , respectively. This relation follows from the standard Hadamard test identity.

The quantum circuit used to measure the imaginary part $[eqn]$ uses the same modified Hadamard test structure as in Figure 5, with an additional phase gate $[eqn]$ applied to the ancilla after the first Hadamard. This phase shift rotates the interference pattern on the ancilla such that the measured $[eqn]$ expectation value is proportional to $[eqn]$ , while the underlying state preparation and time evolution blocks remain unchanged.

Each unitary block appearing in Figure 5 has a specific physical meaning. The operators $[eqn]$ and $[eqn]$ denote the state preparation unitaries that map the computational basis state $[eqn]$ to the reactant and product channel wavepackets, $[eqn]$ and $[eqn]$ , respectively. The unitary $[eqn]$ represents time evolution under the full interacting Hamiltonian, while $[eqn]$ are short-time propagators generated by the asymptotic Hamiltonian and appear in the Møller operator construction. Controlled versions of these unitaries are used in the modified Hadamard test to encode the desired overlap into the ancilla qubit.

Plots of the Hadamard test correlation functions are shown in Figure 6 and Figure 7.

5.5. Multi-Fidelity Estimation (MFE) Method

The multi-fidelity estimation (MFE) method provides a resource-efficient way to estimate complex-valued quantum overlaps such as time-dependent correlation functions without relying on controlled unitaries or ancillary qubits. In this approach, instead of using the conventional Hadamard test, which demands extra circuit depth, the desired matrix element $[eqn]$ is inferred from measurable state fidelities [29,30,31].

The method begins by preparing a normalized superposition state

$$[eqn]$$

where the reference state $[eqn]$ and the target state $[eqn]$ belong to different symmetry sectors of the Hamiltonian. A second superposition,

$$[eqn]$$

is then constructed to serve as a projector state.

The overlap between these two states yields the fidelity

$$[eqn]$$

while a direct fidelity measurement

$$[eqn]$$

provides the magnitude $[eqn]$ of the correlation function. The reference amplitude is independently evaluated as

$$[eqn]$$

where $[eqn]$ and $[eqn]$ denote the magnitude and phase of the reference state under time evolution. The phase of the target correlation function is then reconstructed via

$$[eqn]$$

yielding the complete complex overlap

$$[eqn]$$

Thus, the MFE protocol extracts both the amplitude and phase of time-dependent overlaps solely from fidelity measurements between evolved and reference states, avoiding controlled operations and making it especially suitable for near-term quantum hardware applications such as computing scattering matrix elements.

To summarize,

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are measurable fidelities, while $[eqn]$ and $[eqn]$ denote the magnitude and phase of the reference state overlap $[eqn]$ . This formulation explicitly recovers both the amplitude and phase of the complex overlap from fidelity measurements alone, with the total measurement overhead scaling polynomially as $[eqn]$ [10].

Circuit Structure

Advantages

Avoids the use of controlled unitary operations and ancilla qubits;Reduces the circuit depth, enabling execution on Noisy Intermediate-Scale Quantum (NISQ) devices;Allows the efficient and accurate estimation of time correlation functions and scattering amplitudes.

Both the modified Hadamard test and the multi-fidelity estimation (MFE) protocol reproduce the same time-dependent correlation structure within numerical precision. In particular, the short-time scattering peak, the long-time oscillatory behavior, and the relative phase evolution of the correlation function are consistently captured by both estimators. The MFE protocol achieves this without the use of controlled time evolution gates, as shown in Figure 8, significantly reducing the circuit depth while maintaining full complex overlap reconstruction; this reduction reflects an algorithmic advantage in scaling behavior rather than immediate suitability for NISQ hardware. This establishes MFE as a near-term-compatible alternative to the Hadamard test for computing scattering-type time correlation functions on digital quantum hardware.

5.6. Numerical Implementation

The circuit was implemented in Cirq [32] using dense matrix gates constructed from the kinetic plus London potential Hamiltonian on an $[eqn]$ $[eqn]$ bond coordinate grid. For each time step t in the range $[eqn]$  a.u., the unitary propagator $[eqn]$ was computed using scipy.linalg.expm and wrapped as a controlled gate.

In the simulations reported here, we used a uniform bond coordinate grid with $[eqn]$ points along each coordinate, spanning the interval $[eqn]$ . The resulting Hilbert-space dimension is $[eqn]$ so that the grid can be encoded exactly in an $[eqn]$ -qubit register. The channel-resolved correlation functions $[eqn]$ were evaluated on a uniform time grid $[eqn]$ which was sufficient to resolve both the short-time scattering peak near the barrier region and the long-time oscillatory recurrences in the asymptotic product channel.

The reactant and product scattering wavepackets $[eqn]$ and $[eqn]$ are loaded using the Cirq state preparation scheme. For each target state $[eqn]$ on the D-dimensional grid, we construct a unitary $[eqn]$ such that $[eqn]$ . The resulting $[eqn]$ matrices $[eqn]$ and $[eqn]$ are then wrapped as MatrixGates in Cirq and reused at all times t inside the Hadamard and MFE circuits. This dense approach is not scalable to very large D, but it is sufficient for the modest grid sizes considered here and allows us to isolate the algorithmic features of the correlation estimation protocols.

The $[eqn]$ grid states are encoded in a register of $[eqn]$ qubits, such that each computational basis state corresponds to a single grid point $[eqn]$ . For the grid size used here, this gives $[eqn]$ and $[eqn]$ . The modified Hadamard test therefore acts on a total of $[eqn]$ qubits (one ancilla plus the $[eqn]$ -qubit system register), whereas the MFE protocol uses two system registers and no ancilla, for a total of $[eqn]$ qubits. These sizes are modest enough that we can represent the relevant state preparation and time evolution blocks as dense unitaries in the present proof-of-principle study, while still illustrating the scaling behavior of each algorithm.

All results shown in Figure 9, Figure 10, Figure 11 and Figure 12 were obtained using classical simulation of the quantum circuits and were not executed on physical quantum hardware. The purpose of these simulations is to validate the algorithmic construction and estimator equivalence in a controlled, noiseless setting. Hardware noise effects are therefore not included in the reported figures. Incorporating device-specific noise models and error mitigation strategies is an important direction for future work but is beyond the scope of this proof-of-principle algorithmic study.

6. Results Comparison

6.1. Hadamard vs. MFE

Figure 11 shows that the agreement between the modified Hadamard test and the MFE protocol is not only qualitative but quantitative at the level of numerical precision: the pointwise deviation $[eqn]$ stays at the level of $[eqn]$ over the entire propagation window. To quantify the global discrepancy, we computed discrete $[eqn]$ and $[eqn]$ norms of the difference $[eqn]$ over all time samples. For the elastic channel $[eqn]$ , we obtain

$$[eqn]$$

while, for the inelastic channel $[eqn]$ , the corresponding values are

$$[eqn]$$

These norms are at or below double-precision roundoff and many orders of magnitude smaller than the typical correlation magnitude ( $[eqn]$ near the short-time peak). Thus, the ancilla-free MFE construction reproduces the overlap estimated by the deeper, ancilla-controlled Hadamard test to machine precision in both scattering channels. Since MFE eliminates the need for a global controlled $[eqn]$ and replaces it with SWAP tests and simpler time evolutions, this matching accuracy at dramatically reduced circuit overhead is a key algorithmic advantage, consistent with depth reduction arguments in the literature.

6.2. Standard Quantum vs. Quantum Circuit Correlations

As a basic consistency check, we verify that the quantum circuit-evaluated correlation at $[eqn]$ , $[eqn]$ agrees with the direct inner product $[eqn]$ of the prepared reactant and product wavepackets on a classical computer. For both channels, the real and imaginary parts of $[eqn]$ agree with the standard quantum mechanical reference to better than $[eqn]$ , confirming that the state preparation and overlap estimation procedures are implemented correctly.

Figure 11 compares the modified Hadamard test and the multi-fidelity estimation (MFE) protocol for the present problem. The two estimators show excellent numerical agreement, with residual differences at the level of floating-point roundoff, indicating that the ancilla-free MFE protocol reproduces the same correlation values as the Hadamard test within the precision accessible here. Figure 12 places this algorithmic agreement in a physical context by comparing the quantum circuit-computed correlation magnitudes with the standard quantum calculation wavepacket reference.

The pronounced peak in $[eqn]$ near $[eqn]$ arises when the reactant- and product-like wavepackets overlap in the interaction region, providing a direct time-domain signature of the scattering event. At later times, the correlation magnitude exhibits weaker, quasi-periodic oscillations whose structure depends on the product channel. These long-time features originate from vibrational recurrences of the bound H-H pair in the product well and are consistent with established standard quantum mechanical calculations on classical computer wavepacket descriptions of the H + H_2_ exchange reaction. Small quantitative deviations between the quantum circuit and standard quantum calculations curves can be attributed to the finite grid resolution and statistical sampling effects.

Taken together, these results show that (i) the MFE protocol reproduces the complex time-dependent correlation function with accuracy comparable to the Hadamard test for the present simulations and (ii) the quantum-computed correlations capture the essential dynamical features of the scattering process—both the transient reaction peak and the channel-dependent long-time oscillations—while remaining quantitatively close to high-resolution standard quantum calculations on classical computer wavepacket results.

To validate the correctness of the quantum circuit implementation independently of any specific software backend, we performed several internal consistency checks. First, the exact identity $[eqn]$ was verified by comparing the circuit-evaluated overlap at $[eqn]$ with the direct inner product of the prepared wavepackets, with agreement better than $[eqn]$ . Second, the modified Hadamard test and the multi-fidelity estimation (MFE) protocol—two algorithmically distinct estimators—were shown to reproduce near-identical correlation functions. These checks relied on exact algebraic relations and estimator equivalence and therefore provide validation beyond agreement with external standard quantum calculations on a classical computer propagation code.

7. Algorithmic Complexity and Error Sources

In the present work, time evolution under the full Hamiltonian H and the asymptotic Hamiltonian $[eqn]$ is implemented by explicitly constructing dense unitary propagators. The time evolution operator $[eqn]$ is evaluated numerically via matrix exponentiation and then wrapped as a MatrixGate in the Cirq framework. This avoids Trotter–Suzuki or product formula decompositions at the circuit level but shifts the main computational bottleneck to the cost of generating and storing dense matrices on a classical processor.

This choice was made for benchmarking purposes: since the Hilbert space remains classically tractable, we can compute $[eqn]$ exactly using matrix exponentiation and validate the quantum circuit output against high-precision classical simulations of standard quantum calculations. This approach ensures that all observed differences arise from quantum circuit behavior and not from numerical artifacts in the classical simulations of standard quantum calculations.

Looking ahead, this method is expected to provide significant advantages in regimes where direct classical simulations of standard quantum calculations become infeasible. As the system size increases, the cost of storing and applying $[eqn]$ as a dense matrix scales exponentially, while a quantum computer can implement the same operator through sequences of native gates (e.g., via Trotterization or qubitization) with polynomial overhead. In such cases, our correlation function-based algorithm can be used to simulate correlation dynamics intractable for classical simulations of standard quantum calculations, offering a clear path toward a quantum advantage.

We emphasize that these numerical experiments are intended as conceptual demonstrations of the proposed algorithms on small-scale instances using classical simulators of quantum circuits, rather than as full software performance benchmarks.

7.1. Challenges of Controlled Unitaries in Hadamard Test

The Hadamard test requires controlled versions of both the propagation operator $[eqn]$ and the state preparation unitaries $[eqn]$ and $[eqn]$ that map the vacuum state to reactant and product wavepackets. Mathematically,

$$[eqn]$$

so, in circuit form, each U must be promoted to a controlled operation [33]. If a dense matrix gate acts on n system qubits, its controlled version requires $[eqn]$ qubits and a substantially larger number of two-qubit entangling gates. As D grows, this leads to a sharp increase in both the gate count and circuit depth, even when $[eqn]$ has been precomputed classically [34].

Sampling Error in the Hadamard Test

In the Hadamard test implementation, errors arise from statistical sampling during the execution of the quantum circuit. At each time point t in the correlation grid, the estimator $[eqn]$ obtained from N circuit repetitions (shots) fluctuates around the exact value $[eqn]$ due to measurement noise. This sampling error is bounded by a Chernoff bound so that the number of samples required to estimate the correlation function $[eqn]$ with additive precision $[eqn]$ and a success probability of at least $[eqn]$ scales as $[eqn]$ Thus, the measurement overhead grows polynomially in $[eqn]$ but is independent of the Hilbert-space dimension for a fixed circuit depth.

The rapid growth in circuit depth associated with controlled unitaries is the primary motivation for the multi-fidelity estimation (MFE) protocol. By eliminating global controlled time evolution, MFE provides a pathway for scaling correlation function estimation even when controlled versions of large unitaries become prohibitively expensive.

7.2. Matrix Exponentiation and Storage Cost

For a bond coordinate grid of dimension $[eqn]$ , the Hamiltonian is represented as a $[eqn]$ sparse matrix. Computing the exact propagator $[eqn]$ by dense matrix exponentiation requires $[eqn]$ operations using standard classical linear algebra routines, together with $[eqn]$ memory to store the resulting dense matrix. This cost is manageable for the modest grid sizes ( $[eqn]$ ) employed in our demonstrations but becomes prohibitive for larger discretizations. In practice, the present implementation is therefore limited primarily by classical computer precomputation rather than by the quantum circuit depth.

7.3. Circuit Depth Estimate for MFE

We briefly quantify the reduction in circuit depth afforded by the multi-fidelity estimation (MFE) protocol relative to the Hadamard test. Consider a unitary time evolution operator $[eqn]$ acting on N qubits, decomposed into a circuit containing p CNOT gates and q single-qubit gates. In the Hadamard test, implementing the corresponding controlled unitary requires p Toffoli gates and q controlled single-qubit gates. After decomposing each Toffoli gate into CNOTs, the controlled unitary contains between $[eqn]$ and $[eqn]$ CNOTs, in addition to further two-qubit and single-qubit operations.

By contrast, the MFE protocol requires only $[eqn]$ CNOTs, where r denotes the number of CNOT gates needed to prepare the reference superposition state. Consequently, although the multiplicative factor in front of p in the Hadamard test may appear modest, it substantially limits the achievable circuit depth. Under comparable assumptions, the MFE protocol can therefore enable between three and six times more time steps than the Hadamard test on the same hardware.

Table 1 summarizes the circuit resources required by the two estimators and provides a direct, engineering-level comparison beyond asymptotic depth arguments.

8. Discussion and Conclusions

In summary, this work demonstrates that digital quantum circuits can reproduce channel-resolved correlation overlaps for the benchmark H + H_2_ exchange reaction. Both an ancilla-based Hadamard test and the ancilla-free multi-fidelity estimation (MFE) protocol yield numerically indistinguishable correlation functions. The resulting time-dependent overlaps capture the same transient peak and long-time vibrational recurrences observed in classical simulations of standard quantum calculations of wavepackets on an LEPS potential energy surface.

From an algorithmic standpoint, our results show that time-dependent Møller correlation functions can be evaluated on a gate-based quantum computer using two complementary strategies. The modified Hadamard test directly encodes the complex overlap into the expectation value of a single ancilla qubit but requires controlled implementations of the full time evolution operator $[eqn]$ and of the state preparation unitaries. The MFE protocol removes these global controlled operations in favor of a sequence of SWAP tests and separate evolutions of differently prepared registers. The error analysis (Figure 11) demonstrates that the resulting correlation functions agree to machine precision, establishing MFE as a depth-reduced, ancilla-free alternative to the Hadamard test that is particularly attractive for near-term hardware with limited coherence times and noisy two-qubit gates.

From a chemical physics standpoint, the time traces of $[eqn]$ provide a compact dynamical fingerprint of the reaction [2,5,6]. The pronounced peak near $[eqn]$ in both channels corresponds to the regime where the incoming reactant and outgoing product wavepackets overlap in the interaction region, i.e., the scattering event itself. At later times, the correlation magnitude settles into weaker, quasi-periodic oscillations whose structure depends on the product channel. For the elastic channel $[eqn]$ , these oscillations are relatively strong and regular, consistent with coherent vibrational motion of the bound H-H product [35]. For the inelastic channel $[eqn]$ , the overall amplitude is reduced and the pattern more irregular, reflecting the smaller Franck–Condon overlap between the initial reactant packet and the first excited vibrational state of the product. Figure 12 shows that the quantum results reproduce all of these features, indicating that the quantum circuits not only match the classical simulation of standard quantum calculation benchmark numerically but also faithfully encode the underlying reaction dynamics.

More broadly, the present results position time-domain Møller correlation functions as a practical and chemically meaningful intermediate target for the quantum simulation of molecular dynamics. Rather than attempting full state tomography or the direct solution of the Schrödinger equation on a large grid, the algorithm focuses on physically relevant overlaps that can be accessed through relatively shallow circuits. Extending the present approach to higher-dimensional grids, more accurate three-dimensional H_3_ potential energy surfaces, and larger reactive systems where classical simulations of standard quantum calculations of wavepacket propagation become prohibitive is a natural next step [1,17].