Analog quantum simulation of coupled electron-nuclear dynamics in molecules

Abstract

Quantum computing has the potential to reduce the computational cost required for quantum dynamics simulations. However, existing quantum algorithms for coupled electron-nuclear dynamics simulation either require fault-tolerant devices, or involve the Born–Oppenheimer (BO) approximation and pre-calculation of electronic states on classical computers. We present the first quantum simulation approach for molecular vibronic dynamics in a pre-BO framework with an analog mapping of nuclear degrees of freedom, i.e. without the separation of electrons and nuclei, by mapping the molecular Hamiltonian to a device with coupled qubits and bosonic modes. We perform a proof-of-principle emulation of our ansatz using a single-mode model system which represents vibronic dynamics of chemical systems, such as nonadiabatic charge transfer involving polarization of the medium, and propose an implementation of our approach on a trapped-ion device. We show that our approach has exponential savings in resource and computational costs compared to the equivalent classical algorithms. Furthermore, our approach has a much smaller resource and implementation scaling than the existing pre-BO quantum algorithms for chemical dynamics. The low cost of our approach will enable an exact treatment of electron-nuclear dynamics on near-term quantum devices.

---

1 Introduction

Light–matter interactions are the source of many phenomena in molecular systems such as vision,^1–3 photosynthesis,^4–7 photovoltaics,^8–11 and photocatalysis.^12,13 In order to apply light–matter interactions to the development of high-functional molecular devices, it is important to understand their underlying mechanisms. Such interactions often involve nuclear quantum effects that can dramatically alter reaction dynamics relative to classical predictions. In principle, the molecular time-dependent Schrödinger equation can be exactly solved for systems of electrons and nuclei to simulate molecular dynamics. However, such a solution is practically impossible for molecules with more than a few atoms due to the exponential scaling of the computational cost with respect to system size.

The Born–Oppenheimer (BO) approximation neglects the coupling between electronic and nuclear degrees of freedom, known as nonadiabatic coupling (NAC).¹⁴ While it is useful for properties of the ground electronic state, it fails to describe light-induced chemical reactions where the molecule can reach strong NAC regions after photo-excitation of the molecular electronic state.¹⁴ The BO approximation can be extended to the group BO approximation (GBOA), which takes into account NACs between a group of BO states considered to be relevant to the dynamics, while NACs to states outside of the group are neglected.¹⁴ The molecular dynamics are then described as nuclear dynamics on multiple BO potential energy surfaces coupled via NACs. The majority of existing simulation methods make use of the GBOA since it extends the scope of molecular dynamics simulation while using conventional electronic structure methods.^15–17 Ideally, the couplings neglected by the GBOA are minimal, but choosing the number of BO electronic states requires some intuition. Furthermore, since existing simulation methods with the GBOA require an accurate calculation of BO states and their gradients (including NACs), the computational cost for accurate simulations can become intractable.

Pre-Born–Oppenheimer (pre-BO) methods, on the other hand, naturally include nonadiabatic effects since they treat nuclei and electrons without separation.^18–24 There have been many studies on pre-BO theory with both first quantization^18–21 and second quantization^22–24 of the molecular Hamiltonian. One approach using the pre-BO framework is the nuclear-electronic orbital (NEO) method, which expands certain nuclei in a basis of nuclear orbitals and those the nuclei as quantum particles.²² Recently, an approach was developed based on the multi-configurational time-dependent Hartree method using a second quantization representation of electronic degrees of freedom, which enables the description of coupled electron-nuclear dynamics without electronic potential energy surfaces.²⁴ Nevertheless, the development of methods in the pre-BO framework is still far behind the methods in the BO framework since the computational cost of a pre-BO treatment of molecular dynamics is much greater.

Quantum computing can significantly reduce the computational cost of the simulation of quantum mechanical systems by exploiting the intrinsic quantum nature of the computational device.^25–28 Most near-term quantum computing research in the field of quantum chemistry has focused on obtaining electronic properties at fixed nuclear configurations based on the variational quantum eigensolver (VQE) method^29–33 using the BO approximation, with proposed extensions to pre-BO eigenvalues using NEO.³⁴ For chemical dynamics, although there are several proposed quantum algorithms without the BO approximation for fault-tolerant devices,^35–37 most methods still adopt the GBOA for near-term applications with variational quantum algorithms^38–40 or analog mappings.^41–44 The intractability of quantum dynamics simulations with a pre-BO wavefunction on classical computers suggests that a pre-BO quantum simulation could show an earlier quantum advantage than methods using the BO framework.

In this work, we propose the first quantum simulation method for molecular vibronic dynamics in the pre-BO framework with an analog mapping of nuclear degrees of freedom. Our approach maps the second quantized pre-BO representation of the molecular vibronic Hamiltonian onto an analog quantum device, and thus treats the electron-nuclear interactions exactly. Specifically, we map the nuclear vibrational motions to bosonic modes of a device and use a fermion-qubit mapping for electronic degrees of freedom. We show that our method can efficiently simulate the exact molecular vibronic dynamics within a given single-particle basis set, and therefore it is suitable for accurate near-term quantum simulations of coupled electron-nuclear dynamics. As an example, we show how a single-mode vibronic model system of nonadiabatic charge transfer^45,46 could be simulated on currently available trapped ion quantum computers with existing experimental techniques and noise. Finally, we show how our approach has the potential to out-perform all previous approaches in terms of scaling of hardware resources, implementation cost, and classical pre-calculation.

2 Theory

Our simulation approach is restricted to vibronic (vibrational + electronic) internal degrees of freedom of the molecule, meaning the translational and rotational degrees of freedom of the atomic nuclear coordinates are removed. In general, this can be achieved by any unitary transformation of the Cartesian atomic coordinates that separates the 6 collective translations and rotations of the molecule, given by vectors Q_trans and Q_rot, from the 3N_at − 6 vibrational internal coordinates, Q ≡ Q_vib. For the remainder of the manuscript, we assume that mass-weighted normal mode coordinates are used, whereby the translations and rotations are easily identified as the zero-frequency modes. The electronic coordinates are defined as positions of electrons relative to the nuclear center of mass with a fixed orientation, i.e. in the Eckart frame.⁴⁷ The Coriolis coupling (between nuclear vibrations and rotations) is excluded. While the removal of molecular rotations is an approximation,⁴⁷ it is appropriate for ultrafast (fs–ps) chemical dynamics simulations, as is further justified in the Discussion.

The full molecular vibronic Hamiltonian is given by

where r_i is the position of the electron i in the Eckart frame, ∇_i = ∂/∂r_i is its gradient, Q_ν is the normal mode coordinate for mode ν, P_ν = −i∂/∂Q_ν is the corresponding momentum, N_mode is the number of vibrational modes, N_e is the number of electrons, and v_en, v_ee, and V_nn are electron-nuclear, electron–electron, and nuclear–nuclear interaction potentials, respectively. We use atomic units (ℏ = m_e = e = 4πε₀ = 1) here and throughout this paper.

To reduce the electronic basis size, we adopt a basis of orthonormal spin orbitals that depend parametrically on the nuclear positions, ϕ_p(x; Q) = φ_p(r; Q)σ_p(s), where x is a vector of electronic spatial coordinates r and a spin coordinate s, i.e. x = {r, s}, and φ_p and σ_p are spatial and spin functions of the spin orbital, respectively.

We express the molecular vibronic wavefunction using Slater determinants of the position-dependent spin orbitals for the electronic degrees of freedom, and Hartree products of harmonic oscillator eigenstates of each normal mode for the vibrational degrees of freedom. In the second quantization representation, the Slater determinants and the Hartree product can be expressed as occupation number vectors (ONVs). An electronic ONV for a Slater determinant is written as |n〉_e = |n₁⋯n_No〉_e, where n_j is the occupation of the j-th spin orbital and N_o is the number of spin orbitals, while a vibrational ONV for a Hartree product is given as |v〉_n = |v₁⋯v_Nmode〉_n, where v_ν is the occupation (Fock state) of mode ν. Therefore, our wavefunction ansatz becomes

where C_vn(t) is the coefficient for the collective occupation v ∪ n.

Sasmal and Vendrell showed that the ansatz above transforms the molecular vibronic Hamiltonian to the second quantization representation:²⁴

where h_pq, d_ν,pq, g_ν,pq, and v_pqrs are electron integrals whose definitions are given in the SI. The Q-dependence of spin-orbitals and Slater determinants are absorbed into the electronic integrals, and the couplings between nuclear and electronic degrees of freedom appear as products of the electron integrals and fermionic ladder operators in the Hamiltonian. The orbital vibronic coupling terms (depending on d_ν,pq and g_ν,pq) appear because the nuclear kinetic energy operator acts on the spin orbital basis functions.²⁴ The mass-weighted normal mode position and the conjugate momentum operators for mode ν can be expressed in terms of the bosonic ladder operators {b̂_ν} and {b̂^†_ν} where ω_ν is a reference harmonic frequency of normal mode ν. Therefore, the terms depending on the normal mode coordinates in eqn (3) can be expressed in terms of the bosonic ladder operators by using a Taylor series expansion about a reference geometry Q₀ up to a reasonable order.

For fixed Taylor expansion orders, ω_ν and Q₀ only affect the accuracy of our approach by the number of harmonic eigenstates required to describe the wavefunction. In other words, a quantum simulator with lower oscillator noise can tolerate lower-accuracy values of ω_ν and Q₀, including approximations derived from the electronic integral Taylor expansions. For practical cases in the near term, we expect that density functional theory will provide an ideal balance between accuracy and classical computational cost when finding Q₀ and ω_ν, and Kohn–Sham orbitals provide a compact spin orbital basis.⁴⁸ We can reduce the number of terms in eqn (3) using an active space. The active space defines a subset of spin orbitals for which all possible electron configurations are included, while the remaining (inactive) orbitals have fixed occupations of 0 or 1.⁴⁹ The terms corresponding to the inactive orbitals can be included in the electron integrals of active orbitals, and an external potential V_inact(Q̂) that depends only on the nuclear configuration (i.e. V_nn(Q̂) → V_nn(Q̂) + V_inact(Q̂)).

In general, the orbital vibronic coupling functions d_ν,pq and g_ν,pq can be extremely localized and non-analytical at some nuclear geometries. This requires a high-order expansion in terms of the bosonic ladder operators in eqn (3), which greatly increases the number of terms in the Hamiltonian. One way to reduce the expansion order would be a unitary transformation of orbitals for which the derivative couplings vanish, which must satisfy ∂c/∂Q_ν + d_νc = 0, where c is the transformation matrix to the resulting set of “diabatic” orbitals, . This transformation is identical in form to the many-body BO state diabatization condition.¹⁴ However, much like multi-electron states in multi-mode systems,⁵⁰ a set of strictly diabatic orbitals does not exist for an incomplete basis. Instead, because our goal is simply to achieve a low-order expansion in terms of Q, we can use “quasi-diabatic” orbitals where the above strict diabatization condition is not satisfied but all integral coefficients are smooth functions of Q. These orbitals can likewise borrow techniques from BO state diabatization.⁵¹

3 Analog quantum simulation

We now introduce our approach to map the Hamiltonian Ĥ_mol onto an analog quantum device to perform a pre-BO simulation. Such a mapping must satisfy the symmetry requirements of fermions (electrons) and bosons (vibrations), given by their (anti-)commutation relations. Ideally, we would like to map both types of degrees of freedom onto a device with native fermionic^52–54 and bosonic degrees of freedom. Recently, fermionic quantum processors have been proposed using existing tools such as optical tweezers and neutral fermionic atoms.⁵⁴ However, there is currently no quantum architecture which couples fermionic gates to continuous bosonic modes. Therefore, we will build on the wealth of literature for fermion-qubit mappings^55,56 to propose a digital-analog quantum simulation^57,58 technique. Our technique takes advantage of existing quantum architectures that have qubit levels and bosonic modes with controllable coupling between them, such as ion traps and circuit quantum electrodynamics (cQED). We will refer to such architectures as “coupled multi-qubit-boson” (cMQB) devices. We note that one could achieve the same type of simulation in a fully analog way if coupled fermionic-bosonic operations are developed in the future.

An illustration of our pre-BO analog simulation method in comparison with the BO framework is summarized in Fig. 1. We start by finding a nuclear position dependent spin orbital basis (Fig. 1a) and deriving the expansion coefficients in eqn (3). We then encode the pre-BO wavefunction (Fig. 1b) on a cMQB device, e.g. a trapped-ion device with ion electronic states representing electronic ONVs, and ion motional modes representing the nuclear component of the wavefunction (Fig. 1c). In contrast, existing classical/quantum algorithms in the BO framework require the pre-calculation of BO electronic states (Fig. 1d) with the spin orbital basis. The vibronic wavefunction is then propagated on a truncated BO basis (Fig. 1e).

Fig. 1 Illustration of the pre-BO cMQB analog quantum simulation approach for molecular vibronic dynamics using trapped-ion architectures in comparison with the BO framework. (a) The electronic wavefunction is given in terms of spin orbitals. Each diabatic spatial orbital, which depends on the nuclear positions, yields two spin orbitals with opposite spins: |ϕ₁〉 and |ϕ₃〉 (green), and |ϕ₂〉 and |ϕ₄〉 (yellow). Only the motion of a single ion (grey) is shown. (b) Molecular vibronic dynamics are simulated based on the pre-BO representation with interactions among electronic (cyan) and nuclear (grey) degrees of freedom. (c) In this work, we propose a coupled multi-qubit boson approach which maps occupation numbers {n_p} of spin orbitals {|ϕ_p〉} to qubits and nuclear motions to bosonic modes coupled to the qubits, and thus can perform pre-BO dynamics quantum simulations with a complete electronic basis within a given orbital basis set. (d) Existing quantum/classical algorithms rely on BO states {|Φ_i〉} with potential energy surfaces {|E_i〉} found by solving the electronic time-independent Schrödinger equation using the orbital basis at fixed nuclear positions prior to the simulation. (e) Molecular vibronic dynamics are solved in the BO framework by propagating the nuclear wavefunction (red and blue densities) on a truncated set of coupled BO states.

For the electronic degrees of freedom, we use the Jordan–Wigner transformation which directly maps the occupation number of a spin orbital to the qubit state.⁵⁵ However, other fermion-qubit mappings could be equivalently employed. The fermionic creation operator for spin orbital p is mapped to a tensor product of Pauli operators (X̂_k, Ŷ_k, Ẑ_k) for qubits by â_p = Ẑ₁ ⊗⋯⊗Ẑ_p−1 ⊗ (X̂_p + iŶ_p), and the annihilation operator is its Hermitian conjugate. The general expression for the cMQB Hamiltonian becomes a sum of tensor products of multiple bosonic and qubit operators, i.e.

where N_q is the number of qubits, P̂_k^I is one of the identity or Pauli operators for the k-th qubit (P̂_k ∈ {Î_k, X̂_k, Ŷ_k, Ẑ_k}) in the I-th Pauli string, and f_I is a function of bosonic ladder operators b̂_ν and b̂^†_ν coupled to the I-th Pauli string. Correspondingly, the electronic part of the molecular vibronic state |n〉_e in our ansatz (eqn (2)) is mapped to the multi-qubit state |q〉_q = |q₁⋯q_Nq〉_q of the device, while the multi-mode vibrational wavefunction is directly mapped to the bosonic degrees of freedom (Fig. 1c). The molecular cMQB Hamiltonian may be scaled with a factor F to a simulation scale (Ĥ_sim = FĤ_mol) which depends on the experimental parameters of the simulator.

Our electronic basis spans the entire electronic Fock space, including electron configurations with different number of electrons and spins. This implies that our approach can describe intersystem crossing if the corresponding coefficients derived from integrals of the spatial orbitals and electron angular momentum are non-zero. Our approach is also capable of describing ionization processes by adding interactions with terms that remove electrons, which has potential uses for predicting photoelectron spectra. In both cases, the number of quantum resources is unchanged for the simulation.

A schematic illustration of a circuit for pre-BO dynamics simulation with the cMQB mapping consists of initialization, time-evolution, and measurement, as shown in Fig. 2. We provide details for each step in what follows.

Fig. 2 A schematic quantum circuit diagram for pre-BO vibronic dynamics simulation with our cMQB approach. Straight lines represent qubits and wavy lines represent the motional modes. (a) A simulation consists of initialization, time-evolution, and measurement stages. An initial vibronic state can be prepared by a unitary operator Û_init. The initial state is then propagated by the time-evolution operator Û(t). Finally, an expectation value can be measured with Hadamard tests with the corresponding operator Ô controlled by an ancilla qubit state (represented by a black circle) using Hadamard (Ĥ) and phase gates (Ŝ). Real and imaginary parts of the expectation value are obtained with and without the phase gate, respectively. (b) Trotterization of the time-evolution operator Û(t) and the digital-analog circuit for an example Trotterized operator Û_I(Δt) = exp(−iθΔt(b̂₁ + b̂^†₁)Ŷ₂Ẑ₃Ŷ₄). CNOT gates are represented by symbols with a black circle for the control qubit, connected to an open circle for the target qubit.

As with any quantum dynamics simulation, our approach requires a robust preparation of the initial state which can be generated by a unitary operator Û_init (Fig. 2a). For most applications in photochemical dynamics we can assume that the initial state is well represented as a product state of electronic and nuclear components of the wavefunction which can be prepared by Û_init = Û_e ⊗ Û_n, with separate operations on qubits with Û_e and motional states with Û_n. According to the Franck–Condon approximation, the equilibrium nuclear state is unperturbed by photoexcitation and the excitation is instantaneous. This means the initial state can be prepared with a ground (or coherent) nuclear state, and with qubit states corresponding to the photoexcited spin orbital occupation.

Future devices with longer coherence times and more sophisticated quantum control will allow for increasingly accurate initial state preparation. For example, controlled dissipation of the simulator could be used to prepare the vibronic ground state wavefunction, with the fidelity approaching 1 for longer experimental times. The dynamics could then be initialized by a purely electronic (Franck–Condon) or vibronic transition dipole operator.

The direct implementation of our approach outlined above requires the simultaneous control of all degrees of freedom in the simulation. However, our approach uses a fermion-qubit mapping for the electronic degrees of freedom, which involves a change of basis between non-commuting electronic terms. Therefore, we must use Trotterization⁵⁹ of the time-evolution operator. Trotterization approximates the time-evolution operator as a product of operators, each of which corresponds to a term Ĥ_I in eqn (4). The Trotterized cMQB time-evolution operator is given by where N_op is the number of divided operators and N_t is the number of the Trotter steps (Fig. 2b).

The Trotterized cMQB time-evolution operator for each term with an electronic component in our Hamiltonian (eqn (4)) is given by

where [capital Lambda, Greek, circumflex]^I is an N-qubit Clifford operation which can be composed out of Hadamard, phase, and CNOT gates, to propagate the coupling between bosonic mode(s) and a single qubit q₀^I generated by a laser–ion interaction to multiple qubits for the desired Pauli strings for Ĥ_I. In Fig. 2b we show an example time-evolution operator Û_I(Δt) = exp(−iθΔt(b̂₁ + b̂^†₁)Ŷ₂Ẑ₃Ŷ₄), which can be achieved by a sequence of digital and analog quantum operations. The nuclear-only terms of the Hamiltonian (the final two terms of eqn (3)) are given up to second order by a constant detuning from the bosonic mode frequency. Higher-order nuclear-only terms can be achieved with operations on only the bosonic modes, or qubit-boson operations with an ancilla qubit. The largest coefficients in our Hamiltonian are roughly equal to the orbital energies, thus the maximum Trotter step size Δt is inversely proportional to the (active) orbital energy range to achieve the desired Trotter error.

This approach can be considered a digital-analog quantum simulation^57,58 because the gates involved in generating qubit entanglement take a digital form. Using a coupled fermionic-bosonic operation would eliminate the need for digital gates [capital Lambda, Greek, circumflex]^I for a fully analog simulation. Even with the fermion-qubit mapping, a more compact form of the time-evolution operator could be achieved using time- and spin-dependent squeezing/displacement operators on the collective ionic motions.⁶⁰ For the remainder of this paper we consider the digital-analog approach.

Trapped ions and cQED typically use bosonic modes to achieve entangling gates in digital quantum algorithms. The quantum gates necessary to implement [capital Lambda, Greek, circumflex]^I are thus readily achieved with these architectures. However, since the entangling gates require a bosonic mode, some bosonic modes (the “bus modes”) must be reserved from those used for simulating nuclear vibrational degrees of freedom. The total number of bosonic modes in the simulator must therefore be at least one greater than the total number of molecular vibrational modes.

Our approach encodes the electronic degrees of freedom of the vibronic wavefunction as an ONV, and nuclear degrees of freedom are mapped to the bosonic modes. This difference in encoding means that all nuclear observables can be measured directly from the simulator, whereas some electronic and vibronic observables will depend on the electronic orbital basis functions. In general, the expectation values of observables can be obtained with Hadamard tests⁶¹ for the corresponding unitary operator Ô as shown in Fig. 2a.

For example, a one-electron property 〈Ô_1e〉 ≡ 〈Ψ|Ô_1e|Ψ〉 is obtained by the sum of the measured elements of the one-electron reduced density matrix (1RDM), 〈â^†_pâ_q〉, weighted by the integrals 〈ϕ_p|Ô_1e|ϕ_q〉 for spin orbitals p and q. The 1RDM is calculated from the quantum simulator by measuring expectation values of Pauli strings, whereas the integrals are obtained on a classical computer. Conversely, the expectation value of a unitary nuclear operator that can be expressed as an exponential, Ô_n = exp(Â_n), can be implemented directly on the analog simulator with an effective Hamiltonian Ĥ^eff_n = iÂ_n/τ by evolving for a time τ.

Real-space density functions are important observables for characterizing the dynamics of a molecule. They are expressed in terms of the real-space projection operator, |x̲, Q〉〈x̲, Q| where x̲ = {x₁,…, x_Ne}. For example, the projection operator for the joint nuclear normal modes and electronic density is , which yields

where is a frequency weighted momentum coordinate for mode ν. The variable k = {k₁,…, k_Nmode} is the momentum space variable associated with the normal mode Q, and [scr F, script letter F]_k{f(k)}(Q) represents a multi-dimensional Fourier transform of the function f(k) to nuclear normal mode coordinate (Q) space. Eqn (6) exploits the tomography of the nuclear characteristic function using displacement operators D̂_ν,⁶² where the Fourier transform is performed numerically on the expectation values measured on grid points for displacement operators with different values of ξ_ν. The displacement operator in eqn (6) can be implemented on a cMQB device, and the range and resolution of the density is controlled by the separation of grid points ξ_ν.

Electron and nuclear densities can be obtained by integrating out other degrees of freedom, given by and respectively. Measurement of the nuclear density, which provides the basis for the analysis of chemical reactions such as isomerization and bond dissociation, simplifies to

due to the orthonormality of spin orbitals at all Q. Reduced vibrational densities are likewise found by only measuring the characteristic functions of the desired vibrational modes.

As an analog simulation approach, our approach is subject to environmental noise over the course of the simulation, which is a source of error in addition to the Trotter error. The noise effect on an open quantum system can be described by the Lindblad master equation:⁶³

where [small rho, Greek, circumflex] = |Ψ〉〈Ψ| represents the density operator of the system, Ĥ is the Hamiltonian of the system, is the Lindblad superoperator, L̂_i is the jump operator for the noise source i, and γ_i is the corresponding dissipation rate. The noise effect on a cMQB device translates into the molecular state during a cMQB simulation, where the native rates of the cMQB simulator are scaled to the molecular system depending on the scaling factor F (Ĥ_sim = FĤ_mol), encoding, device parameters, and simulation setup.

The noise does not equally affect electronic and vibrational degrees of freedom in cMQB simulations, since they are encoded in different resources of the simulator. In addition, because cMQB simulators use motional modes to achieve entangling gates, the electronic degrees of freedom are subject to noise on both qubit states and motional modes.

In ion-trap cMQB devices, the dominant source of noise is the decoherence of the motional modes, with a rate orders of magnitude greater than the other sources of noise.⁶⁴ Therefore, in what follows, we only consider the noise effect from motional decoherence (L̂_ν = b̂^†_νb̂_ν) in an ion-trap cMQB device. While the motional modes representing the nuclear degrees of freedom are subject to the decoherence for the entire experiment time, each qubit is subject to the motional noise only when an entanglement operator is applied on the qubit. Thus, the native rates for the motional decoherence and the indirect qubit noise scale differently for vibrational and electronic degrees of freedom of the molecular state. The native decoherence rate γ^d_mot of motional modes for a cMQB device is scaled to the molecular vibrational decoherence rate γ^vib_mol as

while the average indirect spin noise (L̂_q ≈ Ŝ_q) for each qubit in the molecular scale γ^q_mol can be connected to the native motional decoherence rate of the bus mode:where Δt_mol is the Trotter step size in the molecular scale, t_CNOT is the experimental time for a two-qubit entangling gate, and N_CNOT is the number of two-qubit entanglement gates per single Trotter step (see SI). Although a small Trotter step Δt_mol reduces the theoretical Trotter error of the cMQB time-evolution operator, the above equations imply that trade-off exists due to the number of digital qubit entanglement operations.

4 Connection to the Born–Oppenheimer framework

If we define the electronic Hamiltonian Ĥ^e(Q) as the first two terms of eqn (3) at a fixed nuclear position Q, then we can find an electronic Hamiltonian matrix H^e(Q) with elements H^e_mn(Q) = 〈m|Ĥ^e(Q)|n〉_e, where only electronic degrees of freedom are integrated. The Hamiltonian matrix can be transformed by a unitary matrix W(Q) to give E^e(Q) = W(Q)H^e(Q)W^†(Q). When E^e(Q) is a fully diagonal matrix, we arrive at the exact (full configuration interaction, FCI) solution to the electronic structure problem. In practice, full diagonalization is intractable on classical computers, so the eigenvalues of a submatrix of H^e(Q) may be found instead. For example, full diagonalization within an active space is known as complete active space configuration interaction (CASCI). Iterative diagonalization is typically used to further reduce the cost, and the number of eigenvalues found, N_BO, can be less than the rank of the (sub)matrix. The BO electronic states are the corresponding eigenfunctions (Fig. 1d), and molecular dynamics can be expressed by the Born–Huang expansion as nuclear dynamics on multiple BO states:where n(x̲; Q) = 〈x̲, Q|n〉 and v(Q) = 〈Q|v〉 are the real space wavefunctions of the electronic and nuclear ONVs, i.e. the Slater determinant and Hartree product of single particle basis functions, respectively. W_jn is an element of W, is the j-th BO state, and is the corresponding time-dependent BO-projected nuclear wavefunction with an expansion coefficient T_jv(t) (Fig. 1e).

Without spin–orbit coupling, the maximum number of BO electronic states is equal to the number of configuration state functions (CSFs, i.e. spin-adapted linear combinations of Slater determinants) with a fixed electron count and spin, N_CSF, which we will discuss in more detail in Section 6. The summation over the BO states in eqn (11) is exact within the given orbital basis set only when N_BO = N_CSF. The relation between the Born–Huang expansion and our ansatz (eqn (2)) becomes clear via the relation , where . However, N_CSF scales rapidly with the orbital basis set size whereas N_BO must be kept small for practical simulations. Therefore, the GBOA is employed (N_BO < N_CSF) to reduce the computational cost of the BO state calculations and time-propagation of the molecular vibronic state (Fig. 1e) based on prior knowledge on the system and its dynamics.

5 Numerical test

We show a proof-of-principle demonstration of our digital-analog simulation approach using the one-dimensional, two-electron Shin–Metiu model^45,46 and its implementation on a trapped ion quantum simulator, which has been used for the quantum simulation of chemical dynamics for several systems.^41,43,44,65 The model consists of two electrons, two fixed ions, and one moving ion between the fixed ions displaced by a distance L = 5.4 a.u. (Fig. 3a), where the origin is set to the middle point of the two fixed ions. The model is an extension of the original model system inspired by nonadiabatic charge transfer in processes with polarization of the medium induced by electron movement,^45,66 such as proton-coupled electron transfer. The computational details for all simulations performed in this section, including the parameters for the model Hamiltonian, are given in the SI.

Fig. 3 Results for theoretical simulations of the one-dimensional, two-electron Shin–Metiu model with exact time-evolution, our cMQB approach, and the GBOA. (a) The model system consists of two fixed ions (white circles) displaced by distance L, a moving ion (grey circle) at a position R, and two electrons (green circles) at positions r₁ and r₂. The origin is the average of the fixed ion positions. The diabatic orbitals are (b) η_a(r; R) localized around the left fixed ion, and (c) η_b(r; R) localized around the right fixed ion. The density functions ρ(r, R, t) at, (d) t = 0 a.u., (e)–(g) t = 56.1 a.u., and (h)–(j) t = 1514.4 a.u., with exact time-evolution ((e) and (h)), Trotterized cMQB time-evolution (Δt = 5.6 a.u.) ((f) and (i)), and GBOA ((g) and (j)). Spatial functions in (b)–(j) are normalized to their maximum values.

We construct the electronic Fock space with four spin orbitals {|ϕ_i〉} using the spin-up and spin-down configurations of two diabatic spatial orbitals, η_a and η_b in Fig. 3b and c respectively, i.e. φ₁ = φ₃ = η_a, and φ₂ = φ₄ = η_b. As a result, a total of four ions are needed in the experimental trapped ion setup for the simulation, with only two of the 12 available motional modes needed for the simulation. The diabatic orbitals were obtained by numerically integrating two adiabatic orbitals from the one-electron Shin–Metiu Hamiltonian⁴⁵ subject to the diabatization condition. The resulting orbitals are localized around the left and right fixed ions, but delocalize slightly as the moving ion approaches the fixed ion (Fig. 3b and c). We replaced the Coulomb potential for V_nn with a harmonic potential to simplify the Hamiltonian and to confine the spatial extent of the wavefunction. The one- and two-electron integrals are truncated at first order in the nuclear position: v(Q̂) ≈ v₀ + v₁Q̂, where with ladder operators b̂^† and b̂ of our single mode with frequency ω of the harmonic nuclear repulsion potential. With these approximations, the cMQB Hamiltonian can be written as

where V_k^I represents the effective k-th order coefficient of the nuclear position coupled to the I-th Pauli string. The model system has only a single nuclear degree of freedom, meaning the normal mode coordinate is equal to the position of the moving ion multiplied by the square root of the mass of the ion M, i.e. . We show results in terms of the unweighted coordinate R rather than Q to show the ion and electrons in the same position space.

We consider the initial molecular state (Fig. 3d), where the ground state wavefunction of the harmonic oscillator for bosonic degrees of freedom is displaced by R₀ = 0.1 a.u. in real space, and electrons are in the closed shell configuration of the orbital η_a (Fig. 3b). To achieve the initial state experimentally on an ion trap, we would first prepare the ground state of two bosonic modes and of each of the qubits, |0000〉_q. Then, we would apply a digital quantum operation X̂₁X̂₃, and subsequently apply the displacement operator using a spin-dependent force on the ions with a laser.⁶⁷ During the simulation, an additional bosonic mode is required to implement qubit entanglement operations.

For the time-evolution, we assume that the base Hamiltonian of the ion trap (the first term in eqn (12)) is always present during the simulation and rescale the Trotter steps to compensate for it.⁴¹ To achieve the digital-analog time evolution, we would apply a series of digital quantum gates on the trapped ions with single-ion addressing. The second bosonic mode would be used to apply the entangling gates.⁶⁷

In Fig. 3, we show the density functions ρ(r, R, t) obtained from exact time-evolution, Trotterized cMQB time-evolution, and the GBOA. Experimental determination of the density function would require an ancilla qubit which controls the operator for the density measurement. Controlled Pauli strings correspond to CNOT and single-qubit gates, and controlled displacement operators can be achieved with a spin-dependent force on the trapped ions.^62,67 The sum of all orbital pairs with Fourier transforms over a grid of displacements yields the density, as previously shown in eqn (6). For this example, only 12 sets of measurements are required: two Pauli strings with four on-diagonal and two non-zero off-diagonal pairs. This number of measurements can be halved because the densities of the two spin components are equal. We can obtain a nuclear density resolution of 0.02 a.u. using a grid spacing of 1.26 a.u. in the momentum space with 250 points.

N_CSF = 3 is the size of complete electronic basis of singlet electronic states with two electrons in two spatial orbitals. We have excluded the ground BO state from the dynamics with the GBOA (N_BO = 2), which has relatively small NACs with the first excited state (see SI for densities of electronic states in the BO framework). Within a short simulation time (56.1 a.u.), the GBOA density shows a spurious electron transfer from the left fixed ion to the right (Fig. 3g) where the subsequent density deviates even further (Fig. 3j). The electron transfer at t = 56.1 a.u. does not occur in the exact simulation (Fig. 3e) due to vibronic coupling to the ground BO state. The inaccuracy of the GBOA in this example is much greater than one would expect for realistic molecules, due to the minimal nature of the model. On the other hand, the Trotterized cMQB time-evolution reproduces the exact density with a converged Trotter step of Δt = 5.6 a.u (Fig. 3f and i) with the fidelity |〈Ψ|Ψ_exact〉|² > 0.95, where the convergence of fidelity with respect to the Trotter step is reported in the SI. A more detailed discussion on the comparison of our approach with the BO framework is given in Discussion.

We perform open quantum system dynamics simulations using eqn (8) with Trotterization to estimate the noise effect on the dynamics of our model system, using a Rabi rate for ion–laser interaction of Ω = 2π × 1.0 MHz (ref. 68) and the native vibrational decoherence rate of γ^d_mot = 30 s⁻¹ (ref. 64) based on existing ion trap cMQB devices. We report the time evolution of the fractional occupation numbers (FONs) of ϕ₁ and the fidelities in Fig. 4 with vibrational decoherence and the indirect qubit noise effect on the molecular dynamics from the motional decoherence of an ion-trap, in comparison with the result without noise. The indirect qubit noise effect is negligible compared to the direct vibrational decoherence effect on nuclear degrees of freedom in FON, while the fidelity is more sensitive (see SI).

Fig. 4 Comparison of errors originating from noise effects due to motional decoherence and Trotterization, with different Trotter step sizes for cMQB simulations. (a) Time-evolution of the fractional occupation number of ϕ₁. (b) Time-evolution of the fidelity of the molecular wavefunctions: |〈Ψ|Ψ_exact〉|². The blue, red, and green lines represent the noise simulation results with Δt = 5.6, 9.6, and 16.8 a.u., respectively, where the dashed lines represent the corresponding closed system simulation results and the black line represents the exact closed system simulation result without Trotterization.

Although the noise effect leads to a deviation from the exact dynamics, the FONs are qualitatively well reproduced (Fig. 4a). Furthermore, the actual noise effect can be minimized by circuit optimization and experimental techniques to reduce the experiment time required for cMQB simulations. For example, the cMQB Hamiltonian of our model system consists of many two-qubit Pauli strings without bosonic terms, where their time-evolution operators can each be realized by a single native two-qubit entanglement operator (e.g. using ) instead of nesting a single qubit exponential operator with CNOT gates.

Fig. 4b confirms the trade-off between the Trotter error and the noise effect. The fidelity with Δt = 9.6 a.u. is greater than the fidelity with Δt = 16.8 a.u. due to the smaller Trotter error. However, the fidelity decreases again when the Trotter step is further reduced to Δt = 5.6 a.u. due to the contribution from the digital part in the molecular scaling factor (eqn (9) and (10)) for the noise. Therefore, the Trotter step should be chosen carefully considering both the Trotter error and the noise.

6 Discussion

We first compare the resource and computational costs of our pre-BO cMQB simulation approach with those of existing classical and quantum algorithms, where the advantages of our approach are highlighted in Fig. 5. Then, we discuss considerations for practical applications including noise effects and Trotter errors, which depend on the device parameters of cMQB simulators.

Fig. 5 Resource and computational cost comparison of our cMQB approach with other classical and quantum algorithms. (a) Resource scaling of electronic degrees of freedom with respect to N_o in a logarithmic scale for classical (log₂(N_CSF)) and quantum pre-BO simulations (log₂(N_q)) with singlet electronic states. The resource scaling of our method is represented by the red line, wheras those of equivalent classical simulations is represented by the blue lines for the maximum (dashed), minimum (dotted), and N_o/N_e = 6 case (solid). (b) Scaling in the number of interactions/gates required to implement the time-evolution operator (log(N_op)) with respect to N_o for our approach (red), and existing quantum algorithms in the BO framework (yellow) labeled as “qBO” for the maximum (dashed), minimum (dotted), and N_o/N_e = 6 case (solid). (c) The number of qubits (N_q) required vs. the number of electrons N_e for the pre-BO quantum algorithm proposed in ref. 35 using grids where N_q = 3nN_e with n = 10 (black) and our approach with N_q = 6N_e (red).

The number of qubits required for our approach is equal to the number of spin orbitals, N_o. On the other hand, the equivalent number of electronic states in the BO framework corresponds to the exact Born–Huang expansion limit (N_BO = N_CSF). Even for cases where N_BO can be kept small, the resources required for the classical calculation of the BO states equivalent to our approach (i.e. FCI, or CASCI for a subset of orbitals) scales with N_CSF. Therefore, the classical resources of the BO framework equivalent to our approach scale proportionally to N_CSF for a fixed spin-multiplicity, which can be calculated according to the Weyl's formula for the corresponding total spin angular momentum S:

Here, we discuss singlet states (S = 0). Most ground-state electron configurations of neutral organic molecules are singlets since they have an even number of electrons and a lack of orbital degeneracy. In general, the classical resources scale as O(N_o⁻²(m(m − 1)^1/m−1)^N_o) for N_o = mN_e (m > 1) at the large N_o limit. The number of singlet CSFs is minimum (O(N_o²)) when N_e = 2 or N_e = N_o − 2, and maximum (O(2^N_o/N_o²)) when 2N_e = N_o for a fixed N_o. Practical simulations fall in between these two ranges with a number of orbitals roughly proportional to the number of electrons, with the ratio depending on the basis set. For example, conjugated polyene chains with a double-zeta basis set have N_o/N_e = 6. In Fig. 5a we show classical resource scaling for singlet states in comparison with our approach. The linear scaling of our approach has a clear advantage over classical methods, even over the minimum N_CSF case for N_o > 6. Other electronic structure methods may use a different wavefunction ansatz with lower resource costs, but they are not equivalent to our approach, making their comparison difficult due to the many factors involved. The classical simulation has an exponential resource scaling for nuclear degrees of freedom of where N_ν,bas is the size of basis for mode ν, whereas our encoding has a linear scaling O(N_mode).

The simulation cost of our approach, which can be measured by the number of interaction terms required, scales as O(N_o⁴N_mode^k) where k is the maximum order of the Taylor expansion for the functions of bosonic degrees of freedom f_I. This cost also remains far smaller than the exponential cost of quantum dynamics simulations on classical computers:

Previous quantum algorithms, including the MQB approach, can potentially have a smaller quantum resource scaling for electronic degrees of freedom than our method: log₂(N_CSF) < N_o for the complete CSF basis^38,39,41 since the CSF space is a subspace of the spin orbital Fock space. However, the simulation cost N_op would scale as O(N_CSF²) for electronic degrees of freedom in existing quantum simulation approaches, whereas the cMQB mapping scales as O(N_o⁴) (Fig. 5b). For the nuclear degrees of freedom, our approach has an advantage over other digital or hybrid quantum algorithms for molecular quantum dynamics^38,39 since our approach has the same resource scaling as the MQB mapping⁴¹ due to the direct mapping on the bosonic modes.

The GBOA reduces N_op to O(N_BO²), and its use is justified where a limited number of BO electronic states are accessible at different nuclear geometries. Nevertheless, even if the number of BO states can be reduced, the BO states may be obtained approximately since an exact pre-calculation would require fully diagonalizing the FCI or CASCI Hamiltonian on a classical computer at a prohibitive cost (O(N_CSF³)), whereas the cMQB approach requires only the pre-calculation of electron integrals that scales as O(N_o⁴). Highly accurate electronic structure methods have a cost scaling greater than our approach (O(N_o^6–8)),⁶⁹ whereas methods with a lower accuracy require careful characterization. Furthermore, the scalar coupling terms are impractical to calculate for most electronic structure methods and are often neglected, whereas the analogous terms are easily included in our pre-BO approach.

The first proposal of a quantum algorithm for pre-BO molecular dynamics used a real-space, first quantization representation which has a linear resource scaling with respect to the number of grid points per each degree of freedom: 3n(N_n + N_e), where n is the number of qubits required for each degree of freedom and N_n represents the number of nuclei. The authors suggested a minimum grid size of n = 10.³⁵ Although their approach also scales linearly with the number of electrons, our approach has a lower resource cost for most molecules with a reasonable basis set size, where N_mode < N_e and N_o < 30N_e (Fig. 5c).

The noise effect on the cMQB device depends primarily on the Trotter step size and the number of digital gates used to generate multi-qubit entanglement (eqn. (9) and (10)), which has a trade-off with the Trotter error. Relative to the direct effect of motional decoherence on modes representing vibrations, the indirect qubit noise is negligible. The effect of noise is highly specific to the target problem and device; however, the scaling factors for the noise effect suggest desirable improvements in both computational and experimental aspects of the cMQB approach for a near-term quantum advantage: reducing N_CNOT via circuit optimization and reducing t_CNOT in cMQB devices.

From a different perspective, noise can be leveraged as an advantage in analog simulators by enabling the simulation of open quantum systems at minimal cost as demonstrated in previous work,^41,64,70,71 rather than treating it as an effect that must be suppressed. For example, our model problem can be interpreted as a nonadiabatic charge transfer process occurring in solution under the influence of solvent interactions, or within a crystal at an interface with a substrate.

Considering the quantum resource requirements and noise effects, our approach can be implemented on current and near-term hardware. Early applications are well suited to molecules with a small active space, particularly rigid molecules whose nuclear dependence is well-described by low-order Taylor expansions. These molecules are also the target of efficient classical algorithms, but our approach has advantages for the accuracy by including all terms in the vibronic Hamiltonian for the choice of orbital basis. However, practical advantage should be examined carefully by assessing the achievable efficiency and accuracy in the BO framework and the effect of errors in near-term cMQB devices.

For long term applications, we note that the range of orbital energies in a molecule are typically much greater than the range of energetically accessible BO electronic states. The Trotter step size Δt_mol decreases proportional to the orbital energy range. This in turn decreases the scaling factor F due to the minimum experimental time step, leading to longer simulation times and thus proportionally greater noise. To reduce the energy range, we can choose a set of active orbitals at the cost of reducing simulation accuracy. Because our approach is more strongly limited by the orbital energies rather than the number of orbitals, we expect it to accommodate larger active spaces than classical computing approaches. For this reason, our approach will be particularly advantageous for strongly correlated systems, which can have large numbers of nearly degenerate orbitals, leading to extremely large numbers of CSFs and intractable BO-based calculations.

Nuclear quantum effects, such as proton tunneling, are naturally included in our approach. Such effects are essential for describing many chemical processes, especially in biological systems.⁷² Our choice of model system represents two static moieties (fixed ions) coupled to two electrons and a moving ion (e.g. a proton). The time evolution resulted in a superposition of the moving ion between left and right positions, demonstrating the importance of the quantum description. Furthermore, the results for the model system demonstrate how small changes in the electronic description can result in large differences in the dynamics. These results demonstrate that our approach will enable accurate and efficient simulations of the quantum effects in photochemical dynamics and proton-coupled electron transfer.

A key approximation in our approach is the exclusion of the Coriolis (rotational-vibrational) coupling from the Hamiltonian. On the timescale of ultrafast processes (fs–ps), we expect that the effect of Coriolis coupling would be a negligible perturbation on the vibronic states since the rotational frequencies are relatively small. This approximation also comes with an advantage for simulating physical observables. Any observable quantity that depends on the molecular (Eckart) frame is averaged out over the rotational wavefunction, including electronic and nuclear densities. Experimental measurement of these observables in a lab frame results from spontaneous symmetry breaking between the molecule and its environment.⁷³ This gives a potential advantage for interpretability of our wavefunction compared to other quantum simulation proposals.^35,37 Future work could explore how rotations could be included as an environmental effect without compromising the wavefunction interpretability.

As previously mentioned, our approach can describe intersystem crossing by considering spin–orbit coupling electron integrals and ionization by including additional terms in the Hamiltonian (eqn (3)) and incorporating the corresponding operations in the cMQB simulator without any additional quantum resources. Such simulations have the potential to make our approach even more advantageous since the classical and quantum simulations in the BO framework require a greater number of basis states and the corresponding electronic structure calculations for different spin multiplicities. Based on the form of the vibronic Hamiltonian, our Hamiltonian can be considered as an extended version of electron–phonon coupling models used in the solid state physics field such as the Hubbard–Holstein model.⁷⁴ Therefore, our approach can be readily extended to dynamics involving intersystem crossing and ionization, and translated onto electron-phonon coupling dynamics in solids.

7 Conclusions

In this paper, we proposed a quantum simulation approach with an analog mapping of nuclear degrees of freedom to simulate the quantum dynamics of vibrations and electrons in molecules. We presented an emulation of the digital-analog implementation of our approach on a trapped-ion device, using a simple vibronic dynamics model system resembling realistic chemical processes. Our approach uses a pre-BO wavefunction ansatz, which converges to the exact solution of the non-relativistic time-dependent Schrödinger equation with increasing electronic basis set size and nuclear Taylor expansion order. Furthermore, our approach can be extended to describe intersystem crossing and ionization without additional resources, and can be applied to quantum simulations of electron–phonon coupling models for solids. In contrast, most previously proposed approaches for the simulation of molecular dynamics employ the BO framework, which often involves truncation in the electronic basis in dynamics simulation and/or electronic structure calculations for the BO states and couplings. Others use a real-space, first quantization approach, which includes molecular rotations that complicate the interpretability of the wavefunction.

For equivalent descriptions of vibronic dynamics, our approach shows clear scaling advantages over existing quantum algorithms when the accuracy, quantum resources, number of operations, and number of pre-calculations are considered. However, true advantages of our approach for realistic molecules require careful consideration of potential sources of error such as Trotterization, noise, and Taylor expansion truncation, as well as comparison to classical computing approximations. Model systems such as the example presented in this work can be readily implemented on existing quantum hardware such as trapped ion quantum computers with realistic experimental noise. With improvements in quantum hardware and quantum control, we expect that our approach will demonstrate an early quantum advantage for the simulation of quantum chemistry.

Author contributions

Conceptualization: J.-K. H. and R. J. M. method actualization: J.-K. H. simulations: J.-K. H. writing: J.-K. H. and R. J. M. visualization: J.-K. H. and R. J. M.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data supporting this article have been included as part of the SI. Supplementary information: detailed equations, computational details, and supplementary simulation results. See DOI: https://doi.org/10.1039/d5sc04076k.

Acknowledgements

We would like to thank Ting Rei Tan, Ivan Kassal, Seung Kyu Min, and Michael Schuurman for insightful discussion. J.-K. H. and R. J. M. were supported by start-up funding from Dalhousie University, and by National Research Council Canada through the Applied Quantum Challenge program (AQC-100). J.-K. H. was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (RS-2023-00237886).

References

1. D. Polli, P. Altoé, O. Weingart, K. M. Spillane, C. Manzoni, D. Brida, G. Tomasello, G. Orlandi, P. Kukura, R. A. Mathies, M. Garavelli and G. Cerullo, Nature, 2010, 467, 440–443.
2. S. Gozem, H. L. Luk, I. Schapiro and M. Olivucci, Chem. Rev., 2017, 117, 13502–13565.
3. C. Schnedermann, X. Yang, M. Liebel, K. M. Spillane, J. Lugtenburg, I. Fernandez, A. Valentini, I. Schapiro, M. Olivucci, P. Kukura and R. A. Mathies, Nat. Chem., 2018, 10, 449.
4. G. D. Scholes, G. R. Fleming, L. X. Chen, A. Aspuru-Guzik, A. Buchleitner, D. F. Coker, G. S. Engel, R. van Grondelle, A. Ishizaki, D. M. Jonas, J. S. Lundeen, J. K. McCusker, S. Mukamel, J. P. Ogilvie, A. Olaya-Castro, M. A. Ratner, F. C. Spano, K. B. Whaley and X. Zhu, Nature, 2017, 543, 647.
5. E. Romero, V. I. Novoderezhkin and R. van Grondelle, Nature, 2017, 543, 355.
6. M. Kaucikas, K. Maghlaoui, J. Barber, T. Renger and J. J. Van Thor, Nat. Commun., 2016, 7, 13977.
7. D. J. Nürnberg, J. Morton, S. Santabarbara, A. Telfer, P. Joliot, L. A. Antonaru, A. V. Ruban, T. Cardona, E. Krausz, A. Boussac, A. Fantuzzi and A. W. Rutherford, Science, 2018, 360, 1210–1213.
8. J.-C. Blancon, H. Tsai, W. Nie, C. C. Stoumpos, L. Pedesseau, C. Katan, M. Kepenekian, C. M. M. Soe, K. Appavoo, M. Y. Sfeir, S. Tretiak, P. M. Ajayan, M. G. Kanatzidis, J. Even, J. J. Crochet and A. D. Mohite, Science, 2017, 355, 1288–1292.
9. S. Nah, B. Spokoyny, C. Stoumpos, C. Soe, M. Kanatzidis and E. Harel, Nat. Photonics, 2017, 11, 285.
10. J. Huang, Y. Yuan, Y. Shao and Y. Yan, Nat. Rev. Mater., 2017, 2, 17042.
11. L. Qiao, W.-H. Fang, R. Long and O. V. Prezhdo, J. Am. Chem. Soc., 2021, 143, 9982–9990.
12. Z. Wang, C. Li and K. Domen, Chem. Soc. Rev., 2019, 48, 2109–2125.
13. V. K. Singh, C. Yu, S. Badgujar, Y. Kim, Y. Kwon, D. Kim, J. Lee, T. Akhter, G. Thangavel, L. S. Park, J. Lee, P. C. Nandajan, R. Wannemacher, B. n. Milián-Medina, L. Lüer, K. S. Kim, J. Gierschner and M. S. Kwon, Nat. Catal., 2018, 1, 794–804.
14. G. A. Worth and L. S. Cederbaum, Annu. Rev. Phys. Chem., 2004, 55, 127–158.
15. H.-D. Meyer, U. Manthe and L. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78.
16. R. Crespo-Otero and M. Barbatti, Chem. Rev., 2018, 118, 7026–7068.
17. M. Ben-Nun, J. Quenneville and T. J. Martínez, J. Chem. Phys. A, 2000, 104, 5161–5175.
18. E. Mátyus, Mol. Phys., 2019, 117, 590–609.
19. S. Bubin, M. Pavanello, W.-C. Tung, K. L. Sharkey and L. Adamowicz, Chem. Rev., 2013, 113, 36–79.
20. M. Cafiero, S. Bubin and L. Adamowicz, Phys. Chem. Chem. Phys., 2003, 5, 1491–1501.
21. T. Kato and K. Yamanouchi, J. Chem. Phys., 2009, 131, 164118.
22. S. Hammes-Schiffer, J. Chem. Phys., 2021, 155, 030901.
23. M. Sibaev, I. Polyak, F. R. Manby and P. J. Knowles, J. Chem. Phys., 2020, 153, 124102.
24. S. Sasmal and O. Vendrell, J. Chem. Phys., 2020, 153, 154110.
25. S. Lloyd, Science, 1996, 273, 1073–1078.
26. S. Wiesner, arXiv, 1996, preprint, arXiv:quant-ph/9603028,  DOI:10.48550/arXiv.quant-ph/9603028.
27. C. Zalka, Fortschr. Phys., 1998, 46, 877–879.
28. Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis and A. Aspuru-Guzik, Chem. Rev., 2019, 119, 10856–10915.
29. A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik and J. L. O'Brien, Nat. Commun., 2014, 5, 4213.
30. K. M. Nakanishi, K. Mitarai and K. Fujii, Phys. Rev. Res., 2019, 1, 033062.
31. P. J. Ollitrault, A. Kandala, C.-F. Chen, P. K. Barkoutsos, A. Mezzacapo, M. Pistoia, S. Sheldon, S. Woerner, J. M. Gambetta and I. Tavernelli, Phys. Rev. Res., 2020, 2, 043140.
32. S. Yalouz, E. Koridon, B. Senjean, B. Lasorne, F. Buda and L. Visscher, J. Chem. Theory Comput., 2022, 18, 776–794.
33. S. Tamiya, S. Koh and Y. O. Nakagawa, Phys. Rev. Res., 2021, 3, 023244.
34. A. Kovyrshin, M. Skogh, A. Broo, S. Mensa, E. Sahin, J. Crain and I. Tavernelli, J. Chem. Phys., 2023, 158, 214119.
35. I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni and A. Aspuru-Guzik, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 18681–18686.
36. I. D. Kivlichan, N. Wiebe, R. Babbush and A. Aspuru-Guzik, J. Phys. A: Math. Theor., 2017, 50, 305301.
37. Y. Su, D. W. Berry, N. Wiebe, N. Rubin and R. Babbush, PRX Quantum, 2021, 2, 040332.
38. P. J. Ollitrault, G. Mazzola and I. Tavernelli, Phys. Rev. Lett., 2020, 125, 260511.
39. P. J. Ollitrault, A. Miessen and I. Tavernelli, Acc. Chem. Res., 2021, 54, 4229–4238.
40. D. Bultrini and O. Vendrell, Commun. Phys., 2023, 6, 328.
41. R. J. MacDonell, C. E. Dickerson, C. J. T. Birch, A. Kumar, C. L. Edmunds, M. J. Biercuk, C. Hempel and I. Kassal, Chem. Sci., 2021, 12, 9794–9805.
42. C. S. Wang, N. E. Frattini, B. J. Chapman, S. Puri, S. M. Girvin, M. H. Devoret and R. J. Schoelkopf, Phys. Rev. X, 2023, 13, 011008.
43. R. J. MacDonell, T. Navickas, T. F. Wohlers-Reichel, C. H. Valahu, A. D. Rao, M. J. Millican, M. A. Currington, M. J. Biercuk, T. R. Tan, C. Hempel and I. Kassal, Chem. Sci., 2023, 14, 9439–9451.
44. C. H. Valahu, V. C. Olaya-Agudelo, R. J. MacDonell, T. Navickas, A. D. Rao, M. J. Millican, J. B. Pérez-Sánchez, J. Yuen-Zhou, M. J. Biercuk, C. Hempel, T. R. Tan and I. Kassal, Nat. Chem., 2023, 15, 1503–1508.
45. S. Shin and H. Metiu, J. Chem. Phys., 1995, 102, 9285–9295.
46. Y. Suzuki and K. Yamashita, Chem. Phys. Lett., 2012, 531, 216–222.
47. P. Bunker and P. Jensen, Molecular Symmetry and Spectroscopy, NRC Research Press, Ottawa, 2nd edn, 2006.
48. J. Kim, K. Hong, S. Choi, S.-Y. Hwanga and W. Y. Kim, Phys. Chem. Chem. Phys., 2015, 17, 31434–31443.
49. B. O. Roos, P. R. Taylor and P. E. Sigbahn, Chem. Phys., 1980, 48, 157–173.
50. C. A. Mead and D. G. Truhlar, J. Chem. Phys., 1982, 77, 6090–6098.
51. D. R. Yarkony, C. Xie, X. Zhu, Y. Wang, C. L. Malbon and H. Guo, Comput. Theor. Chem., 2019, 1152, 41–52.
52. J. Argüello-Luengo, A. González-Tudela, T. Shi, P. Zoller and J. I. Cirac, Nature, 2019, 574, 215–218.
53. T. Hartke, B. Oreg, N. Jia and M. Zwierlein, Nature, 2022, 601, 537–541.
54. D. González-Cuadra, D. Bluvstein, M. Kalinowski, R. Kaubruegger, N. Maskara, P. Naldesi, T. V. Zache, A. M. Kaufman, M. D. Lukin, H. Pichler, B. Vermersch, J. Ye and P. Zoller, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2304294120.
55. P. Jordan and E. Wigner, Z. Phys., 1928, 47, 631–651.
56. S. B. Bravyi and A. Y. Kitaev, Ann. Phys., 2002, 298, 210–226.
57. L. Lamata, A. Mezzacapo, J. Casanova and E. Solano, EPJ Quantum Technol., 2014, 1, 9.
58. S. Kumar, N. N. Hegade, E. Solano, F. Albarrán-Arriagada and G. A. Barrios, npj Quantum Inf., 2025, 11, 43.
59. N. Hatano and M. Suzuki, in Quantum Annealing and Other Optimization Methods, ed. A. Das and B. K. Chakrabarti, Springer Berlin Heidelberg, Berlin, Heidelberg, 2005, pp. 37–68.
60. O. Katz, M. Cetina and C. Monroe, PRX Quantum, 2023, 4, 030311.
61. R. Cleve, A. Ekert, C. Macchiavello and M. Mosca, Proc. R. Soc. A, 1998, 454, 339–354.
62. C. Flühmann and J. P. Home, Phys. Rev. Lett., 2020, 125, 043602.
63. G. Lindblad, Commun. Math. Phys., 1976, 48, 119.
64. V. C. Olaya-Agudelo, B. Stewart, C. H. Valahu, R. J. MacDonell, M. J. Millican, V. G. Matsos, F. Scuccimarra, T. R. Tan and I. Kassal, Phys. Rev. Res., 2025, 7, 023215.
65. M. Kang, H. Nuomin, S. Chowdhury, J. Yuly, K. Sun, J. Whitlow, J. Valdiviezo, Z. Zhang, P. Zhang, D. Beratan and K. Brown, Nat. Rev. Chem., 2024, 8, 340–358.
66. N. P. Blake, V. Srdanov, G. D. Stucky and H. Metiu, J. Phys. Chem., 1995, 99, 2127–2133.
67. K. Mølmer and A. Sørensen, Phys. Rev. Lett., 1999, 82, 1835–1838.
68. T. R. Tan, J. P. Gaebler, R. Bowler, Y. Lin, J. D. Jost, D. Leibfried and D. J. Wineland, Phys. Rev. Lett., 2013, 110, 263002.
69. D. A. Matthews and J. F. Stanton, J. Chem. Phys., 2016, 145, 124102.
70. A. Lemmer, C. Cormick, D. Tamascelli, T. Schaetz, S. F. Huelga and M. B. Plenio, New J. Phys., 2018, 20, 073002.
71. K. Sun, M. Kang, H. Nuomin, G. Schwartz, D. N. Beratan, K. R. Brown and J. Kim, Nat. Commun., 2025, 16, 4042.
72. J. P. Klinman and A. Kohen, Annu. Rev. Biochem., 2013, 82, 471–496.
73. E. Mátyus and P. Cassam-Chenaï, J. Chem. Phys., 2021, 154, 024114.
74. E. Berger, P. Valášek and W. von der Linden, Phys. Rev. B:Condens. Matter Mater. Phys., 1995, 52, 4806–4814.

---