Quantum algorithm for calculating molecular vibronic spectra

TL;DR

This paper introduces a quantum algorithm for calculating molecular vibronic spectra, effectively handling anharmonicity and preserving state information, which are challenging for classical methods, with potential applications in material design.

Contribution

The authors develop a quantum algorithm that naturally includes vibrational anharmonicity and maintains state information post-measurement, advancing quantum spectral calculations in chemistry.

Findings

Numerical analysis of truncation errors in harmonic approximation for triatomic molecules.

Application of the algorithm to anharmonic spectra of sulfur dioxide.

Potential for future material design and broader spectral calculations.

Abstract

We present a quantum algorithm for calculating the vibronic spectrum of a molecule, a useful but classically hard problem in chemistry. We show several advantages over previous quantum approaches: vibrational anharmonicity is naturally included; after measurement, some state information is preserved for further analysis; and there are potential error-related benefits. Considering four triatomic molecules, we numerically study truncation errors in the harmonic approximation. Further, in order to highlight the fact that our quantum algorithm's primary advantage over classical algorithms is in simulating anharmonic spectra, we consider the anharmonic vibronic spectrum of sulfur dioxide. In the future, our approach could aid in the design of materials with specific light-harvesting and energy transfer properties, and the general strategy is applicable to other spectral calculations in chemistry and condensed matter physics.

Equations

$$\overrightarrow{q^{\prime }}=\text{S}\vec{q}+\vec{d}$$

$$FCP(\omega )=\sum _{\{|f_i⟩\}}|⟨0|f_i⟩{|}^2\delta (\omega -{\omega }_i),$$

$${\vec{q}}_B={\mathbf{\Omega }}_B\mathbf{S}{\mathbf{\Omega }}_A^{-1}{\vec{q}}_A+{\mathbf{\Omega }}_B\vec{d}$$

$${\vec{p}}_B={\mathbf{\Omega }}_B^{-1}\mathbf{S}{\mathbf{\Omega }}_A{\vec{p}}_A,$$

$${\mathbf{\Omega }}_s=diag([{\omega }_{s1},...,{\omega }_{sM}]{)}^{\frac{1}{2}}$$

$$H_B=\frac{1}{2}\sum _k^M{\omega }_{Bk}(q_{Bk}^2+p_{Bk}^2),$$

$${\tilde{a}}_i^{†}=\sum _{l=1}^{L_{max}}\sqrt{l}|l⟩⟨l-1|$$

$$H=H_0+\sum _{ijk}{k}_{ijk}{q}_i{q}_j{q}_k+...$$

$$|0⟩=\sum _i{c}_i|{\psi }_i⟩$$

Full text

TOC Graphic

Calculating the absorption spectrum of molecules is a common and important problem in theoretical chemistry, as it aids both the interpretation of experimental spectra and the a priori design of molecules with particular optical properties prior to performing a costly laboratory synthesis. Further, in many molecular clusters and systems, absorption and emission spectra of molecules are required for calculating energy transfer rates [1]. The widespread use of mature software that solves the vibronic problem is one indication of its relevance to chemistry [2, 3, 4, 5].

Many quantum algorithms have been proposed for practical problems in chemistry, chiefly for solving the fermionic problem of determining the lowest-energy configuration of $N_{e}$ electrons, given the presence of a set of clamped atomic nuclei [6, 7, 8, 9, 10, 11, 12, 13, 14]. However, for many chemical problems of practical interest, solving the ground-state electronic structure problem is insufficient. To calculate exact vibronic spectra, for instance, an often combinatorially scaling classical algorithm must be implemented after the electronic structure problem has been solved for many nuclear positions [15, 16, 17, 18, 19, 20].

In this work, we propose an efficient quantum algorithm for calculating molecular vibronic spectra, within the standard quantum circuit model. A quantum algorithm to solve this problem, for implementation on a boson sampling machine [20, 21], was previously proposed and demonstrated experimentally [22, 23], but to our knowledge no one has previously developed an algorithm for the universal circuit model of quantum computation, nor (more importantly) one that can efficiently include vibrational anharmonic effects while calculating the full spectrum. We are also aware of unpublished work that studies the connection between quantum phase estimation and sampling problems [24].

Other related previous work includes quantum algorithms for calculating single Franck-Condon factors [25, 26] or low-lying vibrational states [27], algorithms for simulating vibrational dynamics [26, 28], and an experimental photonics implementation [28] that simulated several processes related to molecular vibrations in molecules. Though these four works simulate vibrational effects, they do not address the problem of efficiently solving the full vibronic spectrum despite the presence of an exponential number of relevant vibrational states, which is the focus of this work. Note that in this work we use the term “classical” solely to refer to algorithms that run on classical computers for solving the quantum problem of calculating vibronic spectra; we are not referring to methods where nuclear degrees of freedom are approximated with Newtonian physics.

In order to calculate a vibronic spectrum, one needs to consider the transformation between two electronic potential energy surfaces (PESs). The hypersurfaces may be substantially anharmonic, making accurate classical calculations beyond a few atoms impossible [29, 30, 31, 32]. In order to introduce our approach, we begin by assuming that the two PESs are harmonic (i.e. parabolic along all normal coordinates), the relationship between the two PES being defined by the Duschinsky transformation [33],

[TABLE]

where $\vec{q}$ and $\vec{q^{\prime}}$ are the vibrational normal coordinates for the initial (e.g. ground) and final (e.g. excited) PES, respectively, the Duschinsky matrix S is unitary, and $\vec{d}$ is a displacement vector.

A vibronic spectrum calculation consists of determining the Franck-Condon profile (FCP), defined as

[TABLE]

where $\omega$ is the transition energy, $|0\rangle$ is the vibrational vacuum state of the initial PES (i.e. of the initial electronic state), and $|f_{i}\rangle$ is the $i$th eigenstate of the final PES with energy $\omega_{i}$. In practice, the function is desired to some precision $\Delta\omega$. Fig. 1 gives a schematic of the vibronic problem for the (a) one-dimensional and (b) multidimensional case, where each parabola or hypersurface represent an electronic PES.

In the photonics-based vibronic boson sampling (VBS) algorithm [20], a change of basis known as the Doktorov transformation [34] is used to transform between the two PESs, as this harmonic transformation is directly implementable in photonic circuit elements.

Instead of this direct basis change approach, our work is based on constructing a Hamiltonian that encodes the relationship between the two PESs. This provides multiple advantages, outlined below.

We denote dimensionless position and momentum operators as $\hat{q}_{sk}$ and $\hat{p}_{sk}$ respectively, where $s$ labels the potential energy surface ($s\in\{A,B\}$ in this work) and $k$ labels the vibrational mode. These follow standard definitions $\hat{q}_{sk}=(\hat{a}_{sk}+\hat{a}_{sk}^{\dagger})/\sqrt{2}$ and $\hat{p}_{sk}=(\hat{a}_{sk}-\hat{a}_{sk}^{\dagger})/i\sqrt{2}$, where $a_{sk}^{\dagger}$ and $a_{sk}$ are vibrational creation and annihilations operators. The notation $\vec{\cdot}$ denotes standard vectors as well as vectors of operators, such that e.g. $\vec{q}_{A}=[\hat{q}_{A0},...,\hat{q}_{AM}]^{T}$.

The purpose of our classical pre-processing procedure is to express the vibrational Hamiltonian for PES $B$ in terms of $\{q_{Ak},p_{Ak}\}$, by making the following transformations: $\{\vec{q}_{A},\vec{p}_{A}\}$ $\rightarrow\{\vec{Q}_{A},\vec{P}_{A}\}$ $\rightarrow\{\vec{Q}_{B},\vec{P}_{B}\}$ $\rightarrow\{\vec{q}_{B},\vec{p}_{B}\}$, where $\vec{Q}_{s}$ and $\vec{P}_{s}$ are respectively the mass-weighted position and momentum operators of PES $s$ [19]. The full transformations are

[TABLE]

[TABLE]

where

[TABLE]

and $\{\omega_{sk}\}$ are the quantum harmonic oscillator (QHO) frequencies of PES $s$. A more pedagogical explanation as well as an alternate formulation are given in the Supplemental Information. Parameter $\delta$ is often used[35, 19, 20], defined $\delta_{sk}=d_{sk}\sqrt{\omega_{sk}/\hbar}$.

Finally, the vibrational Hamiltonian of PES $B$ is expressed in a standard form as

[TABLE]

after each $q_{Bk}$ and $p_{Bk}$ has been constructed as a function of the ladder operators of PES $A$. Hence the low-level building block of our algorithm is a truncated creation operator,

[TABLE]

where $l$ denotes a vibrational energy level and the imposed cutoff $L_{max}$ denotes the maximum level. Mappings to qubits (i.e. integer-to-bit encodings) are discussed in the SI and errors are analyzed below.

Though this work primarily considers the harmonic case in order to introduce our methodology, the largest quantum advantage will arise from modeling anharmonic effects. In fact, because the harmonic approximation is amenable to clever classical techniques that cannot be applied to anharmonic PESs [19], it is expected that quantum advantage would be more easily demonstrated for the anharmonic problems than in the harmonic ones.

FCPs from anharmonicity are vastly more costly to approximate than the harmonic case using classical algorithms—for calculations that include Duschinsky and anharmonic effects, we are not aware of molecules larger than six atoms that have been accurately simulated [29, 17, 30, 31, 32]. Arbitrary anharmonicity can be straightforwardly included in our quantum algorithm by adding higher-order potential energy terms to the unperturbed (e.g. Eq. 6) vibrational Hamiltonian $H_{0}$:

[TABLE]

The ease with which one includes anharmonic effects is an advantage over the VBS algorithm [20, 21].

Now that we have outlined the required classical steps, we describe our quantum algorithm for determining the Franck-Condon profile. Unlike most quantum computational approaches to Hamiltonian simulation [36, 37, 38, 12], which aim to find the energy of a particular quantum state, the purpose of our algorithm is to construct a full spectrum from many measurements.

As the procedure makes use of the quantum phase estimation (QPE) algorithm [39, 40, 41], we use two quantum registers. QPE is a quantum algorithm that calculates the eigenvalues of a superposition of states, acting on a quantum state as $\sum c_{i}|\psi_{i}\rangle\otimes|0\rangle$ $\rightarrow$ $\sum c_{i}|\psi_{i}\rangle\otimes|\tilde{\phi}_{i}\rangle$, where $\tilde{\phi}_{i}$ and $|\psi_{i}\rangle$ are eigenpairs with respect to an implemented operator. The first register $S$ stores a representation of the vibrational state, and the second register $E$ is used to read out the energy (strictly speaking, it outputs the phase, from which the energy is trivially obtained).

$S$ is initialized to $|0\rangle$, the ground state of $H_{A}$. A simple but key observation is that $|0\rangle$ can be written in the eigenbasis of $H_{B}$, such that

[TABLE]

where $\{|\psi_{i}\rangle\}$ are eigenstates of $H_{B}$ and coefficients $c_{i}$ are not a priori known.

One then runs QPE using the Hamiltonian $H_{B}$ (i.e. implementing $U=e^{-i\tau H_{B}}$ for some arbitrary value $\tau$), with register $E$ storing the eigenvalues. Many quantum algorithms have been developed for Hamiltonian simulation [36, 42, 43, 44, 45, 46, 47, 48, 49], any of which can be used in conjunction with the algorithm’s QPE step. Convincing numerical evidence suggests that Trotterization [36, 42] is likely to be the most viable option for early quantum devices [48]. Computational scaling is briefly discussed in the Supplemental Information.

We define $\varepsilon_{i}$ as the eigenenergy of $|\psi_{i}\rangle$, and $\tilde{\varepsilon}_{i}$ as its approximation, where an arbitrarily high precision can be achieved by increasing the number of qubits in register $E$. Degeneracies in $\tilde{\varepsilon}$ will be ubiquitous, and we define the subspace of states with approximate energy $\tilde{\varepsilon}_{j}$ as $\mathcal{D}_{j}=\{|\psi_{j1}\rangle,...,|\psi_{jK_{j}}\rangle\}$, where $K_{j}$ is the degeneracy in $\tilde{\varepsilon}_{j}$. Measuring register $E$ yields $\tilde{\varepsilon}_{j}$ with probability $\sum_{k\in\mathcal{D}_{j}}|c_{k}|^{2}$. Hence—and this is the key insight—values $\tilde{\varepsilon}_{j}$ are outputted with a probability exactly in proportion to the Franck-Condon factors of Eq. 2. The measurements then produce a histogram that yields the vibronic spectrum. The procedure is depicted in Fig. 2, where for the zero-temperature case one may disregard register $I$ and gate $\hat{V}(\beta)$. See the Supplemental Information for a step-by-step outline of the algorithm.

Note that this is a different approach from how QPE is usually used. Normally one attempts to prepare a state that is as close as possible to a desired eigenstate, whereas here we deliberately begin with a broad mix of eigenstates that corresponds to the particular spectrum we wish to calculate.

We highlight four potential benefits of this algorithm over the VBS algorithm [20, 21]. First, the quantum state in register $S$ is preserved for further analysis, while in VBS the final state is destroyed. After measurement, the state stored in $S$ is a superposition of states with energy $\tilde{\varepsilon}_{j}$. From several runs of the circuit, one may estimate this stored state’s overlap with another quantum state [50], estimate its expectation value with respect to an arbitrary operator, or calculate the transition energy to another PES (i.e. simulate excited-state absorption), though methods for analyzing this preserved information are beyond the scope of this work. Our QPE-based approach has a similar benefit over the canonical quantum circuit method for calculating correlation functions [37], which does not provide this kind of interpretable post-measurement state.

The second potential benefit is that, as stated above, anharmonic effects are easily included in our framework. Third, accurate photon number detection for higher photon counts is a major difficulty in experimental quantum optics [22, 23]; it may be that a scaled-up universal quantum computer is built before quantum optical detectors improve satisfactorally, though this is difficult to predict. Fourth, while there are error correction methods for universal quantum computers, we do not know of such methods for boson sampling devices.

Even at room temperature, the optical spectrum of a molecule can be substantially different from its zero temperature spectrum [51], necessitating methods for including finite temperature effects. These effects can be elegantly included by appending additional steps before and after the zero temperature algorithm, following previous work [52, 21]. Briefly, one appends an additional quantum register, labeled $I$, with the same size as register $S$. An operator called $\hat{V}(\beta)$ is applied to state $|0\rangle_{S}|0\rangle_{I}$, which entangles the $E$ and $I$ registers to produce a thermofield double state. The remainder of the algorithm proceeds as before, except that both $E$ and $I$ are measured and the contribution to the histogram is modified. We elaborate on this procedure in the Supplemental Information.

The least-studied source of error in our algorithm is due to an insufficiently large QHO cutoff $L_{max}$. It is especially important to study this source of error, both qualitatively and quantitatively, because the standard classical algorithms for calculating FCFs [15, 17, 3, 29, 30, 31, 32] do not directly simulate the vibrational Hamiltonian in the Fock basis of PES $A$, and hence do not suffer from this type of truncation error. An analysis of Suzuki-Trotter errors will be dependent on the QHO mapping chosen and is left to future work.

We chose transitions in four triatomic molecules—sulfur dioxide (SO2${}^{-}\rightarrow$ SO2) [53], water (H2O$\rightarrow$H2O++e-)[54], deuterated water (D2O; D $\equiv$ 2H)[54], and nitrogen dioxide (NO2[${}^{2}A_{1}$] $\rightarrow$ NO2[${}^{2}B_{2}$])[55]—and simulated their vibronic spectra using one electronic transition from each (the SI includes additional information, including all physical parameters used). We must note that there is a well-studied conical intersection near the bottom of PES ${}^{2}B_{2}$ in NO2[57, 58, 59, 60] which greatly alters the vibronic properties—in the approximation used here, this feature is ignored. The latter three molecules were chosen explicitly because they have unusally high phonon occupation numbers for a vibronic transition, making them good candidates for a study on $L_{max}$ requirements.

Fig. 3 shows both the theoretically exact vibronic spectra (solid line) and an arbitrarily chosen approximate spectrum (dotted line) for each molecule. These plots show the qualitative error behavior: higher-energy peaks are blue-shifted while low-energy peaks converge rapidly. It is useful to consider this error trend when implementing the algorithm on a future quantum computer. Simulation details and additional analysis are given in the SI.

Finally, in order to highlight the simulation of anharmonicity as our quantum algorithm’s primary application, we consider the zero temperature anharmonic spectrum of SO2${}^{-}\rightarrow$ SO2, using the same Duschinsky matrix as before [53] but an anharmonic PES for SO2  taken from Smith et al. [56]. Fig. 4 shows the stark difference between the full anharmonic simulation and the harmonic approximation. We included both third-order and fourth-order Taylor series terms (Eq. 8), which are easily mapped to the quantum computer using the same procedure as before (see SI for additional details). This relative failure of the harmonic approximation is not uncommon [29, 17, 30, 31, 32], indicating a well-defined set of molecules (those with substantially anharmonic PESs) for which our quantum algorithm would outperform classical computers.

We introduced a quantum algorithm for calculating the vibronic spectrum of a molecule to arbitrary precision. We noted several advantages over the previously proposed vibronic boson sampling (VBS) algorithm. First and perhaps most importantly, anharmonic effects (whose inclusion is very costly classically but often chemically relevant) can be easily included in our approach. Second, measuring the eigenenergy in our algorithm leaves the quantum state preserved, allowing for further analysis that would not be possible in VBS. Third, there are error-related advantages to our approach. Aspects of our algorithm may be extended to other chemical processes for which nuclear degrees of freedom are difficult to simulate on a classical computer. This work’s general strategy, of calculating the energy distribution outputted from quantum phase estimation to arbitrary precision, may be applied to other spectral problems in chemistry and condensed matter physics.

Supporting Information.

Elaboration on vibronic Hamiltonian construction, quantum harmonic oscillator to qubit mappings, the finite temperature algorithm, computational scaling, molecular data and parameters used, and additional error analysis.

Acknowledgements

J.H. acknowledges support by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2015R1A6A3A04059773). The authors thank Gian Giacomo Guerreschi and Daniel Tabor for helpful suggestions on the manuscript.