Approximating Vibronic Spectroscopy with Imperfect Quantum Optics

We study the impact of experimental imperfections on a recently proposed protocol for performing quantum simulations of vibronic spectroscopy. Specifically, we propose a method for quantifying the impact of these imperfections, optimizing an experiment to account for them, and benchmarking the results against a classical simulation method. We illustrate our findings using a proof of principle experimental simulation of part of the vibronic spectrum of tropolone. Our findings will inform the design of future experiments aiming to simulate the spectra of large molecules beyond the reach of current classical computers.


INTRODUCTION
The ability of classical computers to simulate quantum mechanical processes is severely limited, because the size of the Hilbert space used to describe a quantum system increases exponentially with the size of the system. Quantum chemistry is thus expected to benefit greatly from the development of quantum simulators and quantum computers. For example, the calculation of molecular energies can in principle efficiently be done on a quantum computer involving relatively few qubits [1].
It has recently been shown that the estimation of molecular vibronic spectra can in principle be done efficiently using a quantum optics simulator [2,3], whereas there is no known efficient classical algorithm for this task. Vibronic spectra, arising from simultaneous electronic and vibrational transitions in molecules, play an important role in determining the optical and chemical properties of those molecules. In addition to being useful for fundamental research in molecular physics and chemistry, calculating these spectra helps in assessing the performance of different molecules for applications in photovoltaics [4], biology [5], and other forms of industry [6].
The protocol for estimating vibronic spectra is considerably simpler than other quantum simulation protocols which require particle interactions or even full quantum computing. Since both the vibrational modes of a molecule and the modes of light can be approximated as quantum harmonic oscillators, the physics describing vibrational transitions has a direct analogy in quantum optics [7]. An optics experiment that implements simple operations such as squeezing and interferometry, and performs photon number measurements on the resulting state, is sufficient to simulate this physics and therefore reconstruct the vibronic spectrum of a molecule. This protocol scales linearly in the number of vibrational modes to be simulated, does not require any postselection, and can be implemented with readily available tools in quantum optics. It is inspired by the boson sampling protocol [8] which has been demonstrated experimentally [9][10][11].
Other experimental platforms than quantum optics that make use of quantum harmonic oscillators can also be used to simulate vibronic spectra, as shown by recent experiments using superconducting devices [12] and trapped ions [13]. The choice of platform imposes practical limitations on what can be achieved. Quantum optics is promising due to the availability of good sources of squeezing [14] and the large number of modes which can be manipulated, interfered, and measured [15,16].
However, any experimental implementation of a quantum algorithm on a platform that does not have faulttolerant architecture is necessarily degraded by imperfections in the system operations. This is a potential limitation to the performance of all specialized quantum processors. In optical platforms, scattering and absorption losses, mode-mismatches, detector noise and unanticipated correlations may all contribute to less than ideal operation. The presence of experimental imperfections in any platform was not considered in the original proposal. These imperfections can be expected to affect the simulation, possibly reducing both the accuracy and precision of the results.
In this work, we show that imperfect experiments can still yield approximations of molecular vibronic spectra, and we illustrate this result with a quantum optics simulation of the partial vibronic spectrum of tropolone. We Depending on the energy of the absorbed photon, different vibrational states are excited, leading to complex spectra (bottom). The heights of the peaks depend on the overlaps between the ground state of the molecule and the excited vibrational states of the excited molecule. B) We model vibronic transitions using a harmonic approximation of the vibrational modes. The harmonic oscillators describing the excited state (in blue) are squeezed by S and displaced by D with respect to the ground state (in red). The overlaps between the different Fock states determine the heights of the spectral peaks. C) We simulate this process using a quantum optics experiment with squeezing S and displacement D (top). Each optical mode is mapped onto a vibrational mode of the molecule; in our case we consider two coupled vibrational modes of tropolone. The probabilities of measuring photon number outcomes using photon number resolving detectors (PNRD) are mapped onto the heights of the peaks in the spectrum (bottom).
first describe in more detail the analogy between vibronic spectra and quantum optics using a specific transition in tropolone as an example. We then introduce our experiment, which we use to illustrate the three key findings that support our main result. First, we show that optimal experimental parameters can be efficiently determined for any experimental setup. Second, we show that the error on experimental estimates of vibronic spectra can be efficiently bounded. Third, we show that such experiments can outperform a classicality criterion. Finally, we examine the results of our experiment in light of these findings.

THE 370 NM TRANSITION IN TROPOLONE
We first illustrate the connection between vibronic spectrocoscopy and quantum optics using the example of tropolone (C 7 H 6 O 2 ), which is a molecule contributing to the taste and color of black tea [17] (see Fig. 1). In the 370 nm electronic transition in tropolone, the change in molecular configuration caused by the electronic excitation distorts the vibrational modes and couples them to each other. In the following, we focus on two of these modes which couple only to each other due to selection rules [18], and use a harmonic approximation of the potential wells corresponding to the vibrational degrees of freedom. We can write the change in mass-weighted normal coordinates (q 1 , q 2 ) of the two modes under study as [18,19]: According to the Franck-Condon principle [20,21], the intensity of a given vibrational transition is proportional to the overlap between the wave function of its initial vibrational state and that of its final vibrational state. If the initial vibrational state is the ground state and the final state has m 1 energy quanta in mode 1 and m 2 energy quanta in mode 2, then this overlap can be written as: whereÛ Dok is the operator implementing mode transformation 15, known as the Doktorov operator [22]. P (m 1 , m 2 ) is then the normalized intensity of the transition at frequency m 1 ω 1 + m 2 ω 2 , where ω 1 = 176 cm −1 and ω 2 = 110 cm −1 are the excited state vibrational frequencies of modes 1 and 2. P (m 1 , m 2 ) is known as the Franck-Condon factor for this transition.
We now consider the quantum optics analogy to the vibrational transition described above. Equation 15 can be interpreted in quantum optics as a Bogoliubov transformation between two optical modes. If the initial state of these two modes is vacuum, then this transformation can be achieved in quantum optics via two single-mode squeezing operations and a beam splitter [23]. P (m 1 , m 2 ) becomes the probability of detecting m 1 photons in mode 1 and m 2 photons in mode 2. An ideal quantum optics experiment that prepares two appropriate single mode squeezed vacuums (SMSV), interferes them on a beam splitter with the appropriate reflectivity, and measures the resulting photon number distribution using photonnumber resolving detectors can therefore be used to estimate the Franck-Condon factors associated with the 370 nm transition in tropolone.
We choose tropolone to illustrate our findings for the following reasons. First, we note that the vibronic transition described by equation 15 does not include a displacement term, which is present in many other molecules. Displacements can be implemented in quantum optics using classical laser light, and do not affect whether an optical state is classical, as defined by having a regular P-function in phase space [24], or not. Furthermore, the squeezing parameters for the two modes in the ideal experiment are 0.19 and 0.72, which is quite large compared to most other molecules. The absence of displacement and these large squeezing parameters allow us to focus our analysis on the quantum mechanical aspects of the experiment. These factors also allow us to highlight the impact of imperfections such as loss, which squeezing is strongly affected by, on an experiment.

EXPERIMENTAL SETUP
Our setup ( Fig. 2A) consists of a 100 kHz pulsed laser at 780 nm that pumps a periodically poled potassium titanyl phosphate (KTP) waveguide to produce a degenerate two mode squeezed vacuum (TMSV) at 1560 nm [25], the two modes of which we separate and couple into the two inputs of a fiber beam splitter. We then use two transition-edge sensors (TES) [26] to sample from the photon number statistics from the two outputs of the beam splitter. The frequency of the occurrence of photon numbers m 1 and m 2 in modes 1 and 2 gives a direct estimate of the joint probability in equation 2 and thus of the Franck Condon factors for the transition towards the states with m 1 and m 2 vibrational quanta.
Our setup deviates in the following ways from the ideal setup described above. Firstly, we approximate the two independent SMSVs that are required in an ideal setup using a TMSV, since TMSVs are experimentally simpler to generate and mode-match. The TMSV can be converted into two identical SMSVs using the beam splitter in our setup. Furthermore, the two modes are not exactly identical and so do not interfere perfectly on the beam splitter; the Hong-Ou-Mandel dip [27] between the two modes has a 94% visibility. Another significant imperfection is the limited system efficiency from the photon source to the detectors; approximately 60% of the generated photons are not detected. This loss mostly comes from the low coupling efficiency from the photon source into single mode fiber, and is expected to degrade the squeezing. Finally, our TES detectors are noisy, with a 0.2% dark count probability and a probability of 0.1% of detecting pump photons that leak through our setup. Since photons from the pump have twice the energy of the downconverted photons and TESs are energy-resolving detectors, these events register as two-photon events.

DESIGNING THE EXPERIMENT
In the absence of imperfections, equation 15 can be used to determine which squeezing parameter and beam splitter reflectivity to use in an experiment [23]. With the above-mentioned imperfections, however, these experimental parameters will not yield the target optical state. In the following, we demonstrate a method for determining, in a scalable way and in the presence of imperfections, experimental parameters that yield a state that is close to the target optical state.
We first require a metric to compare an experimentally generated state to the target state. We recall that, in the general case, we cannot directly compare the experimentally generated state to the target state in the photon number basis because there is no known efficient classical algorithm for calculating these photon number statistics. We therefore use a phase space description of these states and in particular the Gaussian state formalism [28,29]. This formalism applies to states that have a Gaussian quasi-probability distribution in phase space and to operations that maintain the Gaussian nature of these states. Since the initial state (vacuum) and all the optical operations are Gaussian, as well as the sources of imperfection, both the target state and the experimentally generated optical state can be described as multimode Gaussian states. A Gaussian state can be efficiently characterized experimentally in the phase space description, and although the photon number statistics of Gaussian states cannot be efficiently calculated, the fidelity between two such states can be [30].
To apply this fidelity metric to an experiment, we need a good phase space description of the optical states that are produced, which requires accurate characterization of the experimental setup. This characterization can be done efficiently; losses can be measured using classical light [31], mode overlaps can be determined by their pairwise Hong-Ou-Mandel dip visibility, and detectors and squeezers can be characterized independently. With a good description of the experimental setup, the experiment can be modeled using the Gaussian state formalism and the fidelity to the target state can be efficiently calculated. Furthermore, numerical optimization can then be used to find experimental parameters that maximize the fidelity.
We apply this characterization and optimization procedure to our experiment to determine which squeezing parameter r and beam splitter reflectivity t BS produce the state that maximizes the fidelity to the target state. terization results, and the values of r and t BS yielded by numerical optimization that maximize the fidelity within our model. We estimate the fidelity to the target state to be 0.890(1).

BOUNDING THE ERROR
Due to the best experimentally generated state having non-unit fidelity to the target state, we can expect the resulting estimate of the vibronic spectra to contain some degree of error. For our experimental results to yield a useful approximation of vibronic spectra, we require a method for bounding this error.
The fidelity can be used to provide such a bound. It can be shown (see Supplementary Material) that the fidelity F between any two states ρ 1 and ρ 2 and their trace distance in photon number statistics P 1 and P 2 are related by: Equation 3 now gives us an efficiently calculable measure for determining how close experimental photon number statistics are to those of the ideal state, which we can use to bound the error on an experimental estimate of Franck-Condon factors.
We apply equation 3 to our experiment. By directly using the value for the fidelity that we determined earlier for our model of the experiment, we find a bound for the trace distance of 0.455.

CLASSICALITY CRITERION
Quantum optical experiments aiming to estimate vibronic spectra are worthwhile if they outperform known classical algorithms. While there is no known efficient exact classical algorithm for calculating vibronic spectra, some approximation strategies exist. One case-by-case strategy involves guessing which transitions are likely to contribute the most to the spectrum and only calculating the corresponding Franck-Condon factors [32].
Alternatively, we consider the following more general efficient approximate classical algorithm. It can be shown that there is a unique classical optical state, defined as having a regular P-function in phase space [24], which maximizes the fidelity to the ideal target state within the space of all classical states (see Supplementary Materials). The photon number statistics of this state can efficiently be sampled from [33], and these statistics can then be used to estimate the target vibronic spectrum to within the error bounded by equation 3. Due to the resemblance between this classical algorithm and the quan- We simulate an ideal setup that is only affected by loss (black solid line), to which we add our noisy detectors (blue dashed line), then replace the SMSVs by a TMSV (red dotted line), then add our measured partial overlap between the modes (orange dashed-dotted line). The cross and the circle respectively indicate the best SMSV for our level of loss and the best SMSV without loss but with noisy detectors, yielding the estimated Franck-Condon factors shown in Table 1. The orange square indicates the experimental parameters used for our experiment. For this analysis, we assume equal loss in both modes, hence the discrepancy between the orange square and the position of the orange dashed-dotted line.
tum simulation protocol which also involves sampling from optical states, we propose to use it as a benchmark against which experiments may be compared.
This classical approximation strategy can also be used as a classicality witness: any experimental state with a higher fidelity must have a non-regular P-function. The best classical state can therefore be used to define a criterion for experimentally demonstrating a quantum advantage. Any experimental optical state that beats the witness is both a non-classical state and produces a better approximation of vibronic spectra than would be possible with any classical state.
We apply this analysis to our experiment. In our case, we found that the best classical state to approximate the spectrum of tropolone is vacuum, due to the absence of displacement in the transition under study. Vacuum has a fidelity of 0.879 to the target state. Our experimentally generated state, with a fidelity of 0.890(1) to the target state, beats this witness by roughly 10 standard deviations and is therefore a quantum state. Furthermore, using equation 3, vacuum can be used to estimate the correct Franck-Condon factors to within a trace distance of 0.476, compared to 0.455 for our experiment. Our quantum state therefore does yield a smaller error bound on these Franck-Condon factors than the best classical state. Figure 3 indicates the effect of different sources of imperfection on the fidelity in our experimental setup. An ideal experiment consists of SMSVs and ideal detectors (black solid line). The additional dashed and dotted lines indicate the effect of additional imperfections that must be accounted for in our experiment. The orange dashed and dotted line indicates the effect of all the parameters that we account for in our analysis, and the orange square corresponds to our experiment. The flat green solid line indicates the maximum fidelity that can be achieved with the best classical state.
We note that our classicality criterion is relatively tolerant of experimental imperfections to simulate transitions that involve large amounts of squeezing. An enhancement over the best classical state can be achieved with values of loss up to 90%, with noisy detectors, and with the use of a TMSV instead of two SMSVs. For loss exceeding 90%, the noise on our detectors degrades the fidelity below that of the best classical state. For transitions in other molecules that involve less squeezing, such as the S 0 ( 1 A 1g ) → S 1 ( 1 B 2u ) transition in benzene [34] or the SO − 2 → SO 2 transition [35], we find that our level of detector noise prevents us from outperforming this criterion for any level of loss (see Supplementary Materials).

EXPERIMENTAL RESULTS AND DISCUSSION
We now discuss our experimental results in light of the findings described above. We collected 1638370 samples from the photon number statistics of our state from which we estimate the Franck-Condon factors of the transition under study.   Table 1. The loss leads to an overestimate of the Franck-Condon factor corresponding to vacuum, both for our experiment and for a theoretical lossy SMSV. Furthermore, whereas the Franck-Condon factors for odd numbers of excitations should be 0 due to the difference in symmetry between the ground state and the odd-numbered excited states, our experiment finds these to be non-zero due to photons which would correspond to higher order Franck-Condon factors being lost. The use of a TMSV instead of two independent SMSVs causes us to experimentally find photon numbers that are roughly symmetric in both modes, to within the imbalance in the loss in the two arms.
We note that although the high squeezing and absence of displacement in a simulation of tropolone has allowed us to highlight the issue of imperfections in an experiment, by the same token our simulations result in a large error in estimating Franck-Condon factors. Many transitions in other molecules, such as the SO − 2 → SO 2 transition that was simulated using trapped ions [13], require less squeezing and some amount of displacement, so that the resulting Franck-Condon factors arise mostly from the displacement term. Squeezing is vulnerable to loss and displacement is not, since if the loss in an experimental setup is known then the displacement can be adjusted to exactly compensate for it. Such molecules can therefore be simulated with quantum optics with a much higher fidelity; for the SO − 2 → SO 2 transition we find that with our level of loss, noisy detectors, and with the use of a TMSV instead of two SMSVs, a fidelity of 0.993 could in principle be achieved in an experimental setup that also implements displacements. CONCLUSION We have demonstrated a method for quantum optical approximations of molecular vibronic spectra, which applies to realistic, imperfect experimental setups. We have shown that the optimal experimental parameters for a given setup can be efficiently determined, that the resulting error on the estimate of Franck-Condon factors can be bounded, and that such experiments can be shown to beat a classicality criterion. We illustrated these results using a proof of principle experiment that simulated a partial vibronic spectrum of tropolone.
Our results provide a roadmap for scaling simulations of vibronic spectra to much larger molecules with current quantum optics technology, for example using recent advances in fiber-loop based experiments [15,16] and integrated photonics [10]. We also note that our approach for dealing with imperfections can be applied to the other platforms which have been proposed for vibronic spectra estimation [12,13]. We envisage that our work will be useful for efficiently approximating the spectra of large molecules that are outside the reach of classical computers.

Competing financial interests
The authors declare no competing financial interests.

Supplementary Materials
Accompanies this paper. The underlying data used in this work is available upon request from the authors.

SUPPLEMENTARY MATERIALS
Bounding the estimate on vibronic spectra using the fidelity In this section, we show how equation 3 in the main text is derived.
Both the target state and an experimentally generated optical state can be described as multimode Gaussian states, since they are generated from vacuum using only Gaussian operations such as squeezing, displacements, and rotations, as well as, in an experiment, loss and noise from detectors and imperfect overlaps between optical modes. Within the Gaussian state formalism [28], these states can be completely described by their covariance matrices and their vector of first order moments, which can efficiently be calculated at every step of the state preparation process (even though their photon number statistics cannot be efficiently calculated). Furthermore, the fidelity between two multimode Gaussian states can be efficiently calculated as a function of their covariance matrices and displacement vectors [30]. The fidelity F between two states described by density matrices ρ 1 and ρ 2 is related to the trace distance D by [36]: D is related to the maximum classical l1 distance between different possible measurement outcomes by: where the maximisation is over all sets of detector POVMs {E m } at the output of the network, and p m = tr(ρ 1 E m ) and q m = tr(ρ 2 E m ). If we consider the POVMs projecting onto photon numbers, we then have that: where P 1 and P 2 correspond to the photon number statistics associated with states ρ 1 and ρ 2 . We note that this bound is not a tight bound, and we envisage that future improvements on the theory may yield tighter bounds.

Molecular Simulations
In this section, we explain in more detail the analogy between the calculation of vibronic spectra and boson sampling described by Huh et al. [2], and show how we applied their theory to determine the ideal optical states for estimating the Franck-Condon factors corresponding to the transition under study in tropolone.

Analogy between vibronic spectra and Gaussian Boson Sampling
Upon absorption or emission of light, the change in configuration of a molecule causes a change in its vibrational modes. According to the Franck-Condon principle [20,21], the amplitude of a transition from one vibrational state to another is proportional to the overlap integral between the initial wave function and the final wave function. To calculate these overlap integrals within a harmonic approximation of the potential wells corresponding to the vibrational modes, Duschinsky [19] proposed that the change in mass-weighted normal coordinates associated with the electronic transitions under study can be written as: where U is a real rotation matrix and d is a displacement vector. This is the Duschinsky relation [19]. This canonical transformation maps onto the creation and annihilation operatorsˆ a † ,ˆ a in the following way: and where ω i and ω i are the frequencies of the initial and final electronic state and diag denotes the diagonal matrix [2,22]. Using the Bloch-Messiah decomposition [23], we can write: where U and V are unitary matrices and A D and B D are diagonal matrices with A 2 D = B 2 D + 1. Transformation 8 can then be interpreted as a first multimode rotation described by V , a set of single-mode squeezing operations with squeezing parameters which are the diagonal elements of asinh(A D ), a second rotation described by U , and a displacement described by δ. If the initial physical state is the ground state of the harmonic oscillator, then the transformation can be interpreted as a set of single mode squeezing operations, followed by rotation U and displacement δ.
Since both the vibrational modes of a molecule in this model and modes of the electromagnetic field are described by quantum harmonic oscillators, the operations described above have an analogy in quantum optics. The physical state of the vibrational modes of a molecule after an electronic transition occurring at 0 K can be mapped onto the state of a multimode optical field which starts in vacuum, undergoes single mode squeezing operations described by asinh(A D ), goes through a multiport interferometer implementing transformation U , and is then displaced. Since the Franck-Condon factor corresponding to a specific vibronic transition is proportional to the overlap integral between the ground state and the excited state of the harmonic oscillators, then these Franck-Condon factors can be mapped onto the probabilities of different photon number distributions measured at the output of the interferometer.
Specifically, ifÛ Dok is the operator corresponding to the application of squeezing and rotation as described above, then the Franck-Condon factor (the relative height of a spectral line) corresponding to a specific frequency ω vib is: where |m = |m 1 , m 2 , ... describes the state with m 1 photons in mode 1, m 2 photons in mode 2, etc., and the summation over k is over all the modes. Sampling from these output photon statistics can therefore be used to estimate m|Û Dok |0 and therefore estimate the spectrum of the molecule.

Application to Tropolone
In the 370 nm electronic transition in tropolone, there are two vibrational degrees of freedom that couple only to each other due to selection rules [18]. The Duschinsky rotation and the frequencies of the initial and final electronic states (in cm −1 ) are: According to the theory above, this transition in tropolone can be simulated using two single-mode squeezers and a variable beam splitter. Inserting the Duschinsky rotations and the frequencies for the ground and excited states of the molecules into equations 4-9, we get the following parameters for producing the ideal state: Tropolone Squeezing parameter 1 0.190 Squeezing parameter 2 0.720 Beam splitter transmission 0.10 These parameters allow us to describe the state at any step in this state generation process using its covariance matrix and displacement vector [28,29], and also to determine the corresponding exact photon number statistics (albeit with exponential overhead in the number of modes) [37].

Experimental Details
To produce an optical state that can be used to estimate the vibronic spectrum of tropolone, we built the experiment shown in Fig. 2 in the main text. Here, we provide some additional experimental details.
Our photon source consists of a 8 mm long periodically poled potassium titanyl phosphate (KTP) crystal [25], pumped with a mode-locked titanium-sapphire laser the output of which is stretched, picked off and amplified in a regenerative amplifier, and compressed. The resulting pulses have a center wavelength of about 780 nm, a bandwidth of about 2-4 nm, and occur at a repetition rate of 100 kHz. After the KTP crystal, the pump light is filtered out. The output horizontal and vertical polarization modes have a center wavelength at about 1560 nm. These modes are spatially separated using a Wollaston prism and individually filtered using 3 nm bandpass filters to remove any excess noise at the wrong wavelengths. The two modes are then coupled into polarization-maintaining fiber at the two inputs of a commercial variable fiber beam splitter.
Our detectors are photon-number resolving transitionedge sensors (TES) [26], which have an estimated detection efficiency of 0.98 +0.02 −0.08 [38]. These detectors consist of small 25 µm ×25 µm ×25 nm slabs of superconducting tungsten which are inside an optical cavity optimized for absorption at telecom wavelengths and aligned to an optical fiber [39]. They are kept inside a dilution refrigerator with a base temperature of 40 mK, and maintained at their transition temperature by Joule heating using a bias current via an electro-thermal feedback effect [40]. At this transition temperature, a small increase in temperature causes a large increase in resistance, and this effect is used as the photon detection mechanism. When they absorb one or more photons, their temperature increases, which causes their resistance to increase.
This increase in resistance causes the bias current through the TES to change, and this signal can be read off and amplified using a Superconducting Quantum Interference Device (SQUID) module. After further amplification at room temperature, the resulting analogue signal consists of approximately 5 µs long pulses, which have different characteristics depending on how many photons were absorbed.
These pulses are read off by an analogue-to-digital converter, which compares the incoming pulses with a reference signal and then bins them according to different photon numbers. The processing speed of our computer limits the data taking rate to about 50 kHz; we thus have to discard about half the incoming pulses.

Characterizing the experimental setup
Estimating and optimizing the fidelity of the experimentally generated state to the ideal state calculated above requires characterizing the experimental setup. In this section, we explain how this characterization was performed.
Detector noise: The TESs have a very low intrinsic dark count rate. However, two noise mechanisms are present.
Firstly, the binning procedure that we use to extract photon numbers from the analogue output signal is a bit noisy. Therefore, with the signal and idler fields blocked before the incoupling fibers in our setup, there still is a 0.2% probability of wrongly measuring a one photon event. We approximate this noise as a dark count mechanism, which can be accounted for within our model by considering that our detector is sensitive to both a noiseless input signal and to an additional mode containing a thermal state with a 0.2% single photon component.
Secondly, some pump photons leak through our setup and make it to the detectors. Since TESs resolve the energy of incoming photons, these pump photons are counted as two-photon events. Our TESs have a 0.1% probability of detecting these pump photons. This noise mechanism can be included within our model by considering that our detector is sensitive to an additional pump mode containing a weak coherent state with a 0.1% single photon component.
This noise can be included in our estimate of the fidelities. Both additional noise modes have a fidelity of about 0.998 to vacuum, and since the fidelity for a product state is the product of the fidelities, the total fidelity of these noise modes to vacuum is 0.9958. This fidelity must then be multiplied by the estimated fidelity of the optical state before the detectors in order to determine the total fidelity.
Our TESs also have some degree of inefficiency. This inefficiency is accounted for in our estimate of the total system efficiency described in the following paragraph. Squeezing and loss: We characterize the total loss (including detector inefficiency) and squeezing in our system using a tomography technique similar to that described in [41]. We model our photon source as a perfect TMSV with additional loss in the two modes, use our model for the detectors described above, and numerically find the squeezing parameter and the distribution of the loss in both arms which yield the photon number statistics that most closely match the experimental photon number statistics. We follow this procedure for the beam splitter set first to 100:0, and then to 0:100, in order to determine how the losses in the system are distributed in both modes both before and after the beam splitter. We note that since balanced losses can mathematically be commuted through the beam splitter, we need to consider losses in only one of the output ports of the beam splitter in our model.
To verify that our estimate of the losses is reliable, we perform this tomographic procedure for one of the two beam splitter settings for several different pump powers ranging from 10 µW to 300 µW. We find that at powers exceeding 100 µW our results are skewed as higher order nonlinearities start affecting the pump and the downconverted modes. We therefore use the values for the loss that are found in the low power region in Fig. 4, and estimate the error on these values to be ±2%. These values are also consistent with the heralding efficiencies estimated from the photon number statistics in this plateau region.
We also use the results of our tomography procedure at several pump powers to set the squeezing parameter in experiment and to estimate the error on the squeezing. The squeezing parameter r is related to the pump power P by the following relation: To set the desired squeezing parameter in our experiment, we use the squeezing parameter numerically determined at a given pump power and then use equation 14 to set the correct power. To estimate the error on our estimate of the squeezing parameter, we fit the squeezing parameters determined from our optimization procedure to a curve given by equation 14 in the plateau region. This fit is shown in Fig. 4. We estimate the error on our estimate of the squeezing parameter at low powers to be 0.01. Fig 4 also shows the deviation, quantified by l1 distance, between the measured photon number statistics and those of the closest lossy TMSV measured with noisy detectors determined by our tomography technique. At low powers, this error is less than 10 −3 , so we consider that our description of the optical state and of our detectors is satisfactory. FIG. 4. Characterization data for our experiment. A) Total losses for the two modes in our setup, estimated from our tomography procedure for different pump powers, with the beam splitter set to full transmission. We use the average and the standard deviation of the values in the shaded area for our estimate of the experimental parameters. B) Squeezing parameter, estimated from our tomography procedure, as a function of pump power. Within the shaded area, we find that the relation between squeezing parameter and pump power is very close to the theoretical relation P ∝ sinh(r) 2 . C) Deviation between our experimental photon number statistics and the theoretical photon number statistics given by our model for the optical state, quantified by the l1 distance. D) Coincidence counts on our detectors measured as a function of the angle of the half wave plate (HWP) in our setup. We observe a clear Hong-Ou-Mandel dip.
Distinguishability: We characterized the distinguishability δ between the optical modes by using the depth of the Hong-Ou-Mandel interference (shown in Fig.  4), measured by setting the beam splitter to 50:50 and measuring the number of coincidences as we rotated the HWP in our experimental setup. We used superconducting nanowires as our photon detectors for this procedure. Considering that non-overlapping parts of the two modes contribute to noise photons, the ratio of noise photons to signal photons in the system is then δ. We chose to model this noise as virtual beam splitters of reflectivity δ, placed just after the squeezing operation for both modes, between each mode and a virtual thermal state containing the same average number of photons as the TMSV. The ±2% error on our estimate of δ comes from the error on the estimate of the depth of the measured Hom-Ou-Mandel dip. We note that our model for the noise given by the distinguishability is only an approximation of the full description of this noise, which would require taking several additional modes into account. However, we find that our model is sufficient to provide a rough estimate of the contribution of this noise towards the degradation of the fidelity.
Beam splitter reflectivity: The beam splitter reflectivity was set by blocking one mode in our experiment, setting the beam splitter to be fully transmissive so that the maximum photon number at the detector for that mode could be determined for a given pump power, and then adjusting the beam splitter until the average photon number at the detector was the desired fraction of the maximum. The error on our estimate of the reflectivity was estimated at ±1%.
Fidelity estimation: The squeezing, distinguishability, loss, and beam splitter reflectivity can all be included in our model using the covariance matrix and displacement vector of the state. Characterizing all of these parameters allows us in principle to determine the state that is experimentally generated and to calculate its fidelity to the target state using the relation in [30].
Since there is some error on our estimates of these experimental parameters, we use a bootstrapping method to determine the most likely value of the fidelity that we achieve as well as the error on this value. We simulate 100 states produced by our experiment, where we randomly select the experimental parameters from a Gaussian probability distribution determined by the estimated mean and standard deviations given by our analysis. The mean and standard deviation in the fidelity of these samples are used as our fidelity estimate and as the error on this estimate.

Optimization
After characterizing and modeling the setup, we performed numerical optimization over the squeezing and the beam splitter reflectivity to maximize the fidelity of the experimentally generated states to the target states. Our optimization results are shown in the main text.
In the two mode case, the parameter space is small and finding the optimal parameters is quite straightforward. However, this optimization problem can be difficult for a larger molecule, and finding the global maximum for the fidelity may be a hard problem. However, for many applications one can imagine that finding a good local maximum for the fidelity can yield an estimate of Franck-Condon factors that is good enough for the chosen purpose.

Benzene
To demonstrate the generality of our approach, we also simulated the S 0 ( 1 A 1g ) → S 1 ( 1 B 2u ) transition in benzene [34]. We once again only considered two modes involved in this transition that do not couple to any of the other modes, and for which there is no displacement term.
For this transition, the Duschinsky rotations and the frequencies of the initial and final electronic states (in cm −1 ) are: Benzene: To experimentally simulate benzene, we followed the procedure described above and found the following values of the fidelity and of the optimal parameters, given the experimental parameters and the model shown in figure  2 in the main text. With our estimate of the fidelity for benzene and our estimate for the deviation of our state from a TMSV in photon number statistics (see Fig. 4C), we can bound the error on the l1 distance of our Franck-Condon factors to the actual Franck-Condon factors to 0.1663. By comparing our experimental results to theory, we find an actual error of 0.0198, which is well within the bound.

Benzene, Imperfect Experiment
We find that the loss in our setup and the use of a TMSV instead of two SMSVs have similar effects on our results for benzene as they have for tropolone, discussed in the main text. Due to the loss, we overestimate the Frank-Condon factor corresponding to the vibrational ground state, and find non-zero Franck-Condon factors for transitions which should not be allowed. Our use of a TMSV also yields almost symmetric photon number statistics in the two modes.

Classical state with the highest fidelity
In this section, we show how to find the unique classical optical state which has the highest fidelity to the ideal target state.
We consider an ideal experiment, in which vibronic spectra are estimated using squeezed and displaced optical states sent into an interferometer (we recall that the displacement operation can be performed at any point in the experiment). The fidelity being invariant under unitary transformations, finding the closest classical state to the ideal target state is equivalent to finding the closest classical state to the initial displaced single-mode squeezed vacuums sent into the interferometer. Therefore, since the closest classical state to a single mode displaced Gaussian state is a coherent state with the same displacement [42], the closest classical state to the final target state is a coherent state with the same displacement. Since this state is Gaussian, its fidelity to the target state can be efficiently determined as explained above, and we can then use this fidelity as a classicality witness.
For the two molecules that we simulated, the closest classical state in terms of fidelity is therefore vacuum. We note that since vacuum is a Gaussian state, then in the specific case where the closest classical state is vacuum the difference in photon number statistics between the target optical state and the closest classical state can be efficiently calculated, as opposed to simply bounded using the fidelity.
We can compare our results for benzene to those yielded by the classical approximation procedure discussed in the main text. Vacuum has a fidelity of 0.9893 to the ideal state. We note that this is higher than the fidelity of 0.9861 for our experimental state; the experimental imperfections therefore cause our experiment to perform less well than an efficient classical estimation technique. The error bound using the fidelity for the classical state and equation 3 from the main text is 0.1459, which is smaller than the 0.1663 bound for our experimental state. The actual error in the corresponding estimation of the Franck-Condon factors for benzene is 0.0213 for the classical simulation, which is close to the 0.0198 error that we achieved experimentally.

Sulfur dioxide
The transitions under consideration in tropolone and benzene can be simulated in quantum optics without displacement. To study the case of displacement, and compare our work to that done using a trapped ion simulator [13], we now theoretically study the case of the SO − 2 → SO 2 transition that involves both squeezing and displacement. The Duschinsky transformation between the two normal modes of this molecule is as follows [35]: The ground and excited state frequencies are (in cm −1 ): The ideal optical state to simulate this transition, with the displacement occurring after the squeezing and beam splitter operations, is produced with the following parameters: These low squeezing parameters compared to the displacement cause the Franck-Condon factors to be dominated by the displacement term. This fact is reflected in the fidelity of the closest classical state to this ideal optical state, which is 0.9965. This classical state is the coherent state with the same displacement as the ideal state.
In an experimental setup, displacements can be induced using classical laser light and can be adjusted to account for the loss. An experimental setup that implements these displacements and otherwise contains the same imperfections as our setup (in terms of loss, detector noise, imperfect overlaps between the modes and the use of a TSMV instead of an SMSV) would require the following parameters to maximize the fidelity to the target state: The fidelity of the resulting experimental state to the target state would be 0.993, leading to an error bound in the estimate of the vibronic spectrum following equation 3 in the main text of 0.118. We note that this fidelity is smaller than that of the closest classical state; since the required squeezing is quite low, noisy detectors are enough to degrade the fidelity of the experimental state below that of the closest classical state. We can calculate the expected vibronic spectrum yielded from performing this experiment and compare it to that produced by the ideal optical state. This comparison is shown in figure 5. There is good qualitative agreement between the expected experimental results and theory, with a total l1 distance between the corresponding probability distributions of 0.0394. This simulation of the SO − 2 → SO 2 transition highlights the different roles played by displacement and by squeezing. This transition is dominated by the displacement term. A quantum optical experiment can in principle accurately implement displacements even in the presence of loss, which is why there is good qualitative agreement between the predicted results and the ideal case. Squeezing plays only a marginal role in this transition, so that outperforming the closest classical state requires an experimental setup with a very small amount of imperfections.

Nonclassicality Analysis
Here, we analyze in detail how the classicality of the state that we sample from, determined by the regularity of its P function, is affected by loss and mode distinguishability within our model. We assume balanced loss in our setup and no detector noise for this analysis.
In order for the P function to be regular, the covariance matrix V out of the output state must satisfy [43,44]: where I 4 is the 4 × 4 identity matrix. Since all the optical tranformations in our model are Gaussian, V out can be readily obtained from the symplectic transformations [29] corresponding to squeezing with squeezing parameter r and the virtual beam splitters with transmittivities t A = t B = t AB , t D and 1 − δ which describe the photon losses in both modes, detector inefficiency and the noise photons, respectively. Using some results for positivedefinite matrices, one may show that condition (16) is equivalent to the following expression: t AB t D sinh(r) 2 − sinh(r) cosh(r)(1 − δ) ≥ 0. (17) Inequality (17) is a sufficient condition for the P function of the state to be regular. Inequality (17) indicates that the main source of classicality in our experiment is the distinguishability of the input optical modes. The expression in the parentheses gives a lower value for the distinguishability δ that renders the output state classical, i.e. 0.76 δ for r ≈ 0.25.
Interestingly, this lower value decreases when the squeezing r increases, which means that the output state becomes classical in our model at squeezing r ≥ 1.1 (and δ ≈ 0.2). This is due to the noise photons produced by our model for the distinguishability. Larger values of r yield higher numbers of these noise photons, which washes out the non-classicality caused by the squeezing. Moreover, we see that photon losses and detector inefficiencies contribute equally to the classicality of our experiment, when t AB t D goes to zero.

Computational Complexity
Our protocol for estimating the vibronic spectra of molecules using quantum optics is closely related to boson sampling [8] in general, and to Gaussian boson sampling in particular [45][46][47]. In this section, we outline this relationship and discuss the computational complexity of our protocol.
Aaronson et al. originally demonstrated that sampling from the output distribution of N identical single pho-tons sent into N input ports of a M × M lossless multiport interferometer implementing a Haar-random unitary transformation is probably a hard computational problem even in the approximate case, if M N . This hardness proof was extended to the problem of scattershot boson sampling [45,46], in which every input port of the interferometer is connected to one mode of a two-mode squeezed vacuum, with the other mode connected to a photon detector, in such a way that N photons on average are sent into random inputs of the interferometer. Scattershot boson sampling is a special case of Gaussian boson sampling, which is the general problem of sampling from the output of an interferometer with single mode squeezed state inputs (since single mode squeezed states can be converted into two mode squeezed vacuums on a beam splitter). Since scattershot boson sampling is computationally hard, then the general problem of Gaussian boson sampling must also be hard.
Our experimental protocol is a special instance of Gaussian boson sampling. However, that the general problem of Gaussian boson sampling is hard does not imply that a specific instance of the problem is hard if that specific instance does not map onto the framework described in the previous paragraph. This can clearly be seen, for example, if the input Gaussian states are thermal states, in which case efficient sampling algorithms exist [33]. Our protocol for the estimation of vibronic spectra differs from the framework in which a hardness proof exists in two respects. Firstly, the multiport interferometer does not implement a random unitary transformation but a specific transformation determined by the molecule under study. Secondly, the number of photons that is generated is not necessarily smaller than the number of modes. The computational complexity of our protocol is therefore an open question.
It is worth mentioning, however, than if an experimentally generated state beats the classicality criterion outlined in this work, then currently known simulation algorithms based on the phase space description of the state are generally not applicable [33,48].