Predicting molecular vibronic spectra using time-domain analog quantum simulation

Abstract

Spectroscopy is one of the most accurate probes of the molecular world. However, predicting molecular spectra accurately is computationally difficult because of the presence of entanglement between electronic and nuclear degrees of freedom. Although quantum computers promise to reduce this computational cost, existing quantum approaches rely on combining signals from individual eigenstates, an approach whose cost grows exponentially with molecule size. Here, we introduce a method for scalable analog quantum simulation of molecular spectroscopy: by performing simulations in the time domain, the number of required measurements depends on the desired spectral range and resolution, not molecular size. Our approach can treat more complicated molecular models than previous ones, requires fewer approximations, and can be extended to open quantum systems with minimal overhead. We present a direct mapping of the underlying problem of time-domain simulation of molecular spectra to the degrees of freedom and control fields available in a trapped-ion quantum simulator. We experimentally demonstrate our algorithm on a trapped-ion device, exploiting both intrinsic electronic and motional degrees of freedom, showing excellent quantitative agreement for a single-mode vibronic photoelectron spectrum of SO₂.

---

Spectroscopy—the measurement of light's interaction with matter—is one of the most important and precise experimental techniques for probing microscopic phenomena. The prediction of spectra via computational techniques serves as a benchmark for theoretical models of molecules, and good agreement between theory and experiment is essential if we are to truly understand chemical phenomena.

However, predicting molecular spectra remains difficult, especially for large molecules, those with significant entanglement between degrees of freedom, those that are open to an environment, or when high precision is required.¹ In particular, there are regimes where all common approximations break down, leaving critical cases without practical computational solutions. For example, the Franck–Condon approximation of vibronic (vibrational + electronic) spectroscopy states that a transition is proportional to the overlap of initial and final vibrational wavefunctions;^2,3 while often a good approximation, it can fail when the dipole moment depends on nuclear displacements. More generally, strong vibronic coupling between electronic states can lead to failures of the Born–Oppenheimer approximation and substantial nuclear-electronic entanglement.¹ Methods that include vibronic coupling are generally limited in molecule size or accuracy; for example, surface hopping uses an approximate form of the wavefunction and its evolution to reduce computational cost,⁴ whereas multiconfigurational time-dependent Hartree is numerically exact but heuristic, with an unpredictable computational cost that can be exponential in system size.^1,5

Quantum computers promise to reduce the computational cost associated with predicting molecular spectra by offering a new means of computational simulation. As in other applications of quantum computers to chemistry,^6–15 the advantage of quantum simulation lies in the ability to naturally represent complicated, entangled wavefunctions using quantum coherent degrees of freedom. Indeed, recent proposals and experiments have shown that quantum computers can predict vibronic spectra,^16–22 starting with the simulation of Franck–Condon spectra by encoding Duschinsky transformations in the displacement, squeezing, and unitary rotations of optical modes.¹⁶ This approach has been extended to thermal initial states¹⁷ and non-Condon transitions,²² and demonstrated on other quantum platforms with experimentally accessible bosonic modes, including trapped ions¹⁹ and circuit quantum electrodynamics (cQED).²¹ These approaches are examples of analog quantum simulations, where the Hamiltonian of a system of interest is mapped onto a controllable quantum system in a laboratory. In contrast, a digital quantum simulation (i.e., comprised of qubits and quantum gates) has also been proposed with a straightforward extension to include anharmonic vibrational modes.²⁰

All existing approaches to quantum simulation for molecular spectroscopy suffer the same drawback, requiring exponential resources as the number of molecular vibrational modes increases. Most methods directly simulate the intensity of every spectral line,^{16,17,19–22} whose number can reach astronomical sizes even in small molecules. For example, if 10 states in each quantised vibrational mode are optically accessible, the number of states in an N-atom molecule is 10^3N−6. Most approaches also rely on the Franck–Condon principle^{16,17,19–21} and exclude vibronic coupling, which limits their use to describing single electronic states that are energetically separated from other electronic states. The approach of Hu et al.¹⁸ predicts the spectrum of uncoupled, displaced harmonic oscillators by measuring overlaps between initial and final states; however, extending the method to more general chemical Hamiltonians requires exponential classical resources to predict the final state.

Systems with coupled electronic states and vibrational modes are well-suited for simulation on a class of analog quantum devices known as mixed qudit-boson (MQB) simulators.²⁴ These are comprised of a qudit, i.e., a d-level system with controllable transitions between all d levels, and a set of quantum oscillators or resonators making up the bosonic modes. Example architectures include trapped ions and cQED. In addition to spectroscopic prediction,^14,18,19,21 MQB devices have been proposed for the analog quantum simulation of vibrationally assisted electron transfer,^11,25 spin-boson models,²⁶ and dynamics under vibronic-coupling Hamiltonians.²⁴

Here, we describe a general approach to the quantum simulation of molecular spectra that avoids the exponential measurement requirements, instead requiring a number of measurements that is independent of molecular size. Our approach (Fig. 1e and f) uses analog quantum simulation techniques to predict the spectrum based on molecular dynamics in the time domain, unlike previous approaches that compute transition probabilities in the frequency domain (Fig. 1c and d). We show how our scheme can be efficiently implemented on MQB simulators. The core of our method is a one-to-one mapping between molecular vibrational modes and bosonic degrees of freedom, and between molecular electronic states and qudit states, both available in several quantum-computing architectures. Notably, the time-domain MQB approach allows for the inclusion of effects such as vibronic coupling, mixed states, and open quantum systems, all of which can have dramatic effects on the final spectrum, and are otherwise inaccessible in existing techniques. We validate this approach by performing a proof-of-principle experimental demonstration by simulating a vibronic spectrum of SO₂ using a trapped-ion quantum simulator, and achieving strong agreement between our experimental measurements, theory, and molecular spectroscopy measurements from the literature.

Fig. 1 Different approaches to obtain a molecular vibronic spectrum. (a) Absorption of light by a molecule (shown: SO₂) can be measured experimentally to give a spectrum (grey line in (b)²³). (c) Several quantum simulation techniques can find Franck–Condon factors^{16,17,19,21,22} (illustrated with a boson-sampling optical circuit). (d) These techniques measure the overlaps of the initial state with final eigenstates, here given by vibrational quanta n₁ and n₂ for two vibrational modes. The corresponding intensities give peak heights in the frequency domain (blue sticks in (b)). (e) We show that the molecule can instead be mapped to a time-domain MQB simulation (illustrated with a trapped-ion simulator), with the ability to include vibronic coupling, mixed initial states, and open quantum systems.²⁴ The time-domain procedure is scalable with the number of vibrational modes, because it reconstructs the spectrum directly, not via exponentially many vibronic eigenstates. (f) Measurements of the MQB simulation give the autocorrelation function a(t), whose Fourier transform, , is the vibronic spectrum (red line in (b)), including lineshapes in the presence of noise.

1. Time-domain spectroscopy on an analog quantum simulator

Time-domain simulation avoids the individual measurement of the exponentially growing number of spectroscopically relevant states.^27–30 The desired spectrum contains much less information than all the eigenstates of the molecule, and time-domain approaches can obtain it more directly.

In the frequency-domain approach,^{16,17,19–22} the spectrum is a sum of individually calculated (or sampled) contributions from each eigenstate of the molecule, in proportion to how strongly they interact with light (Fig. 1c). For example, the first-order absorption spectrum of a molecule initially in eigenstate |Ψ₀〉 with frequency ω₀ is

where ω is the frequency, |Ψ_n〉 is the eigenstate of the molecular Hamiltonian Ĥ with frequency ω_n, [small mu, Greek, circumflex] is the dipole moment operator, and ε is the polarisation of the light, both of which are vectors in three dimensions. For simplicity, in what follows we write [small mu, Greek, circumflex] = ε·[small mu, Greek, circumflex]. The computational cost is exponential in the number of modes because of the need to calculate the exponentially many contributions in the sum, even if many peaks overlap or have zero intensity. Although approaches that involve sampling of the wavefunction can reduce the number of measurements required,²¹ the number of relevant eigenstates in the spectrum (i.e., the number of peaks) still grows exponentially with the number of modes. Therefore, obtaining an accurate spectrum requires an exponential number of samples to ensure that relevant features are not missed.

By contrast, the well-established time-domain view of spectroscopy (Fig. 1e) was developed to avoid having to calculate (originally on classical computers) all the eigenstates.^27–30Eqn (1) can be rewritten³⁰ using the Fourier definition of the delta function as

where is the Fourier transform, and |Ψ_μ(t)〉 = e^−iĤt/ℏ[small mu, Greek, circumflex] e^iĤt/ℏ|Ψ₀〉, i.e., the initial state evolved backwards in time for a time t, perturbed by the dipole operator, and evolved forwards in time t. The (dipole) autocorrelation function a(t) is the overlap of the initial and time-evolved wavefunctions. Eqn (2) is a more general expression for the absorption spectrum than eqn (1), since a(t) is the dipole linear response function, 〈μ^†(0)μ(t)〉, which can also describe the evolution of mixed and non-stationary initial states, as well as open quantum systems. It is a complex, Hermitian function, meaning that its Fourier transform is real and can be calculated with only t > 0. Due to the Fourier relation between a(t) and the spectrum, the spectral resolution of the simulation is determined by the maximum propagation time and the frequency range by the size of the time steps. Therefore, the cost of the simulation (i.e., the number of samples of a(t) needed) is determined by the desired properties of the spectrum (resolution and range) and not by the number of eigenstates.

Therefore, a(t) is the only quantity whose measurement is required to generate a vibronic spectrum. It can be measured on a quantum simulator by keeping a copy of the initial wavefunction in superposition with the time-evolved wavefunction, so that their overlap can be determined at the end of the simulation via an interference measurement. Fig. 2a shows the conceptually simplest approach to do so on a quantum simulator (either digital or analog) using the phase kickback technique with an ancilla qubit,^31,32 with a circuit that is similar to more general approaches to correlation-function measurement.^33–35 An operation Â_init prepares the initial state |Ψ₀〉 on the quantum simulator (or [small rho, Greek, circumflex]₀ for a mixed initial state), which is then evolved by e^iĤt/ℏ|Ψ₀〉 (i.e., evolving with −Ĥ for a time t). A Hadamard gate places the ancilla qubit in the superposition . The dipole operator is then controlled by the ancilla to give . Forward time evolution followed by a controlled [small mu, Greek, circumflex]^† returns the |0〉 state to |Ψ₀〉, with a total state given by . The real part of a(t) is then given by the expectation value of [small sigma, Greek, circumflex]_z measurements on the ancilla, while the imaginary part can be obtained by inserting an additional R̂_x(−π/2) = e^{i[small sigma, Greek, circumflex]_xπ/4} rotation before measurement, either before or after the two Hadamard gates. Equivalently, the ancilla can be used to control both of the time evolutions instead of the dipole operators.

Fig. 2 Quantum circuits for simulating vibronic spectra. (a) General circuit to obtain the autocorrelation function from an MQB device with d qudit states (thick line) and N bosonic modes (wavy line). By adding an ancilla qubit (top row) to the MQB simulator, the real and imaginary components of a(t) can be measured as the [small sigma, Greek, circumflex]_z expectation value of the ancilla, depending on whether the dashed-line gate R̂_x(−π/2) is applied or not. Â_init prepares the initial state, and R̂_H is a Hadamard gate. (b) Simplified circuit on an MQB simulator with an additional reference qudit state. R̂^(0,ref)_H and R̂^(0,ref)_x represent gates acting on qudit states |0〉 and |ref〉, and the Hamiltonian Ĥ′ and state preparation [small mu, Greek, circumflex]′ are modified to include |ref〉. The initialisation Â⁽⁰⁾_init and backwards time propagation e^iĤ(0)t/ℏ act only on the initial electronic state. The expectation value of [small sigma, Greek, circumflex]^(0,ref)_z = |0〉〈0| − |ref〉〈ref| leads to the same measurement outputs as in (a).

The potentially impractical controlled unitary gates can, in most cases, be avoided if using an MQB simulator. We assume that the initial state |Ψ₀〉 is well described by a single electronic state, labelled |0〉, although extension to more electronic states is straightforward. To measure a(t) without an ancilla and controlled unitary gates, we add an additional electronic state to the simulator, i.e., for a simulation with d electronic states, we require a (d + 1)-level qudit for measurement of a(t). We call this additional electronic state the reference state, |ref〉, and use it to keep a copy of the initial wavefunction |Ψ₀〉 in superposition with the time-evolving wavefunction (Fig. 2b). This is achieved using the R̂^(0,ref)_H gate, which prepares the state , where |0〉 is the electronic state of the initial wavefunction. |Ψ₀〉 is prepared on the bosonic modes by a single-electronic-state operation Â⁽⁰⁾_init, after which its initial evolution is simulated with e^iĤ(0)t/ℏ, where Ĥ⁽⁰⁾ = 〈0|Ĥ|0〉 is the Hamiltonian describing evolution on only the initial electronic state. The modified operator [small mu, Greek, circumflex]′ = [small mu, Greek, circumflex] + |ref〉〈ref| acts on the original d qudit states, giving . This state then undergoes time evolution under the expanded Hamiltonian

Ĥ′ = Ĥ + Ĥ⁽⁰⁾ ⊗ |ref〉〈ref|, (3)

so that the |ref〉 component of the wavefunction returns to |Ψ₀〉 while the rest of the wavefunction propagates to |Ψ_μ(t)〉. After the final [small mu, Greek, circumflex]^†, a(t) is measured as the expectation value of [small sigma, Greek, circumflex]^(0,ref)_z = |0〉〈0| − |ref〉〈ref| (with the R̂^(0,ref)_x gate differentiating between the real and imaginary parts).

In either scheme, the simulation needs to be repeated sufficiently many times to numerically converge both the real and imaginary parts of a(t) to the required precision for a discrete number of times t. It is possible to halve the number of measurements because the hermiticity of a(t) implies that the spectrum can be reconstructed from only Re a(t) (see Appendix A).

An MQB device such as a trapped-ion or cQED system can simulate a wide range of realistic molecular Hamiltonians.²⁴ Specifically, an MQB simulator with second-order light–matter interactions³⁶ can simulate a quadratic vibronic-coupling (QVC) Hamiltonian,²⁴

which includes N free harmonic oscillators, , and expansion coefficientsdescribing perturbations of individual electronic potential energy surfaces (n = m) and vibronic couplings between them (n ≠ m). In these equations, are the nuclear normal modes q̂_j weighted by reduced mass M_j and frequency ω_j, n̂j = â^†_jâ_j are the number operators, â^†_j and â_j are the bosonic creation and annihilation operators, |n〉 are the electronic states, and c are the vibronic expansion terms.

For the simulation of a QVC model spectrum, Ĥ⁽⁰⁾ can usually be chosen to be equal to Ĥ₀. Mathematically, there are two requirements for this to be met. First, as is usually the case in molecular systems, the ground electronic state should have negligible coupling to excited ones, Ĉ_n,0 = Ĉ_0,n = 0. Second, without loss of generality, we can choose Ĉ_n,m such that Ĉ_0,0 = 0.

2. Example: one-mode model of SO₂

The simplest example of our approach involves dynamics on a single electronic state with a single vibrational mode. This type of model can be used to describe the photoelectron spectrum of the S₀ → D₀ transition of SO₂ (Fig. 3a), which has a bending mode with frequency ω_b = 2π × 12.44 THz, along which the D₀ electronic potential energy surface has a displacement relative to S₀ given by α = 1.716 (ref. 37) (in unitless, mass- and frequency-scaled coordinates). This model is described by the Hamiltonianwhere E_S0 and E_D0 are the potential energies of S₀ and D₀ at Q = 0, and we removed the (constant) zero-point energy ℏω_b/2. For approximations involved, see Sec. 5.

Fig. 3 The trapped-ion simulation of the SO₂ photoelectron spectrum. (a) Photoexcitation of SO₂ from the neutral ground electronic state (S₀) leads to the ground cationic electronic state (D₀), at a geometry that is from the D₀ minimum along the vibrational mode. (b) Time evolution on D₀ (green arrow) corresponds to bending of the SO₂ molecule. (c) In an ion-trap MQB simulator, the initial state is prepared by displacing the ion position by with a laser–ion interaction (blue arrow). (d) Time evolution of the simulator (green arrow) happens under a laser-induced spin-motion interaction that causes the effective potentials to be shifted by in the |±〉 basis (purple arrows).

Under the Condon approximation, we assume that the electronic (dipole) transition completely transfers the population from S₀ to D₀ with no effect on the nuclear coordinates, i.e., [small mu, Greek, circumflex] = |D₀〉〈S₀| + h.c. After transitioning to D₀, the wavefunction is no longer stationary and the molecule begins to vibrate (Fig. 3b).

To measure the autocorrelation function of SO₂ on an MQB simulator, we add a reference state to the model, with corresponding Hamiltonian Ĥ⁽⁰⁾ ⊗ |ref〉〈ref|, where Ĥ⁽⁰⁾ = Ĥ₀ = ℏω_bn̂. This results in a Hamiltonian with three electronic states: |S₀〉, |D₀〉, and |ref〉. However, because [small mu, Greek, circumflex] causes 100% population transfer to D₀ and Ĥ_SO2 contains no terms that couple the two electronic states, we can remove S₀ from the model completely, replacing E_D0 with ΔE = E_D0 − E_S0 to conserve the excitation energy. We can also remove the initial time evolution e^iĤ(0)t/ℏ, since |Ψ₀〉 is a stationary state of Ĥ⁽⁰⁾. The two remaining electronic states can be represented by a qubit, |0〉 = |D₀〉 and |1〉 = |ref〉, giving

which corresponds to eqn (4) with ω₁ = ω_b, , and all other coefficients equal to zero. In this representation, the corresponding dipole operator is [small mu, Greek, circumflex]′ = , because the |D₀〉 electronic state corresponds to the initial qubit state |0〉.

For experimental implementation, can be transformed into a more convenient, but equivalent, form (even though eqn (7) is already in the general form suitable for MQB simulation). First, we remove the constant term ΔE on |0〉, which is the initial excitation energy of the wavefunction from S₀ to D₀; doing so leads to a constant frequency shift of the entire spectrum, which can be restored by adding ΔE/ℏ to the frequencies after the spectrum is predicted. Next, we transform the Hamiltonian into a form that is symmetric about Q = 0. In , the minimum of the D₀ potential energy surface is at and that of the reference state is at Q = 0. These minima can be made symmetric in a displaced coordinate system obtained using the displacement operator D̂(−α/2) = e^{−α(â†+â)/2}. Finally, the two Hadamard gates can be incorporated into the time evolution. Altogether, this gives

which is a Jaynes–Cummings-type interaction that we implement experimentally below. As a result of the Hadamard transformation, |D₀〉 and |ref〉 are now |+〉 and |−〉 (Fig. 3d), where .

The overall circuit for this simulation is shown in Fig. 4a. The initialisation consists of a Â⁽⁰⁾_init = D̂(−α/2) operator, which displaces the initial vibrational ground state into the same displaced coordinates as the Hamiltonian. The time evolution consists of the unitary . The measurement of the real part of a(t) proceeds directly from the qubit state, using the operator [small sigma, Greek, circumflex]^(0,ref)_z = [small sigma, Greek, circumflex]_z. As before, the imaginary part requires the additional R̂_x(−π/2) gate.

Fig. 4 Experimental time-domain simulation of the single-mode SO₂ vibronic spectrum. (a) Quantum circuit diagram for the simulation to extract the real and imaginary components of a(t), using one qubit and one bosonic mode. (b) Experimental pulse sequence implementing the quantum circuit. (c and d) Simulations and measurements of a(t). “Theory + noise” indicates a simulation accounting for known sources of noise. (e) Comparison of the Fourier transformed data from (c) and (d) with theoretical predictions. Dots indicate peak maxima. Inset: comparison of the spectroscopically observed spectrum at 320 K (ref. 23) with frequencies shifted to give E_D0 = 0 and the ion-trap experiment. The decreasing peak spacing in the SO₂ spectrum is caused by a weak anharmonicity that is neglected in the single-mode model.

3. Experimental quantum simulation

We experimentally demonstrate the one-mode SO₂ simulation using a trapped-ion MQB quantum simulator.³⁸ Our system confines a single ¹⁷¹Yb⁺ ion in a linear Paul trap, and we encode a qubit in the ion's ²S_1/2 hyperfine ground-state manifold, |0〉 ≡ |F = 0, m_F = 0〉 and |1〉 ≡ |F = 1, m_F = 0〉. The ion's vibration in the transverse x direction encodes the molecular vibration.

The key tool for manipulating the qubit and motional wavepacket of the ion is a pair of Raman laser beams. As detailed in Sec. 5, a bichromatic laser pulse can apply a state-dependent displacement force (SDF) to the ion, described in the interaction picture by the Hamiltonian

Ĥ^I_SDF = ℏΩ_S[small sigma, Greek, circumflex]_x(â^† e^i(δt+φ) + h.c.), (9)

where the three adjustable parameters are the motional sideband interaction strength Ω_S, the detuning of the bichromatic components δ, and the motional phase φ. In the Schrödinger picture, this Hamiltonian takes the time-independent formwhere P̂ is the conjugate momentum of Q̂. This equation is identical to eqn (8) when φ = 0, δ = ω_b and Ω_S = ω_bα/2. In practice, these values need to be scaled from molecular frequencies to experimental frequencies by a constant scaling factor F whose value depends on the type of MQB simulator.²⁴

The experimental pulse sequence, Fig. 4b, describes the four stages of the simulation protocol: (i) cooling, (ii) initialisation, (iii) time evolution, and (iv) measurement.

Cooling prepares the ion in the ground state, from which further operations can be executed. First, Doppler cooling and sideband cooling are used to cool the motional degree of freedom as close as possible to the ground state (we obtained n̄ ≈ 0.05). Second, optical pumping on the internal electronic state is used to prepare the qubit state |0〉. For details, see Sec. 5.

Initialising the SO₂ simulation requires preparing the state |0〉 ⊗ |−α/2〉, where the first ket refers to the qubit and the second to the displaced motional ground state. In the three-pulse initialisation sequence, the first R̂_y(π/2) pulse rotates the qubit to the |+〉 state while the second R̂_y(−π/2) pulse returns the qubit to the |0〉 state. Between the two rotations, Ĥ_SDF is applied on resonance (δ = 0, φ = −π/2) for 0.093 ms, implementing the operation |+〉 ⊗ |0〉 → |+〉 ⊗ |−α/2〉 with α/2 = 0.858. The overall sequence produces the desired |0〉 ⊗ |−α/2〉. For details, see Sec. 5.

Time evolution is accomplished using with φ = 0. We use the scaling factor F = 1.37 × 10⁻¹⁰ to convert molecular timescales and frequencies (fs, THz) to trapped-ion timescales and frequencies (ms, kHz). Doing so gives δ = Fω_b = 2π × 1.705 kHz and Ω_S = Fω_bα/2 = 2π × 1.463 kHz. To obtain the time trace of a(t), we measure its value at 200 different times t by repeating the experiment with the duration of unitary evolution under Ĥ_SDF varying between 0 and 2 ms, corresponding to molecular durations of up to 274 fs.

The final step in the simulation is the measurement of a(t), which is carried out by measuring the qubit in the computational basis (for details, see Sec. 5). Reading out the imaginary part of a(t) requires the additional R̂_x(−π/2) gate on the qubit following the displacement in the initialisation step.

The full experimental sequence above is repeated 500 times for each duration t of the simulated time evolution in order to converge the measurement observables.

Fig. 4c–e shows the agreement between our experimental results and theoretical predictions. The predicted and measured a(t) are shown in Fig. 4c and d. The theoretical curve is calculated as shown by the circuit in Fig. 4a, using eqn (8) for the time evolution. To give a non-zero linewidth in the theoretical spectrum, the predicted a(t) was multiplied by an exponential decay of 6 ms (corresponding to 822 fs at the molecular timescale). Fig. 4e shows the agreement between the predicted and measured spectra, i.e., the Fourier transforms of the theoretical and experimental a(t).

Despite the good overall agreement between theory and experiment, there are minor differences between the two, most of which can be explained by the presence of noise in the trapped-ion simulator. Points of difference include a drift in the a(t) signal with simulation time, discrepancies in peak heights, and asymmetric lineshapes. Most of the discrepancies can be accounted for by adding a model of experimental noise to our theory. This simulation includes a linear frequency drift with the Hamiltonian and uses an initial thermal state with an average motional state population n̄. The evolution of the density operator [small rho, Greek, circumflex] obeys the master equation

where is a Lindblad superoperator acting on [small rho, Greek, circumflex] for the jump operator L̂. The dissipation is described by a motional heating rate γ_h and a pure motional dephasing lifetime τ_d. The noise parameters are fitted to the measured a(t) using non-linear least squares, giving d_δ = 2π × 52 Hz ms⁻¹, n̄ = 0.061, γ_h = 43 s⁻¹ and τ_d = 110 ms. This yields an effective motional coherence time of 33 ms, in agreement with the experimentally measured value (see Sec. 5). The simulation that includes noise agrees better with the experiment, accounting for the signal drift in a(t) using the frequency drift (Fig. 4c and d), as well as linewidths and peak heights using n̄, γ_h and τ_d (Fig. 4e).

4. Discussion

Our approach has two types of advantages over existing proposals for the quantum simulation of spectroscopy: those which result from our theoretical framing of the simulation, and those that come from our choice of experimental platform.

The advantages of our theoretical framework stem from framing molecular spectroscopy in the time domain. Our work mirrors the development of classical computing methods in the time domain, which greatly simplified the calculation of spectra for high-dimensional systems without the need to resolve eigenvalues.^29,30 In the time domain, spectroscopy is an initial-value problem, rather than an eigenvalue problem where the number of solutions grows exponentially with system size. This reframing leads to two distinct advantages over competing proposals for the analog quantum simulation of spectroscopy: scalability and generality.

The scalability of our approach stems from the exponentially reduced number of measurements needed to predict the spectrum. In frequency-domain approaches, the number of eigenvalues (i.e., peaks) grows exponentially with the number of modes (i.e., with molecule size), each of which needs to be sampled to determine its intensity. Even if the number of eigenvalues is truncated on an ad hoc basis, the number of significant eigenvalues grows rapidly. For example, the technique employed by Shen et al.¹⁹ involves a sequence of laser pulses to project the population of each multimode motional state |n₁, n₂, …, n_N〉 onto the qubit state population; if each mode occupation is truncated at n_max, the computational cost scales exponentially as n_max^N. By contrast, in time-domain approaches such as ours, the number of measurements required is independent of system size. Instead, the number of measurements is determined by the desired frequency range and resolution of the spectrum, which are the inverses of the time step and the total simulation time, respectively. Perfect spectral resolution is not necessary for characterising a spectrum, since measured spectra of even modestly sized molecules have broad features of overlapping peaks, especially when environmental effects, strong coupling, and limited measurement resolution are considered. Therefore, the cost of a time-domain simulation is determined by experimentally relevant parameters (spectral range and resolution), rather than the size of the underlying Hamiltonian.

As for generality, our method can be used to predict the spectroscopy of any chemical system due to the fully general relationship between σ(ω) and shown in eqn (2). The observable we measure, a(t) = 〈μ^†(0)μ(t)〉,^28,30 is defined without an eigenstate expansion and can, in principle, be efficiently measured on any quantum simulator, including those simulating open quantum systems, vibronic couplings, nonlinearities, or non-Condon effects.

The simulation of open quantum system is the most striking example of the generality of the time-domain approach. Introducing controlled noise into a simulation allows an MQB simulator to simulate environmental effects such as peak broadening²⁴ with the same number of measurements of a(t). By contrast, frequency-domain approaches typically fail on open systems, which no longer have discrete eigenstates that can be measured one by one. Furthermore, in the time-domain approach, initial conditions can likewise include mixed states such as thermal states, allowing simulations of spectroscopy of molecules at finite temperature without additional experimental resources. By contrast, doing so in the frequency domain requires multiple experiments¹⁶ or doubling the simulator size.¹⁷

The generality of our approach also extends to the ability to include vibronic couplings (eqn (4), n ≠ m) and non-Condon effects. Vibronic coupling is ubiquitous in UV-visible spectroscopy, and, in the time domain, any approach able to simulate dynamics with vibronic coupling²⁴ can also predict the spectrum. In contrast, all frequency-domain analog simulation approaches use Duschinsky transformations to prepare the initial state in the vibrational basis of the final electronic state.^{16,17,19,21,22} This is a powerful technique for energetically separated electronic states, but cannot describe vibronic coupling. In addition, our approach has the potential to simulate non-Condon effects (the dependence of [small mu, Greek, circumflex] on nuclear coordinates) with no additional experimental resources. In comparison, non-Condon effects require multiple simulations for current frequency-domain approaches.²²

Turning to the advantages of our experimental demonstration, the use of MQB simulators offers a significant reduction in required quantum resources over digital quantum simulation approaches. Both analog MQB and digital simulations of dynamics—either of which could be used for our time-domain simulation of spectroscopy—require resources that scale linearly with system size;²⁴ however, the resource cost per mode is considerably higher in digital approaches, where many qubits would be needed to adequately represent a single vibrational mode. For example, our demonstration of SO₂ with a single trapped ion is equivalent to a digital encoding with at least 6 qubits (assuming 32 Fock states per mode). MQB simulation also comes with a time advantage, since the harmonic motion native to an MQB simulator would require many gates to implement digitally. Finally, the relative frequencies of qudit and bosonic levels on a trapped-ion MQB simulator lead to relative electronic and vibrational noise strengths that are similar to those in molecules, which can be exploited for simulating open quantum systems and which would not be the case for digital simulation, where all qubits are ordinarily assumed to experience comparable noise.

Our experimental results demonstrate the availability of the essential MQB building blocks in existing trapped-ion technology. Our demonstration is a proof of principle, and further work is necessary to reach a scale where it could outperform classical computers by simulating larger, more complicated molecules. Indeed, eqn (6) captures the absorption of a single displaced harmonic oscillator, and can be solved analytically.²⁸ Nevertheless, we see a clear path towards integrating complicated, non-linear vibronic-coupling Hamiltonians into analog spectroscopy simulations of molecules large enough to be intractable on classical computers. All of the necessary components for a more general, QVC simulation²⁴ have already been demonstrated in trapped-ion systems, including those with qudits and more vibrational modes. In particular, our approach can incorporate more modes using an additional Raman interaction for every mode, which can be efficiently implemented with the same experimental setup by interleaving different interactions using trotterisation,^24,39 an established technique in trapped ions.⁴⁰ Furthermore, higher-order terms in the vibronic-coupling Hamiltonian—responsible for anharmonicities and nonlinearities that are particularly difficult to simulate classically—can be incorporated into both the simulated Hamiltonian and the initial state preparation using techniques such as motional squeezing and mode-mixing^19,36,41 or using ancillary ions.⁴² In addition, more electronic states could be simulated using recent experimental advances in trapped-ion qudits.^43,44

As in any analog simulation—quantum or classical—the absence of error correction means that excessive noise can lead to inadequate results. However, our demonstration shows that, despite the lack of error correction, existing trapped-ion technology can provide remarkable agreement with theoretical predictions. Moreover, our analysis of experimental noise sources shows that most of the imperfections in our simulation can be accounted for, making it clear which experimental improvements are necessary if greater accuracy is desired. When simulating larger molecules or those open to the environment—where classical chemical simulations struggle the most—the presence of noise in the simulator becomes a feature that can be controlled. For instance, the inset in Fig. 4e shows that our simulation gives narrower peaks than are measured in a high-precision spectroscopic experiment at 320 K, meaning that we would have to inject additional noise to fully reproduce the spectroscopic observations.

Our method's most likely path to quantum advantage is by simulating a combination of effects that make classical simulation challenging, including vibronic coupling and a finite-temperature bath. We previously outlined²⁴ the favourable resource scaling that our approach can achieve for such systems. For example, a full-dimensional quadratic vibronic coupling model of pyrazine is a challenging system for classical computers.⁴⁵ It involves 24 modes and two electronic states, meaning that our technique could simulate its spectrum with 8 trapped ions.²⁴ For comparison, previous ion-trap experiments have controlled interactions of as many as 20 ions,⁴⁶ putting our example within existing experimental feasibility.

Overall, our approach shows the remarkable advantages—both theoretical and experimental—of using the time-domain representation of spectroscopy in analog quantum simulation. By using a one-to-one mapping between simulated molecular vibrations and bosonic modes in a quantum simulator, our scheme provides an exponential improvement in resource requirements compared to existing quantum methods using frequency-domain simulations. In addition, it straightforwardly generalises to simulations of more complicated chemical systems or those open to the environment. Our proof-of-principle demonstration of the simplest example of our approach showcases all of the necessary experimental building blocks, giving us confidence that foreseeable developments in quantum technology will allow larger simulations of molecular spectroscopy to occur in the near term, including of molecules that could not be simulated on any classical computer.

5. Methods

5.1 Hamiltonian for SO₂

We selected SO₂ for our proof-of-principle demonstration because it is well described by the Hamiltonian in eqn (6). This Hamiltonian includes several simplifications to the physics of SO₂, all of which could be relaxed in a more detailed model. First, the model assumes equal frequency ω_b = 2π × 12.44 THz for both electronic states, although the ground electronic state has a slightly higher vibrational frequency of 2π × 15.55 THz.³⁷ Second, it neglects higher-lying electronic states (and vibronic couplings to them) due to their large energetic separation from S₀ and D₀. Finally, only the bending mode is included because the displacements between S₀ and D₀ are small along the other two vibrational modes, the symmetric and asymmetric stretches, being α = −0.026 and α = 0, respectively.³⁷

5.2 Numerical methods

Time evolution for all theoretical simulations was simulated using the master equation solver in QuTiP.⁴⁷ Curve fitting was performed with the Levenberg–Marquardt algorithm in SciPy.⁴⁸ All Fourier transforms were calculated using the Fourier–Padé approximation (FPA). Unlike a discrete Fourier transform (DFT), the FPA finds coefficients of continuous rational functions that approximate the Fourier transform from discrete data. Relative to DFT, it has a faster convergence for features of a spectrum generated from a finite-length time series.⁴⁹ To avoid poles due to the rational function expansion, a(t) was multiplied by a frequency shift e^iθt with θ = 2π × 8 kHz before the Fourier transform, which was later corrected by shifting the spectrum by −θ. Testing values of θ in the range 2π × 7–9 kHz showed no change in the spectrum, indicating no poles or spurious peaks.

5.3 Ion trap characteristics

The motional mode frequencies of our ion-trap MQB simulator are {ω_x, ω_y, ω_z} = 2π × {1.31, 1.45, 0.5} MHz. Using Ramsey-type measurements, we find native (uncorrected) coherence times of for the qubit and an effective 35 ms for the transverse motional mode along x. We measured a motional mode rethermalisation (heating) rate of 0.2 quanta s⁻¹ in the absence of laser light.

5.4 Coherent operations

The qubit states and motional modes are manipulated by stimulated Raman transitions driven by a 355 nm pulsed laser.^50,51 Two separately controllable 355 nm laser beams are focused on the ion's location; they form an orthogonal geometry such that only the x and y transverse modes of motion can be driven. The applied laser light is controlled using acousto-optic modulators (AOM), driven by radio-frequency (RF) signals. Changing the RF signal amplitude, frequency and phase allows the tuning of Ω_S, δ and φ in eqn (9), respectively. These parameters are controlled using RF signal generators as part of the experiment control system.⁵²

5.5 Cooling and qubit state preparation

A laser, red-detuned from the ²S_1/2 → ²P_1/2 transition near 369.5 nm, is used to Doppler cool the motional modes to a thermal state. The Doppler-cooled ion temperature is further reduced with pulsed resolved-sideband cooling.^53,54 The qubit state is prepared by optically pumping to |0〉 with a beam resonant with the ²S_1/2 → ²P_1/2 transition. The ²P_1/2 state has a non-zero probability of decaying which results in population in ²D_3/2 and ²F_7/2 states; this leaves the ion state off-resonant from the Doppler cooling laser. Therefore, additional light fields at 935 nm and 760 nm are used to depopulate these states and return the ion to the cooling cycle.^55,56

5.6 SDF interaction

Using a two-tone RF signal to drive one of the 355 nm laser beams, a bichromatic light field is generated, which can simultaneously drive the red and blue motional sideband transitions adjacent to the qubit resonance. This bichromatic light creates the state-dependent force described by Ĥ_SDF, which acts in the [small sigma, Greek, circumflex]_x eigenbasis of the qubit.⁵⁷ Through a prior rotation of the internal qubit state into the eigenbasis, the SDF enacts a coherent displacement of the motional wavepacket in a particular direction in phase space.⁵⁸ For details of the calibration procedure, see Appendix B.

5.7 Measurement

The qubit state is measured using state-dependent fluorescence, induced by a laser beam resonant with the ²S_1/2|1〉 → ²P_1/2|F = 0〉 transition.^55,59 The |0〉 state is off-resonant from the detection beam by roughly 2π × 14.75 GHz and only scatters an average of 0.1 photon counts during a detection window. Measurement of the |1〉 state produces 16.7 photon counts. Using a threshold of 4 counts we are able to infer the qubit state with a detection error of 1.2%.⁵⁶

Appendices

A. Simulation of positive-frequency spectra

The most general approach to measure a(t) with an MQB simulator requires two experiments at each time t: one to measure Re a(t), and another to measure Im a(t). However, the total number of measurements can be halved by exploiting the hermiticity of a(t) to only measure Re a(t).

The autocorrelation is a complex, Hermitian function, meaning a(t) = a*(−t) and thus Re a(t) = Re a(−t) is even and Im a(t) = −Im a(−t) is odd. The Fourier transform maintains parity and multiplies odd functions by i. Thus, is real and even, is real and odd, and the spectrum is the sum of the two, as shown in Fig. 5a. If all features of the spectrum appeared at strictly positive frequencies, it would mean that for ω > 0, and the two components would cancel for ω < 0. Therefore, only one of the two would be necessary to produce the spectrum.

Fig. 5 Fourier transforms of different components of the autocorrelation function a(t), using the D₀ photoelectron spectrum of SO₂ as an example.³⁷ (a) Fourier transforms of the real and imaginary components of a(t) and their sum. (b) Fourier transforms of real and imaginary components of a(t) simulated with a frequency shift of ε = 2π × 30 THz and with frequencies corrected by −ε.

In cases where features of the spectrum (i.e., peaks and their linewidths) appear at ω ≤ 0, we can introduce a frequency shift ε such that all features appear at ω + ε > 0. This corresponds to subtracting a constant frequency term from the reference state, or

Fig. 5b shows the spectrum and its components when a frequency shift ε is applied. After the Fourier transform, the correct frequencies are restored by subtracting ε, i.e., translating the spectrum by ε in frequency. Because the real and imaginary components are equal to for ω > −ε, experimental measurement of Re a(t) is sufficient to obtain the spectrum.

Choosing ε requires an estimate of the frequency and linewidth of the lowest-frequency peak in the spectrum. For some of Hamiltonians, such as QVC Hamiltonians with weak vibronic coupling near Q = 0, the lowest frequency can be estimated as the zero-point energy of the lowest excited electronic state. The linewidth can be estimated from the noise conditions of the simulator using previous experimental measurements. In the absence of a good estimate of the lowest frequency or its linewidth, low-resolution experiments (i.e., with a short propagation time) can be used to increase ε until all features of the spectrum are well separated from ω = 0.

B. Experimental calibration

Operations with the trapped ion's degrees of freedom are driven by coherent laser interactions, and in this appendix we describe how the relevant laser parameters were calibrated.

In our experiment, three main laser interactions are used: carrier, and red- and blue-sideband transitions. Their Hamiltonians in the interaction picture are

where the Rabi frequency Ω quantifies the coupling strength between the qubit states and the applied laser light, and [small sigma, Greek, circumflex]_± = ([small sigma, Greek, circumflex]_x ∓ i[small sigma, Greek, circumflex]_y)/2. The light–ion interaction imprints a phase relationship, which can be controlled by the parameter ϕ_C, allowing for qubit state rotations around the x or y axis of the Bloch sphere. The sideband interactions in both eqn (B2) and (B3) are similar to the carrier interaction, but they also contain the bosonic ladder operators â and â^† for a single motional mode. Their interaction strength is scaled by the Lamb–Dicke parameter , where λ is the laser wavelength and m is the ion's mass. η is included in eqn (9)via the sideband Rabi frequency Ω_S = ηΩ/2. Similar to the carrier interaction, the sideband interactions have associated phases ϕ_R and ϕ_B. The detuning δ is a symmetric frequency offset from resonant motional sidebands.

Simultaneously driving the red and blue sidebands implements a bichromatic pulse described by

which depends on the spin phase ϕ_S = (ϕ_R + ϕ_B)/2 and the motional phase ϕ_M = (ϕ_B − ϕ_R)/2. Experimentally, ϕ_S and ϕ_M can have an arbitrary, constant offset, and they are adjusted so that eqn (B4) has equivalent phase relationships to eqn (9). The bichromatic pulse is used to both initialise the ion's motional wavepacket into a displaced coherent state and to drive the simulated molecular time evolution. In the following, we describe how the physical parameters that make up the bichromatic fields are calibrated, namely their frequencies, amplitudes, phases, and durations.

Laser frequencies are calibrated as described in Fig. 6, allowing the control of the detuning δ in eqn (B4). The sequence consists of two SDF pulses of equal duration, with a π phase shift between them.⁵⁸ If the bichromatic fields are resonant with the motional sidebands, the spin and motion are disentangled after the SDF evolutions (see Fig. 6b). However, in the presence of motional frequency offsets, the spin and motion remain entangled (see Fig. 6c), leading to measurements of a partially mixed qubit state. Static frequency offsets are therefore calibrated by symmetrically varying the bichromatic fields' frequencies and measuring qubit population in the [small sigma, Greek, circumflex]_z basis. Using this method, we are able to calibrate the frequencies of the bichromatic fields to within 50 Hz. Furthermore, calibrations of the motional frequency are scheduled every 5 minutes to mitigate the effects of drift.

Fig. 6 Calibration of symmetric frequency detuning. (a) Pulse sequence for calibrating the symmetric detuning δ with respect to the motional sidebands. (b) At resonance, the spin states move in straight trajectories and return to the origin, leading to final qubit state |0〉. (c) In the presence of a frequency offset δ ≠ 0 during the first SDF pulse, the spin states |+〉 and |−〉 follow a circular trajectory, completing one loop in time t = 2π/δ. For a pulse time τ < 2π/δ, the spin states will be displaced by a small amount from the origin. The phase offset of ϕ_M = π for the second SDF pulse again displaces along a similar trajectory, but with the centre of a circle displaced due to the first pulse. The residual displacement along the imaginary axis results in measuring significant population in the |1〉 state. (d) The resonance condition is found by sweeping the offset detuning and observing where the population measurement is closest to 0. This data was collected by repeating the pulse sequence 300 times with an SDF pulse time of 300 μs. The error bars correspond to uncertainty due to quantum projection noise.

The amplitudes of the two tones of the bichromatic field, corresponding to the red and blue sidebands, are calibrated independently. An imbalance leads to an unwanted AC-Stark shift affecting the qubit frequency and residual coupling between the qubit state and the bosonic modes. We calibrate the amplitudes by measuring their Rabi frequencies from Rabi oscillations. The Rabi frequencies are equalised by adjusting the respective RF signal's amplitude. In practice, we find a difference of <3% between Rabi frequencies of two tones and negligible variation on the timescale of an experiment.

We use the pulse sequence depicted in Fig. 7a to calibrate the phases of the bichromatic fields such that the displacement operator enacted by the SDF interaction acts in the [small sigma, Greek, circumflex]_x eigenbasis of the qubit. In the calibration, we set ϕ_R = ϕ_B and vary them simultaneously. We find the phase for which the spin state is unchanged, indicating that the overall operation is acting on the [small sigma, Greek, circumflex]_x eigenbasis.

Fig. 7 Bichromatic phase calibration. (a) The pulse sequence used to calibrate the phase of the bichromatic field ϕ_S with respect to the carrier phase ϕ_C = 0 set by the first π/2 pulse. (b) An example of the sine wave traced out by a full cycle of ϕ_S. The phase is varied and the operation of the SDF on its eigenstate [small sigma, Greek, circumflex]_x is indicated by measurement of the qubit state having a near-zero population in |1〉. In this case, the phases ϕ_R and ϕ_B are set to 1.59 radians. The pulse sequence was repeated 100 times for each phase setting.

Once the frequency, amplitude, and phase are calibrated with sufficient precision, the duration of the pulse can be calibrated to set the displacement to the desired value. The pulse sequence for the calibration is shown in Fig. 8a. The magnitude |β| of the displacement operator D̂(β) is determined by the Rabi frequency and the SDF pulse duration; for simplicity, we do not change the Rabi frequency and instead adjust only the pulse duration. After the displacement operation, the magnitude is estimated by observing the change in qubit state population after driving a blue-sideband transition.⁶⁰ The resulting spin probability follows

where is the Rabi frequency of the blue sideband for Fock state k, L_k¹(x) is the Laguerre polynomial in x of order k, and the fitting parameter ζ introduces amplitude damping that might be present due to motional state decoherence. The observed oscillations are fitted to extract |β| (see Fig. 8b–d), allowing the bichromatic fields' duration to be varied to correct the displacement magnitude.

Fig. 8 Calibrating the duration of the bichromatic pulses. (a) The pulse sequence used to calibrate the displacement operation. Longer SDF pulses create larger displacements D̂(β). The displacement distance can be inferred from fitting the time evolution of a blue-sideband-driven population to eqn (B5). Examples in (b)–(d) correspond to applying D̂(β) for t = 0.05, 0.15 and 0.4 ms with Ω_S = 2π × 0.850 kHz. Each data point is an average of 200 repetitions.

Author contributions

RJM designed the quantum algorithm and performed theoretical simulations. TN performed the experiment. TN, TFWR, CHV, ADR, MJM, and TRT developed the experimental system. MAC simulated experimental noise. MJB, TRT, and CH supervised the experimental work. IK supervised the theoretical work. RJM and TN drafted the manuscript. All authors refined the project design and edited the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We were supported by a Westpac Scholars Trust Research Fellowship (I. K.), 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), by the Sydney Quantum Academy (T. R. T.), by the University of Sydney Nano Institute, and by the Australian National Computational Infrastructure.

References

1. Conical Intersections: Electronic Structure, Dynamics & Spectroscopy, ed. W. Domcke, D. R. Yarkony, and H. Köppel, Adv. Ser. Phys. Chem., World Scientific, Singapore, 2004, vol. 15.
2. J. Franck, Trans. Faraday Soc., 1926, 21, 536–542.
3. E. Condon, Phys. Rev., 1926, 28, 1182–1201.
4. J. C. Tully, J. Chem. Phys., 1990, 93, 1061.
5. G. A. Worth, H.-D. Meyer, H. Köppel, L. S. Cederbaum and I. Burghardt, Int. Rev. Phys. Chem., 2008, 27, 569.
6. L. Lamata, A. Mezzacapo, J. Casanova and E. Solano, EPJ Quantum Technol., 2014, 1, 9.
7. 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.
8. J. R. McClean, J. Romero, R. Babbush and A. Aspuru-Guzik, New J. Phys., 2016, 18, 023023.
9. J. Argüello-Luengo, A. González-Tudela, T. Shi, P. Zoller and J. I. Cirac, Nature, 2019, 574, 215.
10. A. Di Paolo, P. K. Barkoutsos, I. Tavernelli and A. Blais, Phys. Rev. Res., 2020, 2, 033364.
11. F. Schlawin, M. Gessner, A. Buchleitner, T. Schaetz and S. S. Skourtis, PRX Quantum, 2021, 2, 010314.
12. K. Seetharam, D. Biswas, C. Noel, A. Risinger, D. Zhu, O. Katz, S. Chattopadhyay, M. Cetina, C. Monroe, E. Demler and D. Sels, arXiv, 2021, preprint, arXiv:2109.13298,  DOI:10.48550/arXiv.2109.13298.
13. S. M. Young and M. Sarovar, Phys. Rev. Res., 2023, 5, 013027.
14. 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.
15. P. Richerme, M. C. Revelle, D. Saha, S. A. Norrell, C. G. Yale, D. Lobser, A. D. Burch, S. M. Clark, J. M. Smith, A. Sabry and S. S. Iyengar, J. Phys. Chem. Lett., 2023, 14, 7256.
16. J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean and A. Aspuru-Guzik, Nat. Photonics, 2015, 9, 615.
17. J. Huh and M.-H. Yung, Sci. Rep., 2017, 7, 7462.
18. L. Hu, Y.-C. Ma, Y. Xu, W.-T. Wang, Y.-W. Ma, K. Liu, H.-Y. Wang, Y.-P. Song, M.-H. Yung and L.-Y. Sun, Sci. Bull., 2018, 63, 293.
19. Y. Shen, Y. Lu, K. Zhang, J. Zhang, S. Zhang, J. Huh and K. Kim, Chem. Sci., 2018, 9, 836.
20. N. P. D. Sawaya and J. Huh, J. Phys. Chem. Lett., 2019, 10, 3586.
21. 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.
22. H. Jnane, N. P. D. Sawaya, B. Peropadre, A. Aspuru-Guzik, R. Garcia-Patron and J. Huh, ACS Photonics, 2021, 8, 2007.
23. D. M. P. Holland, M. A. MacDonald, M. A. Hayes, P. Baltzer, L. Karlsson, M. Lundqvist, B. Wannberg and W. von Niessen, Chem. Phys., 1994, 118, 317.
24. 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.
25. D. J. Gorman, B. Hemmerling, E. Megidish, S. A. Moeller, P. Schindler, M. Sarovar and H. Haeffner, Phys. Rev. X, 2018, 8, 011038.
26. A. Lemmer, C. Cormick, D. Tamascelli, T. Schaetz, S. F. Huelga and M. B. Plenio, New J. Phys., 2018, 20, 073002.
27. R. G. Gordon, J. Chem. Phys., 1965, 43, 1307.
28. L. S. Cederbaum and W. Domcke, J. Chem. Phys., 1976, 64, 603.
29. E. J. Heller, J. Chem. Phys., 1978, 68, 2066.
30. E. J. Heller, Acc. Chem. Res., 1981, 14, 368.
31. D. S. Abrams and S. Lloyd, Phys. Rev. Lett., 1999, 83, 5162.
32. A. Aspuru-Guzik, A. D. Dutoi, P. J. Love and M. Head-Gordon, Science, 2005, 309, 1704.
33. B. M. Terhal and D. P. DiVincenzo, Phys. Rev. A, 2000, 61, 022301.
34. R. Somma, G. Ortiz, E. Knill and J. Gubernatis, Int. J. Quant. Inf., 2003, 1, 189.
35. J. S. Pedernales, R. D. Candia, I. L. Egusquiza, J. Casanova and E. Solano, Phys. Rev. Lett., 2014, 113, 020505.
36. O. Katz, M. Cetina and C. Monroe, PRX Quantum, 2023, 4, 030311.
37. C.-L. Lee, S.-H. Yang, S.-Y. Kuo and J.-L. Chang, J. Mol. Spectrosc., 2009, 256, 279.
38. A. R. Milne, C. L. Edmunds, C. Hempel, F. Roy, S. Mavadia and M. J. Biercuk, Phys. Rev. Appl., 2020, 13, 024022.
39. S. Lloyd, Phys. Rev. Lett., 1995, 75, 346.
40. B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Gerritsma, F. Zähringer, P. Schindler, J. T. Barreiro, M. Rambach, G. Kirchmair, M. Hennrich, P. Zoller, R. Blatt and C. F. Roos, Science, 2011, 334, 57.
41. K. Marshall and D. F. V. James, Appl. Phys. B, 2016, 123, 26.
42. 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.
43. P. J. Low, B. M. White, A. A. Cox, M. L. Day and C. Senko, Phys. Rev. Res., 2020, 2, 033128.
44. M. Ringbauer, M. Meth, L. Postler, R. Stricker, R. Blatt, P. Schindler and T. Monz, Nat. Phys., 2022, 18, 1053.
45. A. Raab, G. A. Worth, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1999, 110, 936.
46. 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.
47. J. R. Johansson, P. D. Nation and F. Nori, Comput. Phys. Commun., 2013, 184, 1234.
48. 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.
49. A. Bruner, D. LaMaster and K. Lopata, J. Chem. Theory Comput., 2016, 12, 3741.
50. D. Hayes, D. N. Matsukevich, P. Maunz, D. Hucul, Q. Quraishi, S. Olmschenk, W. Campbell, J. Mizrahi, C. Senko and C. Monroe, Phys. Rev. Lett., 2010, 104, 140501.
51. R. Islam, W. C. Campbell, T. Choi, S. M. Clark, C. W. S. Conover, S. Debnath, E. E. Edwards, B. Fields, D. Hayes, D. Hucul, I. V. Inlek, K. G. Johnson, S. Korenblit, A. Lee, K. W. Lee, T. A. Manning, D. N. Matsukevich, J. Mizrahi, Q. Quraishi, C. Senko, J. Smith and C. Monroe, Opt. Lett., 2014, 39, 3238.
52. S. Bourdeauducq, Whitequark, R. Jördens, D. Nadlinger, Y. Sionneau and F. Kermarrec, Zenodo, 2021,  DOI:10.5281/zenodo.6619071.
53. F. Diedrich, J. C. Bergquist, W. M. Itano and D. J. Wineland, Phys. Rev. Lett., 1989, 62, 403.
54. C. Monroe, D. M. Meekhof, B. E. King, S. R. Jefferts, W. M. Itano, D. J. Wineland and P. Gould, Phys. Rev. Lett., 1995, 75, 4011.
55. S. Olmschenk, K. C. Younge, D. L. Moehring, D. N. Matsukevich, P. Maunz and C. Monroe, Phys. Rev. A, 2007, 76, 052314.
56. C. L. Edmunds, T. R. Tan, A. R. Milne, A. Singh, M. J. Biercuk and C. Hempel, Phys. Rev. A, 2021, 104, 012606.
57. K. Mølmer and A. Sørensen, Phys. Rev. Lett., 1999, 82, 1835.
58. A. R. Milne, C. Hempel, L. Li, C. L. Edmunds, H. J. Slatyer, H. Ball, M. R. Hush and M. J. Biercuk, Phys. Rev. Lett., 2021, 126, 250506.
59. B. B. Blinov, D. Leibfried, C. Monroe and D. J. Wineland, Quantum Inf. Process., 2004, 3, 45.
60. D. Leibfried, D. M. Meekhof, B. E. King, C. Monroe, W. M. Itano and D. J. Wineland, Phys. Rev. Lett., 1996, 77, 4281.

---

Footnotes

† These authors contributed equally to this work.

‡ Present address: Department of Chemistry, Dalhousie University, Halifax, NS B3H 4R2, Canada.

---