Efficient Multiphoton Sampling of Molecular Vibronic Spectra on a Superconducting Bosonic Processor

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 [J. Huh et al. , Nat. Photonics 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 reprogrammability 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. DOI


I. INTRODUCTION
Simulation of quantum systems with quantum hardware offers a promising route toward understanding the complex properties of those systems that lie beyond the computational power of classical computers [1].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 [2][3][4].Boson sampling is another example of a computationally challenging task that can be performed by manipulating bosonic excitations [5,6].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 [7,8].The ability to perform QND photon number measurements not only enables the high-fidelity measurement of bosonic qubits via repeated detections [9,10], 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 beam splitter operations between separate modes [11][12][13].These capabilities motivate an alternative approach to performing both boson sampling [14] and bosonic quantum simulation protocols [15] 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 [16].The algorithm requires, in addition to beam splitters, 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 have made a linear optical implementation of this algorithm challenging [17].Recent work leveraging the bosonic nature of two phonon modes of a single trapped ion has been successful in generating more accurate spectra [18]; however, an efficient detection scheme capable of directly sampling from the exponentially growing bosonic Hilbert space remains elusive.Alternative approaches to this problem that directly emulate an electronic transition in a superconducting architecture via quenching dynamics have been proposed [19], but such quenching may be experimentally challenging to implement with high precision on a large scale.
Previous experimental work on simulating molecular vibronic spectra in superconducting systems has been restricted to a single vibrational mode [20,21].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 the first n max ¼ 16 Fock states per mode.This detection scheme, when operated without errors on a multiphoton distribution bounded by n max states 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 vibration 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 â ¼ ðâ 1 ; …; âN Þ for N vibrational modes may be expressed using the Duschinsky transformation [22] which can be decomposed into Gaussian operations via the Doktorov operator [23]: where correspond to a tensor product of single-mode displacement and squeezing operations across all N modes, respectively.RðUÞ is an N-mode rotation operator corresponding to a N × N rotation matrix U which can be decomposed into a product of two-mode beam splitter operations (see Appendix A) [24].For N ¼ 2, U is a twodimensional rotation matrix parametrized by a single angle θ.The set of dimensionless Doktorov parameters and U originate from molecular structural information in the different electronic configurations (see Appendix A).Applying ÛDok to an initial state jψ 0 i of the bosonic processor directly emulates the physical process of a pretransition molecular vibrational state experiencing a sudden change in the electronic 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 pretransition vibrational eigenstate, j⃗ ni ¼ jn; m; …i, with a final posttransition vibrational eigenstate, j⃗ n 0 i ¼ jn 0 ; m 0 ; …i: In spectroscopic experiments, the validity of this "sudden approximation" depends on the vastly different energy scales that typically characterize electronic and nuclear motions (≫10 000 cm −1 versus ∼100-1000 cm −1 , 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, ĉA and ĉB (Fig. 1).In our design, a coupler transmon tC dispersively couples to both cavities, enabling beam splitter [13] and squeezing operations through driven four-wave mixing processes.The coupler transmon is also dispersively coupled to a readout resonator rC .Displacement operations on the cavity modes are performed via resonant drives through local coupling ports.Additional ancillary transmon-readout systems f tA ; rA g and f tB ; rB g 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 jψ 0 i corresponding to the pretransition molecular vibrational state of interest.A vacuum state of both bosonic modes is prepared through feedback cooling protocols, and Fock states jn; mi are initialized using optimal control techniques [12] (see Appendix D).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 beam splitter.Both the single-mode squeezing and beam splitter operations utilize the four-wave mixing capabilities of the coupler transmon.Two pump tones that fulfill the appropriate frequency matching condition (ω 1 þ ω 2 ¼ 2ω A=B for the squeezing operation and ω 2 − ω 1 ¼ ω B − ω A for the beam splitter operation, where ω 1=2 are the two pump frequencies and ω A=B are the cavity frequencies) are sent through a port that primarily couples to the coupler transmon, which enacts the desired Hamiltonians [13]: where i ∈ fA; Bg.Importantly, the phase of the operations fφ; ϕ A ; ϕ B g implemented by these Hamiltonians is controlled by the phase of the pump tones, allowing the microwave control system to generate the correct family of beam splitter and squeezing operations for performing ÛDok .A set of transmon measurements is then carried out for the purpose of postselecting 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 ancillas otherwise dephases the cavities while coupler heating effectively halts the pumped operations by shifting the requisite frequency matching conditions.This postselection also serves to ensure that the ancillas 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, fn 0 ; m 0 g, gives the desired FCFs.The full set of controls is detailed in Table I; the relatively small error rates due to photon loss provide a sense of scale for attainable circuit depths while maintaining high-fidelity performance (see Sec. VI C).

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.2(a)]: , where the dispersive interaction strengths in our experiment are χ A ¼ 2π × 748 kHz and χ B ¼ 2π × 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.

A. Single-bit extraction
The first scheme [Fig.2(b)] maps a given joint cavity photon number population fn 0 ; m 0 g onto the joint state of the two transmons via state-selective π pulses.These pulses have frequencies , where the small second-order dispersive shift is also taken into account: The pulses are applied with a Gaussian envelope truncated at AE2σ t such that the bandwidth is approximately σ f ¼ 1=ð2πσ t Þ, where σ 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 σ t ¼ 1 μ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 jψi ¼ P i;j c ij ji; ji for a chosen set of photon numbers fn 0 ; m 0 g to probe: c n 0 j jn 0 ;ji ⊗ je;giþ c n 0 m 0 jn 0 ;m 0 i ⊗ je;ei; ð8Þ where jgi and jei 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.
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 max , where n max is the maximum number of Fock states considered for each mode (so for instance, n max ¼ 16 corresponds to j0i − j15i).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 [18].

B. Photon number resolved sampling
In optical systems, detectors that can simultaneously resolve the photon number and do so with a high quantum efficiency [25,26] are a necessary ingredient for protocols TABLE I. Full set of controls of the bosonic processor.The subscripts i ∈ fA; Bg correspond to operations on each cavity and their respective ancillary modules.κ A and κ 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Þ and ϵ i ðtÞ, respectively, generate a desired initial Fock state (see Appendix D). jn 0 ; m 0 i denotes a particular joint photon number to be measured in the cavities, and fb i g represent the bits associated with the binary decomposition of the photon number.For example, j5i ¼ j0101i.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 Ĥdrive ðtÞ ¼ such as quantum computing [27] and communication based on repeater networks [28], among others.The boson sampling algorithm implemented in this work also requires photon number resolved detection, as the output distribution usually has multiple photons per mode.
Our second measurement scheme implements singleshot photon number resolving detection [Fig.2(c)] that overcomes the scalability issue associated with the singlebit extraction method.At the heart of this technique is the dispersive Hamiltonian between each cavity mode ĉ and their respective transmon ancilla t: It can be shown that driving the transmon with numerically optimized waveforms ϵðtÞ, enables the QND mapping of any binary valued operator of the cavity Hilbert space onto the state of the transmon [11,12].Our strategy for measuring the number of photons for a given cavity state is to represent photon number in binary jni ¼ jΠ k max k¼0 b k i, where n ¼ P k max k¼0 2 b k and fb k g ∈ f0; 1g, and then to sequentially measure each bit b k on a given run of the experiment.This amounts to identifying a set of parity operators Pk in the cavity Hilbert space whose eigenvalues λ k;AE ¼ AE1 correspond to b k ¼ f0; 1g, respectively, with the following matrix elements ij: where b•c denotes the floor function.We synthesize a set of optimal control pulses ϵ k ðtÞ, k ∈ f0; 1; 2; 3g, that excite the transmon from its ground state jgi to the excited state jei conditioned on the cavity state's projection onto the parity operators Pk .In our experiment, these pulses are numerically optimized for the Hilbert space of each cavity up to n max ¼ 16 and have a duration between 800 and 1200 ns.
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 jψi ¼ P 15 n¼0 c n jni with bounded support within n max ¼ 16 into a definite Fock state jni with a probability jc n j 2 .In order to minimize errors due to decoherence when the transmon is excited, after the kth bit is mapped onto the transmon and measured, the transmon is reset to its ground state using real-time feedforward control to prepare for the measurement of the Here, instead of having a binary output, each detection module serves as a photon number resolving detector.Sequential QND measurements of the operators associated with the first four bits of each photon number's binary decomposition resolve up to 15 photons per mode.For a general state of the cavity, each measurement projects the state into the eigenspace of the outcome, ultimately projecting out a single Fock state with its corresponding probability.In the schematic, a sequence sampling the 6 photon component (j6i ¼ j0110i) from a displaced Fock state is shown.The state is first probabilistically measured to be in the b 0 ¼ 0 (even parity) subspace, thus projecting out only even photon numbers.The next measurement further projects this state into the b 1 ¼ 1 subspace, and so on, until the measurement converges on a single photon number.
(k þ 1)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 simply call sampling.These binary detectors optimally resolve the photon number in the cavities in N log 2 ðn max Þ measurements, but are prone to more errors due to transmon decoherence that are, in principle, correlated bit to bit.Errors in the measurement of individual bits can lead to unintuitive misassignments of the photon number-for instance, j8 ¼ 1000i would be misassigned as j0 ¼ 0000i if a single error misassigns b 3 .
We leave the task of characterizing these errors in full detail, optimizing the technique, and potentially applying deconvolution methods [29] to improve the accuracy of the sampled distribution as the subject of future work.
The advantage of the sampling protocol compared to the single-bit extraction method can thus be seen by comparing the required number of simulator runs for a desired statistical error on all output probabilities (see Supplemental Material [30]).Specifically, the latter method requires a factor of n N max more runs, directly reflecting the need to independently query each state in the entire Hilbert space.

V. SIMULATED PHOTOELECTRON SPECTRA
The full set of FCFs for a given electronic transition can provide, within the Condon approximation, the relative intensities of vibronic progressions appearing in corresponding photoelectron spectra.The algorithm exploited in this work extracts FCFs under the commonly deployed harmonic approximation.As such, the resulting relative intensities will deviate from those measured in actual molecular spectra owing to the effects of vibrational anharmonicity in the true PES.We simulate four different photoelectron processes starting in various vibrational initial states jψ 0 i ¼ jn; mi: Experimental results for photoionization of water starting in jψ 0 i ¼ j0; 0i and photodetachment of the ozone anion starting in jψ 0 i ¼ j1; 2i, respectively, are presented in Fig. 3 (see Supplemental Material for additional data [30]).We highlight the results for these two processes due to the long vibronic progressions they possess, thus verifying the capability of our simulator to generate and measure large photon numbers.The particular electronic state of the water cation considered here ð B2 B 2 Þ is the second excited state of doublet spin multiplicity, with the attendant electronic wave function having B 2 symmetry.The other 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 pretransition 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 ¼  .The abscissa scale corresponds to vibrational term values (energies in cm −1 ) within the final (posttransition) electronic PES calculated from the harmonic frequencies for symmetric-stretching and bending degrees of freedom: Solid lines depict theoretical FCFs, artificially broadened with Lorentzian profiles (10 cm −1 FWHM).Circles represent experimental data using the single-bit extraction (purple) and sampling (red) measurement schemes explained in the main text; statistical error bars for the latter measurement are not visible on this scale (see Supplemental Material [30]).Systematic errors associated with transmon decoherence during the selective π pulses are corrected for (see Appendix E).Additional errors are present in the sampled values, owing to decoherence effects during the binary decomposition measurement chain.
distances for the two simulated processes in Fig. 3 are D ¼ 0.049 (H 2 O) and 0.105 (O 3 ) for the single-bit extraction scheme and D ¼ 0.152 (H 2 O) and 0.148 (O 3 ) for the sampling scheme.The sampling distance will ultimately be the relevant figure of merit for evaluating the practical performance of this approach as it is scaled up.The other relevant metric is run time.As previously discussed in Sec.IV B, the single-bit extraction scheme requires a factor of n N max more runs compared to the sampling scheme for a desired statistical error.Our experiment operates at an effective repetition rate of roughly ∼300 Hz, thus requiring ∼7 h versus ∼100 s of data acquisition, respectively, for ∼3 × 10 4 samples of the photoionization of water.
Furthermore, these distance metrics are accompanied with a success probability due to postselection 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 multidimensional optimization problem that seeks to maximize the desired interaction rates while minimizing induced decoherence and dissipation rates [31].Self-Kerr Hamiltonian terms of the form ĤKerr =ℏ ¼ − P i∈fA;Bg ðK i =2Þĉ † i ĉ † i ĉi ĉi , photon loss, and imperfect state preparation account for errors that are undetected by postselection.The first two effects are captured through full time-domain master equation simulations (see Supplemental Material [30]).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 A. Classical methods
To properly contextualize the efficiency of the bosonic algorithm implemented in our work, it is useful to consider state-of-the-art computational methods for calculating Franck-Condon factors using classical computers.A pervasive and representative approach is one that utilizes generating functions to compute a single Franck-Condon overlap integral hn 1 ; n 2 ; …; n N jψ 0 i in terms of other integrals with fewer occupations fhn 1 − i 1 ; n 2 − i 2 ; …; n N − i N jψ 0 ig.This can be compactly represented through a recursion relation [32,33], which means that the total number of necessary integrals to compute a Franck-Condon factor with a total number of quanta Therefore, the formal task of obtaining the full Franck-Condon profile, which by itself contains an exponential number of entries corresponding to the size of the Hilbert space n N max , can be seen to be at least exponentially challenging.
In practice, however, electronic spectra of many polyatomic molecules initially cooled to near their rovibronic ground states typically exhibit vibronic progressions confined to a small region in the vibrational Hilbert space.This stems primarily from two effects.First, a single electronic transition in large molecules typically does not induce large shifts of the equilibrium coordinates for a significant fraction of all the normal modes, which keeps n max for most modes relatively low.Second, following symmetry considerations, there typically is structure in the Duschinsky rotation matrix which allows it to be approximated as block diagonal [34].This is significant because each block-diagonal subspace can be treated independently as they do not couple to the remaining vibrational degrees of freedom of the molecule.This effectively reduces the dimensionality of the problem from N to several subproblems with dimensionality of each block.In the limit where no mode mixing occurs, i.e., U ¼ 1, the full multimode Franck-Condon integrals can be easily and efficiently computed via the product of all single-mode Franck-Condon integrals.
Consequently, classical methods can implement convergence criteria that can provide a trade-off between accuracy and computational resources [35].In cases where the Duschinsky rotation matrix cannot be well approximated as block diagonal, however, recursive methods inevitably need to explore the entire Hilbert space and therefore become intractable for large system sizes.The advantage of using a quantum simulator that samples from the full distribution would then be to quickly identify the relevant FCFs with large weights, for which a classical computation can then be done for the specific FCFs of interest.Thus, the utility of the quantum simulator for obtaining Franck-Condon factors in large polyatomic molecules will depend on the configurational details for each electronic transition on a case-by-case basis.

B. Resource comparison: Bosonic versus qubit processor
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-1=2 quantum computer [36] needs to map the problem onto qubits and a univeral gate set.This first requires encoding the Hilbert space of size n N max onto n q ¼ Nlog 2 ðn max Þ qubits.The choice of n 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 [37], the approximate number of gates n g then needed to implement ÛDok using a universal qubit gate set to within an error ε is n g ¼ O½N 2 n 2 max log 3 ð1=εÞ.For our experiment with N ¼ 2 modes, taking n max ¼ 16 and desiring an error ε ¼ 5 × 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, where we do not yet have fault tolerant quantum processors.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 beam splitter operations [24,38].This translates to a total of OðN 2 Þ operations and a corresponding circuit depth of OðNÞ when nonoverlapping beam splitters are applied simultaneously.A possible 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.

C. Scalability and error budget
As shown in Fig. 1, a linear array of bosonic memories with nearest-neighbor coupling is sufficient for scaling to larger system sizes.In our architecture of fixed frequency cavity modes, the bilinear Gaussian operations are enacted via robust frequency converting four-wave mixing processes that obviate the need for any in situ frequency tuning.
As the hardware is scaled up, the fidelity of the individual operations will determine the fidelity of the overall simulation.Though different photoelectron processes will have different errors, we can consider a simplified model for quantifying performance with system size.We associate a success probability for each operation fρ i g where the index i encompasses state preparation (SP), displacements (D), squeezing (sq), beam splitters (BS), and measurements (M).For simplicity, we assume a uniform probability for each operation across all modes.We then specify a target success probability threshold ρ th that reflects the accuracy of the full simulation.The number of modes that can be accurately simulated for a given ρ th , therefore, can be determined by Each of these probabilities can be taken to be the average fidelity of each operation across a representative set of Doktorov transformations.Taking the expected bounds on the error rates of the operations due to photon loss in our experiment (Table I), while assuming a measurement error rate of 10 −2 and targeting ρ th ¼ 0.5, we get N ≈ 5. Modest improvements in cavity lifetimes [39,40] and further circuit optimization for engineering the bilinear interactions [31] can reduce the error rate of the squeezing and beam splitter operations to 10 −3 , which increases the number of modes to N ≈ 25.Beyond this, further reduction of the error rates or implementing bosonic error correction protocols that preserve bosonic statistics at the logical level [41] will be required for maintaining performance with increasing system size.

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 [42], 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.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 [43].
To simulate a broader class of molecular systems, the quantum simulator needs the capability of systematically incorporating inevitable deviations from harmonic behavior in the potential energy surfaces [19].This is certainly an interesting and exciting direction to pursue in future experiments, especially given that classical computational methods have difficulty incorporating such anharmonic effects even for relatively modest system sizes.Combining Josephson nonlinearities with external flux control, circuits with higher-order nonlinearities can be designed.These nonlinear modes can then represent distinct vibrational modes such that, when combined with bilinear [13] or higher-order interaction terms [44], a broader class of anharmonic molecular Hamiltonians can be represented accurately.More generally, for condensed matter manybody systems, the tools utilized here enable the simulation of tight-binding lattice Hamiltonians, and the inclusion of controllable self-Kerr interactions widens the scope to extended Bose-Hubbard models [45,46].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.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 No. CHE-1900160 (V. S. B.), No. CHE-1464957 (P.H. V.), and No. DMR-1609326 (S.M. G.), the NSF Center for Ultracold Atoms (I.L. C.), and the Packard Foundation (L.J.).Fabrication facilities use was supported by the Yale Institute for Nanoscience and Quantum Engineering (YINQE) and the Yale SEAS clean room.

APPENDIX A: OBTAINING DOKTOROV PARAMETERS
The Doktorov parameters originate from the physical properties of a given molecule in the two electronic states of interest.Specifically, it is the structural information of the molecular configurations and the relationship between the two that fully parametrize the problem.

Description of quantum-chemical analyses
Theoretical predictions of optimized equilibrium geometries (with imposed C 2v symmetry constraints), harmonic (normal-mode) vibrational displacements, and Franck-Condon parameters (Duschinsky rotation matrices and associated shift vectors) exploited the commercial (G16 rev.A.03) version of the GAUSSIAN quantum-chemical suite (Table II) [47], with canonical Franck-Condon matrix elements for specific vibronic bands being evaluated through use of the open-source ezSpectrum (version 3.0) package [48].All analyses relied on the CCSD(T) coupledcluster paradigm, which includes single and double excitations along with noniterative correction for triples.Dunning's correlation-consistent basis sets [49][50][51] of triple-ζ quality augmented by supplementary diffuse functions (aug-cc-pVTZ ≡ apVTZ) were deployed for all targeted molecules except water, where a larger doubly augmented, quadruple-ζ basis was employed (daugcc-pVQZ ≡ dapVQZ).
The Duschinsky rotation matrices and associated shift vectors provided by the commercial package GAUSSIAN are defined via where Q 0 and Q 00 are mass-weighted normal coordinates of the pre-and post-transition molecular configurations, respectively.Because our simulation considers the transformation from a vibrational state in the pretransition configuration to the post-transition configuration, we must redefine the Duschinsky rotation matrices and associated shift vectors accordingly:

Conversion from molecular parameters to Doktorov parameters
The Doktorov transformation as given in Eq. ( 2) of the main text is where for N ¼ 2 modes, the squeezing and displacement operations are defined as is the vibrational frequency of mode i in the pretransition (post-transition) TABLE II.Theoretically optimized molecular parameters.Vibrational frequencies for the symmetric-stretching and bending modes of each molecule in pre-(ν) and post-transition ( ν0 ) states are provided in wave numbers (cm −1 ), which is related to angular frequency ω via ν ¼ ω=2πc, where c is the speed of light.The rotation angle corresponding to the Duschinsky rotation matrix is defined in Eq. (A2).The shift vector K ¼ ðk 1 ; k 2 Þ is provided in mass weighted normal coordinates (where a 0 is the Bohr radius and m e is the electron mass) and reflects the relative displacement of equilibrium geometries between the two molecular configurations.The Duschinsky rotation matrix U generates the N-mode rotation operator RðUÞ.The multimode mixing elements implemented in this experiment are two-mode beam splitters, necessitating a decomposition of U into nearest-neighbor rotations, and thus R into nearest-neighbor beam splitters.RðUÞ becomes a product of two mode beam splitters parametrized by fθ k g and fi k ; j k g, a sequence of angles and rotation axes derived from the decomposition of

Molecular photoelectron process
We can then write The decomposition of U is analogous to generalizing Euler angles to SO(N); any rotation in R N can be written as a product of rotations in a plane R i k ;j k ðθ k Þ, known as rotations.Following an algorithm similar to that in Refs.[24,38], but simplified to real orthogonal matrices, produces a decomposition of U as a product of nearestneighbor rotation matrices.The Duschinsky matrix for N ¼ 2 is a single Givens rotation parametrized by an angle θ which is enacted with one beam splitter:

Optimization of squeezing parameters
The modification of the creation and annihilation operators under the mode transformation is given in [16] where The structure of L allows for a free scaling parameter η which leaves L invarant; namely, Given that the squeezing operations of the Doktorov transformation take ζ ¼ lnð ⃗ ΩÞ as inputs, an optimization may be performed, as done in Ref. [18], that minimizes the total amount of squeezing while leaving the unitary invariant.This is desirable as less squeezing corresponds to shorter gate times in the simulation, which reduces the overall error rate.Table III lists the final set of dimensionless Doktorov parameters used in the experiment.

APPENDIX B: THEORETICALLY PREDICTED
HAMILTONIAN TERMS

Derivation of ancilla-mediated operations
In this Appendix, we derive the ancilla-mediated beam splitter and single-mode squeezing interactions as shown in the main text as well as the associated ancilla-induced cavity frequency shifts.Derivations based on the perturbative four-wave frequency mixing enabled by a weak ancilla nonlinearity have been presented previously in Ref. [13].Here, we follow the formalism used in Ref. [31] and sketch the general results without assuming weak ancilla nonlinearity or weak pumps.We also give explicit expressions TABLE III.Dimensionless Doktorov parameters.All values are truncated to the precision that the operations are able to be implemented experimentally.for the strength of the engineered interactions in the case of weak pumps.We start from the Hamiltonian of the two bare cavity modes A and B coupled to the coupler transmon in module C: We emphasize that here the operators ĉA and ĉB are the annihilation operators for the bare cavity modes whereas in the main text the operators correspond to the dressed cavity modes that are weakly hybridized with the ancilla transmons.
ĤC is the Hamiltonian of the bare coupler transmon in module C.After expanding the transmon potential energy to fourth order in the phase across the Josephson junction and neglecting counterrotating terms, we obtain [52] where tð †Þ C again is the bare annihilation (creation) operator for the coupler transmon with frequency ω C and anharmonicity K C .
ĤI is the interaction energy between the coupler transmon and the two cavity modes.Neglecting counterrotating terms, it may be written as Ĥpump ðtÞ represents two pumps on the coupler transmon: Of primary interest to us is the dispersive regime where the cavity-transmon coupling strengths are much smaller than their detuning: jg A;B j ≪ jω A;B − ω C j [53].In this regime, we can treat the cavity-transmon interaction as a perturbation (while treating the remaining parts of the Hamiltonian exactly), and to second order in the interaction strength, we obtain an effective Hamiltonian after performing a unitary on where jΨ m i is the mth Floquet state that quasiadiabatically connects to the mth Fock state jmi of the bare transmon as the pumps are ramped up or down [31].At zero pump amplitudes, jΨ m i ¼ jmi.The first term in Ĥeff thus represents transmon-induced cavity frequency shifts δω A;m and δω B;m when the transmon is in jΨ m i.
The difference between δω A;m or δω B;m with different m leads to cavity-photon-number-dependent transmon transition frequencies.In particular, at zero pump amplitudes, the transmon's transition frequency from the ground to the first excited state decreases linearly with the cavity photon number with a proportionality constant: Physically, the factor jg i =δ i j 2 quantifies the participation ratio of cavity A or B in the coupler transmon.In the experiment, this factor is 0.3% for cavity A and 0.2% for cavity B.
Pumps on the coupler transmon can induce effective inter-or intracavity interactions (single-mode squeezing or beam splitter) denoted as V in Ĥeff .For the case of the beam splitter interaction For the case of single-mode squeezing ( Similar to the transmon-induced frequency shifts on the cavities, here both the strength and phase of the transmonmediated interactions depend on the state of the transmon.Of primary interest to us is the strengths of the transmonmediated interactions when the transmon is in the Floquet state jΨ 0 i.For weak drives, these strengths are where δ 1;2 ¼ ω 1;2 − ω C .In the case where the transmon anharmonicity K C is small compared to the detunings jδ 1;2;A;B j, the expressions above reduce to those obtained based on perturbative multiwave frequency mixing [13]. EFFICIENT MULTIPHOTON SAMPLING OF MOLECULAR … PHYS.REV.X 10, 021060 (2020) 021060-11 Note that the interaction strengths presented in the main text and the rest of the appendixes refer to the values of g BS;0 and g sq;i;0 .For weak drives, Eqs.(B9) and (B10) show that the strengths of the engineered beam splitter and single-mode squeezing increase linearly with both drive amplitudes Ω 1;2 .For strong drives, this dependence becomes nonlinear in Ω 1;2 and can be accurately captured using Floquet theory for the driven transmon [31].We have verified that the experimentally measured beam splitter and single-mode squeezing rates match the expressions (B9) and (B10) for weak drives and the full Floquet analysis at strong drives.

Transmon-induced cavity Kerr
Another important effect and a source of infidelity is the cavity nonlinearity induced by the transmons.To fourth order in the cavity-transmon coupling, this nonlinearity is a Kerr nonlinearity and has the following form: where K A;m and K B;m are the self-Kerr of cavities A and B and K AB;m is the cross-Kerr between cavities A and B when the transmon is in the state jΨ m i.First, we consider the case in the absence of pumps.Of interest to us is the cavity Kerr when the transmon is in the ground state j0i: Also of interest to us is the difference between K i;0 and K i;1 .This difference leads to a nonlinear dependence of the transmon transition frequency on the cavity photon number.This difference is usually denoted as We note that there is also a contribution to χ 0 iC from a term in the sixth-order expansion of the transmon cosine potential, but for ω C ≫ jδ i j this correction is negligible.
Here we have only considered the cavity Kerr induced by the coupler transmon.In general, the transmon ancillas in modules A and B also induce Kerr in their respective cavities.The total Kerr of each cavity will then be the sum of all contributions.
In the presence of pumps on the coupler transmon, the cavity Kerr can be strongly modified due to a relatively strong hybridization between cavity photons and excitations of the coupler transmon.To illustrate this effect, we consider as an example the pumps used in generating the beam splitter interaction between the two cavities.For the choice of pumps used in the experiment, the sum of the frequency of cavity A and the higher-frequency pump is close to the frequency of transition from transmon ground to the second excited state: ω A þ ω 2 ≈ ω 02 .As a result, the cavity photons become relatively strongly hybridized with the second excited state of the transmon, thus modifying their nonlinearity.Using a sixth-order perturbation theory (fourth order in g A and second order in Ω 2 ), we find that the modification to the cavity Kerr is , and ω 02 is the Stark shifted transmon transition frequency from the ground to the second excited state.The above expression, which applies for small jΔj, qualitatively captures the observed enhanced self-Kerr of cavity A in the experiment during the beam splitter operation.Comparing this expression with that of the bare cavity Kerr K A;0 without pumps, we see that δK A;0 becomes comparable to K A;0 when K C jΩ 2 =δ 2 j ∼ jΔj.
We note that such dependence of the cavity Kerr on the drive parameters also potentially provides a knob to control the cavity Kerr for the purpose of simulating nonlinear bosonic modes.

APPENDIX C: SYSTEM CHARACTERIZATION 1. Calibration of Gaussian operations
In the dispersive regime, the transition frequency ω l t i of ancilla ti depends on the photon number l in the respective cavity: where ω 0 t i is the ancilla frequency when there are no photons in its respective cavity and χ i and χ 0 i are the dispersive shifts originating from fourth-and sixth-order Hamiltonian terms, respectively, as introduced in the main text.Using this, the photon number population of each cavity can be extracted via π pulses selective on each photon number after applying various strengths of each operation (Fig. 4).These populations are then fit to the corresponding expected models, including an overall offset and scaling factor to take into account errors due to ancilla relaxation and readout imperfections (Table IV).For the beam splitter, we assume an effective detuning between cavities A and B in a frame where δ BS ¼ 0 if the beam splitter resonance condition is satisfied.
Transmon heating leads to fluctuations in δ BS , which dephases the beam splitter operation with a dephasing rate: Thus, the oscillating populations of a single photon in each cavity P 10=01 is given to leading order in κ A;B =g BS and κ BS ph =g BS by the expression in Table IV, where κ ¼ ðκ BS A þ κ BS B Þ=2 and κ BS A;B are the effective linewidths of cavities A and B during the beam splitter operation.

Static and pump-induced self-Kerr Hamiltonian terms, −ðK
B , are estimated using the protocol detailed in Ref. [54] (Fig. 5).For the pumpinduced cases, one of the two pumps is detuned by δ ¼ 20 and 50 kHz for squeezing and beam splitter operations, respectively, to make the engineered interaction off resonant.We assume that the induced self-Kerr is not a strong function of this detuning.
Static and pump-induced cavity decay rates are measured via T 1 experiments (Fig. 6).A single photon is prepared in each cavity, followed by either a delay or an offresonant pumped operation (with the same detunings as above).Again, we assume that the pump-induced decay rates are not a strong function of the pump detuning.The ancillas are then flipped via selective π rotations conditioned on n ¼ 1 photon.In both cases, the data are postselected on the ancilla being in the ground state before the selective π rotation.We attribute the higher decay rate to beginning with a single photon in cavities A (top, jψ 0 i ¼ j1; 0i) and B (bottom, jψ 0 i ¼ j0; 1i).The length of two beam splitter pump tones is varied, and the probability that the photon remains in the cavity that it started in is plotted over time.
the hybridization of the cavities with the shorter-lived coupler transmon.Measured cavity Kerr and T 1 values are given in Table V.

APPENDIX D: CIRCUIT IMPLEMENTATION
The full quantum circuit implemented in our experiment is shown in Fig. 7. State preparation in our experiment (Fig. 8) consists of first performing measurement-based feedback cooling of all modes to their ground state (this protocol is described in full detail in the Supplemental Material of Ref. [10]).For preparing Fock states, optimal control pulses are then played that perform the following state transfers: These state transfers, however, suffer a finite error probability on the order of a few percent due to decoherence during the operation.This error is suppressed by performing a series of QND measurements of each cavity photon number and postselecting on outcomes that verify that the correct state was prepared.This is done via k selective π rotations on the ancilla transmons conditioned on the desired photon numbers in the cavities, followed by transmon measurements-even if the desired state is joint vacuum.The final data are postselected on the ancilla measurement outcomes being ð"e; ""g"Þ ⊗k=2 for both modules, where k is chosen to be even.In our experiment, we choose k ¼ 2 for the single-bit extraction measurement scheme and k ¼ 6 for the sampling measurement scheme.

APPENDIX E: CORRECTING SYSTEMATIC ERRORS DUE TO TRANSMON DECOHERENCE DURING SINGLE-BIT EXTRACTION
Errors due to ancilla decoherence during the single-bit extraction measurement scheme may be systematically calibrated out.Specifically, decay and heating events during selective π rotations and readout errors result in a systematic bias in the final estimate of the photon number population.For the case of a single ancilla qubit coupled to a single cavity, these effects result in a reduction of contrast for a Rabi experiment when both the ancilla and the cavity are prepared in their ground state (Fig. 9).
When using this pulse to infer cavity photon number populations, we assume that there is no photon number dependence to either the Rabi or decoherence rates of the ancilla.Under this model, we can relate the measured probabilities ⃗ Q to the true probabilities ⃗ P via (c) Single-bit extraction.Selective π pulses ðR π Þ flip each ancilla transmon conditioned on having n 0 and m 0 photons in cavities A and B, respectively, for a given run of the experiment.The ancillas are then simultaneously read out using standard dispersive techniques.Subsequent runs of the experiment thus need to scan n 0 and m 0 over the photon number range of interest up to the desired n max .(d) Sampling.Optimal control pulses are designed to excite each ancilla transmon from jgi to jei conditioned on the value of the binary bits b i of each cavity state, followed by dispersive readouts.Here, we measure the first 4 bits on a given run of the experiment, thus resolving the first 16 Fock states for each cavity.Real-time feed-forward control is used to dynamically reset the state of the ancilla in between bit measurements to minimize errors due to ancilla relaxation.Measurement-based feedback cooling techniques prepare the full system in its ground state (i.e., both cavities in j0i and all transmons in jgi).Optimal control pulses of the form listed in Table I of the main text are played simultaneously on each module to prepare a desired photon number state, followed by a set of k check measurements.

⃗
where f and t are the probabilities of assigning the ancilla measurement to the excited state when it is prepared in the ground and excited states, respectively.Thus, inferring the true probabilities from the measured probabilities is a relatively straightforward task.
For two modes, however, the problem becomes more complicated as a measurement of a joint probability relies on shot-by-shot correlations of the individual ancilla outcomes.Thus, false positive counts due to heating and readout errors lead to misassignment in a nonlinear fashion.We can again write what a given joint measured probability Q nm is in terms of the true distribution P nm : Q nm ¼ t A t B P nm þt A f B P n m þf A t B P nm þf A f B P nm : Equation (E2) may be solved for P nm by noting that It is worth noting that this requires Q nm to be a square matrix, which translates to measuring both n 0 and m 0 up to a prespecified n max .

t f
Calibrated amplitude Excited state probability FIG. 9. Calibration of systematic measurement errors using selective ancilla pulses.A standard Rabi calibration experiment of a selective pulse used for measurement (here shown for ancilla B).The maximum probability t is limited by decoherence of the transmon during both the pulse and the readout.The floor f is set by the probability of heating out of the ground state during both the pulse and the readout.

3 FIG. 2 .
FIG.2.QND measurements of cavity photon number.(a) Depiction of two 3D λ=4 coaxial cavities dispersively coupled to individual ancillary modules (outlined in black), each consisting of an ancilla qubit and a readout resonator.From the point of view of the measurement, the ancillary modules serve as reconfigurable black boxes used to detect the photon number in each cavity.The dotted trace serves to illustrate that the state of the two cavities can generally be entangled.(b) Single-bit extraction.In this scheme, on a given run of the experiment, each ancilla is excited conditioned on a predetermined photon number n in its respective cavity.(c) Sampling.Here, instead of having a binary output, each detection module serves as a photon number resolving detector.Sequential QND measurements of the operators associated with the first four bits of each photon number's binary decomposition resolve up to 15 photons per mode.For a general state of the cavity, each measurement projects the state into the eigenspace of the outcome, ultimately projecting out a single Fock state with its corresponding probability.In the schematic, a sequence sampling the 6 photon component (j6i ¼ j0110i) from a displaced Fock state is shown.The state is first probabilistically measured to be in the b 0 ¼ 0 (even parity) subspace, thus projecting out only even photon numbers.The next measurement further projects this state into the b 1 ¼ 1 subspace, and so on, until the measurement converges on a single photon number.

1 2 PFIG. 3 .
FIG. 3. Experimental Franck-Condon factors.Measured data for (a) photoionization of water to the ð B2 B 2 Þ excited state of the starting in the vacuum (vibrationless, n ¼ 0, m ¼ 0) state and (b) the photodetachment of the ozone anion to the ground state of the neutral species O − 3 → hν O 3 þ e − starting from a vibrational eigenstate possessing one quantum of symmetric-stretching and two quanta of bending excitation (n ¼ 1, m ¼ 2)

FIG. 4 .
FIG. 4. Calibrations for Gaussian operations.(a) Displacement calibrations for cavities A (top) and B (bottom) starting in vacuum.Here, the amplitude of a resonant pulse is varied.(b) Squeezing calibration for cavities A (top) and B (bottom) starting in vacuum.Here, the length of two squeezing pump tones is varied.Legends indicate population in photon number n. (c) Beam splitter calibration for beginning with a single photon in cavities A (top, jψ 0 i ¼ j1; 0i) and B (bottom, jψ 0 i ¼ j0; 1i).The length of two beam splitter pump tones is varied, and the probability that the photon remains in the cavity that it started in is plotted over time.

FIG. 5 .FIG. 6 .
FIG.5.Estimation of intrinsic and pump-induced self-Kerr.In all plots, the y coordinate corresponds to the effective frequency of a coherent state with average photon number n on the x coordinate.The slope determines the self-Kerr, and the offsets reflect pumpinduced Stark shifts.Experiments for estimating the in the (a) of pumps, (b) presence of off-resonant squeezing pumps, and (c) presence of off-resonant beam splitter pumps for cavity A (top panels) and B (bottom panels).

FIG. 7 .
FIG. 7. Circuit implementation of the Franck-Condon simulation.(a) Overview of the quantum simulation algorithm, consisting of state preparation, unitary Doktorov transformation, and measurement.A set of verification measurements is performed after the unitary Doktorov transformation for the purpose of postselecting the final data on measuring the transmons in their ground state.(b) The twomode circuit decomposition of the Doktorov transformation used in this experiment.The nonlinearity of the coupler transmon is primarily utilized to perform all three pumped operations, though in principle that of the ancilla transmons could have been used as well.(c)Single-bit extraction.Selective π pulses ðR π Þ flip each ancilla transmon conditioned on having n 0 and m 0 photons in cavities A and B, respectively, for a given run of the experiment.The ancillas are then simultaneously read out using standard dispersive techniques.Subsequent runs of the experiment thus need to scan n 0 and m 0 over the photon number range of interest up to the desired n max .(d) Sampling.Optimal control pulses are designed to excite each ancilla transmon from jgi to jei conditioned on the value of the binary bits b i of each cavity state, followed by dispersive readouts.Here, we measure the first 4 bits on a given run of the experiment, thus resolving the first 16 Fock states for each cavity.Real-time feed-forward control is used to dynamically reset the state of the ancilla in between bit measurements to minimize errors due to ancilla relaxation.

FIG. 8 .
FIG. 8. Circuit implementation of heralded state preparation.Measurement-based feedback cooling techniques prepare the full system in its ground state (i.e., both cavities in j0i and all transmons in jgi).Optimal control pulses of the form listed in TableIof the main text are played simultaneously on each module to prepare a desired photon number state, followed by a set of k check measurements.
FIG.1.Circuit schematic of the superconducting bosonic processor.The device consists of two microwave cavity modes (blue ĉA , and green ĉB ) which represent the symmetric-stretching and bending modes of a triatomic molecule in the C 2v point group, respectively.Ancilla measurement and control modules (shaded in blue and green) consisting of a transmon qubit ð tA ; tB Þ and readout resonator ðr A ; rB Þ couple to each cavity mode for state preparation and measurement.The coupler module (shaded in red) consists of a coupler transmon tC and readout resonator rC .The coupler transmon is used to facilitate bilinear Gaussian operations on the two cavity modes through four-wave mixing, and the readout resonator is used for both characterization and postselection.This configuration is extensible to a general linear array of N cavity modes with nearest-neighbor coupler modules and N ancillary modules, as depicted in light gray.
where fd i g are the vector elements of d in Eq. (A3).

TABLE IV .
Calibrated rates of Gaussian operations.The amplitude of a displacement operation of fixed length τ ¼ 72 ns is calibrated for generating a coherent state with α ¼ 1 [Fig.4(a)].The rates for the squeezing and beam splitter operations are extracted from the fits [Figs.4(b) and 4(c)].

TABLE V .
Estimated cavity Kerr and T 1 values.The beam splitter decay rates are extracted from the fit performed in the calibration of the operation in Fig.4(c) assuming that κ BS A ¼ κ BS B ¼ κ.