Analog quantum simulation of chemical dynamics

Abstract

Ultrafast chemical reactions are difficult to simulate because they involve entangled, many-body wavefunctions whose computational complexity grows rapidly with molecular size. In photochemistry, the breakdown of the Born–Oppenheimer approximation further complicates the problem by entangling nuclear and electronic degrees of freedom. Here, we show that analog quantum simulators can efficiently simulate molecular dynamics using commonly available bosonic modes to represent molecular vibrations. Our approach can be implemented in any device with a qudit controllably coupled to bosonic oscillators and with quantum hardware resources that scale linearly with molecular size, and offers significant resource savings compared to digital quantum simulation algorithms. Advantages of our approach include a time resolution orders of magnitude better than ultrafast spectroscopy, the ability to simulate large molecules with limited hardware using a Suzuki–Trotter expansion, and the ability to implement realistic system-bath interactions with only one additional interaction per mode. Our approach can be implemented with current technology; e.g., the conical intersection in pyrazine can be simulated using a single trapped ion. Therefore, we expect our method will enable classically intractable chemical dynamics simulations in the near term.

---

I. Introduction

Computational chemistry aims to predict energies, structures, reactivity, and other properties of molecules. Although molecular dynamics is, in principle, best simulated with a fully quantum-mechanical treatment of coupled electrons and nuclei, the computational cost of doing so scales exponentially with molecular size, making it intractable for most chemical systems. The cornerstone of the vast majority of quantum-chemistry methods is the Born–Oppenheimer (adiabatic) approximation, which neglects the coupling between electronic and nuclear degrees of freedom.¹ This approximation fails in regions of strong non-adiabatic coupling, particularly near degeneracies between electronic states, known as conical intersections.¹ Non-adiabatic couplings are essential to photochemistry, where conical intersections act as funnels from one electronic state to another on ultrafast timescales (those comparable to nuclear vibrational periods, fs–ps).^1,2 State-of-the-art algorithms, such as (multilayer) multiconfigurational time-dependent Hartree (MCTDH),^3,4 can significantly reduce the simulation cost using a careful choice of basis-set contractions optimized for each system, enabling simulations of particular systems with tens to hundreds of modes.^5–7 Nevertheless, the optimal form of MCTDH wavefunctions cannot be predicted a priori, and the method—like all fully quantum-mechanical treatments of chemical dynamics—has an exponential worst-case scaling with system size.

Quantum computing promises to reduce the steep computational cost of simulating quantum systems by representing the problem on a controllable quantum device.^8–16 Most chemical applications of quantum computing have focused on time-independent observables,^17–27 in an effort to reduce the cost of electronic-structure methods. Such quantum methods optimize trial wavefunctions and map observables onto registers of qubits, and can predict molecular properties with linearly many qubits and polynomially many quantum gates as a function of molecule size.^17,18 Due to limitations in qubit count and coherence times of current quantum computers, many recent methods use hybrid approaches, such as variational quantum eigensolvers, to divide the calculation into classical and quantum tasks.^20–24,26 By contrast, few quantum algorithms focus on the simulation of molecular dynamics. Proposed methods include the simulation of electronic²¹ and vibrational^28,29 dynamics and the coupling between them.^30,31 Although each of these methods scales polynomially with system size, their qubit requirements restrict them to small model systems.

Analog quantum simulators provide an alternative approach to quantum simulation by mapping a desired Hamiltonian onto a purpose-built quantum system.^11–15 As a result, the dynamics of the simulator directly corresponds to the dynamics of the simulated system, and allowing any observables to be read out, in principle, using suitable measurements. Qubit-based analog simulators have already been used to demonstrate quantum advantage by carrying out classically-intractable simulations of spin models.^32,33 As a chemical application, ultracold atoms in optical lattices were proposed to simulate the electronic degrees of freedom of a molecule in a grid-based adiabatic picture.³⁴

Particularly promising recent approaches to analog quantum simulation have taken advantage of bosonic modes natively present in certain architectures, including photonic chips, microwave resonators, and trapped ions. Using native bosonic degrees of freedom gives analog simulators an advantage over digital ones, which would require many qubits to accurately represent a single harmonic oscillator. Chemical applications of bosonic simulation have mainly focused on the vibrational degrees of freedom of molecules, such as vibrational dynamics^35,36 and Franck–Condon vibronic spectra.^19,36–39 Approaches that include the entanglement of internal states to bosonic modes could simulate more complex processes, such as the dynamics of fermionic lattice models,^40,41 with possible extensions to electron-nuclear dynamics for quantum chemistry,^41,42 as well as models of charge and energy transfer.^43–46

Here, we show that analog quantum simulators can efficiently simulate non-adiabatic chemical dynamics. Our approach can be implemented using any quantum system containing a qudit with controllable couplings to a set of bosonic modes, a device we call a mixed qudit-boson (MQB) simulator. We achieve linear scaling in molecular size by mapping a molecule with d electronic states and N vibrational modes onto an MQB with a d-level qudit and N bosonic modes. We analyze the non-adiabatic dynamics of pyrazine as an example, showing it can be simulated using a single trapped-ion MQB, and we detail how to extend the procedure to large systems using a Trotterization scheme and to open systems using laser cooling.

II. Hamiltonian mapping

Our approach uses vibronic coupling (VC) Hamiltonians, which are widely used to model chemical dynamics. They express molecular electronic matrix elements as analytical functions of the nuclear coordinates, often as a power series,¹where |n〉 is the nth electronic state and Q̂_j = (μ_jω_j)^1/2q̂_j and P̂_j = (μ_jω_j)^−1/2p̂_j are the dimensionless position and momentum of mode j (which has reduced mass μ_j, frequency ω_j, position Q̂_j, and momentum P̂_j). We set ℏ = 1 throughout. The coefficients c₀^(n,m), c_j^(n,m), etc. are the (real) expansion coefficients of the molecular potential energy about the reference geometry, typically the minimum of the ground electronic state. While adiabatic potential energy surfaces obtained from electronic structure methods are ill-behaved, with features such as conical intersections and singular kinetic-energy couplings, they can be transformed to the diabatic surfaces and couplings of VC Hamiltonians given by eqn (2) (Fig. 1).¹ For most applications, linear vibronic coupling (LVC) or quadratic vibronic coupling (QVC) models accurately represent photochemical observables.^1,47 However, the expansion order of the Hamiltonian has little effect on the classical memory required to store the wavefunction.

Fig. 1 Representations of molecular potential-energy surfaces and their link to chemical dynamics. (a) Adiabatic surfaces correspond to the eigenvalues of the Born–Oppenheimer Hamiltonian at all nuclear displacements. The upper and lower surfaces shown correspond to higher- and lower-energy electronic states. Two inset molecules illustrate different molecular geometries at two sets of nuclear coordinates Q₁ and Q₂. An initial wavepacket (blue, top left) can slide down the upper surface and, upon reaching the conical intersection (the cusp where the two surfaces are degenerate) it can split into two entangled branches (blue and cyan), one on each electronic surface. (b) Diabatic surfaces (red, orange) of a vibronic coupling Hamiltonian (eqn (1)) formed from the adiabatic surfaces in a. The diabatic surfaces are analytic, avoiding numerical divergences near conical intersections. The couplings between the diabatic surfaces are not shown.

To map VC Hamiltonians onto MQB simulators, we exploit the similarities between terms in Ĥ_VC and common entangling gates in MQB architectures. We start by expressing an LVC model in the interaction picture,

where and â_j are the bosonic creation and annihilation operators such that and c₀^(n,m) is the potential energy of state n at the reference geometry. Ĥ^LVC_I is a generalized Jaynes–Cummings Hamiltonian of a d-level system coupled to N bosonic modes, the same interaction that is used for digital quantum computation with existing architectures, including ion traps and cQED. Typically, ion-trap quantum computers encode qubits in the electronic states of the ions, and lasers are used to entangle them by transient coupling through the ions' vibrations in the trapping potential.^48–50 Likewise, circuit quantum electrodynamics (cQED) uses microwave cavities coupled to superconducting qubits to implement multi-qubit gates.⁵¹ Any MQB architecture can be described, to first order, with the Hamiltonian

Each of the parameters in Ĥ^MQB_I can be tuned using well-developed light-matter interactions. For example, for an ion trap, corresponds to a time-dependent AC Stark shift from a pair of non-copropagating lasers, is proportional to the Rabi frequency of a Mølmer–Sørensen interaction, χ_n is a time-independent AC Stark shift and δ_j is the detuning from (sideband) resonance for mode j (see Appendix A for details). Therefore, because eqn (4) is of the same form as eqn (3), the terms of the MQB Hamiltonian can be tuned to match those of an LVC Hamiltonian.

The MQB approach can also be extended to higher-order expansions; the second-order terms of a QVC Hamiltonian, in the interaction picture, become

Both the dispersive (j = k) and mode-mixing (j ≠ k) terms have been demonstrated experimentally in trapped ions^52,53 and cQED.^54,55

Overall, the numbers of qudit states and bosonic modes required for the MQB mapping scale linearly with the corresponding numbers of molecular degrees of freedom, a significant improvement over the exponential classical scaling of fully quantum-mechanical chemical dynamics methods. Fig. 2 shows that even optimised classical methods such as MCTDH are limited to tens of modes for many systems, meaning that small to medium MQBs could achieve quantum advantage. Fig. 2 also shows that the MQB approach can simulate chemical dynamics with a significant improvement in quantum hardware requirements compared to digital approaches.^30,31 Our advantage comes from the larger Hilbert space in a bosonic mode, which would otherwise require many qubits to simulate, possibly with complicated interactions between them being required to simulate simple bosonic processes. In MQB architectures, this advantage can be decisive, because a single well-controlled bosonic mode might be comparably difficult to implement experimentally as a single well-controlled qubit.

Fig. 2 Resource (space) scaling of different methods for simulating vibronic dynamics. (a) Classical memory required for conventional (conv.) and MCTDH approaches in bits (b) and bytes (B = 8b), and quantum hardware resources for digital and MQB approaches. For MCTDH, an approximate range is shown, because the performance of the method depends on the Hamiltonian and cannot be accurately predicted a priori. The quantum resources are qubits (for digital algorithms, assuming 8 qubits per simulated mode), ions (for ion-trap MQB, assuming 3 vibrations per ion), and microwave resonators (for cQED MQB, assuming 20 modes per resonator). (b) Linear-scale comparison of quantum approaches (range as the grey box in a).

The advantage of the MQB approach is due to the natural mapping between molecular vibrations and MQB bosonic modes, which obviates the need, present on digital computers, to represent the large bosonic Hilbert space using qubits. The precise quantum requirements depend on the nature of the MQB simulator: if trapped ions are used, each ion provides three vibrational modes, while in cQED, each resonator could yield dozens of microwave modes.^56,57 The cost of the simulation can also be measured in the total number of interaction terms required, which scales as N^kd(d + 1)/2 for an order-k expansion of eqn (2). The vast majority of VC Hamiltonians are LVC models, which scale linearly with the number of vibrational modes, or QVC models, which scale quadratically.¹ The number of vibrations is a measure of the size of the molecule because it equals N = 3N_n − 6 in a molecule with N_n atoms. Furthermore, the simulation may be simplified by excluding certain modes or interactions that do not participate in the dynamics, either exactly (due to symmetry) or approximately (due to weak coupling).

III. Analog MQB simulation

An MQB simulation of molecular dynamics comprises three steps: initialization, evolution, and measurement.

Many initial states can readily be prepared on an MQB simulator using established techniques. In particular, in most photochemical studies, the Franck–Condon principle implies that photoexcitation only promotes the electronic degree of freedom to an excited state, leaving the vibrations unchanged.¹ On the MQB simulator, this corresponds to a change of the qudit state with no change to the bosonic state. In more general cases, superpositions of simulator qudits states could be prepared to simulate molecules with nearly degenerate electronic states or undergoing broadband excitation. Sometimes, it is convenient to work with displaced nuclear coordinates, in which case the initial state must be displaced as well, as discussed in Appendix B.

Simulating the time evolution involves implementing qudit-boson interactions according to eqn (3) and (5) to mimic the VC Hamiltonian for the duration of the evolution. Given the correct light-matter interactions, the MQB simulator evolves in real time with the same dynamics as the VC model, except that the Hamiltonian is scaled by a factor F so that the evolution occurs on the natural timescale of the MQB simulator. The scaling freedom is a major advantage for the simulation: whereas molecular vibronic frequencies are typically 10–100 THz, MQB simulators operate at frequencies that are many orders of magnitude less (e.g., kHz for trapped ions, F ∼ 10⁻⁹, and GHz for cQED, F ∼ 10⁻³). The dynamics on the simulator thus occur in extreme slow motion relative to the simulated molecule, meaning the quantum simulator, like classical algorithms, has a much greater time resolution than ultrafast molecular spectroscopy.

Finally, after evolution for the desired amount of time, the MQB simulator can be probed to measure observables. Most chemically important observables do not require knowledge of the full wavefunction (which would require exponentially scaling full state tomography), and are instead easily measurable properties of the qudits or the bosons, such as the electronic-state populations or nuclear positions and momenta. Most MQB architectures have established methods for measuring qudit states; boson observables can also be measured, often by first mapping them onto the qudits.^38,58 Due to the statistical nature of quantum measurement, observables must be averaged over many experiments, while repeating the experiment for different simulation times would yield temporal information. Importantly, the results at any simulation time are independent of the measurements from previous times and the time intervals at which the measurements are made.

IV. Example: 2D LVC model of pyrazine

To demonstrate our approach, we consider the MQB simulation of a two-state, two-mode (or two-dimensional, 2D) LVC model of pyrazine using a single trapped ion (Fig. 3). Pyrazine is a canonical model system for molecular dynamics because the conical intersection between its two excited electronic states, the optically bright ππ* state and the dark nπ* state, can be described with only two modes. The vibrations comprise a fully symmetric (tuning) mode and a symmetry-breaking (coupling) mode.^47,59 The ground electronic state is energetically separated from and weakly coupled to the excited states, meaning it can be excluded from the dynamics simulation.⁴⁷ The 2D LVC Hamiltonian iswhere |0〉 and |1〉 are the nπ* and ππ* states, modes 1 and 2 are the tuning and coupling modes, [small sigma, Greek, circumflex]_z = |0〉〈0| − |1〉〈1| and [small sigma, Greek, circumflex]_x = |0〉〈1| + |1〉〈0| are Pauli matrices, and ΔE = c₀^(1,1) − c₀^(0,0) is the electronic energy difference at the ground-state-minimum geometry.⁴⁷ The initial state of the simulation, corresponding to Franck–Condon excitation, is in the ground vibrational state, with the electronic state promoted to the ππ* state, i.e., |1〉 ⊗ |v₁ = 0〉 ⊗ |v₂ = 0〉 for vibrational eigenstates v_j.

Fig. 3 Example of the mapping between the degrees of freedom of a molecule (pyrazine, eqn (6)) and an MQB simulator (one trapped-ion qubit, eqn (8)). (a) The degrees of freedom of the molecule: electronic, shown as red and orange electronic orbitals; vibrational, shown as cyan and blue atomic displacements; and their couplings shown in purple and magenta. (b) The corresponding degrees of freedom for the trapped ion: qubit states (red and orange levels); bosonic modes (cyan and blue potential curves); and laser-mediated interactions (purple and magenta).

In the absence of a laser field, the Hamiltonian of a single, two-level trapped ion is

where ω₀ is the frequency difference between qubit states and ω^ion_j is the trap frequency of mode j. The diagonal and off-diagonal couplings of eqn (6) can be generated, respectively, by Raman interactions with bichromatic laser fields having frequency differences ω^ion_j − δ_j and ω₀ ± (ω^ion_j − δ_j), where δ_j is a detuning.⁵⁰ The resulting Hamiltonian can be expressed in a rotating frame (i.e., by transforming eqn (4)) into the time-independent formwhere Δχ is the energy difference between qubit states induced by a time-independent AC Stark shift, is proportional to the time-dependent AC Stark shift of state n and Ω′ is proportional to the Rabi frequency of the internal states of the ion (details in Appendix A).

The initialization, evolution, and measurement steps of a pyrazine simulation are illustrated in Fig. 4a–c. This simulation can be simplified by working in a displaced frame, which can transform the tuning term (third term in eqn (8)) to either which are more easily implemented in trapped ions (details in Appendix B). To compensate for the frame displacement, Fig. 4a shows an initial wave function displacement.

Fig. 4 Illustrative molecular-dynamics simulation using an MQB simulator (trapped-ion qubit for concreteness). (a) The initial wavefunction is prepared by an excitation of the qudit state (yellow arrow/laser) and a displacement along the relevant mode(s) (blue arrow/laser). (b) The VC Hamiltonian is implemented through interactions with laser fields (purple and magenta, see eqn (4)), which are equivalent to displacements and couplings between potential energy surfaces. These interactions are maintained for the duration of the simulation. (c) After evolution by the effective Hamiltonian, observables of the wavefunction (inset: excited state population) are measured by a detector (e.g., fluorescence of the qubit state induced by the red laser). Multiple measurements are required to obtain expectation values, and the time dependence of observables is obtained by varying the time delay from initialization to measurement. (d) The evolution step can be discretized using a Suzuki–Trotter expansion into a series of short laser interactions, particularly for the treatment of many modes. (e) Open quantum systems can be simulated by adding sideband laser cooling (green) during the evolution step.

V. Analog Trotterization

When simulating larger systems, it may be impractical to implement all VC interactions simultaneously on an MQB simulator, because doing so would require a separate light-matter interaction for each term in eqn (1). Instead, the evolution of the wavefunction can be discretized with a Suzuki–Trotter expansion⁹ by applying sets of terms in sequence for short times (Fig. 4d), allowing the same hardware to implement multiple interactions. A consequence of Trotterization is a trade-off between the simulation accuracy and the size of Trotter steps; for example, in the first-order Suzuki–Trotter expansion, the error scales with the square of the Trotter time step, Δt². Furthermore, as in many analog simulations, the problem is complicated by the fact that the Hamiltonian of an MQB simulator is non-zero in the absence of interactions. For example, the vibrational degrees of freedom in trapped ions continue to oscillate whether or not they are perturbed by a laser field.

Here, we describe two approaches to analog Trotterization. We assume that the base Hamiltonian Ĥ₀ cannot be turned off, while interaction Hamiltonians Ĥ_k, for 1 ≤ k ≤ M, can be turned on and off on demand.

The rescaling approach relies on the fact that, in an analog simulation, the total Hamiltonian can be scaled by a constant factor, changing only the speed of the simulation. The individual terms can also be rescaled, as long as the total Hamiltonian remains proportional to the molecular Hamiltonian. If we re-express the full Hamiltonian as the first-order evolution (for timestep Δt) is Essentially, the evolution by Ĥ₀ is slowed down relative to the interaction terms to compensate for the continuous evolution under Ĥ₀ while the Ĥ_k are implemented in turn. Higher-order Suzuki–Trotter expansions can be used to reduce the Trotterization error.

The rewinding approach treats Ĥ₀ and the interaction terms on the same timescale (i.e., without rescaling), instead reversing the excessive evolution under Ĥ₀. This approach assumes that time evolution under can be implemented between Trotter steps. At first order, this scheme is given by i.e., by correcting the evolution by exp(iĤ₀Δt) after each Trotter step, which can be thought of as implementing another unitary of the form exp(−i(−2Ĥ₀)Δt) while the base Hamiltonian is always on.

The two approaches have similar performance, as can be seen in the example of a two-mode system, where

where ‘res.’ and ‘rew.’ denote rescaling and rewinding, respectively. For the pyrazine example, we let , and , with Δχ₁ = Δχ₂ = Δχ/2. The effect of Trotterization is shown in Fig. 5: as expected, increasing the Trotter step Δt increases the error in the electronic population (Fig. 5a). The fidelities (Fig. 5b) show larger errors; however, fidelities of both rescaling and rewinding approaches are nearly identical, indicating that either may be used for simulations.

Fig. 5 Evaluation of Trotterization approaches using the pyrazine LVC Hamiltonian. (a) Electronic populations as a function of time with different Trotter steps Δt, using the rescaling approach. (b) The rescaling (‘res.’) and rewinding (‘rew.’) approaches result in comparable fidelities (overlaps of the Trotterized wavefunction, |Ψ_trot(t)〉, with the exact wavefunction, |Ψ(t)〉, with Δt = 0).

A difference between the two approaches arises in the limit of many interaction terms, where one approach may be preferred over another depending on experimental limitations; the rescaling approach requires fewer operations, while the rewinding approach ensures that all components evolve at the same rate.

VI. Open quantum systems

MQB devices also allow the simulation of open quantum systems—such as molecules in the condensed phase—with minimal overhead, using existing dissipation processes in the simulator. In general, open-system simulation is even more computationally expensive on classical computers, as it involves the evolution of a density matrix rather than a state vector. Classical MCTDH simulations have been achieved with weak coupling to baths, but increasing the number of bath modes quickly becomes intractable.^3,4 By contrast, dissipation on an MQB simulator could be implemented using available environmental couplings. For example, analog trapped-ion devices have simulated electronic dissipation with optical pumping⁶⁰ and vibrational dissipation with sympathetic laser cooling.^45,46

Here, we show how ion-trap MQB simulators could simulate system-bath coupling with a desired coupling strength and effective temperature (Fig. 4e). While simulating arbitrary system-bath interactions would be complicated, we focus on the leading contribution to dissipation in molecules in a thermal bath, coupling to vibrational modes.⁶¹ Our approach extends system-bath simulations using laser cooling,⁴⁶ and requires a maximum number of interactions that scales linearly with the number of modes coupled to the environment.

The conditions for a molecular simulation at a finite temperature are commonly given by the coupling of a constant-temperature bath of harmonic oscillators to the molecular modes to represent inter-molecular collisions.^4,61 We consider a bath at temperature T with jump operators (γ_j(n̄_j + 1))^1/2â_j for cooling and for heating of mode j, where γ_j are the system-bath coupling strengths and n̄_j = (exp(ω_j/k_BT) − 1)⁻¹ are the Bose–Einstein occupation numbers of the bath. This gives the master equation

where [small rho, Greek, circumflex] is the density operator of the system and is the Linblad superoperator for the jump operator L̂_j, given by

Similar time evolution can be obtained using laser cooling, which is used in ion-trap quantum computation to cool the ions to their ground motional state. Applied lasers both cool and heat; in the resolved sideband limit, the combined evolution due to laser cooling and heating takes the form^45,62

wherewhere η_j is the Lamb–Dicke parameter, Γ_j is the spontaneous emission rate of the carrier transition (i.e., the |0〉 → |1〉 transition), Δ is the detuning from the carrier transition, α is the angular factor for spontaneous emission (α = 2/5 for dipole transitions⁶²), and Ω₀ ≪ Γ_j is the Rabi frequency of the carrier transition.⁶²

Comparing eqn (11) and (12) shows that laser cooling and heating can simulate molecular thermal dissipation if the experimental parameters are chosen so that γ_j = A⁻_j − A⁺_j and n̄_j = A⁺_j/(A⁻_j − A⁺_j), which can be achieved in two steps. First, because A⁺_j/(A⁻_j − A⁺_j) depends only on Δ and Γ_j, those two parameters can be chosen to simulate a specific n̄_j and, therefore, the temperature. Second, A⁻_j − A⁺_j is proportional to the square of the remaining free parameter, Ω₀, which depends on the tunable laser intensity and can therefore be used to set the system-bath interaction γ_j for each mode j.

The experimental couplings γ_j must be scaled to simulator frequencies, as was done for Ĥ. Such a rescaling changes the timescale in eqn (11), giving Rabi frequencies in Hz to kHz instead of GHz to THz. By contrast, the occupation numbers n̄_j are dimensionless and require no scaling. The difference between molecular and trapped-ion frequencies results in a difference in their vibrational temperatures because the invariable n̄_j depends only on ω_j/k_BT. For example, a simulation at 300 K corresponds to an ion-trap vibrational temperature of ∼40 μK, which is readily achievable in ion traps.⁶³

Fig. 6 shows an example of weak system-bath coupling with the LVC model of pyrazine at 300 K. The simulation assumes the Ohmic system-bath couplings that are widely used in chemical dynamics,⁵⁹γ_j = γ₀ω_je^{−ω_j/ω₀}, where γ₀ is constant and ω₀ is the high-frequency cutoff.

Fig. 6 Electronic populations of pyrazine, using the LVC Hamiltonian and different system-bath couplings γ₀, all at a temperature of 300 K. As γ₀ is increased, the details of the closed-system result (γ₀ = 0) are lost and the coherent oscillations are damped.

Our scheme requires at most one laser-cooling interaction per mode. However, because chemical dynamics often takes place at low temperatures (k_BT ≪ ω_j), it may be possible to simultaneously cool all the vibrational modes with only a single broadband cooling laser. Even at room temperature, n̄_j ≪ 1 for typical molecular vibrations with frequencies of tens of THz. In that limit, only the cooling terms in eqn (11) are significant, . If we also assume that all the γ_j are comparable (which is true for a slowly varying ohmic spectral density and comparable values of ω_j), we arrive at a master equation that can be simulated with just one broadband cooling laser, which can cool all the modes simultaneously.⁶⁴ For example, repeating the simulation of Fig. 6 with a single cooling laser and equal γ_j leads to imperceptible changes in the populations shown.

VII. Limitations

The main limitation of our approach—and of any analog quantum simulation—is the accumulation of simulation errors due to decoherence, dissipation, and the absence of error correction. As a result, hardware limitations place limits on the size of molecular systems that can be simulated and the duration of those simulations. Of course, the same was true of analog classical simulations, which were indispensable tools in the early decades of classical computation. Even with these limitations, early MQBs are likely to be useful, especially if the decoherence is harnessed as a resource.

The principal sources of error for MQB simulation will be decoherence and thermalization, but the subspace most affected may depend on the architecture. For example, while trapped-ion simulators have electronic coherence times of seconds or longer, motional decoherence (typically dephasing) will limit simulations because it sets in within 1–100 ms.^65,66 By contrast, in cQED, the limiting factor is likely to be the qudit coherence time, typically 50–100 μs and shorter than the coherence time of the cavity modes, which is limited by cavity losses and may be 2 ms or more.^54,67

Whether an MQB simulator will be useful depends on the ratio of its decoherence time and the necessary simulation time at realistic experimental frequencies. A simulation requires a scaling factor F that relates molecular and simulator frequencies, ω_sim = Fω_mol. In order to complete a simulation for a time t_max (in molecular timescales) within a decoherence time τ_d, F must exceed t_max/τ_d. If, in addition, the time evolution is broken into M Trotter components, the lower bound of F is F_min = Mt_max/τ_d. There are also two upper limits on F determined by experimental conditions: first, as the ratio of the maximum experimentally achievable frequency or coupling to the corresponding molecular parameter, e.g., for a coupling interaction; and second, in a Trotterized simulation, as the ratio of the maximum acceptable molecular-scale timestep to the minimum timestep achievable experimentally, F⁽²⁾_max = Δt_mol/Δt_sim. The simulation will be completed within the decoherence time if F can be chosen to meet all these constraints, i.e., if both F⁽¹⁾_max and F⁽²⁾_max are larger than F_min. To minimize decoherence, F should ideally be set to the upper bound.

The constraints above can be relaxed if the native dissipation in an MQB simulator can be harnessed to simulate dynamics in open quantum systems, which are often of greater practical relevance than isolated molecules. For example, decoherence and dissipation of bosonic modes might be tuned (e.g., via temperature) to values that simulate elastic and inelastic collisions in condensed phases or at high pressures. For such a simulation, care must be taken to consider the relative dissipation rates of electronic and vibrational modes, particularly for cQED simulators where the relevant internal and bosonic frequencies have similar magnitudes.

The limits of a simulation also depend on the observable of interest, meaning that many simulations will be able to run longer than worst-case estimates from fidelities. In general, errors in analog simulation depend strongly on the observable, with global properties being more robust.^68,69 For example, in Fig. 5 with a Trotter step of 0.5 fs (blue lines), the population deviates by at most 0.05 in 300 fs, whereas the fidelity rapidly decays to less than 0.6 in the same time. Fortunately, most chemically interesting molecular observables are global properties, rather than properties of individual vibronic eigenstates, which would allow for greater error tolerance than fidelities would suggest.

Finally, in very large systems, the number of bosonic modes could be limited by the experimental apparatus. In trapped ions, the frequencies of the modes become increasingly spectrally congested as the number of ions increases, which may eventually lead to cross-talk due to multiple modes interacting with a single laser frequency.^50,70 So far, ion trap experiments have been performed with as many as 40 resolved modes.⁷⁰ In cQED, resonators are physically coupled to qudits, placing spatial limitations in connecting many resonators a single qudit. However, if multimode resonators^56,57 are used, their total number may be suitable for medium-scale simulations. Therefore, the current technological capabilities are more than enough for near-term demonstration and will undoubtedly improve as technology develops.

VIII. Conclusions

The one-to-one correspondence of molecular and simulator degrees of freedom gives MQB simulators a clear advantage in the simulation of non-adiabatic chemical dynamics over both classical and qubit-based schemes. For an increasing number of vibrational modes, our approach requires linearly scaling quantum hardware and polynomially many interaction terms, in either isolated or open molecular systems. For large systems, the number of simultaneous interactions can be significantly reduced with a Suzuki–Trotter expansion, allowing the same driving fields to be reused. Scaling the frequencies of the VC Hamiltonian to the MQB simulator also offers a significant advantage in time resolution relative to ultrafast spectroscopy experiments. Many of the tools necessary for MQB simulation have been developed in the context of digital quantum computation and can be repurposed for analog simulation, making MQB simulators a promising platform for quantum advantage in the simulation of chemical dynamics. Indeed, we expect that dynamical simulations will more readily achieve quantum advantage than solving time-independent electronic-structure problems.

IX. Appendices

Appendix A: Light-matter interactions in trapped-ion simulators

Here, we follow ref. 50 to derive the light-matter interactions required to map an LVC Hamiltonian onto an ion-trap Hamiltonian.

Tuning terms in the LVC Hamiltonian can be generated using two non-copropagating lasers, which can induce Raman transitions that couple the internal states of an N-level trapped ion with its vibrational mode j (frequency ω^ion_j). If the frequency between the two lasers is ω^ion_j − δ_j, the resulting interaction-picture Hamiltonian becomes

where η_j is the Lamb–Dicke parameter, is the Debye-Waller factor, and Θ_n,j is the time-dependent AC Stark shift of state n. The interaction of eqn (15) is well established; in particular, the two-qubit generalisation is widely used to implement the phase gate .⁴⁹ This form of the Hamiltonian assumes that the intermediate state |e〉 can be adiabatically eliminated, i.e., that the global detuning Δ of ω_α and ω_β from |e〉 is much larger than the linewidth γ_e of |e〉, and that the detuning δ_j is small relative to the ion trap frequency . The laser field also causes a time-independent AC Stark shift for each state, χ⁽ⁿ⁾, which modulates the internal state energies and can be experimentally tuned by Δ.⁴⁹ Without loss of generality, the phase ϕ can be dropped as it corresponds to an arbitrary initial phase-space angle.

Coupling terms in the LVC Hamiltonian can be generated if, instead, the difference between the laser frequencies is ΔE^ion_n,m ± (ω^ion_k − δ_k), where ΔE^ion_n,m is the energy difference between internal states n and m, and the plus or minus determining whether the interaction corresponds to a blue- or red-sideband transition. The sum of red and blue sideband terms gives a Hamiltonian of the form

where Ω_n,m,k is the base Rabi frequency of the internal states, ϕ_S is the spin-phase of the ion resulting from the sum of the two transitions, and ϕ_M is the phase difference between the red and blue sidebands. As with the interaction above, this interaction is widely used, especially to generate the Mølmer–Sørensen entangling gate.^48,71 Setting ϕ_S = 0 results in an interstate coupling in the same form as in the LVC Hamiltonian. As with ϕ before, the phase ϕ_M can be dropped.

The total interaction Hamiltonian is obtained by adding an instance of eqn (15) for each tuning mode and an instance of eqn (16) for each coupling mode. Once the time-independent AC Stark shifts χ_n are also included, the total Hamiltonian becomes

where and the sets t and c are indices of tuning and coupling modes, respectively. This Hamiltonian can be transformed into a rotating frame to give

The results above are general for any number of internal states and any number of vibrational modes. Although a single ion has only three vibrational modes, additional modes can be obtained by adding ancillary, optically inactive ions into the trap. The ancillary ions can be made optically inactive in multiple ways, including the use of different isotopes⁷² or by pumping them into long-lived shelved states.^73,74

For a two-state system (n, m ∈{0, 1} with two modes (t = {1}, c = {2}), eqn (18) simplifies to eqn (8) of the main text.

Appendix B: Displaced vibronic coupling Hamiltonians

For experimental implementation of VC Hamiltonian simulations, it may be convenient to transform the elements of the Hamiltonian into a different form. In particular, the interstate coupling factors of the form can be replaced with more easily implemented operators, such as [small sigma, Greek, circumflex]_x or [small sigma, Greek, circumflex]_z.

One way to do this is by a coherent displacement of the Hamiltonian along a single mode. This corresponds to the action of a displacement operator, D̂_k(β) = exp(−iβP̂_k) for mode k and real β, where

is a change in position by β along mode k. For a general VC Hamiltonian (eqn (1)), the displacement gives

Dropping the constant term ω_kβ²/2 and truncating to linear terms (c_k,j^(n,m) = 0) yields the displaced LVC Hamiltonian,

where δ_n,m is a Kronecker delta, meaning only the diagonal elements (n = m) of first-order terms are affected by displacement.

In the two-state case (n,m ∈ {0,1}), the electronic components of the LVC Hamiltonian can be expressed in terms of Pauli operators to yield

where ΔE = c₀^(1,1) − c₀^(0,0), W₀ = c₀^(0,1), [small kappa, Greek, macron]_j = c_j^(1,1) + c_j^(0,0)/2, Δκ_j = c_j^(1,1) − c_j^(0,0) and λ_j = c_j^(0,1). After displacement along mode k,where the constant energy β[small kappa, Greek, macron]_k is excluded as it only contributes a global phase. This is equivalent to making the replacements ΔE → ΔE − βΔκ_k, W₀ → W₀ − βλ_k and [small kappa, Greek, macron]_k → [small kappa, Greek, macron]_k − βω_k, which could be used to simplify the Hamiltonian. For example, choosing β = [small kappa, Greek, macron]_k/ω_k would remove the term from eqn (22).

The dynamics simulated with the displaced VC Hamiltonian is exactly equivalent (i.e., produces the same observables) to that of the original model as long as the initial wavefunction is displaced by an equal amount, which is achievable for common bosonic simulators.^72,75

Appendix C: Numerical methods

Simulations of the closed-system LVC model of pyrazine⁵⁹ were performed using NumPy⁷⁶ and SciPy⁷⁷ with the number of vibrational eigenstates truncated to 20 for both modes (ν_6a and ν_10a). Open-system dynamics used the Linblad master-equation solver in QuTiP.⁷⁸

Author contributions

CH and IK conceived the original idea; RJM, CED, CJTB, AK, and IK developed the theoretical methods and carried out simulations; RJM, CED, CJTB, CLE, and CH developed mappings to concrete architectures; all authors discussed the results and interpreted the data; RJM wrote the manuscript with feedback from all authors.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge valuable discussion with Alistair Milne, Ramil Nigmatullin, Ting Rei Tan, and Joel Yuen-Zhou. We were supported by a Westpac Scholars Trust Research Fellowship, by the Lockheed Martin Corporation, by the Australian Government's Defence Science and Technology Group, by the United States Office of Naval Research Global (N62909-20-1-2047), and by the University of Sydney Nano Institute Grand Challenge Computational Materials Discovery. We were supported by computational resources and assistance from the National Computational Infrastructure (NCI) and the University of Sydney's computing cluster Artemis.

References

1. Conical Intersections: Electronic Structure, Dynamics & Spectroscopy, Adv. Ser. Phys. Chem., ed. W. Domcke, D. R. Yarkony, and H. Köppel, World Scientific, Singapore, 2004, vol. 15.
2. L. S. Cederbaum, H. Köppel and W. Domcke, Int. J. Quantum Chem., 1981, 20, 251.
3. G. A. Worth, H.-D. Meyer, H. Köppel, L. S. Cederbaum and I. Burghardt, Int. Rev. Phys. Chem., 2008, 27, 569.
4. H. Wang, J. Phys. Chem. A, 2015, 119, 7951.
5. Q. Meng and H.-D. Meyer, J. Chem. Phys., 2013, 138, 014313.
6. Y. Xie, J. Zheng and Z. Lan, J. Chem. Phys., 2015, 142, 084706.
7. J. Schulze, M. F. Shibl, M. J. Al-Marri and O. Kühn, J. Chem. Phys., 2016, 144, 185101.
8. R. Feynman, Int. J. Theor. Phys., 1982, 21, 467.
9. S. Lloyd, Phys. Rev. Lett., 1995, 75, 346.
10. M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, UK, 2000).
11. I. Buluta and F. Nori, Science, 2009, 326, 108.
12. R. Blatt and C. F. Roos, Nat. Phys., 2012, 8, 277.
13. A. Aspuru-Guzik and P. Walther, Nat. Phys., 2012, 8, 285.
14. I. Bloch, J. Dalibard and S. Nascimbène, Nat. Phys., 2012, 8, 267.
15. I. M. Georgescu, S. Ashhab and F. Nori, Rev. Mod. Phys., 2014, 86, 153.
16. F. Schäfer, T. Fukuhara, S. Sugawa, Y. Takasu and Y. Takahashi, Nat. Rev. Phys., 2020, 2, 411.
17. A. Aspuru-Guzik, A. D. Dutoi, P. J. Love and M. Head-Gordon, Science, 2005, 309, 1704.
18. B. P. Lanyon, J. D. Whitfield, G. G. Gillett, M. E. Goggin, M. P. Almeida, I. Kassal, J. D. Biamonte, M. Mohseni, B. J. Powell, M. Barbieri, A. Aspuru-Guzik and A. G. White, Nat. Chem., 2010, 2, 106.
19. M.-H. Yung, J. Casanova, A. Mezzacapo, J. McClean, L. Lamata, A. Aspuru-Guzik and E. Solano, Sci. Rep., 2014, 4, 3589.
20. 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.
21. Y. Li and S. C. Benjamin, Phys. Rev. X, 2017, 7, 021050.
22. A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow and J. M. Gambetta, Nature, 2017, 549, 242.
23. J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. d. Jong and I. Siddiqi, Phys. Rev. X, 2018, 8, 011021.
24. C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt and C. F. Roos, Phys. Rev. X, 2018, 8, 031022.
25. N. P. D. Sawaya and J. Huh, J. Phys. Chem. Lett., 2019, 10, 3586.
26. Y. Nam, J.-S. Chen, N. C. Pisenti, K. Wright, C. Delaney, D. Maslov, K. R. Brown, S. Allen, J. M. Amini, J. Apisdorf, K. M. Beck, A. Blinov, V. Chaplin, M. Chmielewski, C. Collins, S. Debnath, K. M. Hudek, A. M. Ducore, M. Keesan, S. M. Kreikemeier, J. Mizrahi, P. Solomon, M. Williams, J. D. Wong-Campos, D. Moehring, C. Monroe and J. Kim, Npj Quantum Inf., 2020, 6, 33.
27. F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, S. Demura, A. Dunsworth, D. Eppens, E. Farhi, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, S. Habegger, M. P. Harrigan, A. Ho, S. Hong, T. Huang, W. J. Huggins, L. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, J. Kelly, S. Kim, P. V. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, M. Lindmark, E. Lucero, O. Martin, J. M. Martinis, J. R. McClean, M. McEwen, A. Megrant, X. Mi, M. Mohseni, W. Mruczkiewicz, J. Mutus, O. Naaman, M. Neeley, C. Neill, H. Neven, M. Y. Niu, T. E. O'Brien, E. Ostby, A. Petukhov, H. Putterman, C. Quintana, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, D. Strain, K. J. Sung, M. Szalay, T. Y. Takeshita, A. Vainsencher, T. White, N. Wiebe, Z. J. Yao, P. Yeh and A. Zalcman, Science, 2020, 369, 1084.
28. S. McArdle, A. Mayorov, X. Shan, S. Benjamin and X. Yuan, Chem. Sci., 2019, 10, 5725.
29. N. P. D. Sawaya, F. Paesani, and D. P. Tabor, 2020, arXiv 2009.05066.
30. I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni and A. Aspuru-Guzik, Proc. Natl. Acad. Sci. U.S.A., 2008, 105, 18681.
31. P. J. Ollitrault, G. Mazzola and I. Tavernelli, Phys. Rev. Lett., 2020, 125, 260511.
32. J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk and J. J. Bollinger, Nature, 2012, 484, 489.
33. J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong and C. Monroe, Nature, 2017, 551, 601.
34. J. Argüello-Luengo, A. González-Tudela, T. Shi, P. Zoller and J. I. Cirac, Nature, 2019, 574, 215.
35. C. Sparrow, E. Martín-López, N. Maraviglia, A. Neville, C. Harrold, J. Carolan, Y. N. Joglekar, T. Hashimoto, N. Matsuda, J. L. O'Brien, D. P. Tew and A. Laing, Nature, 2018, 557, 660.
36. W. Chen, J. Gan, J.-N. Zhang, D. Matuskevich, and K. Kim, 2021, arXiv 2103.14299.
37. J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean and A. Aspuru-Guzik, Nat. Photonics, 2015, 9, 615.
38. C. S. Wang, J. C. Curtis, B. J. Lester, Y. Zhang, Y. Y. Gao, J. Freeze, V. S. Batista, P. H. Vaccaro, I. L. Chuang, L. Frunzio, L. Jiang, S. M. Girvin and R. J. Schoelkopf, Phys. Rev. X, 2020, 10, 021060.
39. H. Jnane, N. P. D. Sawaya, B. Peropadre, A. Aspuru-Guzik, R. Garcia-Patron, and J. Huh, 2020, arXiv 2011.05553.
40. J. Casanova, A. Mezzacapo, L. Lamata and E. Solano, Phys. Rev. Lett., 2012, 108, 190502.
41. L. Lamata, A. Mezzacapo, J. Casanova and E. Solano, EPJ Quantum Technol., 2014, 1, 9.
42. L. Lamata, A. Parra-Rodriguez, M. Sanz and E. Solano, Adv. Phys.: X, 2018, 3, 1457981.
43. L. García-Álvarez, U. Las Heras, A. Mezzacapo, M. Sanz, E. Solano and L. Lamata, Sci. Rep., 2015, 6, 27836.
44. D. J. Gorman, B. Hemmerling, E. Megidish, S. A. Moeller, P. Schindler, M. Sarovar and H. Haeffner, Phys. Rev. X, 2018, 8, 011038.
45. A. Lemmer, C. Cormick, D. Tamascelli, T. Schaetz, S. F. Huelga and M. B. Plenio, New J. Phys., 2018, 20, 073002.
46. F. Schlawin, M. Gessner, A. Buchleitner, T. Schaetz, and S. S. Skourtis, 2020, arXiv 2004.02925.
47. L. Seidner, G. Stock, A. L. Sobolewski and W. Domcke, J. Chem. Phys., 1992, 96, 5298.
48. K. Mølmer and A. Sørensen, Phys. Rev. Lett., 1999, 82, 1835.
49. D. Leibfried, B. DeMarco, V. Meyer, D. Lucas, M. Barrett, J. Britton, W. M. Itano, B. Jelenković, C. Langer, T. Rosenband and D. J. Wineland, Nature, 2003, 422, 412.
50. P. J. Lee, K.-A. Brickman, L. Deslauriers, P. C. Haljan, L.-M. Duan and C. Monroe, J. Opt. B: Quantum Semiclassical Opt., 2005, 7, S371.
51. A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin and R. J. Schoelkopf, Phys. Rev. A, 2004, 69, 062320.
52. J. S. Pedernales, I. Lizuain, S. Felicetti, G. Romero, L. Lamata and E. Solano, Sci. Rep., 2015, 5, 15472.
53. K. Marshall and D. F. V. James, Appl. Phys. B, 2016, 123, 26.
54. A. Blais, S. M. Girvin and W. D. Oliver, Nat. Phys., 2020, 16, 247.
55. Y. Y. Gao, B. J. Lester, Y. Zhang, C. Wang, S. Rosenblum, L. Frunzio, L. Jiang, S. M. Girvin and R. J. Schoelkopf, Phys. Rev. X, 2018, 8, 021073.
56. N. M. Sundaresan, Y. Liu, D. Sadri, L. J. Szőcs, D. L. Underwood, M. Malekakhlagh, H. E. Türeci and A. A. Houck, Phys. Rev. X, 2015, 5, 021035.
57. A. J. Kollár, M. Fitzpatrick and A. A. Houck, Nature, 2019, 571, 45.
58. R. Gerritsma, B. P. Lanyon, G. Kirchmair, F. Zähringer, C. Hempel, J. Casanova, J. J. García-Ripoll, E. Solano, R. Blatt and C. F. Roos, Phys. Rev. Lett., 2011, 106, 060503.
59. A. Kühl and W. Domcke, J. Chem. Phys., 2002, 116, 263.
60. J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller and R. Blatt, Nature, 2011, 470, 486.
61. A. Nitzan, Chemical Dynamics in Condensed Phases, Oxford University Press, 2006.
62. S. Stenholm, Rev. Mod. Phys., 1986, 58, 699.
63. F. Diedrich, J. C. Bergquist, W. M. Itano and D. J. Wineland, Phys. Rev. Lett., 1989, 62, 403.
64. R. Lechner, C. Maier, C. Hempel, P. Jurcevic, B. P. Lanyon, T. Monz, M. Brownnutt, R. Blatt and C. F. Roos, Phys. Rev. A, 2016, 93, 053401.
65. D. M. Lucas, B. C. Keitch, J. P. Home, G. Imreh, M. J. McDonnell, D. N. Stacey, D. J. Szwer, and A. M. Steane, 2007, arXiv 0710.4421.
66. I. Talukdar, D. J. Gorman, N. Daniilidis, P. Schindler, S. Ebadi, H. Kaufmann, T. Zhang and H. Häffner, Phys. Rev. A, 2016, 93, 043415.
67. H. Wang, M. Hofheinz, M. Ansmann, R. C. Bialczak, E. Lucero, M. Neeley, A. D. O'Connell, D. Sank, M. Weides, J. Wenner, A. N. Cleland and J. M. Martinis, Phys. Rev. Lett., 2009, 103, 200404.
68. M. Sarovar, J. Zhang and L. Zeng, EPJ Quantum Technol., 2017, 4, 1.
69. P. M. Poggi, N. K. Lysne, K. W. Kuper, I. H. Deutsch and P. S. Jessen, PRX Quantum, 2020, 1, 020308.
70. N. Friis, O. Marty, C. Maier, C. Hempel, M. Holzäpfel, P. Jurcevic, M. B. Plenio, M. Huber, C. Roos, R. Blatt and B. Lanyon, Phys. Rev. X, 2018, 8, 021012.
71. A. Sørensen and K. Mølmer, Phys. Rev. Lett., 1999, 82, 1971.
72. C. Hempel, B. P. Lanyon, P. Jurcevic, R. Gerritsma, R. Blatt and C. F. Roos, Nat. Photonics, 2013, 7, 630.
73. K. Hosaka, S. A. Webster, A. Stannard, B. R. Walton, H. S. Margolis and P. Gill, Phys. Rev. A, 2009, 79, 033403.
74. C. L. Edmunds, T. R. Tan, A. Milne, A. Singh, M. Biercuk, and C. Hempel, 2020, arXiv 2012.14606.
75. B. Vlastakis, G. Kirchmair, Z. Leghtas, S. E. Nigg, L. Frunzio, S. M. Girvin, M. Mirrahimi, M. H. Devoret and R. J. Schoelkopf, Science, 2013, 342, 607.
76. C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke and T. E. Oliphant, Nature, 2020, 585, 357.
77. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt and SciPy 1.0 Contributors, Nat. Methods, 2020, 17, 261.
78. J. R. Johansson, P. D. Nation and F. Nori, Comput. Phys. Commun., 2013, 184, 1234.

---

Footnote

† Present address: Paul Scherrer Institute, 5232 Villigen PSI, Switzerland

---