Efficient multiphoton sampling of molecular vibronic spectra on a superconducting bosonic processor

TL;DR

This paper demonstrates a superconducting bosonic processor capable of efficiently sampling molecular vibronic spectra using multiphoton states, high-fidelity photon detection, and reprogrammability for various molecules, advancing quantum simulation capabilities.

Contribution

It introduces the first high-fidelity single-shot photon number resolving detection in a superconducting platform and showcases reprogrammable bosonic simulations of molecular spectra.

Findings

Successful implementation of boson sampling for molecular spectra

High-fidelity detection of up to 15 photons per mode

Reprogrammable simulation of multiple molecular vibronic spectra

Abstract

The efficient simulation of quantum systems is a primary motivating factor for developing controllable quantum machines. For addressing systems with underlying bosonic structure, it is advantageous to utilize a naturally bosonic platform. Optical photons passing through linear networks may be configured to perform quantum simulation tasks, but the efficient preparation and detection of multiphoton quantum states of light in linear optical systems are challenging. Here, we experimentally implement a boson sampling protocol for simulating molecular vibronic spectra [Nature Photonics $\textbf{9}$, 615 (2015)] in a two-mode superconducting device. In addition to enacting the requisite set of Gaussian operations across both modes, we fulfill the scalability requirement by demonstrating, for the first time in any platform, a high-fidelity single-shot photon number resolving detection scheme capable of resolving up to 15 photons per mode. Furthermore, we exercise the capability of synthesizing non-Gaussian input states to simulate spectra of molecular ensembles in vibrational excited states. We show the re-programmability of our implementation by extracting the spectra of photoelectron processes in H$_2$O, O$_3$, NO$_2$, and SO$_2$. The capabilities highlighted in this work establish the superconducting architecture as a promising platform for bosonic simulations, and by combining them with tools such as Kerr interactions and engineered dissipation, enable the simulation of a wider class of bosonic systems.

Tables

Table 1. Table 1: Full set of controls of the bosonic processor. The subscripts i ∈ { A , B } 𝑖 𝐴 𝐵 i\in\{A,B\} correspond to operations on each cavity and their respective ancillary modules. κ A subscript 𝜅 𝐴 \kappa_{A} and κ B subscript 𝜅 𝐵 \kappa_{B} are the intrinsic linewidths of the cavity modes. Optimal control pulses that utilize both resonant cavity and ancilla drives with complex envelopes, ε i ​ ( t ) subscript 𝜀 𝑖 𝑡 \varepsilon_{i}(t) and ϵ i ​ ( t ) subscript italic-ϵ 𝑖 𝑡 \epsilon_{i}(t) , respectively, generate a desired initial Fock state (see supplementary text IV). | n ′ , m ′ ⟩ ket superscript 𝑛 ′ superscript 𝑚 ′ \ket{n^{\prime},m^{\prime}} denotes a particular joint photon number to be measured in the cavities, and { b i } subscript 𝑏 𝑖 \{b_{i}\} represent the bits associated with the binary decomposition of the photon number. For example, | 5 ⟩ = | 0101 ⟩ ket 5 ket 0101 \ket{5}=\ket{0101} . The third column presents an estimate of the expected limits on fidelity due to photon loss for the hardware in this current implementation of the simulator.

	Interaction	Error Rate
State preparation
Optimal Control	${\hat{H}}_{\text{drive}}\left(t\right)={\epsilon }_i^{\ast }\left(t\right){\hat{c}}_i+{\epsilon }_i\left(t\right){\hat{c}}_i^{†}+{ϵ}_i^{\ast }\left(t\right){\hat{t}}_i+{ϵ}_i\left(t\right){\hat{t}}_i^{†}$	${\kappa }_i{\tau }_{\text{prep},i}\sim {10}^{-3}-{10}^{-2}$
Operations
Displacement	${\tilde{\epsilon }}_i^{\ast }\left(t\right){\hat{c}}_i+{\tilde{\epsilon }}_i\left(t\right){\hat{c}}_i^{†}$	${\kappa }_i{\tau }_{\text{disp},i}\sim {10}^{-4}$
Squeezing	$g_{\text{sq},i}\left(t\right)\left(e^{-i{\varphi }_i}{\hat{c}}_i^2+e^{i{\varphi }_i}{\hat{c}}_i^{†2}\right)$	${\kappa }_i/g_{\text{sq},i}\sim 5\times {10}^{-2}$
Beamsplitter	$g_{\text{BS}}\left(t\right)\left(e^{i\phi }{\hat{c}}_A{\hat{c}}_B^{†}+e^{-i\phi }{\hat{c}}_A^{†}{\hat{c}}_B\right)$	${\kappa }_i{\tau }_{\text{BS}}\sim {10}^{-2}$
Measurement
Single-bit extraction	${\hat{M}}_0=|n^{\prime },m^{\prime }⟩⟨n^{\prime },m^{\prime }|{\hat{M}}_1=\widehat{\mathrm{𝟙}}-{\hat{M}}_0$	${\kappa }_i{\tau }_{\text{meas}}\sim {10}^{-3}-{10}^{-2}$
Sampling	$|n^{\prime }⟩=|b_3,b_2,b_1,b_0⟩$	${\kappa }_i{\tau }_{\text{meas}}\sim {10}^{-2}-{10}^{-1}$

Equations

$$\hat{\text{a}}\to {\hat{U}}_{\text{Dok}}\hat{\text{a}}{\hat{U}}_{\text{Dok}}^{†}$$

$${\hat{U}}_{\text{Dok}}=\hat{\mathbf{D}}(\mathbf{\alpha }){\hat{\mathbf{S}}}^{†}({\mathbf{\zeta }}^{\prime })\hat{\mathbf{R}}(U)\hat{\mathbf{S}}(\mathbf{\zeta })$$

$$\hat{\mathbf{D}}(\mathbf{\alpha })$$

$${\hat{\mathbf{S}}}^{(†)}({\mathbf{\zeta }}^{{(}^{\prime })})$$

$${\text{FCF}}_{\vec{n},{\vec{n}}^{\prime }}=|⟨{\vec{n}}^{\prime }|{\hat{U}}_{\text{Dok}}|\vec{n}⟩{|}^2$$

$${\hat{H}}_{\text{BS}}/\hslash$$

$${\hat{H}}_{\text{sq},i}/\hslash$$

$$\sum _{i,j}{c}_{ij}|i,j⟩\otimes |g,g⟩\to \sum _{i\ne n^{\prime },j\ne m^{\prime }}{c}_{ij}|i,j⟩\otimes |g,g⟩$$

$$+$$

$$+c_{n^{\prime }{m}^{\prime }}|n^{\prime },m^{\prime }⟩\otimes |e,e⟩$$

$$\hat{H}=-\chi {\hat{c}}^{†}\hat{c}{\hat{t}}^{†}\hat{t}+{ϵ}^{\ast }(t)\hat{t}+ϵ(t){\hat{t}}^{†}$$

$$(\hat{P}_{k})_{ij}=\left\{\begin{array}[]{ll}0&\textrm{if}\quad i\neq j\\ 1-2\bigg{(}\lfloor\frac{i}{2^{k}}\rfloor\quad(\textrm{mod}\,2)\bigg{)}&\textrm{if}\quad i=j\\ \end{array}\right.$$

Full text

I Introduction

Simulation of quantum systems with quantum hardware offers a promising route towards understanding the complex properties of those systems that lie beyond the computational power of classical computers Feynman1982 . A particularly efficient approach is one that utilizes the natural properties of the quantum hardware to simulate physical systems that share those properties. This approach serves as the foundation for bosonic quantum simulation, where the natural statistics and interference between bosonic excitations are directly exploited. A prominent example of this is the manipulation of bosonic atoms in an optical lattice to explore various types of many-body physics Bakr2009 ; Tai2017 ; Bernien2017 .

Boson sampling is another example of a computationally challenging task that can be performed by manipulating bosonic excitations Aaronson2013 ; Hamilton2017 . Conventionally described in the context of linear optics, boson sampling, in its many forms, involves the single photon detection of nonclassical states of light passing through a linear interferometric network. Current technologies for optical waveguides allow the creation of complex interferometers across many modes. An outstanding challenge in the optical domain, however, is generating and detecting nonclassical states of light with high efficiencies.

In the microwave domain, nonlinearities provided by Josephson junctions enable the powerful preparation and flexible quantum nondemolition (QND) measurement of quantum states of light in superconducting circuits Schuster2007 ; Johnson2010 . The ability to performing QND photon number measurements not only enables the high fidelity measurement of bosonic qubits via repeated detections Hann2018 ; Elder2019 , but also enables the construction of more complex measurement operators. The circuit QED platform has also demonstrated universal control over individual bosonic modes as well as robust beamsplitter operations between separate modes Heeres2017 ; Gao2018 . These capabilities motivate an alternative approach to performing both boson sampling Peropadre2016 and bosonic quantum simulation protocols using superconducting circuits.

An instance of bosonic quantum simulation that maps onto a generalized boson sampling problem is obtaining molecular vibronic spectra associated with electronic transitions Huh2015 . The algorithm requires, in addition to beamsplitters, single-mode displacement and squeezing operations. Furthermore, the output will generally be a multimode multiphoton state, thus requiring a series of photon number resolving detectors at the output. Experimental imperfections in the controls, photon loss, and detector inefficiency has made a linear optical implementation of this algorithm challenging Clements2018 . Recent work leveraging the bosonic nature of two phonon modes of a single trapped ion has been successful in generating more accurate spectra Shen2018 ; However, an efficient detection scheme capable of directly sampling from the exponentially growing bosonic Hilbert space remains elusive.

Here, we experimentally implement a two-mode superconducting bosonic processor with a full set of controls that enables the scalable simulation of molecular vibronic spectra. The processor combines arbitrary (Gaussian and non-Gaussian) state preparation and a universal set of Gaussian operations enabled by four-wave mixing of a Josephson potential. Most importantly, we implement a single-shot QND photon number resolving detection scheme capable of resolving up to $n_{\textrm{max}}$ = 15 photons per mode. This detection scheme, when operated without errors on a multiphoton distribution bounded within $n_{\textrm{max}}$ per mode, extracts the maximum possible amount of information from the underlying spectra per run of the experiment.

II Bosonic algorithm for Franck Condon factors

The mapping of molecular vibronic spectra onto a bosonic simulation framework can be understood by considering the nature of vibrational dynamics that accompanies an electronic transition. In keeping with the tenets of the adiabatic Born-Oppenheimer approximation, the presumed separability between electronic and nuclear degrees of freedom results in distinct electronic states, each of which forms a potential-energy surface (PES) that supports a distinct manifold of vibrational eigenstates. The normal modes of vibrations for a given electronic state are obtained by expanding the PES in powers of displacement coordinates referenced to the minimum-energy (equilibrium) configuration and retaining only up to quadratic terms. Under this harmonic approximation, the corresponding transformation of the set of creation and annihilation operators $\hat{\textbf{a}}=(\hat{a}_{1},...,\hat{a}_{N})$ for $N$ vibrational modes, may be expressed using the Doktorov transformation Doktorov1977 :

[TABLE]

where

[TABLE]

correspond to $N$-dimensional vectors of single-mode displacement and squeezing operations, respectively. $\hat{\boldsymbol{R}}(U)$ is an $N$-mode rotation operator corresponding to a $N\times N$ rotation matrix $U$ which can be decomposed into a product of two-mode beamsplitter operations (see supplementary text I) Reck1994 . For $N$ = 2, $U$ is a two-dimensional rotation matrix parameterized by a single angle $\theta$. The set of dimensionless Doktorov parameters $\boldsymbol{\alpha}=(\alpha_{1},...,\alpha_{N}),\boldsymbol{\zeta}=(\zeta_{1},...,\zeta_{N}),\boldsymbol{\zeta^{\prime}}=(\zeta^{\prime}_{1},...,\zeta^{\prime}_{N})$ and $U$ originate from molecular structural information in the different electronic configurations (see supplementary text I). Applying $\hat{U}_{\textrm{Dok}}$ to an initial state $\ket{\psi_{0}}$ of the bosonic processor directly emulates the physical process of a pre-transition molecular vibrational state experiencing a sudden change in the vibrational PES and being expressed in the post-transition vibrational basis. This projection gives rise to the Franck-Condon factors (FCFs) defined as the vibrational overlap integrals of an initial pre-transition vibrational eigenstate, $\ket{\vec{n}}=\ket{n,m,...}$, with a final post-transition vibrational eigenstate, $\ket{\vec{n}^{\prime}}=\ket{n^{\prime},m^{\prime},...}$:

[TABLE]

In spectroscopic experiments, the validity of this “sudden approximation” depends on the vastly different energy and time scales that typically characterize electronic and nuclear motions (e.g., $\gg$10000 cm*-1* vs. $\sim$1000 cm*-1* and $\ll 10^{-18}$ s vs. $\geq 10^{-15}$ s, respectively). By further assuming that the transition electric dipole moment does not depend on nuclear coordinates (i.e., the Condon approximation), the relative intensities of features appearing in vibrationally resolved absorption and emission spectra will be directly proportional to the corresponding FCFs. Indeed, the practical importance of these quantities often stems from the structural and dynamical insights that they can provide about excited electronic states, information that can be challenging to obtain through other conventional spectroscopic means.

III Experimental implementation

Our superconducting processor is designed to manipulate the bosonic modes of two microwave cavities, $\hat{c}_{A}$ and $\hat{c}_{B}$ (Fig. 1). In our design, a coupler transmon $\hat{t}_{C}$ dispersively couples to both cavities, enabling beamsplitter Gao2018 and squeezing operations through driven four-wave mixing processes. The coupler transmon is also dispersively coupled to a readout resonator $\hat{r}_{C}$. Displacement operations on the cavity modes are performed via resonant drives through local coupling ports. Additional ancillary transmon-readout systems $\{\hat{t}_{A},\hat{r}_{A}\}$ and $\{\hat{t}_{B},\hat{r}_{B}\}$ are inserted into the device and couple to each cavity for the purposes of state preparation and tomography.

The goal of the present simulation is to emulate the transformation of a molecular vibrational state due to an electronic transition using the photonic state of the quantum processor. The processor is first initialized in a state $\ket{\psi_{0}}$ corresponding to the pre-transition molecular vibrational state of interest. A vacuum state of both bosonic modes is prepared through feedback cooling protocols, and Fock states $\ket{n,m}$ are initialized using optimal control techniques Heeres2017 (see supplementary text IV). The ability to reliably synthesize arbitrary Fock states translates to the powerful capability of simulating FCFs starting from vibrationally excited states, a task which is challenging in most other bosonic simulators. The Doktorov transformation is then applied, producing a basis change to that of the post-transition vibrational Hamiltonian. In our case where $N$ = 2, the rotation operator corresponds directly to enacting a single beamsplitter. Both the single-mode squeezing and beamsplitter operations utilize the four-wave mixing capabilities of the coupler transmon. Two pump tones that fulfill the appropriate frequency matching condition ($\omega_{1}+\omega_{2}=2\omega_{A/B}$ for the squeezing operation, and $\omega_{2}-\omega_{1}=\omega_{B}-\omega_{A}$ for the beamsplitter operation, where $\omega_{1/2}$ are the two pump frequencies and $\omega_{A/B}$ are the cavity frequencies) are sent through a port that primarily couples to the coupler transmon, which enacts the desired Hamiltonians Gao2018 :

[TABLE]

where $i\in\{A,B\}$. Importantly, the phase of the operations $\{\varphi,\phi_{A},\phi_{B}\}$ implemented by these Hamiltonians is controlled by the phase of the pump tones, allowing the microwave control system to generate the correct family of beamsplitter and squeezing operations for performing $\hat{U}_{\textrm{Dok}}$. A set of transmon measurements is then carried out for the purpose of post-selecting the final data on measuring all transmons in their ground states. This verification step primarily aims to reject heating events of the transmons out of their ground state; the heating of these ancillae otherwise dephases the cavities while coupler heating effectively halts the pumped operations by shifting the requisite frequency matching conditions. This post-selection also serves to ensure that the ancillae begin in their ground states for the subsequent measurement of the cavities. In our experiment, we reject $5-10\%$ of the data depending on which transformation is simulated. Finally, averaging many measurements of the cavities in a photon number basis, $\{n^{\prime},m^{\prime}\}$, give the desired FCFs. The full set of controls is detailed in Table 1; the relatively small error rates due to photon loss provide a sense of scale for attainable circuit depths while maintaining high fidelity performance.

IV Measurement Protocols

Two complementary measurement schemes are used to extract FCFs from the final state of the processor, both of which utilize the dispersive coupling of each microwave cavity to its ancillary transmon-readout system (Fig. 2a): $\hat{H}_{\textrm{int}}/\hbar=-\sum_{i\in\{A,B\}}\chi_{i}\hat{c}^{\dagger}_{i}\hat{c}_{i}\hat{t}^{\dagger}_{i}\hat{t}_{i}$, where the dispersive interaction strengths in our experiment are $\chi_{A}=2\pi\times 748$ kHz and $\chi_{B}=2\pi\times 1240$ kHz. Fundamentally, the difference between these two measurement schemes arises from the ability of the latter to extract more than one bit of information on a given single shot of the experiment.

IV.1 Single-bit extraction

The first scheme (Fig. 2b) maps a given joint cavity photon number population $\{n^{\prime},m^{\prime}\}$ onto the joint state of the two transmons via state-selective $\pi$ pulses. These pulses have frequencies $\omega_{t_{A}}=\omega^{0}_{t_{A}}-n^{\prime}\chi_{A}+(n^{\prime 2}-n^{\prime})\frac{\chi^{\prime}_{A}}{2}$ and $\omega_{t_{B}}=\omega^{0}_{t_{B}}-m^{\prime}\chi_{B}+(m^{\prime 2}-m^{\prime})\frac{\chi^{\prime}_{B}}{2}$ where the small second-order dispersive shift is also taken into account: $\hat{H}^{\prime}_{\textrm{int}}/\hbar=\sum_{i\in\{A,B\}}\frac{\chi^{\prime}_{i}}{2}\hat{c}^{\dagger}_{i}\hat{c}^{\dagger}_{i}\hat{c}_{i}\hat{c}_{i}\hat{t}^{\dagger}_{i}\hat{t}_{i}$. In our experiment, $\chi^{\prime}_{A}=2\pi\times 1.31$ kHz and $\chi^{\prime}_{B}=2\pi\times 1.35$ kHz. The pulses are applied with a Gaussian envelope truncated at $\pm 2\sigma_{t}$ such that the bandwidth is approximately $\sigma_{f}=1/(2\pi\sigma_{t})$, where $\sigma_{t}$ is the standard deviation of the pulse in time. The selectivity is defined as the probability of exciting the ancilla given occupation in an adjacent photon number state of interest. Pulses with $\sigma_{t}$ = 1 $\mu$s for both ancillas are used in the experiment, which give selectivities above 99.9$\%$ and implement the following mapping for a general state of the two cavities $\ket{\psi}=\sum_{i,j}c_{ij}\ket{i,j}$ for a chosen set of photon numbers $\{n^{\prime},m^{\prime}\}$ to probe:

[TABLE]

where $\ket{g}$ and $\ket{e}$ are the ground and first excited state of each ancilla transmon. The transmons are then individually read out using standard dispersive techniques, and the results are correlated on a shot-by-shot basis to extract a single bit of information for each joint photon number state probed. We thus call this measurement scheme “single-bit extraction.”

IV.2 Photon number resolved sampling

Extracting FCFs using the “single-bit extraction” scheme, however, is not scalable. The bosonic Hilbert space grows exponentially with the number of modes $N$ as $n^{N}_{\textrm{max}}$, where $n_{\textrm{max}}$ is the maximum number of excitations considered for each mode. Thus, any measurement protocol that only extracts a single bit of information from the underlying distribution at a time must necessarily query the exponentially growing number of final states. This is the case for a recent implementation of this simulation using two phonon modes of a single trapped ion Shen2018 .

The second measurement scheme circumvents this problem and implements single-shot photon number resolving detection (Fig. 2c). At the heart of this technique is the fact that the following Hamiltonian:

[TABLE]

enables the QND mapping of any binary valued operator of the cavity Hilbert space onto the state of the transmon by numerically optimizing an appropriate waveform $\epsilon(t)$ that is applied to only the transmon. Through this capability, we synthesize pulses $\epsilon_{k}(t)$, $k\in\{0,1,2,3\}$ that excite the transmon from its ground state $\ket{g}$ to the excited state $\ket{e}$ conditioned on the cavity state’s projection onto a series of parity subspaces that represent a binary decomposition $b_{k}$ of the photon number $\ket{n}=\ket{b_{3}b_{2}b_{1}b_{0}}$:

[TABLE]

with eigenvalues $\pm$1. In our experiment, we synthesize pulses each with a duration between $800-1200$ ns that explicitly designate the above parity operators over the Hilbert space of each cavity up to $n_{max}=15$.

The QND nature of each of these mapping pulses on the cavity state ensures that the pulses can be applied sequentially (with transmon measurements following each pulse) to project an initial cavity state $\ket{\psi}=\sum_{n=0}^{15}c_{n}\ket{n}$ with bounded support within $n_{max}=15$ into a definite Fock state $\ket{n}$ with a probability $|c_{n}|^{2}$. In order to minimize errors due to decoherence when the transmon is excited, after the $k^{\textrm{th}}$ bit is mapped onto the transmon and measured, the transmon is reset to its ground state using real-time feed-forward control to prepare for the measurement of the $(k+1)^{\textrm{th}}$ bit of the same cavity state. This protocol is performed simultaneously on each cavity, which returns a sample from the underlying joint photon number distribution on a shot-by-shot basis. This measurement scheme thus implements bosonic sampling of the final state distribution, which we will simply call “sampling.” These binary detectors optimally resolve the photon number in the cavities in $N\textrm{log}_{2}(n_{\textrm{max}})$ measurements, but are prone to more errors due to transmon decoherence that are, in principle, correlated bit-to-bit. We leave the task of characterizing these errors in full detail, optimizing the technique, and potentially applying de-convolution methods Nachman2019 to improve the accuracy of the sampled distribution as the subject of future work.

V Simulated photoelectron spectra

The full set of FCFs for a given electronic transition can provide, under the Condon approximation, the relative intensities of vibronic progressions appearing in corresponding photoelectron spectra. We simulate four different photoelectron processes starting in various vibrational initial states $\ket{\psi_{0}}=\ket{n,m}$: H2O $\xrightarrow{h\nu}$ H2O${}^{+}(\tilde{B}^{2}$B${}_{2})$ + e-, O${}^{-}_{3}\xrightarrow{h\nu}$ O3 + e-, NO${}^{-}_{2}\xrightarrow{h\nu}$ NO2 + e-, and SO${}_{2}\xrightarrow{h\nu}$ SO${}^{+}_{2}$ + e-. Experimental results for the first two processes starting in $\ket{\psi_{0}}=\ket{0,0}$ and $\ket{\psi_{0}}=\ket{1,2}$, respectively, are presented in Fig. 3 (see supplementary text VI for additional data). The particular electronic state of the water cation considered here $(\tilde{B}^{2}\textrm{B}_{2})$ is the second excited state of doublet spin multiplicity, with the attendant electronic wave function having B2 symmetry. The other four processes consider transitions beginning and ending in electronic ground states. All of the triatomic molecules targeted here retain $C_{2v}$ point-group symmetry for their equilibrium configurations in both the pre-transition and post-transition electronic states, thereby ensuring that analyses can be restricted to the two-dimensional subspace of the symmetric-stretching and bending modes.

A figure of merit for quantifying the quality of the quantum simulation is the distance $D=\frac{1}{2}\sum_{i=0}^{n_{\textrm{max}}}\sum_{j=0}^{n_{\textrm{max}}}|p^{\textrm{meas}}_{ij}-p^{\textrm{ideal}}_{ij}|$ between the measured probabilities $\{p^{\textrm{meas}}_{ij}\}$ and the ideal distribution $\{p^{\textrm{ideal}}_{ij}\}$. The distances for the two simulated processes in Fig. 3 are D = 0.049 (H2O) and 0.105 (O3) for the “single-bit extraction” scheme. This distance metric is accompanied with a success probability due to post-selection of transmon heating events, which are 95$\%$ and 93$\%$ for the aforementioned simulations. The heating events are dominated by the coupler transmon; the dynamics of a driven Josephson element for engineered bilinear operations presents a multi-dimensional optimization problem that seeks to maximize the desired interaction rates while minimizing induced decoherence and dissipation rates Zhang2019 . Self-Kerr Hamiltonian terms of the form $\hat{H}_{\textrm{Kerr}}/\hbar=-\sum_{i\in\{A,B\}}\frac{K_{i}}{2}\hat{c}^{\dagger}_{i}\hat{c}^{\dagger}_{i}\hat{c}_{i}\hat{c}_{i}$, photon loss, and imperfect state preparation account for errors that are undetected by post-selection. The first two effects are captured through full time-domain master equation simulations (see supplementary text VI). The magnitudes of these two errors depends on the molecular process; each corresponding Doktorov transformation will have different squeezing and rotation parameters thus leading to varying lengths of the pumped operations. Simulations of shorter length circuits will therefore have lower error rates. Additionally, errors due to self-Kerr interactions of the cavities are larger for higher photon number states.

VI Resource requirements and Scalability

The central advantages of simulating the transformation of a bosonic Hamiltonian using a bosonic system lie in both the native encoding and the efficient decomposition of the Doktorov transformation into Gaussian operations. A recent proposal for obtaining Franck-Condon factors on a conventionally envisioned spin-$\frac{1}{2}$ quantum computer Sawaya2019 needs to map the problem onto qubits and a univeral gate set. This first requires encoding the Hilbert space of size $n^{N}_{\textrm{max}}$ onto $n_{q}=N\textrm{log}_{2}(n_{\textrm{max}})$ qubits. The choice of $n_{\textrm{max}}$ is dependent on the initial state as well as the magnitude of the displacement and squeezing; both operations can produce states with large photon numbers. Using quantum signal processing Low2017 , the approximate number of gates $n_{g}$ then needed to implement $\hat{U}_{\textrm{Dok}}$ using a universal qubit gate set to within an error $\varepsilon$ is $n_{g}=O(N^{2}n^{2}_{\textrm{max}}\textrm{log}^{3}(1/\varepsilon))$. For our experiment with $N$ = 2 modes, taking $n_{\textrm{max}}$ = 16 and desiring an error $\varepsilon=5\times 10^{-2}$, this translates to $n_{q}$ = 8 qubits and $n_{g}=O(10^{3})$ gates. The coherence requirements for performing such a computation this way is thus relatively demanding and exceeds the capabilities of current technologies. The advantage of the qubit-based algorithm, however, is the ability to systematically incorporate anharmonicities in the PES, a task which still needs to be theoretically investigated for the bosonic implementation. By comparison, our native bosonic simulator containing $N$ modes simply requires $2N$ squeezing operations, $N$ displacement operations, and a maximum of $N(N-1)/2$ nearest-neighbor beamsplitter operations Reck1994 ; Cybenko2001 . This translates to a total of $O(N^{2})$ operations and a corresponding circuit depth of $O(N)$ when non-overlapping beamsplitters are applied simultaneously.

VII Discussion and Outlook

The superconducting platform demonstrated here is capable of successfully integrating all of the necessary components for performing a high-fidelity, scalable implementation of a practical computational task of interest. Looking ahead, there exist concrete steps toward scaling up, improving performance, and mitigating sources of error. High-Q superconducting modules controlling up to three modes have been experimentally demonstrated Chou2018 , which would allow simulations to encompass nonlinear triatomic molecules of $C_{s}$ symmetry. In general, a molecule composed of $M$ atoms will have $3M-6$ vibrational degrees of freedom ($3M-5$ for linear species), which sets the requirement for the number of modes needed in the simulator. A linear array of cavity modes (Fig. 1) is sufficient for preserving the efficiency of the implementation as discussed in Section VI. The photon loss rates of the microwave cavities in the device presented here can also be improved Reagor2016 ; Romanenko2018 . Finally, the self-Kerr Hamiltonian terms of the cavity modes may be cancelled with extra off-resonant pump tones, or three-wave mixing methods that avoid self-Kerr nonlinearities altogether may be used Frattini2017 .

The capabilities of the bosonic processor extend beyond the simulation shown for the estimation of FCFs. For applications in quantum chemistry, the full set of Gaussian operations enables the simulation of time-domain vibrational dynamics, potentially probing coherent evolution after a vibronic transition. Beyond the Gaussian framework, however, there also exist possibilities to utilize programmable self-Kerr nonlinearities to probe anharmonic effects in molecular dynamics Lloyd1999 . More generally, for condensed matter many-body systems, the tools utilized here enable the simulation of tight-binding lattice Hamiltonians, and the inclusion of controllable self-Kerr interactions extends the scope to extended Bose-Hubbard models Ma2019 ; Yanay2019 . The ability to use Josephson nonlinearities to manipulate microwave photons in these rich and diverse ways opens up promising avenues for the simulation of bosonic quantum systems in a superconducting platform.

Acknowledgements.

We thank Salvatore Elder, Hannah Lant, Zachary N. Vealey, Lidor Foguel, Kenneth A. Jung, and Shruti Puri for helpful discussions. This research was supported by the U.S. Army Research Office (W911NF-18-1-0212, W911NF-16-1-0349). We gratefully acknowledge support of the U.S. National Science Foundation under grants CHE-1900160 (VSB), CHE-1464957 (PHV), and DMR-1609326 (SMG), the NSF Center for Ultracold Atoms (ILC), and the Packard Foundation (LJ). Fabrication facilities use was supported by the Yale Institute for Nanoscience and Quantum Engineering (YINQE) and the Yale SEAS clean room.