Classical Algorithm for Generating Multi-Mode Bosonic Transition Spectra by Phase Space Sampling

Abstract

This disclosure relates to a method and system for efficiently and analytically generating spectra associated with transitions between thermal states in a multi-mode bosonic system using a classical algorithm to compute Fourier components of the transition spectra of the multi-mode bosonic system based on a representation space sampling method inspired by quantum Gaussian boson sampling followed by an inverse Fourier transform. The disclosed method and system particularly apply to efficient estimation of molecular vibronic spectra. For example, an exact solution of the Fourier components of molecular vibronic spectra at zero temperature may be analytically computed in a representation space by using, for example, a positive P-representation of the multi-mode bosonic vibronic quantum states of the molecule. Such a method may be further applied to more general vibronic spectroscopy, such as computing molecular vibronic spectra at finite temperatures by introducing additional auxiliary bosonic modes and computing vibronic spectra associated with non-thermal states, such as Fock states.

Description

FIELD OF THE INVENTION

This disclosure relates to a method and system for efficiently generating spectra associated with transitions between thermal states in a multi-mode bosonic system using a classical algorithm based on a phase space sampling method inspired by a quantum Gaussian boson sampling methodology. The disclosed method and system particularly apply to efficient estimation of molecular vibronic spectra.

BACKGROUND OF THE INVENTION

Vibronic spectra represent an essential and characteristic property of a molecule composition and structure. As such, molecular vibronic spectra may be relied upon to analyze chemical composition of molecules and to study molecular structures. Molecular vibronic spectra of a plurality of reference molecular structures may be computationally generated and used as a library of molecular fingerprints for experimentally identifying and studying molecular compounds. However, the computation of vibronic spectra of a particular molecular structure has been known to be a challenging task, with its complexity rapidly growing with the size and/or the number of vibrational modes of the molecular structure. For example, the best-known classical algorithm scales combinatorially with the size of the molecular system. While a quantum simulator, particularly the ones based on Gaussian boson sampling of multi-mode photons, may be implemented to efficiently generate the characteristic molecular vibronic spectra, such a simulator does require a quantum computing machine, such as a linear optical circuits constructed to process delicate non-classical photonic states generated by non-classical photon sources.

SUMMARY OF THE INVENTION

This disclosure relates to a method and system for efficiently and analytically generating spectra associated with transitions between thermal states in a multi-mode bosonic system using a classical algorithm to compute Fourier components of the transition spectra of the multi-mode bosonic system based on a representation space sampling method inspired by quantum Gaussian boson sampling followed by an inverse Fourier transform. The disclosed method and system particularly apply to efficient estimation of molecular vibronic spectra. For example, an exact solution of the Fourier components of molecular vibronic spectra at zero temperature may be analytically computed in a phase space by using, for example, a positive P-representation of the multi-mode bosonic vibronic quantum states of the molecule. Such a method may be further applied to more general vibronic spectroscopy, such as computing molecular vibronic spectra at finite temperatures by introducing additional auxiliary bosonic modes and computing vibronic spectra associated with non-thermal states.

In some example implementations, a method for estimating a transition spectra in a bosonic system is disclosed. The bosonic system may be associated with M bosonic modes, M being an integer equal to or larger than 2. The method may include determining a unitary operation modified from a transition operator associated with the bosonic system; processing each of N sampling bosonic states of the bosonic system in a representation space of the bosonic system to generate N sets of samples using at least the unitary operation, wherein each set of samples correspond to one of the N sampling bosonic states and are sampled from at least one variable in the representation space, N being a positive integer; analytically generating Fourier components the transition spectra based on the N sets of samples; and inverse-transforming the Fourier components to estimate the transition spectra.

In any one of the example implementations above, the unitary operation may include a complex mixed position and momentum operator.

In any one of the example implementations above, each of the N sampling bosonic states may be a non-Gaussian state of the M bosonic modes.

BRIEF DESCRIPTION OF THE DRAWINGS

For a more complete understanding of the invention, reference is made to the following description and accompanying drawings, in which:

FIG. 1 illustrates molecular vibronic spectra generation process based on a quantum Gaussian algorithm and a quantum-inspired classical algorithm;

FIG. 2 shows an example data and logic flow for the quantum-inspired classical algorithm of FIG. 1 for generating molecular vibronic spectra;

FIG. 3 shows comparison of molecular vibronic spectra for CH₂O₂, 1¹A′→1² A′ as generated based on direct computation (ideal bars) and based on the quantum-inspired classical and analytical algorithm of FIG. 1 (“Analytical Solution” points and curve);

FIG. 4 illustrates transformation of molecular vibronic spectra at finite temperature (upper) and treating the problem by introducing auxiliary modes and two-mode squeezing operation and using the zero-temperature algorithm (lower);

FIG. 5 illustrates transformation of molecular vibronic spectra with a Fock state as an input;

FIG. 6 illustrates comparison of molecular vibronic spectra for CH₂O₂, 1¹ A′→1² A′ as generated based on direct computation (ideal bars) and based on Gaussian boson sampling simulation (“GBS” points and curve);

FIG. 7 illustrates comparison of total variance distance (TVD) between ideal direct computation (upper curve with circular points) and the quantum-inspired classical algorithm (lower curve with triangular points) of FIG. 3 (the insert is in logarithm scale);

FIG. 8 illustrates virtual ideal vibronic spectra and empirical vibronic spectra generated using the classical algorithm with N=10⁵ samples and with a total variation distance is of 0.00222;

FIG. 9 illustrates virtual ideal vibronic spectra and empirical vibronic spectra generated by Gaussian boson sampling experimental raw data with N=10⁶ samples and with a total variation distance of 0.01215; and

FIG. 10 illustrates Virtual FCP with a randomly generated weights for all 144 vibronic modes of a molecular system. The unit of wave number is cat. The generation of each instance using classical algorithm costs about 90 secs using a single core computing system.

DETAILED DESCRIPTION

Bosonic Systems and Molecular Vibronic Spectra

Atomic/nuclear components of a compound molecular structure interact with one another and thus may vibrate in multiple number (e.g., integer M) of vibrational modes. Each vibrational mode, approximately, may be modeled as a harmonic oscillator of a particular frequency. The multiple vibrational modes thus may correspond to a multitude of characteristic vibrational frequencies cv, within the molecular structure, where i∈[1, 2, . . . , M]. Quantum mechanically, vibrations in each of the modes may be described by quantum states |m_i (where m_i∈[0, 1,]) with ladder-like energy levels each associated with a quantum number m_i and an energy determined by ω_i(m_i−½) (for m_i=0, 1, 2, . . . , ∞), where the vacuum level is represented by ω_i/2. An arbitrary molecular vibrational state thus may be described by a quantum superposition of these multi-mode base states. The molecular vibrations may thus be considered a bosonic system, similar to a system of multi-mode photons. For convenience, a molecular vibrational quanta may be referred to as a photon in this disclosure.

The various vibrational modes are essential characteristics of a molecular structure. Molecular vibrations may be further coupled to the electronic states or transitions in the molecule. For example, electronic transitions may change or perturb the nuclear structure of the molecule, thereby introducing new sets of vibrational modes. The vibrational modes of a molecular may be investigated spectroscopically, for example, as fine spectral lines or structures within an electronic transition spectra. Such spectra may be referred to as molecular vibronic spectra. Vibronic spectra thus represent a fundamental property of a molecule. As such, molecular vibronic spectra may be relied upon to analyze chemical compositions and study molecular structures. In many circumstances, computation of such spectra may be critical in order to derive a molecular structure based on measured vibronic spectra. For example, molecular vibronic spectra of a plurality of reference molecular structures may be computationally generated and used as a library of molecular fingerprints for experimentally identifying and studying molecular compounds.

However, traditional classical computation of vibronic spectra of a particular molecular structure has been known to be a challenging task, with its complexity rapidly growing with the size and/or the number of the vibrational modes of the molecular structure. For example, as shown below, the best-known classical algorithm scales combinatorially with the size of the molecular system. While a quantum simulator, particularly the ones based on Gaussian boson sampling of multi-mode photons, may be implemented to efficiently generate the characteristic molecular vibronic spectra, such a simulator does require a quantum computing machine, such as a linear optical circuits constructed to process delicate non-classical photonic states generated from a non-classical multi-mode photon source, which is usually challenging from experimental standpoint.

In this disclosure, a method and system for efficiently generating molecular vibrational spectra associated with transitions between vibrational thermal states in a molecular structure are described. The disclosed method and system are based on a quantum inspired classical algorithm. In particular, the disclosed classical algorithm analytically obtain Fourier components of the molecular vibrational spectra at zero temperature by sampling phase-space variables in a positive P-representation of a thermal multi-mode molecular vibrational state. Such a sampling technique is inspired by a quantum approach involving Gaussian boson sampling. While this algorithm relies on a quantum-inspired simulation, it nevertheless provides an exact and efficient solution of the Fourier components of molecular vibronic spectra in a classical manner without needing a quantum simulator. The molecular vibrational spectra may then be generated by performing inverse Fourier transform of the Fourier components. Such a method is further extended to generate molecular vibronic spectra at finite temperatures by introducing additional auxiliary modes and two-mode squeezing operations.

In more general cases, a Fock state rather than a thermal state as an initial state may be considered. The Fock state may correspond to single vibronic levels and may be written as a loop halfnian of a matrix. A Monte Carlo method using the positive P-representation may be used to estimate the molecular vibronic spectra in that case. However, the Monte Carlo method becomes inefficient when the number of quanta in the initial Fock states grows. In some other implementations of a similar method setup in single-photon boson rather than multi-photon sampling to obtain corresponding virtual molecular vibronic spectra, the corresponding Fourier components may become hard to compute (#P-hard), in stark contrast to Gaussian boson sampling which has an exact solution, rendering a computational advantage for Gaussian boson sampling over single-photon boson sampling. Nonetheless, a practically relevant question is whether the Fourier components can be additively approximated to within the same 1/poly(N) accuracy that can be obtained using N number of samples for a boson sampling experiment. Interestingly, when squeezing parameters are zero, the quantity reduces to a permanent, which can be approximated by Gurvits's algorithm within an additive error. As such, an approximation method may be implemented to tackle the problem using the positive P-representation with post-selection and another method that approximates a loop hafnian. While there may be a limitation of the method of positive P-representation with post-selection as the number of quanta in the initial Fock state grows, a regime may be identified where the approximation error can be efficiently controlled even for large system size.

While the example implementations below are provided in the context of molecular vibronic spectral generation, the underlying principles and procedure may be applicable to generating transition spectra in other multi-mode bosonic systems.

Example Implementations of Quantum-Inspired Classical Algorithm for Estimating Molecular Vibronic Spectra

An example procedure for analytically generating molecular vibronic spectra using an efficient classical algorithm is shown in FIG. 1 in comparison to a typical quantum simulator for the same purposes. In particular, the analytical and classical algorithm is shown by the lower data-flow path 110 whereas the quantum simulation procedure is illustrated by the upper data-flow path 120. For both approaches, the computational output includes the vibronic spectra 130 for the zero temperature molecule (e.g., absorption vibronic spectra), although the example method/system disclosed herein can be expanded in non-zero temperatures, as described in further detail below.

For the quantum simulation data path 120, an example Gaussian boson sampling procedure may be performed. For example, a product quantum state of single-mode Gaussian quantum states of multi-mode (e.g., M-mode) photons may be initially determined, prepared, and generated from appropriately configured and controlled photon sources. Such a product quantum state may represents a dressed and displaced initial state (e.g., vacuum state for zero temperature and generated as an input dressed and displaced M-mode photon state, as shown by 122, to a quantum circuit 124. The quantum circuit 124 may include a plurality of beam splitters and other linear optical components to simulate a unitary quantum operation. A quantum unitary operator associated with the quantum circuit 124 may be determined in accordance with example implementations described in further detail below.

For example, the unitary operator corresponding to the unitary quantum operation of the quantum circuit 124, in combination with the dressing operator and the displacement operator for generating the dressed and displaced initial quantum state of the input product quantum state of the M-mode photon state, may correspond to a transition operator of the molecular system between the ground state and vibronic thermal states. The linear optical components of the quantum processing circuit may thus be configured to preform the operation corresponding to the unitary operator. The output of the linear operation may then be projected to the various photon modes and detected as photon numbers in the M photon modes, {m}. The array of linear optical components and the detection in the M-photon modes are thus equivalent to a Gaussian Boson sampler, as shown by 124 and 126.

A set of N samples, represented by S={m⁽ⁱ⁾}_i=1^N and as shown by 128 in FIG. 1, may be selected. The samples may then be post processed by first being sorted according to total energy of each sample, where the total energy of each sample corresponds to the energy of all vibronic quanta (e.g., m⁽ⁱ⁾·ω′, where ω′ represents the M-mode frequencies). In some practical implementations, for a particular spectra frequency, ω_vib, a window Δω_vib, may be introduced for the post processing procedure. In other words, when a samples total energy-equivalent frequency m⁽ⁱ⁾·ω′ lies between [ω_vib−Δω_vib/2, ω_vib+Δω_vib/2] for a certain ω_vib, the count for a bin corresponding to ω_vib, is incremented. The distribution of final counts in all ω_vib bins forms the simulated molecular vibronic spectra FCP(ω_vib), as long as sufficient N number of samples are collected from the quantum linear optical circuit and are post processed.

In some example implementations, a classical algorithm may be implemented as shown by the procedure 110 in FIG. 1. While such an algorithm as described in further detail below is inspired by the quantum Gaussian boson sampling procedure described above (procedure 120 of FIG. 1), it nevertheless would not need to rely on an actually quantum simulation machine. In order to generate molecular vibronic spectra due to transition to vibronic states at zero temperature, a particular quantum state of the multi-mode molecular vibronic state may be determined as an input state to the classical algorithm. Such an input vibronic state may be a general M-mode vibronic state determined and generated for an initial state of the molecule (e.g., a vacuum state at zero temperature, a Gaussian state at a finite initial temperature, or Fock state, and the like). Such an input vibronic state corresponding to an initial vacuum state at zero temperature, for example, may be a Gaussian state associated with the M vibronic modes. In some example implementations, the input M-mode vibronic quantum state of the molecule to the disclosed classical algorithm may be represented in a phase space as a quasi-probability distribution of a set of phase space variables associated with the phase space. Such a phase-space representation of a quantum state may be referred to as a P-representation.

An important property of a generalized P-representation may be that the underlying phase-space distribution may always be chosen non-negatively such that the expectation value of a normal-ordered operator can be readily estimated by sampling from the distribution. Such non-negative phase-space distribution may be referred to as positive P-representation.

Unlike traditional classical algorithm for the calculation of molecular vibronic spectra, the quantum-inspired classical algorithm can provide fast and efficient estimation of the molecular vibronic spectra. The gain of such a classical algorithm over traditional classical calculation with complexity that grows combinatorially with respect to the size of the problem or the number of the vibronic modes of the molecular system mostly relies on the principle that information being sought or of interest (e.g., the molecular vibronic spectra) represents a coarse information that need not require detailed knowledge whose calculation cannot be avoided when the traditional classical calculation method is being used. The example quantum-inspired algorithm described below preforms sampling and analytics that more resemble the quantum simulation in that it is not capable of revealing or calculating these knowledge items pertinent to the molecular vibronic system, and fortunately, those knowledge items are not necessary and thus need not be calculated for the estimation of the molecular vibronic spectra of interest.

As shown in 110 of FIG. 1, this classical algorithm starts by preparing an input state of the M-mode molecular vibronic system based on an initial vibronic state of interest (e.g., a vacuum state at zero temperature, a Gaussian state at a finite temperature, a Fock state, and the like). For an initial vacuum state, the input state may correspond to an Gaussian state, as described in further detail below. Merely as an example, the input state may be expressed in the P-representation in a phase space. Again, the P-representation may be characterized by a quasi-probability distribution function of a bosonic state over the phase space that can always be chosen non-negative. For example the phase space variables for the quasi-probability distribution function of the P-representation may include a set of M-dimensional vectors. For a particular example, the phase space variables may include a set of two M-dimensional vectors. In FIG. 1, such example quasi-probability distribution function of the example input Gaussian state is shown by 112 and represented by P(α, β), where α and β represent the set of two M-dimensional vectors.

As further shown in 110 of FIG. 1, the input state of the M-mode molecular vibronic state as represented by P(α, β) may then be transformed to generate a transformed P-representation P(α′,β′), as shown by 114, by applying an operator U, shown by 113, to the set of two phase space variables α and β to transform them into α′ and β′. The phase space of α and β may be further sampled. Specifically, 2M complex numbers of (α, β) may be sampled. The 2M samples may be transformed by the operator U, as expressed by (α, β)→(α′, β′)=(Uα, Uβ). The operator U is implemented as a matrix and represents the transmission/transition matrix for the M-mode molecular system in the phase space. The operator U, for example, may comprise a mix of position and momentum operators.

As further shown by 115 of FIG. 1, Fourier components of the molecular vibronic spectra may be estimated classically from the phase space sampling above. Such estimate may be exact or may be approximated to a certain predetermined error level. For example, the Fourier components of the molecular vibronic spectra may by represented by {tilde over (G)}(k), where k represents the Fourier space of the molecular vibronic spectra. The classical process for estimating {tilde over (G)}(k) from the transformed sampled phase space variables (α′, β′) is provided in further detail below for various initial and input states of the M vibronic modes of the molecule. Such estimate may be based on a correspondence between the Fourier components of the molecular vibronic spectra and an effective photon number operator {circumflex over (n)}, as shown by 115 of FIG. 1. The estimated Fourier components {tilde over (G)}(k) are exemplarily illustrated in 116 of FIG. 1. Thereafter, the estimated Fourier components {tilde over (G)}(k) may then be used to generate the molecular vibronic spectra 130 via a classical inverse Fourier transform process, as indicated by 118 in FIG. 1.

Although the example of FIG. 1 illustrates sampling in the P representation phase space, the underlying principles apply to implementations where the sampling is performed in other spaces. Other representations or spaces may be used under different circumstances, depending on, for example, the initial state of the molecular system. For example, the sampling space may be generally a mix of position and momentum.

The example procedure above for estimating the molecular vibronic spectra in the quantum-inspired classical manner is further described in the data and logic flow 200 of FIG. 2. As shown in FIG. 2, the data and logic flow may include step 210 for determining a unitary operation configured to project an input superposition state of the bosonic system into number states corresponding to the M bosonic modes. Such an input state may be determined form an initial state of moledule system. For example, such an input state may be a Gaussian state of the M bosonic modes corresponding to an vacuum initial state at zero temperature, for calculating zero-temperature vibronic spectra of the molecule.

In step 220, each of N sampling bosonic states of the bosonic system may be processed in a representation space of the bosonic system e.g., a phase space to generate N sets of samples based on at least the unitary operation, wherein each set of samples correspond to one of the N sampling bosonic states and are sampled from at least one variable in the representation space, N being a positive integer.

In step 230, Fourier components the transition spectra may be analytically generated based on the N sets of samples. Such generation may be exact or may be approximate to a predetermined error level corresponding to the size (number of vibronic modes) of the molecule.

In step 240, inverse-transformation may be performed on the Fourier components to estimate the transition spectra. The correspondence between the Fourier components and the final transition spectra may be exact or may be approximate to a predetermined error level corresponding to the size (number of vibronic modes) of the molecule for different initial and input sates, as described in further detail below.

The example implementations above may be applied to generate vibronic spectral for a molecular system at zero temperature. In particular, the initial vibronic state of the molecular system would be the vacuum state of all M modes. The corresponding input state to the sampling process above may be an Gaussian state of the M vibronic modes, and the correspondence between the Fourier components and the final transition spectra may be exact, as described in further detail below. The representation space in such a case may be the P-representation space.

In some other specific example implementations, vibronic spectra may be estimated for a molecular system at a finite non-zero temperature following the zero-temperature implementation. For example, as described in further detail below, M auxiliary modes may be added in comparison to the zero-temperature case. Two-mode squeezed states may be prepared to simulate thermal states and detect the additional auxiliary modes as well, as depicted in FIG. 4 and described in further detail below. As such, the procedure described above and in the data and logic flow of FIG. 2 for estimating zero-temperature vibronic spectra of the molecular system applies to the finite temperature situation, except that, as a cost, the additional M auxiliary modes are included in the preparation of the input Gaussian state and in the sampling in the phase space.

In yet some other example implementations, the sampling method above may be applied to an estimation of molecular vibronic spectra for non-Gaussian initial state, such as a Fock state. In order to employ a mapping from molecular vibronic spectra to Gaussian boson sampling, additional M modes and two-mode squeezed states of some squeezing parameters may be introduced. In comparison to the finite temperature Gaussian initial state case, a post selection may be further performed for the additional modes. As described in further detail below, the estimation of the Fourier components may not be exact in such a situation. The error may exponentially increase with the photon number in the initial Fock state, suggesting an exponential increase in samples needed according to the sampling process of FIGS. 1 and 2. However, as described in further detail below, the above sampling process may be combined with, for example, the Gurvits' algorithm generalized to include multiphotons to approximate the Fourier components of the molecular spectra with reduced costs in such situation within a reasonable error. The sampling may be performed in the particularly determined representation space as described in further detail below.

The present disclosure implies that Gaussian boson sampling itself is a hard task, but when it comes to the coarse-grained version of sampling, an efficient classical counterpart may exist. Therefore, on the contrary to the common understanding, the disclosure herein suggests that molecular vibronic spectra may not be the candidate for which only the power of a quantum sampler can boost the computational performance beyond classical means. Rather, classical algorithm as described above may diminish or reduce the quantum advantage.

On the other hand and as shown in further detail below, beyond Gaussian framework where a direct integral does not work, the disclosure herein shows that a Monte Carlo method may be inefficient when the number of quanta of the input state is large. When a similar problem in single-photon boson sampling is considered, computing the Fourier components is a #P-hard problem, which is in stark contrast to Gaussian boson sampling case. Nevertheless, a quantum-inspired classical efficient algorithm has been identified herein to approximate Fourier components of molecular vibronic spectra, which may be further applied to other application scenarios.

Non-Sampling Classical Molecular Vibronic Spectra Computation

The transition probability between an initial vibronic mode and a certain final vibronic mode is referred to as the Franck-Condon (FC) factor. FC profiles (FCP) may be obtained by computing many FC factors corresponding to a given vibrational transition frequency (ω_vib). For example, the FC factor at 0 K may be obtained with the initial vacuum state as

$\begin{array}{cc}\mathrm{FCP}\left({\omega }_{\mathrm{vib}}\right)=\sum _{m=0}^{\infty }{❘\langle m❘{\hat{U}}_{\mathrm{Dok}}❘0\rangle ❘}^2\delta \left({\omega }_{\mathrm{vib}}-{\omega }^{\prime }·m\right),& \left(1\right)\end{array}$

where the Doktorov transformation Û_Dok is given by

Û _Dok ={circumflex over (D)}(δ/√{square root over (2)})S_Ω^\,{circumflex over (R)}_UŜ_Ω,  (2)

where δ(·) represents the delta function and m=(m₁, . . . , m_M) is the final vibronic modes excitation vector with M vibronic modes. Here, Ŝ, {circumflex over (D)}, and {circumflex over (R)}_u represent squeezing, displacement, and rotation operators, respectively, and Ω and Ω′ are given by

Ω≡diag(√{square root over (ω₁)}, . . . ,√{square root over (ω_M)}),Ω′≡diag(√{square root over (ω′₁)}, . . . ,√{square root over (ω′_M)}),  (3)

where ω_i's and ω_i's account for initial and final harmonic angular frequencies, respectively.

The FCP above represents the molecular vibronic spectra. Eq. (1) shows that the direct computation of all FC factors are inevitably computational hard. In particular, not only is the number of combinations of m satisfying the delta function combinatorially large, the FC factor |m|Û_Dok|0|² corresponds to the hafnian of a matrix, the computation of which costs exponentially in the size of matrix or a number of vibronic modes. For a fixed m=Σ_i=1^Mm_i, the number of FC factors may be given by

$\begin{array}{cc}\left(\begin{array}{c}M+m-1\\ m\end{array}\right).& \left(4\right)\end{array}$

Gaussian Boson Sampling for Molecular Vibronic Spectra Generation via Quantum Simulation

Gaussian boson sampling represents one of the sampling problems that are proved hard for classical computers under some plausible assumptions. To implement Gaussian boson sampling, a product of M single-mode Gaussian states, |ψ_in, may first be prepared, and the product state may be fed into an M-mode linear-optical circuit Û_LON, composed of beam splitters. The output state |ψ_out=Û_LON|ψ_in may be detected or measured by photon-number detectors.

The output photon patterns in may follow the probability distribution below

ρ(m)=|m|Û_LON|ψ_in|²,  (5)

which may be related to the (loop) hafnian of a certain matrix.

Due to Gaussian boson samplings computational power beyond classical computers, a Gaussian boson sampler may be employed to generate molecular vibronic spectra. For example, the correspondence between the FC factor |m|Û_Dok|0|² in Eq. (1) for zero temperature and the output probability ρ(m) of detection photon number pattern in in Gaussian boson sampling may be established with an appropriate choice of input state and the linear-optical operation.

In some implementations, an appropriate Gaussian boson sampling may always be chosen by exploiting the fact that the Doktorov operation Û_Dok is a Gaussian unitary, and can be decomposed as Û_LON{circumflex over (D)}(α)Ŝ(r){circumflex over (V)}_LON. Since {circumflex over (V)}_LON|0=|0, the decomposed unitary operation can be implemented by preparing displaced squeezed states and applying a linear-optical circuit, which is exactly the procedure of Gaussian boson sampling. In other words, for the Doktorov operation Û_Dok of a molecular structure, a Û_LON can be identified and used to build the linear optical circuits. The linear operators above, may be mix of both position and momentum operators.

Therefore, by implementing Gaussian boson sampling and post-processing the output samples using the relation (1), FCP(ω_vib) may be efficiently generated via simulation by the linear optical quantum circuit, as shown by 120 of FIG. 1 and as described above.

More specifically, sampling may be first taken from ρ(m) using Gaussian boson sampler. A number of samples may be collected, represented by S={m⁽ⁱ⁾}_i=1^N. The samples may then be sorted according to the delta function above with respect to the frequency (or energy). In some implementations, a window around each frequency Δω_vib may be introduced. The window may be provided with a predetermined and/or configurable size Δω_vib.

In other words, when a samples frequency ω′·m lies on [ω_vib Δω_vib/2,ω_vib+Δω_vib/2] for a certain ω_vib, a statistical bin corresponding to ω_vib may be incremented. Once sufficient number of samples are post-processed, an estimate of eFCP(ω_vib) (denoting empirical FCP, referring to the experimental quantum simulation aspect of the methodology) may be generated.

Exact Quantum-Inspired Classical Solution for Molecular Vibronic Spectra at Zero Temperature

Inspired by the mapping from a molecular vibronic spectra problem to a Gaussian boson sampling and the phase-space method in quantum optics, the quantum-inspired solution for classically obtaining the molecular vibronic spectral without replying on a quantum mechanical simulation machine may be identified.

For example, it is noted that the Fourier components of the molecular vibronic spectra may be analytically derived and simplified as:

$\begin{array}{cc}\tilde{G}\left(k\right)\equiv \langle \underset{i=1}{\overset{M}{\otimes }}\sum _{m_i=0}^{\infty }❘m_i\rangle \langle m_i{❘}_i{e}^{-{\mathrm{ikm}}_i{\omega }_i^{\prime }\theta }\rangle & \left(6\right)\end{array}$ $\begin{array}{cc}=\langle e^{-\mathrm{ik}\hat{n}·{\omega }_i^{\prime }\theta }\rangle & \left(7\right)\end{array}$ $\begin{array}{cc}G\left({\omega }_{\mathrm{vib}}\right)=\frac{1}{K+1}\sum _k\tilde{G}\left(k\right)e^{\mathrm{ik}\theta {\omega }_{\mathrm{vib}}}& \left(8\right)\end{array}$ $\begin{array}{cc}=\sum _{m=0}^{\infty }p\left(m\right)\delta \left({\omega }_{\mathrm{vib}}-{\omega }^{\prime }·m\right),& \left(9\right)\end{array}$

where θ≡2π/(K+1), k∈{0, . . . , K}, K is the number of bins in discretized distribution of ω_vib, and |m_im_i|_i, is the projector on m_i, photon state for the ith mode.

The final expression above shows the exact correspondence between G(ω_vib) and FCP (ω_vib). Hence, in general, a Fourier component is an overlap between an output state |ψ_out and a phase-shifted state e^{−ikñ·w′θ}|ψ_out. Once such overlap is computed efficiently (e.g., in a classical manner), the Fourier components may then be obtained and used to recover the spectra by an inverse Fourier transform.

Especially when the relevant quantum states are Gaussian states, which is the case for molecular vibronic spectra at zero temperature to transition to thermal states, the Fourier components can be analytically computed. One example way to compute the Fourier components may involve using the positive P-representation of the Gaussian state in a phase space as described above.

More generally, a generalized P-representation of an M-mode quantum state may represent one of the quasi-probability distributions of a corresponding bosonic state of the molecular vibronic modes:

$\begin{array}{cc}\hat{\rho }={\int }_{{ℂ}^{2M}}P\left(\alpha ,\beta \right)\hat{\Lambda }\left(\alpha ,\beta \right)d\alpha d\beta ,\hat{\Lambda }\left(\alpha ,\beta \right)\equiv \frac{❘\alpha \rangle \langle {\beta }^{*}❘}{\langle {\beta }^{*}|\alpha \rangle },& \left(10\right)\end{array}$

where |α=|α₁ ⊗ . . . ⊗|α_M is an M-mode coherent state.

An important property of the generalized P-representation is that the distribution P(α, β) representing the bosonic state can always be chosen non-negative (hence positive P-representation), and the expectation value of a normal-ordered operator can be readily computed. Specifically, for molecular vibronic spectra generation, 2M complex numbers (α,β) may be sampled from the generalized P-representation of an input Gaussian state and the samples may be transformed by the linear-optical unitary operation U as (α,β)→(α′,β′)=(Uα, U*β), where U represents the trans-mission/transition matrix in the phase space. For the measurement portion of the analytics, the photon number operator {circumflex over (n)}_i may be replaced with α′β′. For example, for the measurement outcome (m₁, m₂) for two modes, the corresponding observable may be written as |m₁m₁|₁⊗|m₂m₂|₂=: {circumflex over (n)}₁^m1e^{−{circumflex over (n)}1}/m₁!⊗{circumflex over (n)}₂^m2e^{−{circumflex over (n)}2}/m₂!:, so that the phase-space variable represented by (α′₁β′₁)^m1(α′₂β′₂)^m2e^{−α′1β′1−α′2β′2}/m₁!m₂! may be obtained. Here, |jj|_i represents the projector on j photon state for the ith mode, and :Ô: represents the normally ordered form of an operator Ô. For a sufficiently large number of samples, it may be guaranteed that for normal-ordered observables, the average converges to the expectation value for the output quantum state

For example, the positive P-representation for a single-mode squeezed state may be chosen as:

$\begin{array}{cc}P\left(\alpha ,\beta \right)=\frac{\sqrt{1+\gamma }}{\pi \gamma }{e}^{-\left({\alpha }^2+{\beta }^2\right)\left({\gamma }^{-1}+1/2\right)+\mathrm{\alpha \beta }},& \left(11\right)\end{array}$

where the domain is the real plane and γ≡e^2r−1 for a squeezing parameter r>0.

Using the positive P-representation of single-mode Gaussian states and the relation for a normal-ordered operator such as e^{−γ{circumflex over (n)}}=: e^{{circumflex over (n)}(e−γ−1)}:, the Fourier components may be rewritten as:

$\begin{array}{cc}\tilde{G}\left(k\right)=\langle :\exp \left[\sum _{i=1}^M{\hat{n}}_i\left(e^{-\mathrm{ik}{\omega }_i^{\prime }\theta }-1\right)\right]:\rangle & \left(12\right)\end{array}$ $\begin{array}{cc}={\int }_{{ℝ}^{2M}}d\alpha d\beta P_{\mathrm{in}}\left(\alpha ,\beta \right)\exp \left[\sum _{i=1}^M{\alpha }_i^{\prime }{\beta }_i^{\prime }\left(e^{-\mathrm{ik}{\omega }_i^{\prime }\theta }-1\right)\right],& \left(13\right)\end{array}$

where P_in(α,β) is the positive P-representation of an input squeezed state. Here, :Ô: represents the normal-ordered form of an operator Ô, and (α′,β′)=(Uα,U*β) accounts for the linear-optical unitary operation.

More specifically, when a coherent state is processed through an M-mode linear-optical circuit, it is transformed as Û|α=|Uα, where U is the corresponding M×M unitary matrix for the circuit.

Since {acute over (G)}(k) is now written as a Gaussian integral, it can be analytically obtained, which is given by

$\begin{array}{cc}\tilde{G}\left(k\right)=𝒩\frac{{\left(2\pi \right)}^M}{\sqrt{\det \left(Q\right)}}\exp \left(\frac{1}{2}{c}^T{Q}^{-1}c+c_0\right),& \left(14\right)\end{array}$ $\mathrm{where}{\varphi }_j=-k\theta {\omega }_j,$ $\begin{array}{cc}Q\equiv \left(\begin{array}{cc}2{\Gamma }^{-1}+M_{}& -U^T{\mathrm{diag}\left(e^{i{\varphi }_j}\right)}_{j=1}^M{U}^{*}\\ -U^{†}{\mathrm{diag}\left(e^{i{\varphi }_j}\right)}_{j=1}^{M}U& 2{\Gamma }^{-1}+M_{}\end{array}\right),\mathrm{and},& \left(15\right)\end{array}$ $\begin{array}{cc}𝒩\equiv \prod _{i=1}^M\frac{\sqrt{1+{\gamma }_i}}{\pi {\gamma }_i},\Gamma \equiv {\mathrm{diag}\left({\gamma }_i\right)}_{i=1}^M,\Phi ={\mathrm{diag}\left(e^{i{\varphi }_j}-1\right)}_{j=1}^M,& \left(16\right)\end{array}$ $c\equiv {\left(a^T,b^T\right)}^T,a\equiv \frac{U^T\Phi \delta }{\sqrt{2}},b\equiv \frac{U^T\Phi \delta }{\sqrt{2}},c_0\equiv \frac{{\delta }^T\Phi \delta }{2}.$

Here, δ represents a real vector accounting for the displacement in the Doktorov transformation (see Eq. (2)). More detailed derivation above and proof of the convergence of the integral are further provided later in this disclosure.

It implies that Fourier components of the molecular vibronic spectra at zero temperature have a solution. As such, the spectra may be obtained by simply taking the inverse Fourier transform, thereby avoiding having to use a quantum simulator to generate the molecular vibronic spectra.

As an example, FIG. 3 shows a comparison between an analytical solution of an example molecular vibronic spectra based on the formalism above and direct classical computation. FIG. 3 shows that the analytical solution based on the Fourier formalism above matches with the spectra obtained by directly computing probabilities. For the analytical solution, the Fourier components are first computed using Eq. (14) followed by the inverse Fourier transform to obtain the spectra. The corresponding Doktorov transformation is determined in accordance with the simulated molecular system being CH₂O₂, 1¹A′→1²A′. The number of samples processed in FIG. 2 is N=10⁴. An increase of number of bins in the frequency domain may be achieved with a polynomial computational cost.

Exact Solution for Molecular Vibronic Spectra at Finite Temperature

The same method above may be applied to initial states of the molecular system at finite temperatures. Specifically, in some example implementations, the quantum-inspired classical approach to the generation of molecular vibronic spectra may be further applied to finite temperatures other than zero temperatures.

At finite temperature T, the molecular vibronic spectra may be written as

$\begin{array}{cc}\mathrm{FCP}\left({\omega }_{\mathrm{vib}}\right)=\sum _{n,m=0}^{\infty }p\left(n,m\right)\delta \left({\omega }_{\mathrm{vib}}-\left({\omega }^{\prime }·m-\omega ·n\right)\right),\mathrm{where}& \left(17\right)\end{array}$ $\begin{array}{cc}p\left(n,m\right)=p_T\left(n\right){❘\langle n❘{\hat{U}}_{\mathrm{Dok}}❘m\rangle ❘}^2.& \left(18\right)\end{array}$

Here, p_T(n)≡n|{circumflex over (ρ)}_T|n is the probability of photon number |n from an M-mode thermal state {circumflex over (ρ)}_T at temperature T.

In some specific implementations, a Gaussian boson sampling quantum procedure described above may be correspondingly constructed for the finite temperature situations. For example, M auxiliary modes may be determined and added to the quantum simulation. Two-mode squeezed states may be further prepared to simulate thermal states. The photon numbers in the additional modes may be detect as well. The process is exemplary depicted in FIG. 4. Inspired by such quantum simulation, a classical sampling algorithm of Fourier components in the P-representation similar to the zero-temperature situation may also be employed with an addition of the auxiliary modes.

In such example implementations, the original M modes output photons correspond to in and the additional M modes' output photons corresponds to n. Therefore, the same framework described above for zero temperature may be employed to obtain the solution for molecular vibronic spectra at finite temperature, as illustrated in FIG. 4, showing a comparison of state transform between zero-temperature scenario and the finite temperature scenario.

For example, the FCP of the molecular system at the finite temperature may be rewritten as:

$\begin{array}{cc}\mathrm{FCP}=\sum _{m,n=0}^{\infty }{❘\langle m,n❘{\hat{U}}_{\mathrm{Dok}}^{\prime }❘0\rangle ❘}^2\delta \left({\omega }_{\mathrm{vib}}-\left({\omega }^{\prime }·m-\omega ·n\right)\right),& \left(19\right)\end{array}$

which is the same form as Eq. (1) except for the negative part in the delta function. Here, Û′_Dok=Û_DokÛ_TMSO is obtained by adding two-mode squeezing operators. The term “TMSO” denotes “Two-Mode-Squeezing Operator”. Each pair of two modes include an original mode and an auxiliary mode.

Again, the Fourier transform of the FCP at finite temperature may be given by

{tilde over (G)}(k)=e^{−ik{circumflex over (n)}·ω″θ},  (20)

where ω″=(ω′^T, −ω^T)^T. Therefore, the analytic solution for Fourier components can thus be found for the molecular system at a finite temperature. The vibronic spectra at the finite temperature thus may be estimated by inverse Fourier transform of the analytically generated Fourier components above.

Molecular Vibronic Spectra for Non-Gaussian Initial State

The above implementations assume a vacuum or otherwise a thermal or Gaussian initial state for the molecular vibronic modes. In a more general case, for example, the molecular system may be in a non-thermal initial state. In such circumstances, the molecular vibronic system may not be described by Gaussian states.

In a particular circumstance, the initial state of the molecular vibronic system may be a Fock state or number state of the vibronic, |n, rather than vacuum |0 a Gaussian state. For optical processes including the single vibronic level fluorescence and the resonance Raman scattering, for example, a spectra from specific vibronic levels needs to be analyzed, instead of a thermal distribution or a ground state. Indeed, a quantum simulation of such a process has also been experimentally implemented for various photoelectron processes such as photodetachment of ozone anion. In these cases, the FCP of the molecular vibronic system for such a non-Gaussian initial state may be expressed as:

$\begin{array}{cc}\mathrm{FCP}\left({\omega }_{\mathrm{vib}}\right)=\sum _{m=0}^{\infty }{❘\langle m❘{\hat{U}}_{\mathrm{Dok}}❘n\rangle ❘}^2\delta \left({\omega }_{\mathrm{vib}}-{\omega }^{\prime }·m\right),& \left(21\right)\end{array}$ $\begin{array}{cc}\tilde{G}\left(k\right)=\langle n❘{\hat{U}}_{\mathrm{Dok}}^{†}{e}^{-\mathrm{ik}\hat{n}·{\omega }^{\prime }\theta }{\hat{U}}_{\mathrm{Dok}}❘n\rangle =\langle n❘\hat{V}❘n\rangle ,& \left(22\right)\end{array}$

where {circumflex over (V)}≡Û_Dok^†e^{−k{circumflex over (n)}·w′θ}Û_Dok is a Gaussian unitary operation. Using Bloch-Messiah decomposition, the Gaussian unitary operation is decomposed as {circumflex over (V)}={circumflex over (D)}(ξ)Û_lin2Ŝ(r)Û_lin1. Because of the non-Gaussian initial state, the spectra cannot be directly mapped to a Gaussian boson sampling and so it may require a non-Gaussian type of boson sampling.

In some implementations, the Fourier components can be written as a loop hafnian of a matrix:

$\begin{array}{cc}\tilde{G}\left(k\right)=\frac{\mathrm{lhaf}\left({\tilde{\Sigma }}_n\right)}{n!Z},\mathrm{where}& \left(23\right)\end{array}$ $\begin{array}{cc}\Sigma \equiv \left(\begin{array}{cc}{U}_{\mathrm{lin}2}\tanh {\mathrm{rU}}_{\mathrm{lin}2}^T& U_{\mathrm{lin}2}\mathrm{sech}{\mathrm{rU}}_{\mathrm{lin}1}\\ U_{\mathrm{lin}1}^T\mathrm{sech}{\mathrm{rU}}_{\mathrm{lin}2}^T& -U_{\mathrm{lin}1}^T\tanh {\mathrm{rU}}_{\mathrm{lin}1}\end{array}\right)& \left(24\right)\end{array}$

is an 2M×2M complex symmetric matrix,

$\begin{array}{cc}Ϛ\equiv \left(\begin{array}{c}-U_{\mathrm{lin}2}\tanh {\mathrm{rU}}_{\mathrm{lin}2}^T·{\xi }^{*}+\xi \\ -U_{\mathrm{lin}1}^T\mathrm{sech}r{U}_{\mathrm{lin}2}^T·{\xi }^{*}\end{array}\right)& \left(25\right)\end{array}$

is an 2M-dimensional vector, and

$\begin{array}{cc}{Z}^{-1}\equiv \langle 0❘\hat{V}❘0\rangle =\frac{e^{\frac{1}{2}{\xi }^T·\Sigma ·\xi }}{\sqrt{{\prod }_{i=1}^M\cosh r_i}},& \left(26\right)\end{array}$

is the normalization factor, which is the same as the Fourier components of molecular vibronic spectra at zero temperature.

Here, {tilde over (Σ)}_n is obtained by first replacing the diagonal elements of Σ by ζ to obtain {tilde over (Σ)} and repeating ith row and column of each block matrix of {circumflex over (Σ)}n_i, times, and is thus an n×n matrix with n≡Σ_i=1^Mn_i. Therefore, computing the Fourier components reduces to computing loop hafnians of n×n complex symmetric matrices.

Loop hafnian is a quantity that is related to perfect matchings of a graph including loops (hafnian does not allow loops). The best-known algorithms computational cost of computing a loop hafnian is O(n³2^n/2) with n being the matrix size. Thus, if FC factors, which is written as a loop hafnian of a matrix whose size is n+m is to be directly computed to obtain FC profiles in Eq. (1), then the complexity would be exponential in (n+m)/2. On the other hand, the complexity of computing the Fourier components in the proposed method relies on the input photons n rather than the output photons m and the system size, i.e., the proposed method can be efficiently implemented as long as n is small enough. More precisely, since the redundancy of rows and columns for n_i≥2 does not increase the complexity of computing a loop hafnian because it does not increase the rank of the matrix, the important factor for complexity is the number of nonzero elements of n.

Meanwhile, computing the loop hafnian of a general complex symmetric matrix in multiplicative error is known to be #P-hard. It implies that computing the Fourier components of general molecular vibronic spectral with Fock state inputs and squeezing is also #P-hard. More formally, an arbitrary complex matrix may be embedded into a submatrix of a unitary matrix with an appropriate normalization so the corresponding spectra generation problem reduces to computing a permanent, which is #P-hard. The unitary matrix for Fourier coefficient is written as Ûe^{i{circumflex over (n)}·ω}Û^†, which is a diagonalized form of arbitrary unitary matrix, as shown further below.

Therefore, the Fourier components for molecular vibronic spectra mapped to Gaussian boson sampling have an analytic solution while the corresponding problem for non-Gaussian cases is a computationally hard problem, which presents a separation between Gaussian boson sampling and non-Gaussian type of boson sampling such as single-photon boson sampling. The Fourier components reduce to hafnian and permanent when there is no displacement or squeezing, respectively.

Even if a boson sampling experiment is performed using a quantum device the spectra is reproduced from the outcomes, estimation of the spectra has an additive error depending on the number of samples. Therefore, reproducing the spectra within a multiplicative error is not expected to be achieved by running a boson sampling circuit, so an additive error is the relevant quantity using a classical simulation to compare with a quantum device.

Fock States: Sampling Using Positive P-Representation

To employ a mapping from molecular vibronic spectra to Gaussian boson sampling, again, additional M modes and two-mode squeezed states of squeezing parameters si may be introduced, as illustrated in FIG. 5, which shows an comparison between Gaussian boson sampling for a thermal state, and mapping to Gaussian sampling for a Fock state. The difference from the finite temperature case above is that now the outcome n for the additional modes need to be post-selected.

Thus, the corresponding FCP can be written as

$\begin{array}{cc}G\left({\omega }_{vib}\right)=\frac{1}{N\left(n\right)}\sum _{m=0}^{\infty }p\left(m,n\right)\delta \left({\omega }_{vib}-{\omega }^{\prime }·m\right),& \left(27\right)\end{array}$

where is the normalization due to the post-selection,

$\begin{array}{cc}𝒩\left(n\right)=\prod _{i=1}^M\frac{{\tanh }^{2n_i}{s}_i}{{\cosh }^2{s}_i},& \left(28\right)\end{array}$

and p(m, n)=(n)m|Û_Dok|n².

Its Fourier components are given by

$\begin{array}{cc}\tilde{G}\left(k\right)=\frac{1}{𝒩}\langle e^{-\mathrm{ik}\theta {\sum }_{j=1}^M{\hat{n}}_j{\omega }_j^{\prime }}\otimes ❘n\rangle \langle n|\rangle .& \left(29\right)\end{array}$

Approximation of Hafnian

In some implementations, another method is to use the expression in Eq. (23). Assume that the displacement is zero, the loop hafnian in Eq. (23) then reduces to a hafnian. In this case, a randomized algorithm is constructed to estimate a hafnian of an n×n complex symmetric matrix Σ using the equality

$\begin{array}{cc}\mathrm{haf}\left(\Sigma \right)={𝔼}_{\upsilon \in {\left\{0,1\right\}}^n}\left[{\left(-1\right)}^{\sum _{i=1}^n{\upsilon }_i}\frac{2^{n/2}}{\left(n/2\right)!}{\left(h^T\Sigma h\right)}^{n/2}\right],& \left(30\right)\end{array}$

where _Σ[·] is the average over the normal distribution of covariance matrix Σ, and h=(½−v₁, . . . , ½−v_M)^T and n is even (otherwise the hafnian is zero).

Therefore, by sampling v∈{0,1}ⁿ uniformly and averaging over the samples, the hafnian may be estimated, which, together with the fact that ∥Σ∥=1 in Eq. (24), leads to the bound for estimate q of a Fourier component as

$\begin{array}{cc}\mathrm{Pr}\left[❘q-\frac{\mathrm{haf}\left({\Sigma }_n\right)}{Zn!}❘\gtrsim ϵ\frac{e^{n/2}}{\sqrt{\pi n}Z}\right]\le 2e^{-N{ϵ}^2/2},& \left(31\right)\end{array}$

where Z=√{square root over (Π_i=1ⁿ cosh r_i)} and N is the number of samples. Thus, if Π_i=1ⁿ cosh r_i>eⁿ, the estimation error is smaller than e with exponentially small failure probability if the number of samples is chosen as N=O(1/ϵ²).

It suggests that if the condition is satisfied, classically estimating the Fourier components and taking the inverse Fourier transformation renders the same scaling of precision to running a boson sampling circuit.

For nonzero displacement and loop hafnian, although there is a similar equality as in Eq. (30), the resultant error bound is more complicated than hafnian.

It is further provided another method of estimating a loop hafnian that uses the integral representation of loop hafnian in Methods. Since loop hafnian corresponds to the high order moment of multivariate normal distribution, loop hafnian of an n×n complex symmetric matrix {tilde over (Σ)} may be written as an integral:

$\begin{array}{cc}\mathrm{lhaf}\left(\tilde{\Sigma }\right)=\frac{1}{{\pi }^{4n}}{\int }_{{ℂ}^{4n}}d\alpha d\gamma d\omega d\beta \underset{i=1}{\overset{n}{\prod }}\left[\left({\alpha }_i+{\xi }_i+{\zeta }_i^{\alpha }\right){\left({\beta }_i+{\zeta }_i^{{\beta }^{*}}\right)}^{*}\right]\times \exp \left(-\frac{1}{2}{\upsilon \left(\alpha ,\beta ,\gamma ,\omega \right)}^T·{\Sigma }^{-1}·\upsilon \left(\alpha ,\beta ,\gamma ,\omega \right)\right).& \left(32\right)\end{array}$

The expression is obtained in the section below entitled “Approximation algorithm of hafnian’,’ where the explicit definition of each term is given. Therefore, the loop hafnian can be estimated by decomposing the exponential term into an imaginary part and a real part, which follows the normal distribution, and, by sampling from the normal distribution and averaging over the random variable including the imaginary part.

Fock States: Approximation of Fourier Components of Molecular Vibronic Spectra

The observation that the additive error is the achievable error from running a quantum device, some classical algorithms that might be potentially able to efficiently estimate the Fourier components within an additive error may be conceived.

For example, Gurvits's algorithm can efficiently approximate the corresponding quantity within an additive error e in single-photon boson sampling, where the FCP is written as permanent of a submatrix of a unitary matrix (see the section entitled “Computing Fourier components of vibronic spectra for single-photon input cases” below). Indeed, if a boson sampling is directly run to obtain the FCP, the Chernoff bound poses a limitation of accuracy, which is given by O(1/poly(N)) when the number of samples is N. It suggests that even though exactly computing the Fourier components might be difficult, it may be possible to efficiently approximate them within a reasonable additive error.

Since inverse Fourier transformation does not amplify the additive error of Fourier components, if there is an efficient algorithm that approximates the Fourier components within an error e in a running time O(poly(1/ϵ)), the algorithm may allow for achieving a similar accuracy as running a boson sampling.

Fock-State: Estimating Fourier Components of Spectra with Boson Sampling with Generalized Gurvits' Algorithm

For Fock-state boson sampling, the Fourier components can be written as

$\begin{array}{cc}\tilde{G}\left(k\right)=\langle n❘{\hat{U}}^{†}{e}^{-\mathrm{ik}\theta \hat{n}·\omega }\hat{U}❘n\rangle \equiv \langle n❘V❘n\rangle =\frac{\mathrm{Per}\left(V_{n,n}\right)}{n!},& \left(33\right)\end{array}$

where a unitary operator {circumflex over (V)} is defined as {circumflex over (V)}≡Û^†e^{−k{circumflex over (n)}·ωθ}Û, consisting of Û and a phase-shift operator e^{−ik{circumflex over (n)}·wθ} and Û^†.

Thus, it is a linear-optical circuit characterized by a unitary matrix V=U^†DU, with D≡diag(e^−ikθω1, . . . , e^−ikθωM) characterizing the phase-shift operator. Here, such a diagonal form suggests that V is a general form of an arbitrary unitary matrix. Together with this fact, since computing the permanent of a general matrix in multiplicative error is #P-hard and one can embed an arbitrary complex matrix into a submatrix of a unitary matrix by normalizing the matrix, computing its Fourier component in multiplicative error is also #P-hard. Therefore, it shows that computing Fourier components in multiplicative error is a #P-hard problem, and consequently, it proves that the exact computation of the spectra is also a #P-hard problem.

In fact, it is noted that even if a boson sampling experiment runs using a quantum device and reproduces the spectra from the sampling outcomes, the resultant spectra has an additive sampling error. Therefore, reproducing the spectra within a multiplicative error is not expected to be achievable by running a boson sampling circuit with polynomially many samples. As such, an additive error is a relevant target using a classical simulation to compare with a quantum device.

Interestingly, there exists a classical algorithm, the so-called Gurvits' algorithm, which can efficiently approximate the permanent of a matrix within an additive error. However, this algorithm is not sufficient for general Fock-state boson sampling cases where n_i contains more than a single-photon for some i's because the error can increase exponentially in the system size in that case (described in further detail below). The Gurvits' algorithm, however, may be modified to estimate the spectra which includes multiphotons n using the following equality (see detailed derivation in the next section):

$\begin{array}{cc}\frac{\mathrm{Per}\left(V_{n,n}\right)}{n!}={𝔼}_{x\in 𝒳}\left[\prod _{i=1}^M{\left(\frac{{{\stackrel{\_}{y}}_i\left(\mathrm{Vy}\right)}_i}{n_i}\right)}^{n_i}\right],& \left(34\right)\end{array}$

where y_i≡√{square root over (n_i)}x_i, x∈≡[n₁+1]× . . . ×[n_M+1], where [j] is the set of jth roots of unity.

Thus, by sampling the random variable with uniform x∈, the randomized algorithm gives an estimate μ of permanent of an n×n matrix such that

$\begin{array}{cc}❘\frac{\mathrm{Per}\left(V_{n,n}\right)}{n!}-\mu ❘<ϵ{V}^n,& \left(35\right)\end{array}$

with high probability 1−δ with running time T=O(poly(n, 1/ϵ, log δ⁻¹)). Here, ∥V∥ is the spectral norm of the matrix V, which is always 1 for this case because V is a unitary matrix.

Thus, the modified (or generalized) Gurvits' algorithm enables an efficiently approximation of Fourier components within a reasonable additive error even though computation of the Fourier components in a multiplicative error is hard (#P-hard). Note that although a modified or generalized Gurvits' algorithm might be used to approximate individual probabilities in Eq. (9), it is nontrivial to approximate the sum of exponentially many probabilities, which is enabled by approximating Fourier components instead of probabilities.

The remaining challenge is the propagation of the error of Fourier coefficients to that of the spectra through the inverse Fourier transformation. Using Parseval's relation, we prove that as long as we estimate the Fourier coefficients with a small error ϵ, the transformed spectra's error is also small as ϵ:

$\begin{array}{cc}\sum _{\Omega =0}^{{\Omega }_{\max }}{❘\Delta G\left(\Omega \right)❘}^2=\frac{1}{{\Omega }_{\max }+1}\sum _{k=0}^{{\Omega }_{\max }}{❘\Delta G\left(k\right)❘}^2\le {ϵ}^2,& \left(36\right)\end{array}$

which proves that |ΔG(Ω)|≤ϵ for any Ω, where ΔG and Δ{tilde over (G)} represent the error of spectra and Fourier component estimation.

Hence, if there is an efficient algorithm that approximates the Fourier components within an error e in a running time T=O(poly(M, 1/ϵ, log δ⁻¹)), the algorithm enables us to achieve the same accuracy as running a boson sampling. For Fock-state boson sampling case, Gurvits's algorithm is evidently such an algorithm estimating Fourier components. Consequently, we have shown that the molecular vibronic spectra problem corresponding to Fock-state boson sampling can be efficiently solved by a classical computer as accurately as running a boson sampler, which indicates that there is no quantum advantage from this problem.

Gurvits' Algorithm with Repeated Rows and Columns

The disclosure below first presents the original Gurvits's algorithm and generalized algorithm to improve the bound when rows and columns are repeated. The Gurvits's algorithm exploits the following equality:

$\begin{array}{cc}\mathrm{Per}\left(X\right)={𝔼}_{x\in {\left\{-1,1\right\}}^n}\left[\prod _{i=1}^n\left(\sum _{j=1}^n{x}_i{X}_{\mathrm{ij}}{x}_j\right)\right].& \left(37\right)\end{array}$

Thus, by sampling the random variable with uniform x∈{−1, 1}ⁿ, the randomized algorithm gives an estimate μ of permanent of an n×n matrix X such that

|Per(X)−μ|<ϵ∥X∥ⁿ,  (38)

with high probability 1−δ with running time T=O(poly(n, 1/ϵ, log δ⁻¹)).

When a submatrix of a unitary matrix is considered, this expression is sufficient because the submatrix's norm is bounded by 1. However, when the rows and columns in X from a unitary matrix are repeated, which is the case when dealing with multiphoton inputs, the norm ∥X∥ is not generally bounded by 1. Therefore, the upper bound of the error can increase as ∥X∥ⁿ, i.e., exponentially in n. To deal with such multiphoton cases with a reduced error, Gurvits's algorithm may need to be generalized. Note that although the Gurvits's algorithm for multiphoton outputs may be generalized, i.e., Per(V_n,m)=n|{circumflex over (V)}|m, where n∈{0, 1}ⁿ and in is a nonnegative integer vector, that is not sufficient for the purpose of handling multiphoton input and output, i.e., n|{circumflex over (V)}|n=Per (V_n,n).

Consider a case where the goal is to approximate a quantity

$\begin{array}{cc}\frac{\mathrm{Per}}{n!},& \left(39\right)\end{array}$

where A∈^n×n which is obtained by repeating ith row and column of a matrix B∈^k×k for n_i times. Here, n∈_≥0^k with Σ_i=1^k n_i=n. Further define a random variable

$\begin{array}{cc}{\mathrm{GenGly}}_x\left(A\right)\equiv \frac{1}{{\prod }_{i=1}^k{n}_i^{n_i}}\prod _{i=1}^k{\stackrel{\_}{y}}_i^{n_i}\prod _{i=1}^k{\left(y_i{B}_{i,1}+\dots +y_k{B}_{i,k}\right)}^{n_i}=\prod _{i=1}^k{\left(\frac{{{\stackrel{\_}{y}}_i\left(\mathrm{By}\right)}_i}{n_i}\right)}^{n_i},& \left(40\right)\end{array}$

where γ_i≡n_i^1/2x_i. Here, x∈[n₁+1]× . . . ×[n_k+1]≡, where [j] is the set of jth roots of unity. Then, its absolute value is given as

$\begin{array}{cc}❘{\mathrm{GenGly}}_x\left(A\right)❘={\left(\prod _{i=1}^k{n}_i^{n_i}\right)}^{-1/2}❘\prod _{i=1}^k{\left(y_i{b}_{i,1}+\dots +y_k{b}_{i,k}\right)}^{n_i}❘={\left(\prod _{i=1}^k{n}_i^{n_i}\right)}^{-1/2}❘\prod _{i=1}^k{c}_i^{n_i}❘,& \left(41\right)\end{array}$

where c_i≡γ₁B_i,1+ . . . +γ_kB_i,k. Now, let us find its upper bound, which determines the error of the corresponding randomized algorithm by Chernoff bound. First, we have a constraint

$\begin{array}{cc}\sum _{i=1}^k{❘c_i❘}^2=\sum _{i=1}^k{❘{\left(\mathrm{By}\right)}_i❘}^2\le {B}^2{y}^2=n{B}^2.& \left(42\right)\end{array}$

Further considering:

$\begin{array}{cc}\log \prod _{i=1}^k{❘c_i❘}^{n_i}=\sum _{i=1}^k{n}_i\log ❘c_i❘.& \left(43\right)\end{array}$

Using lagrange multiplier, with a constraint Σ_i=1^k|c_i|²=c,

$\begin{array}{cc}Z\equiv \sum _{i=1}^k{n}_i\log ❘c_i❘-\lambda \left(\sum _{i=1}^k{❘c_i❘}^2-c\right),& \left(44\right)\end{array}$

It may be obtained that:

$\begin{array}{cc}\frac{\partial Z}{\partial ❘c_i❘}=\frac{n_i}{❘c_i❘}=2\lambda ❘c_i❘=0,\forall i.& \left(45\right)\end{array}$

Using the constraint Eq. (42),

$\begin{array}{cc}\sum _{i=1}^k{❘c_i❘}^2=\sum _{i=1}^k\frac{n_i}{2\lambda }=\frac{n}{2\lambda }=c,& \left(46\right)\end{array}$

the solution is obtained as

$\begin{array}{cc}❘c_i❘=\sqrt{\frac{n_i}{2\lambda }}=\sqrt{\frac{{\mathrm{cn}}_i}{n}}.& \left(47\right)\end{array}$

Therefore, using the solution, we have an upper bound as

$\begin{array}{cc}\prod _{i=1}^k{❘c_i❘}^{n_i}\le \exp \left(\sum _{i=1}^k{n}_i\log \sqrt{\frac{{\mathrm{cn}}_i}{n}}\right)=\prod _{i=1}^k{n_i^{n_i/2}\left(\frac{c}{n}\right)}^{n/2}\le \prod _{o=1}^k{n}_i^{n_i/2}{B}^n.& \left(48\right)\end{array}$

Hence, there is an upper bound of the absolute value of the random variable,

$\begin{array}{cc}❘{\mathrm{GenGly}}_x\left(A\right)❘={\left(\prod _{i=1}^k{n}_i^{n_i}\right)}^{-1/2}❘\prod _{i=1}^k{c}_i^{n_i}❘\le {B}^n.& \left(49\right)\end{array}$

Next, it may be proved that the average of the random variables gives an unbiased estimator of Per(A)/n!, i.e.,

$\begin{array}{cc}{𝔼}_{𝕩\in \chi }\left[{\mathrm{GenGly}}_x\left(A\right)\right]=\frac{\mathrm{Per}\left(A\right)}{n!}.\mathrm{First},& \left(50\right)\\ {𝔼}_{𝕩\in \chi }\left[\frac{1}{\prod _{i=1}^k{n}_i^{n_i}}\prod _{i=1}^k{\stackrel{\_}{y}}_i^{n_i}\prod _{i=1}^k{\left(y_1{B}_{i,1}+\cdots +y_k{B}_{i,k}\right)}^{n_i}\right]& \left(51\right)\\ =\frac{1}{\prod _{i=1}^k{n}_i^{n_i}}{𝔼}_{𝕩\in \chi }\left[\prod _{i=1}^k{{\stackrel{\_}{y}}_i^{n_i}\left(y_1,B_{i,1}+\cdots +y_k{B}_{i,k}\right)}^{n_i}\right]& \left(52\right)\\ =\frac{1}{\prod _{i=1}^k{n}_i^{n_i}}\sum _{{\sigma }_1,\dots ,{\sigma }_n\in \left[k\right]}\left(\prod _{i=1}^k{B}_{i,{\sigma }_i}\right){𝔼}_{𝕩\in \chi }\left[\prod _{i=1}^k{\stackrel{\_}{y}}_i^{n_i}\prod _{i=1}^n{y}_{{\sigma }_i}\right],& \left(53\right)\end{array}$

where by the symmetry over the roots of unity,

$\begin{array}{cc}{𝔼}_{𝕩\in \chi }\left[\prod _{i=1}^k{\stackrel{\_}{y}}_i^{n_i}\prod _{i=1}^n\right]=0,& \left(54\right)\end{array}$

unless the product insides is equal to Π_i=1^k|γ_i|²ⁿⁱ, in which case it may be derived that

$\begin{array}{cc}{𝔼}_{𝕩\in \chi }\left[\prod _{i=1}^k{\stackrel{\_}{y}}_i^{n_i}\prod _{i=1}^n\right]=\prod _{i=1}^k{n}_i^{n_i}.& \left(55\right)\end{array}$

Hence,

$\begin{array}{cc}{𝔼}_{𝕩\in \chi }\left[{\mathrm{GenGly}}_x\left(A\right)\right]=\sum _{t_1\dots ,t_k\in \left[k\right]:❘\left\{i:t_i=j\right\}❘=n_j\forall j\in \left[k\right]}\prod _{i=1}^k{B}_{i,t_i}=\frac{\mathrm{Per}\left(A\right)}{n!},& \left(56\right)\end{array}$

where n! appears to take into account the fact that each produce Π_i=1^kb_i,tk appears n! times.

Therefore, GenGly_x(A) gives us an unbiased estimator of Per(A)/n!, and the random variables are bounded by ∥B∥ⁿ. Due to the Chernoff bound, the randomized algorithm with number of samples N=O(1/ϵ²) provides an estimate μ with high probability 1−δ such that

$\begin{array}{cc}❘\frac{\mathrm{Per}\left(A\right)}{n!}-\mu ❘<ϵ{B}^n.& \left(57\right)\end{array}$

More Details of Exact Solution of Fourier Components of Molecular Vibronic Spectra

A detailed derivation of the analytic solution of Fourier components of the vibronic spectra is described below.

The initial positive P-representation above may be represented by:

$\begin{array}{cc}{P}_{\mathrm{in}}\left(\alpha ,\beta \right)=\prod _{i=1}^M\left[\frac{\sqrt{1+{\gamma }_i}}{{\mathrm{\pi \gamma }}_i}\exp \left(-\left({\alpha }_i^2+{\beta }_i^2\right)\left({\gamma }_i^{-1}+1/2\right)+{\alpha }_i{\beta }_i\right)\right]& \left(58\right)\\ =𝒩\exp \left[-\alpha ·A·\alpha -\beta ·B·\beta -\alpha ·C·\beta -\beta ·C^T·\alpha \right],& \left(59\right)\end{array}$

where

$\begin{array}{cc}A\equiv {\mathrm{diag}\left({\gamma }_i^{-1}+1/2\right)}_{i=1}^M,B\equiv {\mathrm{diag}\left({\gamma }_i^{-1}+1/2\right)}_{i=1}^M,C\equiv -1/2,𝒩\equiv \prod _{i=1}^M\frac{\sqrt{1+{\gamma }_i}}{{\mathrm{\pi \gamma }}_i}.& \left(60\right)\end{array}$

Here, α and β are real M-dimensional vectors. In the disclosure above, the expression of the Fourier components of molecular vibronic spectra has been found as:

{tilde over (G)}(k)=e^{−ikθ{circumflex over (n)}·w′},  (61)

where the expectation value is over the output state of Gaussian boson sampling.

By using the relation e^{−γ{circumflex over (n)}}=: e^{{circumflex over (n)}(e−γ−1)}:, the operators can be written in a normal-ordered form and the Fourier transform may have the following property:

$\begin{array}{cc}\tilde{G}\left(k\right)=\sum _{{\omega }_0=0}^K\sum _{m=0}^{\infty }p\left(m\right)\delta \left({\omega }_0-\omega ·m\right)e^{-\mathrm{ik}{\mathrm{\theta \omega }}_0}& \left(62\right)\\ =\sum _{m=0}^{\infty }p\left(m\right)e^{-\mathrm{ik}\mathrm{\theta \omega }·m}& \left(63\right)\\ =\langle \underset{i=1}{\overset{M}{\otimes }}\sum _{m_i=0}^M❘m_i\rangle \langle m_i❘e^{-\mathrm{ik}{\mathrm{\theta \omega }}_i{m}_i}\rangle & \left(64\right)\\ =\langle \underset{i=1}{\overset{M}{\otimes }}\sum _{m_i=0}^M❘m_i\rangle \langle m_i❘e^{-\mathrm{ik}{\mathrm{\theta \omega }}_i{\hat{n}}_i}\rangle & \left(65\right)\\ =\langle \underset{i=1}{\overset{M}{\otimes }}{e}^{-\mathrm{ik}{\mathrm{\theta \omega }}_i{\hat{n}}_i}\rangle & \left(66\right)\\ =\langle :\exp \left[\sum _{i=1}^M{\hat{n}}_i\left(e^{-{\mathrm{ikw}}_i\theta }-1\right)\right]:\rangle ,& \left(67\right)\end{array}$

where θ=2π/(K+1). For the last equality, we have used the relation e^{−γ{circumflex over (n)}}=: e^{{circumflex over (n)}(e−γ−1)}:.

The expectation value may be computed as:

$\begin{array}{cc}\tilde{G}\left(k\right)=\langle :\exp \left[\sum _{i=1}^M{\hat{n}}_i\left(e^{i{\varphi }_i}-1\right)\right]:\rangle ={\langle \exp \left[\sum _{i=1}^M{\alpha }_i^{\prime }{\beta }_i^{\prime }\left(e^{i{\varphi }_i}-1\right)\right]\rangle }_P& \left(68\right)\end{array}$ $\begin{array}{cc}={\langle \exp \left[\sum _{i=1}^M\left(\sum _{j=1}^M{U}_{\mathrm{ij}}{\alpha }_j+{\alpha }_{0i}\right)\left(\sum _{k=1}^M{U}_{\mathrm{ik}}^{*}{\beta }_k+{\beta }_{0i}\right)\left(e^{i{\varphi }_i}-1\right)\right]\rangle }_P& \left(69\right)\end{array}$ $\begin{array}{cc}={\langle \exp \left[\alpha ·𝒰·\beta +a·\alpha +b·\beta +c_0\right]\rangle }_P& \left(70\right)\end{array}$ $\begin{array}{cc}=𝒩\int d\alpha d\beta \exp \left[-\alpha ·A·\alpha -\beta ·B·\beta -\alpha ·\tilde{C}·\beta -\beta ·{\tilde{C}}^T·\alpha +a·\alpha +b·\beta +c_0\right]& \left(71\right)\end{array}$ $\begin{array}{cc}=𝒩\int \mathrm{dq}\exp \left[-q·Q·q/2+c·q+c_0\right],\mathrm{where}{\varphi }_i\equiv -k{\mathrm{\theta \omega }}_i^{\prime },& \left(72\right)\end{array}$ $\begin{array}{cc}𝒰\equiv U^T\Phi U^{*},\Phi \equiv \mathrm{diag}\left(e^{i{\varphi }_1}-1,\dots ,e^{i{\varphi }_M}-1\right),\tilde{C}\equiv C-𝒰/2,Q=2\left(\begin{array}{cc}A& \tilde{C}\\ {\tilde{C}}^T& B\end{array}\right),& \left(73\right)\end{array}$ $\begin{array}{cc}c\equiv \left(\begin{array}{c}a\\ b\end{array}\right),a=U^T{\mathrm{\Phi \beta }}_0,b=U^{†}{\mathrm{\Phi \alpha }}_0,c_0={\alpha }_0{\mathrm{\Phi \beta }}_0,& \left(74\right)\end{array}$

and represents the average over the positive P-representation.

It is noted that the displacement occurs after the linear-optical unitary, which matches the displacement δ/√{square root over (2)} in the Doktorov transformation, i.e., α₀=β₀=δ/√{square root over (2)}. It is further noted that Q's real part is positive-definite.

To see this, Q may be decomposed as:

$\begin{array}{cc}Q=22_{}\otimes {\Gamma }^{-1}+\left(\begin{array}{cc}{M}_{}& -M_{}\\ -M_{}& M_{}\end{array}\right)-\left(\begin{array}{cc}0& 𝒰\\ {𝒰}^T& 0\end{array}\right)& \left(75\right)\\ =22_{}\otimes {\Gamma }^{-1}+\left(\begin{array}{cc}{M}_{}& -U^T{\mathrm{diag}\left(e^{i{\varphi }_j}\right)}_{j=1}^M{U}^{*}\\ -U^{†}{\mathrm{diag}\left(e^{i{\varphi }_j}\right)}_{j=1}^{M}U& M_{}\end{array}\right),& \left(76\right)\end{array}$

where Γ≡diag(γ₁, . . . , γ_M).

The real part of Q may then be written as:

$\begin{array}{cc}\mathrm{Re}\left(Q\right)=22_{}\otimes {\Gamma }^{-1}+\left(\begin{array}{cc}{M}_{}& L\\ L^T& M_{}\end{array}\right),& \left(77\right)\end{array}$

where L≡Re(−U^T diag(e^iϕj)_j=1^MU*). Since ∥L∥_op≤1, Re(Q) is positive-definite, the integral always converges.

Using the Gaussian integral formula,

$\begin{array}{cc}\int \mathrm{dq}\exp \left[-q·Q·q/2+c·q\right]=\sqrt{\frac{{\left(2\pi \right)}^{2M}}{❘Q❘}}\exp \left(\frac{1}{2}c·Q^{-1}·c\right),& \left(78\right)\end{array}$

the analytical form of the Fourier components may be obtained as:

$\begin{array}{cc}\tilde{G}\left(k\right)=𝒩\sqrt{\frac{{\left(2\pi \right)}^{2M}}{❘Q❘}}\exp \left(\frac{1}{2}c·Q^{-1}·c+c_0\right),\mathrm{where}& \left(79\right)\\ Q=\left(\begin{array}{cc}2{\Gamma }^{-1}+M_{}& -U^T{\mathrm{diag}\left(e^{i{\varphi }_j}\right)}_{j=1}^M{U}^{*}\\ -U^{†}{\mathrm{diag}\left(e^{i{\varphi }_j}\right)}_{j=1}^{M}U& 2{\Gamma }^{-1}+M_{}\end{array}\right).& \left(80\right)\end{array}$

More Detail of the Positive P-Representation

The generalized (positive) P-representation introduced above is described in further detail. Such generalized positive P representation may be written as:

$\begin{array}{cc}\hat{\rho }={\int }_{{ℂ}^{2M}}P\left(\alpha ,\beta \right)\hat{\Lambda }\left(\alpha ,\beta \right)d\alpha d\beta ,\hat{\Lambda }\left(\alpha ,\beta \right)\equiv \frac{❘\alpha \rangle \langle {\beta }^{*}❘}{\langle {\beta }^{*}❘\alpha \rangle },& \left(81\right)\end{array}$

where |α and |β* are coherent states.

The generalized P-representation represents one of the quasi-probability distributions of a bosonic state. An important property of the generalized P-representation is that the distribution P(α,β) can always be chosen non-negative (hence and name positive P-representation) and the expectation value of a normal-ordered operator can be readily computed.

For simplicity, a simpler version may be used, where it is assumed that the Glauber-Sudarshan P-function exists, which is written as

{circumflex over (ρ)}=d²αP_GS(α)|αα|,  (82)

for a single-mode state.

It is noted that although the existence of the Glauber-Sudarshan P-function is assumed for simplicity, such assumption may be true for arbitrary density matrices. The assumption is that the positive P-representation can be written as

$\begin{array}{cc}P\left(\alpha ,\beta \right)=\frac{1}{4{\pi }^2}\exp \left(-\frac{{❘\alpha -{\beta }^{*}❘}^2}{4}\right)\langle \left(\alpha +{\beta }^{*}\right)/2❘\hat{\rho }❘\left(\alpha +{\beta }^{*}\right)/2\rangle .& \left(83\right)\end{array}$

If the assumption is true, then it guarantees that the representation is positive and real.

By directly substituting the representation, it is obtained that:

$\begin{array}{cc}P\left(\alpha ,\beta \right)=\frac{1}{4{\pi }^2}\int P_{\mathrm{GS}}\left(\gamma \right)\exp \left(-\frac{1}{2}{❘\alpha -\gamma ❘}^2-\frac{1}{2}{❘{\beta }^{*}-\gamma ❘}^2\right)d^2\gamma .& \left(84\right)\end{array}$

In the meanwhile, for an analytical function ƒ:

$\begin{array}{cc}f\left(\gamma \right)=\frac{1}{2\pi }{\int }_{ℂ}f\left(\alpha \right)\exp \left(-\frac{1}{2}{❘\alpha -\gamma ❘}^2\right)d^2\alpha .& \left(85\right)\end{array}$

It can be shown that, for example at γ=0,

$\begin{array}{cc}{\int }_{ℂ}f\left(\alpha \right)\exp \left(-\frac{1}{2}{❘\alpha ❘}^2\right)d^2\alpha =2\pi \exp \left(\frac{1}{2}\frac{{\partial }^2}{\partial \alpha \partial {\beta }^{*}}\right)f\left(\alpha \right){❘}_{\alpha =0}=2\pi f\left(0\right),& \left(86\right)\end{array}$

where the last equality holds because the function is analytic.

Using this property, the following may be obtained from Eq. (84)

$\begin{array}{cc}{\int }_{{ℂ}^2}\hat{\Lambda }\left(\alpha ,\beta \right)P\left(\alpha ,\beta \right)d^2\alpha d^2\beta & \left(87\right)\end{array}$ $\begin{array}{cc}=\frac{1}{4{\pi }^2}{\int }_{{ℂ}^3}{d}^2\alpha d^2\beta d^2\gamma \hat{\Lambda }\left(\alpha ,\beta \right)P_{\mathrm{GS}}\left(\gamma \right)\exp \left(-\frac{1}{2}{❘\alpha -\gamma ❘}^2-\frac{1}{2}{❘{\beta }^{*}-\gamma ❘}^2\right)& \left(88\right)\end{array}$ $\begin{array}{cc}={\int }_{ℂ}{d}^2\alpha P_{\mathrm{GS}}\left(\alpha \right)❘\alpha \rangle \langle \alpha ❘=\hat{\rho },& \left(89\right)\end{array}$

which proves the assumption above.

It may be further shown that the expectation value of a normal-ordered operator can be computed by averaging over the phase-space variable. For example, when computing for a single-mode: {circumflex over (α)}{circumflex over (α)}^†{circumflex over (α)}{circumflex over (α)}^†:=α^†2{circumflex over (α)}²,

$\begin{array}{cc}\mathrm{Tr}\left[\hat{\rho }:\hat{a}{\hat{a}}^{†}\hat{a}{\hat{a}}^{†}:\right]=\mathrm{Tr}\left[\hat{\rho }{\hat{a}}^{\mathrm{†2}}{\hat{a}}^2\right]={\int }_{{ℂ}^2}P\left(\alpha ,\beta \right)\mathrm{Tr}\left[\hat{\Lambda }\left(\alpha ,\beta \right){\hat{a}}^{\mathrm{†2}}{\hat{a}}^2\right]d\alpha d\beta & \left(90\right)\end{array}$ $\begin{array}{cc}={\int }_{{ℂ}^2}P\left(\alpha ,\beta \right)\frac{\langle {\beta }^{*}❘{\hat{a}}^{\mathrm{†2}}{\hat{a}}^2❘\alpha \rangle }{\langle {\beta }^{*}❘\alpha \rangle }d\alpha d\beta ={\int }_{{ℂ}^2}P\left(\alpha ,\beta \right){\beta }^2{\alpha }^{2}d\beta d\beta .& \left(91\right)\end{array}$

Therefore, the expectation value may be directly integrated or estimated by sampling (α,β) from P and averaging α²β² over the samples. The same procedure works for any normal-ordered operator ƒ({circumflex over (α)}, {circumflex over (α)}^†) by replacing {circumflex over (α)}→α and {circumflex over (α)}^†→β. Especially when the operator is written as a function of {circumflex over (n)}, g({circumflex over (n)}), then the corresponding phase-variable is given by g(αβ).

When applying a linear-optical circuit to a given input state, the quantum state transforms as

{circumflex over (ρ)}= P(α,β){circumflex over (Λ)}(α,β)dαdβ→{circumflex over (ρ)}=P(α,β){circumflex over (Λ)}(Uα,U*β)dαdβ.  (92)

which is written for the positive P-representation as

{circumflex over (ρ)}=P(α,β){circumflex over (Λ)}(α,β)dαdβ→{circumflex over (ρ)}=P(α,β){circumflex over (Λ)}(Uα,U*β)dαdβ.  (93)

Therefore, sampling from the positive P-representation of the output state of the linear-optical circuit is equivalent to sampling from the positive P-representation of the input state and transform the phase variable as (α,β)→(α′,β′)=(Uα,U*β).

Numerical Simulation

The various algorithms above may be numerically simulated to generate molecular vibronic spectra.

The first example involves photo-electron spectra of formic acid (CH₂O₂, 1¹A′→1²A′). To quantify the performance of the scheme based on Gaussian boson sampling and the proposed quantum-inspired classical algorithm, the ideal FCP distribution, generated by directly computing probability distributions of the outcomes of Gaussian boson sampling, may be compared to the empirical distributions obtained by Gaussian boson sampling simulation and the proposed classical algorithm.

It may be emphasized that the ideal Gaussian boson sampling assumes no errors and losses. FIGS. 3 and 6 show that both schemes predict the ground truth distribution well. In particular, FIG. 3 shows an example data and logic flow for the quantum-inspired classical algorithm of FIG. 1 for generating molecular vibronic spectra, whereas FIG. 6 shows comparison of molecular vibronic spectra for CH₂O₂, 1¹A′→1²A′ as generated based on direct computation (ideal bars) and based on Gaussian boson sampling simulation (“GBS” points and curve).

In FIG. 7, the total variation distances by varying the number of samples is shown. The total variance distance (TVD) is defined as

$\begin{array}{cc}\mathrm{TVD}\equiv \frac{1}{2}\sum _{{\omega }_{\mathrm{vib}}}❘\mathrm{FCP}\left({\omega }_{\mathrm{vib}}\right)-\mathrm{eFCP}\left({\omega }_{\mathrm{vib}}\right)❘.& \left(94\right)\end{array}$

FIG. 7 particular shows comparison of total variance distance (TVD) between ideal direct computation (upper curve with circular points) and the quantum-inspired classical algorithm (lower curve with triangular points) of FIG. 3. FIG. 7 shows that while both schemes accuracy increases as the number of samples increases, the proposed classical algorithm provides a better precision (about 5 times).

The second example is generated from a recent Gaussian boson sampling experiments samples in order to verify that the quantum-inspired classical algorithm operates in the quantum supremacy regime.

Specifically, arbitrarily chosen virtual harmonic angular frequencies w′ are introduced for the output modes to demonstrate that the algorithm can even simulate a quantum supremacy regime for molecular vibronic spectra generation. For FIGS. 8 and 9, nonzero weights are set only for the first 6 modes to compare with its ideal distribution. More details for this method are further disclosed below.

Since it corresponds to a Gaussian boson sampling on the small number of marginal modes, all ideal FC factors may be directly computed. FIGS. 8 and 9 again show that both schemes predict the virtual spectra correctly, while the spectra from the classical algorithm is more accurate (about 4 times) in terms of total variation distance.

Also, it is further observed that the total variance distance from the ideal Gaussian boson sampling is smaller than the one from the experimental Gaussian boson sampling, which is due to the experimental noise.

Finally, a random weight vector w′ has been generated for all 144 modes and compared the spectra from Gaussian boson sampling and classical algorithm, which is shown in FIG. 10. While the spectra generated from Gaussian boson sampling fluctuates, which may be due to experimental noise, the classical algorithms spectra has the very similar shape to Gaussian boson samplings spectra. Therefore, it confirms that the quantum-inspired classical algorithm gives a consistent result as the Gaussian boson sampler in the quantum supremacy regime.

Approximation Algorithm of Hafnian

An example algorithm to approximate the hafnian above is provided. First, it is known that

$\begin{array}{cc}{x}_1^{n_1}{x}_2^{n_2}\cdots x_M^{n_M}=\frac{1}{n!}\sum _{v_1=0}^{n_1}\cdots \sum _{v_M=0}^{n_M}{\left(-1\right)}^{\sum _{i=1}^M{v}_i}\left(\begin{array}{c}{n}_1\\ v_1\end{array}\right)\cdots \left(\begin{array}{c}{n}_M\\ v_M\end{array}\right){\left(\sum _{i=1}^M{h}_i{x}_i\right)}^n,& \left(95\right)\end{array}$

where h_i=n_i/2−v_i.

Also a known proposition states that for z=(z₁, . . . , z_M)^T˜(0,Σ), where Σ is a real symmetric matrix and z can be a complex variable vector,

$\begin{array}{cc}\mathrm{haf}\left({\Sigma }_n\right)={𝔼}_z\left[\prod _{i=1}^M{z}_i^{n_i}\right]=\frac{1}{\left(n/2\right)!}\sum _{v_1=0}^{n_1}\cdots \sum _{v_M=0}^{n_M}{\left(-1\right)}^{\sum _{i=1}^M{v}_i}\left(\begin{array}{c}{n}_1\\ v_1\end{array}\right)\cdots \left(\begin{array}{c}{n}_M\\ v_M\end{array}\right){\left(\frac{h^T\Sigma h}{2}\right)}^{n/2},& \left(96\right)\end{array}$

where n=n₁+ . . . +n_M and h=(n₁/2−v₁, . . . , n_M/2−v_M)^T and n is even.

It is noted that the equality between the hafnian and the sum in Eq. (96) holds even for general complex symmetric matrices. Especially when n_i=1 for all i∈{1, . . . , n},

$\begin{array}{cc}{𝔼}_z\left[\prod _{i=1}^n{z}_i\right]=\frac{1}{\left(n/2\right)!}\sum _{v_1=0}^1\dots \sum _{v_M=0}^1{\left(-1\right)}^{{\sum }_{i=1}^n{v}_i}{\left(\frac{h^T\sum h}{2}\right)}^{n/2}& \left(97\right)\end{array}$ $\begin{array}{cc}=\frac{1}{2^n}\sum _{v_1=0}^1\dots \sum _{v_n=0}^1{\left(-1\right)}^{{\sum }_{i=1}^n{v}_i}\frac{2^n}{\left(n/2\right)!}{\left(\frac{h^T\sum h}{2}\right)}^{n/2}& \left(98\right)\end{array}$ $\begin{array}{cc}={𝔼}_v\left[{\left(-1\right)}^{{\sum }_{i=1}^n{v}_i}\frac{2^{n/2}}{\left(n/2\right)!}{\left(h^T\sum h\right)}^{n/2}\right],& \left(99\right)\end{array}$

where h=(½−v₁, . . . , ½−v_M)^T. In this case, M=n.

Therefore,

$\begin{array}{cc}{\left(-1\right)}^{{\sum }_{i=1}^n{v}_i}\frac{2^{n/2}}{\left(n/2\right)!}{\left(h^T\sum h\right)}^{n/2}& \left(100\right)\end{array}$

is an unbiased estimator for the hafnian in the sense that

$\begin{array}{cc}{𝔼}_v\left[{\left(-1\right)}^{{\sum }_{i=1}^n{v}_i}\frac{2^{n/2}}{\left(n/2\right)!}{\left(h^T\sum h\right)}^{n/2}\right]=\mathrm{haf}(\sum ).& \left(101\right)\end{array}$

Note that

$\begin{array}{cc}❘\frac{2^{n/2}}{\left(n/2\right)!}{\left(h^T\sum h\right)}^{n/2}❘=❘\frac{2^{n/2}}{\left(n/2\right)!}{\left(\frac{n}{4}{\hat{h}}^T\sum \hat{h}\right)}^{n/2}❘=❘\frac{n^{n/2}}{\left(n/2\right)!2^{n/2}}{\left({\hat{h}}^T\sum \hat{h}\right)}^{n/2}❘\le \frac{{\left(n\sum \right)}^{n/2}}{\left(n/2\right)!2^{n/2}}.& \left(102\right)\end{array}$

Therefore, in general, the estimate by sampling v and using the estimator has an error bound given by Chernoff bound as

$\begin{array}{cc}\mathrm{Pr}\left[❘\stackrel{\sim }{p}-\mathrm{haf}(\sum )❘>ϵ\frac{{\left(n\sum \right)}^{n/2}}{\left(n/2\right)!2^{n/2}}\right]\le 2\exp \left(-\frac{1}{2}{ϵ}^{2}T\right).& \left(103\right)\end{array}$

Here, ∥Σ∥ is the spectral norm, i.e., the maximum singular value.

The Fourier components of molecular vibronic spectra, which is given by

$\begin{array}{cc}\stackrel{\sim }{G}\left(k\right)=\frac{\mathrm{haf}\left({\sum }_n\right)}{Zn!}.& \left(104\right)\end{array}$

Note that the operator norm of the matrix Σ_α,β* is 1 (see Eq. (??)). Since a matrix's submatrix's operator norm is smaller than or equal to the matrix's norm, when the vector n is composed 0 or 1, Σ_n's operator norm is smaller than or equal to 1.

In addition, when n has a component larger than 1, hafnian for a matrix obtained may be computed by repeating the row and column. Nevertheless, when repeating the row and column the operator norm does not increase as much as n!, i.e.,

$\begin{array}{cc}\frac{{\sum }_n}{n!}\le \sum .& \left(105\right)\end{array}$

Therefore, the estimate by sampling x and using the estimator has an error bound given by Chernoff bound as

$\begin{array}{cc}\mathrm{Pr}\left[❘\stackrel{\sim }{q}-\frac{\mathrm{haf}(\sum )}{Zn!}❘\gtrsim ϵ\frac{e^{n/2}}{\sqrt{\pi n}Z}\right]\le 2\exp \left(-\frac{1}{2}{ϵ}^{2}T\right),& \left(106\right)\end{array}$

where Z=√{square root over (Π_i=1ⁿ cosh r_i)} and we have used Stirling's formula, N!≈√{square root over (2πN)}(N/e)^N (more precisely, √{square root over (2πN)}(N/e)^Ne^1/12N+1<N!<√{square root over (2πN)}(N/e)^Ne^1/12N). Thus, if Z>e^n/2, the estimation error is smaller than e with exponentially small failure probability if we choose the number of samples as T=O(1/ϵ²).

For nonzero mean, one may again employ a similar formula for z˜N(μ, Σ):

$\begin{array}{cc}\mathrm{lhaf}\left({\sum }_n\right)={𝔼}_z\left[\prod _{i=1}^M{z}_i^{n_i}\right]& \left(107\right)\end{array}$ $\begin{array}{cc}=\sum _{v_1=0}^{n_1}\dots \sum _{v_M=0}^{n_M}{\left(-1\right)}^{{\sum }_{i=1}^M{v}_i}\left(\begin{array}{c}{n}_1\\ v_1\end{array}\right)\dots \left(\begin{array}{c}{n}_M\\ v_M\end{array}\right)\sum _{r=0}^{\left[n/2\right]}\frac{{\left(\frac{h^T\sum h}{2}\right)}^r{\left(h·\mu \right)}^{n-2r}}{r!\left(n-2r\right)!},& \left(108\right)\end{array}$

where n=n₁+ . . . +n_M and h=(n₁/2−v₁, . . . , n_M/2−v_M)^T and Σ_n is obtained by replacing the diagonal elements of Σ by μ and repeating ith row and column for n_i times. Again, let n_i=1 for all i's. In this case, M=n, and

$\begin{array}{cc}\mathrm{lhaf}(\sum )={𝔼}_z\left[\prod _{i=1}^n{z}_i\right]=\sum _{v_1=0}^1\dots \sum _{v_n=0}^1{\left(-1\right)}^{{\sum }_{i=1}^n{v}_i}\sum _{r=0}^{\left[n/2\right]}\frac{{\left(\frac{h^T\sum h}{2}\right)}^r{\left(h·\mu \right)}^{n-2r}}{r!\left(n-2r\right)!}& \left(109\right)\end{array}$ $\begin{array}{cc}=\frac{1}{2^n\left[n/2\right]}{𝔼}_{v,r}\left[{\left(-1\right)}^{{\sum }_{i=1}^n{v}_i}\frac{2^n\left[n/2\right]{\left(\frac{h^T\sum h}{2}\right)}^r{\left(h·\mu \right)}^{n-2r}}{r!\left(n-2r\right)!}\right].\mathrm{Here},& \left(110\right)\end{array}$ $\begin{array}{cc}❘\frac{2^n\left[n/2\right]{\left(\frac{h^T\sum h}{2}\right)}^r{\left(h·\mu \right)}^{n-2r}}{r!\left(n-2r\right)!}❘=\frac{\left[n/2\right]{n^{n/2}\left({\hat{h}}^T\sum \hat{h}\right)}^r{\left(\hat{h}·\mu \right)}^{n-2r}}{2^{r}r!\left(n-2r\right)!}\le \frac{\left[n/2\right]n^{n/2}{\Sigma }^r{\mu }^{n-2r}}{2^{r}r!\left(n-2r\right)!}& \left(111\right)\end{array}$

Therefore, the Chernoff bound provides

$\begin{array}{cc}\mathrm{Pr}\left[❘p-lh\mathrm{af}(\sum )❘>ϵ\frac{\left[n/2\right]n^{n/2}{\Sigma }^r{\mu }^{n-2r}}{2^{r}r!\left(n-2r\right)!}\right]\le 2\exp \left(-\frac{1}{2}{ϵ}^{2}T\right).& \left(112\right)\end{array}$

Especially when ∥Σ∥≤1, which is the case that we are interested in,

$\begin{array}{cc}\mathrm{Pr}\left[❘p-lh\mathrm{af}(\sum )❘>ϵ\frac{\left[n/2\right]n^{n/2}{\mu }^{n-2r}}{2^{r}r!\left(n-2r\right)!}\right]\le 2\exp \left(-\frac{1}{2}{ϵ}^{2}T\right).& \left(113\right)\end{array}$

The error bound clearly depends on the mean vector μ, which makes the bound more nontrivial than the hafnian case.

Other Techniques: Chernoff Bound

The accuracy of molecular vibronic spectra generation using Gaussian boson sampling may be derived. For example, let S={m_i}_i=1^M be the sample set, and for a fixed ω_vib, the random variable may be the indicator I(m, ω_vib) i.e., S→I(S, ω_vib)={I(m_i,ω_vib)}_i=1^N, where

$\begin{array}{cc}I\left(m,{\omega }_{vib}\right)=1\mathrm{if}{w}^{\prime }·m\in \left[{\omega }_{vib}-\frac{{\mathrm{\Delta \omega }}_{vib}}{2},{\omega }_{vib}+\frac{{\mathrm{\Delta \omega }}_{vib}}{2}\right]& \left(114\right)\end{array}$ $\begin{array}{cc}I\left(m,{\omega }_{vib}\right)=0\mathrm{otherwise}.& \left(115\right)\end{array}$

Thus, the Chernoff inequality for a given ω_vib may become

$\begin{array}{cc}\mathrm{Pr}\left[❘\mathrm{FCP}\left({\omega }_{vib}\right)-\frac{1}{N}\sum _{m\in 𝒮}I\left(m,{\omega }_{vib}\right)❘\ge ϵ\right]\le 2e^{-2N{ϵ}^2},& \left(116\right)\end{array}$

where Σ_m∈SI(m, ω_vib)/N is the estimate eFCP(ω_vib).

Then the Chernoff inequality may be expressed as:

Pr[|(X)−X|≥ϵ]≤2e^−2Nϵ2,  (117)

where X=(X₁+ . . . +X_N)/N is the sample average and X_i∈[0, 1].

Methods for Numerical Results Above

The parameters for FIGS. 3 and 6-7 are provided in J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean, and A. Aspuru-Guzik, Boson sampling for molecular vibronic spectra, Nature Photonics 9, 615 (2015), which is herein incorporated by reference.

For FIGS. 8 and 9 above, the first 6 elements are selected as ω′ as [1000, 500, 300, 200, 400, 250] and zeros are used for the remaining elements. In addition, the raw data of the recent Gaussian boson sampling is provided in https://quantum.ustc.edu.cn/web/node/951.

For the simulation, the frequencies in the unit of Δω_vib are approximated depending on the number of bins and domain. As explained above, increasing the resolution requires an additional cost for the fast Fourier transform, which is only poly(K).

The generalized P-representation of a squeezed state is given by

$\begin{array}{cc}P\left(\alpha ,\beta \right)=\frac{\sqrt{1+\gamma }}{\pi \gamma }{e}^{-\left({\alpha }^2+{\beta }^2\right)\left({\gamma }^{-1}+1/2\right)+\alpha \beta },& \left(118\right)\end{array}$

where α and β are real numbers, and γ=e^2r with r>0 being the squeezing parameter.

The various classical algorithm above may be performed by any type of computing system include a set of processors, memories, and other components including but not limited to various user interfaces and network units. At least one memory may computer instructions, when executed by the set of processors, causes the computing system to perform the various algorithms described above.

Finally, it is to be understood that the following claims are intended to cover all of the generic and specific features of the invention herein described and all statements of the scope of the invention which, as a matter of language, might be said to fall therebetween.

Claims

1. A method for estimating a transition spectra in a bosonic system having M bosonic modes, M being an integer equal to or larger than 2, the method comprising: determining a unitary operation modified from a transition operator associated with the bosonic system;processing each of N sampling bosonic states of the bosonic system in a representation space of the bosonic system to generate N sets of samples using at least the unitary operation, wherein each set of samples correspond to one of the N sampling bosonic states and are sampled from at least one variable in the representation space, N being a positive integer;generating Fourier components of the transition spectra based on the N sets of samples; andinverse-transforming the Fourier components to estimate the transition spectra.

2. The method of claim 1, wherein the unitary operation comprises a complex mixed position and momentum operator.

3. The method of claim 1, wherein each of the N sampling bosonic states comprises a non-Gaussian state of the tM bosonic modes.

4. The method of claim 1, wherein each of the N sampling bosonic states comprises a Gaussian state of the M bosonic modes.

5. The method of claim 1, wherein each of the N sampling bosonic states of the M bosonic modes is expressed in a positive P-representation in the representation space of the bosonic system for the processing of generating the corresponding set of samples in the representation space.

6. The method of claim 5, wherein the M bosonic modes comprise M molecular vibronic modes.

7. The method of claim 5, wherein the positive P-representation of the each of the N sampling bosonic states of the M bosonic modes and for the processing therein represents a quasi-probability distribution of the each of the N sampling bosonic states of the M bosonic modes.

8. The method of claim 7, wherein the phase space for the positive P-representation comprises 2M dimensions.

9. The method of claim 7, wherein processing the each of the N sampling bosonic states of the M bosonic modes in the positive P-representation to generate the corresponding set of samples in the representation space comprises: selecting M pairs of complex numbers for two representation space variables of the each of the N sampling bosonic states in the representation space;transforming the M pairs of complex numbers by applying a unitary matrix corresponding to the unitary operation to generate a transformed M pairs of complex numbers; andgenerating the corresponding set of samples from the transformed M pairs of complex numbers.

10. The method of claim 9, wherein the transition operator associated with the bosonic system is decomposed into the unitary operation, a dressing operation, and a displacement operation.

11. The method of claim 9, wherein at least one of the two phase space variables is associated with mean quantum numbers of M-mode coherent states of the bosonic system.

12. The method of claim 1, wherein: the transition operator associated with the bosonic system is decomposed into the unitary operation, a dressing operation, a displacement operation, and a two-mode squeezing operation;M comprises a even integer;the bosonic system comprises a molecular vibronic system; andthe M bosonic modes comprises M/2 vibronic modes of the molecular vibronic system and M/2 auxiliary modes.

13. A computing system for estimating a transition spectra in a bosonic system having M bosonic modes, M being an integer equal to or larger than 2, the system comprising one or more processors and a memory for storing computer instructions, the one or more processors, when executing the computer instructions, are configured to cause the computing system to: determine a unitary operation modified from a transition operator associated with the bosonic system;process each of N sampling bosonic states of the bosonic system in a representation space of the bosonic system to generate N sets of samples using at least the unitary operation, wherein each set of samples correspond to one of the N sampling bosonic states and are sampled from at least one variable in the representation space, N being a positive integer;analytically generate Fourier components the transition spectra based on the N sets of samples; andinverse-transform the Fourier components to estimate the transition spectra.

14. The computing system of claim 13, wherein the unitary operation comprises a complex mixed position and momentum operator.

15. The computing system of claim 13, wherein each of the N sampling bosonic states comprises a non-Gaussian state of the tM bosonic modes.

16. The computing system of claim 13, wherein each of the N sampling bosonic states comprises a Gaussian state of the M bosonic modes.

17. The computing system of claim 16, wherein each of the N sampling bosonic states of the M bosonic modes is expressed in a positive P-representation in the representation space of the bosonic system for the processing of generate the corresponding set of samples in the representation space.

18. The computing system of claim 17, wherein: the M bosonic modes comprise M molecular vibronic modes, and wherein the positive P-representation of the each of the N sampling bosonic states of the M bosonic modes and for the processing therein represents a quasi-probability distribution of the each of the N sampling bosonic states of the M bosonic modes; andthe phase space for the positive P-representation comprises 2M dimensions.

19. The computing system of claim 18, wherein to process the each of the N sampling bosonic states of the M bosonic modes in the positive P-representation to generate the corresponding set of samples in the representation space comprises: select M pairs of complex numbers for two representation space variables of the Gaussian state in the representation space;transform the M pairs of complex numbers by applying s unitary matrix corresponding to the unitary operation to generate a transformed M pairs of complex numbers; andgenerate the corresponding set of samples from the transformed M pairs of complex numbers.

20. The computing system of claim 19, wherein: the transition operator associated with the bosonic system is decomposed into the unitary operation, a dressing operation, and a displacement operation; andat least one of the two phase space variables is associated with mean quantum numbers of M-mode coherent states of the bosonic system.