Introduction

It is widely believed that quantum computers can perform computational tasks exponentially more efficient than classical computers, such as for simulating quantum many body systems1,2, or factoring large integer3. A series of remarkable studies in the last decades however also revealed that quantum computational algorithms, even those that invoke a lot of peculiar quantum effects such as entanglement, do not always prevail. The first important result in this direction is the Gottesman-Knill theorem4,5,6, which states that a class of quantum algorithms which start from a stabilizer state, followed by the operations of Clifford logic gates—which can generate a large amount of entanglement—, can be simulated efficiently on classical computers. This seminal result was further extended to different classes of quantum algorithms7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26. Identifying quantum algorithms that can be efficiently simulated classically, and characterizing the distinctive quantum features that are lacking in them, may offer insight to specify the elusive physical resources responsible for the quantum computational speedup27,28,29,30, and therefore are crucial to better understand the boundary and correspondence between the quantum and classical computations. In turn, deeper grasp of these fundamental problems may provide useful physical intuition to devise new quantum algorithms which outperform their classical counterparts, and to study the minimal requirement for nonuniversal but significantly easier to implement quantum computational schemes, which nevertheless exhibit quantum speedup31,32,33,34,35.

Clearly, the above basic questions in quantum computation are fundamentally linked to the longstanding foundational problem of quantum-classical correspondence and contrast. What are the deep physical concepts which uniquely define quantum mechanics relative to classical mechanics? Or, what is lacking in classical mechanics relative to quantum mechanics? Is it possible to obtain the latter by modifying the former, supplementing it with the necessary conceptual patches? To this end, it has been shown that a significant large set of phenomena traditionally seen as specifically quantum could in fact be explained within classical statistical models with some kinds of “epistemic (statistical) restrictions”36,37,38,39. Partly inspired by these remarkable results, recently, we have developed a novel phase space representation of quantum mechanics, showing that its opaque formalism in the complex Hilbert space can be phrased more intuitively as a classical statistical mechanics of ensemble of trajectories subjected to a specific epistemic restriction parameterized by a global-nonseparable variable fluctuating randomly on the order of Planck constant40,41. The epistemic restriction manifests directly the quantum uncertainty relation between noncommuting (quantum) positon and momentum operators acting on the abstract Hilbert space into a more intuitive statistical constraint on the allowed distributions over the (classical) phase space. It was further argued that the new phase space formalism can be interpreted as a calculus for estimation of momentum given information on the conjugate positions under the epistemic restriction42, respecting a plausible inferential-causality principle of estimation independence43,44.

Here, guided by the intuition offered by the epistemically restricted (ER) phase space representation, we devise a classical algorithm which efficiently computes the quantum expectation values arising in a class of continuous variable (CV) quantum computational circuits or CV quantum processes widely encountered in quantum optical settings45,46. The idea is to employ the ER phase space representation to transform the quantum circuits into unconventional classical stochastic processes and evaluate the resulting multidimensional integral using the Monte Carlo sampling technique. The classical algorithm thus goes along the spirit of those based on the quasiprobability representations47,48,49 reported in Refs.9,10,11,12,13. However, unlike the latter, our classical algorithm applies only when the final quantum observable after the Heisenberg evolution associated with the quantum circuits, is at most second order in momentum. The classical simulations based on the quasiprobability representations require that the quasiprobability distributions associated with the quantum circuits are nonnegative9,10,11,13; or, that the amount of negative value is sufficiently small to guarantee a fast convergence rate12. By contrast, our classical algorithm does not suffer from, and thus is not limited by, such problem of negative ‘probability’. The results therefore show that, for the specific class of computational circuits, negativity of quasiprobability is not a sufficient resource for quantum speedup. A similar suggestion based on different results and method is reported recently in Ref.26.

Results

Epistemically restricted phase space representation of quantum mechanics

We summarise the phase space representation of quantum mechanics based on the ER ensemble of trajectories proposed in Refs.40,41. Consider a system of N spatial degrees of freedom \(q=(q_1,\dots ,q_N)^{{\mathrm{T}}}\), arranged as N-elements vector in \({\mathbb {R}}^N\) where the superscript T denotes transposition, with the corresponding canonical conjugate momentum \(p=(p_1,\dots ,p_N)^{{\mathrm{T}}}\). In classical mechanics, for a system with a classical Hamiltonian H(q, p), the dynamics of the system is deterministic governed by the Hamilton–Jacobi equation, i.e.,

$$p_n(q;t)=\partial _{q_n}S_{{\mathrm{C}}},$$
(1a)
$$\begin{aligned} -\partial _tS_{\text {C}}(q;t)&=H(q,p), \end{aligned}$$
(1b)

\(n=1,\dots ,N\), where \(S_{{\mathrm{C}}}(q;t)\) is a real-valued function of positions q and time t, known as the Hamilton’s principal function. The Hamilton–Jacobi equation is formally equivalent to the Hamilton’s equation50. Moreover, the Hamilton–Jacobi equation offers a geometrical picture of an ensemble of trajectories in configuration space: it describes the dynamics of the whole ensemble of trajectories, all characterized by a single Hamilton’s principal function \(S_{{\mathrm{C}}}(q;t)\). Solving Eq. (1) in terms of \(S_{{\mathrm{C}}}(q;t)\), a single trajectory in configuration space is picked up by choosing a configuration q at any particular time t, and the momentum along the trajectory is obtained by computing Eq. (1a).

In classical statistical mechanics of ensemble of trajectories, given \(S_{{\mathrm{C}}}(q)\), the probability distribution of p conditioned on q thus reads as, noting Eq. (1a),

$$\begin{aligned} {\mathbb P}_{S_{{\mathrm{C}}}}(p|q)=\prod _{n=1}^{N}\delta \big (p_n-\partial _{q_n}S_{{\mathrm{C}}}(q)\big ), \end{aligned}$$
(2)

(trivial time dependence is implicit). Here and below, a subscript F in \({\mathbb P}_{F}\) indicates a dependence of the probability on a function F. For later convenience, we denote the probability distribution of q at time t with a specific notation as \(\rho (q;t)\). The phase space distribution is therefore given by \({\mathbb P}_{\{S_{{\mathrm{C}}},\rho \}}(p,q)={\mathbb P}_{S_{{\mathrm{C}}}}(p|q)\rho (q)\), so that the average of any physical quantity O(q, p) over the phase space distribution can be written as

$$\begin{aligned} \langle O\rangle _{\{S_{{\mathrm{C}}},\rho \}}\doteq & {} \int \mathrm{d}^Nq~\mathrm{d}^Np~O(q,p){\mathbb P}_{S_{{\mathrm{C}}}}(p|q)\rho (q)\nonumber \\= & {} \int \mathrm{d}^Nq~O(q,\partial _qS_{{\mathrm{C}}})\rho (q), \end{aligned}$$
(3)

where \(\mathrm{d}^{N}q\doteq \mathrm{d}q_1\dots \mathrm{d}q_N\), \(\mathrm{d}^{N}p\doteq \mathrm{d}p_1\dots \mathrm{d}p_N\), and we have used Eq. (2) to get the second line.

Note that since the classical dynamics governed by Eq. (1) conserves the energy along each single trajectory, it automatically conserves the average energy over the whole ensemble of trajectories, i.e., \((\mathrm{d}/\mathrm{d}t)\langle H\rangle _{\{S_{{\mathrm{C}}},\rho \}}=0\). Conversely, it has been shown in Ref.40 that assuming a momentum field which can be written as in Eq. (1a), and imposing the conservation of average energy and conservation of trajectories (probability current) manifested by the continuity equation \(\partial _t\rho +\partial _q\cdot (\rho \dot{q})=0\) where \(\dot{q}\doteq \mathrm{d}q/\mathrm{d}t\), singles out the Hamilton–Jacobi equation of Eq. (1b).

Next, suppose that solving the Hamilton–Jacobi equation, or more straightforwardly solving the equivalent Hamilton’s equation, the initial position and momentum are mapped at some later time as: \(q\mapsto q'=f_{{\mathrm{C}}}(q,p)\) and \(p\mapsto p'=g_{{\mathrm{C}}}(q,p)\), where \(f_{{\mathrm{C}}}=(f_{\mathrm{C}_1},\dots ,f_{\mathrm{C}_N})^{{\mathrm{T}}}\) and \(g_{{\mathrm{C}}}=(g_{\mathrm{C}_1},\dots ,g_{\mathrm{C}_N})^{{\mathrm{T}}}\) are vectors of functions on phase space. This deterministic mapping can be expressed in terms of transition probability as \({\mathbb P}_{\{f_{{\mathrm{C}}},g_{{\mathrm{C}}}\}}(p',q'|p,q)=\prod _{n=1}^N\delta \big (q'_n-f_{\mathrm{C}_n}(q,p)\big )\prod _{m=1}^N\delta \big (p'_m-g_{\mathrm{C}_m}(q,p)\big )\). The statistical average of a physical quantity O(q, p) over the phase space distribution after the time evolution can thus be computed as

$$\begin{aligned} \langle O\rangle _{\{(f_{{\mathrm{C}}},g_{{\mathrm{C}}});(S_{{\mathrm{C}}},\rho )\}}\doteq & {} \int \mathrm{d}^Nq'~\mathrm{d}^Np'~\mathrm{d}^Nq~\mathrm{d}^Np~O(q',p'){\mathbb P}_{\{f_{{\mathrm{C}}},g_{{\mathrm{C}}}\}}(p',q'|p,q){\mathbb P}_{\{S_{{\mathrm{C}}},\rho \}}(p,q)\nonumber \\= & {} \int \mathrm{d}^Nq'~\mathrm{d}^Np'~\mathrm{d}^Nq~\mathrm{d}^Np~O(q',p')\nonumber \\&\times \prod _{n=1}^N\delta \big (q'_n-f_{\mathrm{C}_n}(q,p)\big )\prod _{m=1}^N\delta \big (p'_m-g_{\mathrm{C}_m}(q,p)\big )\prod _{k=1}^N\delta (p_k-\partial _{q_k}S_{{\mathrm{C}}}(q))\rho (q). \end{aligned}$$
(4)

We have argued in Refs.40,41 that the abstract formulas of quantum mechanics can be expressed as a specific modification of the above classical statistical mechanics of ensemble of trajectories in phase space. First, we introduce an “ontic extension” in the form of a global-nonseparable variable \(\xi\): it is real-valued with the dimension of action and depends only on time (i.e., spatially uniform). We assume that \(\xi\) fluctuates randomly on a microscopic time scale with a probability density \(\chi (\xi )\) such that the first and second moments are constant in time, given by

$$\begin{aligned} \overline{\xi }\doteq \int \mathrm{d}\xi ~\xi ~\chi (\xi )=0,~~~\overline{\xi ^2}=\hbar ^2. \end{aligned}$$
(5)

We then impose a specific epistemic restriction on the class of phase-space probability distributions that Nature allows us to prepare, as follows. Consider an ensemble of identical preparations (characterized by the same set of controllable macroscopic parameters) resulting in an ensemble of trajectories in configuration space following a momentum field. In conventional classical statistical mechanics, from Eq. (1a), it is possible to prepare an ensemble of trajectories with a given targeted probability distribution of positions \(\rho (q)\) using an arbitrary form of momentum field p(q) by choosing an arbitrary \(S_{{\mathrm{C}}}(q)\). This means that, in classical mechanics, as shown in Eq. (2), the probability distribution of p conditioned on q is independent of \(\rho (q)\), i.e., \({\mathbb P}_{\{S_{{\mathrm{C}}},\rho \}}(p|q)={\mathbb P}_{S_{{\mathrm{C}}}}(p|q)\). Or, equivalently, a given momentum field can be used to prepare an ensemble of trajectories with an arbitrary targeted \(\rho (q)\), i.e., each trajectory in the momentum field can be weighted with an arbitrary choice of \(\rho (q)\).

Let us assume that such an “epistemic (statistical) freedom” in conventional classical statistical mechanics is no longer fully granted in microscopic world. Hence, we assume that in microscopic regime, each trajectory in a given momentum field can no longer be assigned an arbitrary weight \(\rho (q)\), i.e., the allowed distribution of positions is fundamentally restricted by the underlying momentum field42,43. Or, equivalently, given a targeted \(\rho (q)\), it is no longer possible to prepare an ensemble of trajectories realizing \(\rho (q)\) following an arbitrary momentum field. This means that the probability distribution of p conditioned on q must depend on the targeted \(\rho (q)\), i.e.: \({\mathbb P}_{\{\dots ,\rho \}}(p|q,\dots )\ne {\mathbb P}_{\{\dots \}}(p|q,\dots )\). Let us then consider a statistical model with such an epistemic restriction so that an ensemble of identical preparations yields a probability distribution of p conditioned on q which fundamentally depends on \(\rho (q)\) as follows (Cf. Eq. (2))40:

$$\begin{aligned} {\mathbb P}_{\{S ,\rho \}}(p|q,\xi )=\prod _{n=1}^N\delta \Big (p_n-\Big (\partial _{q_n}S +\frac{\xi }{2}\frac{\partial _{q_n}\rho }{\rho }\Big )\Big ), \end{aligned}$$
(6)

parameterized by the global random variable \(\xi\) satisfying Eq. (5). Here, S(q; t) is a real-valued scalar function of (q; t) with the dimension of action, replacing the role of \(S_{{\mathrm{C}}}(q;t)\) in Eq. (2). From Eq. (6), the allowed (conditional) phase space distributions are thus also restricted to have the following specific form:

$$\begin{aligned} {\mathbb P}_{\{S,\rho \}}(p,q|\xi )= & {} {\mathbb P}_{\{S,\rho \}}(p|q,\xi )\rho (q)\nonumber \\= & {} \prod _{n=1}^N\delta \Big (p_n-\Big (\partial _{q_n}S +\frac{\xi }{2}\frac{\partial _{q_n}\rho }{\rho }\Big )\Big )\rho (q). \end{aligned}$$
(7)

One can straightforwardly show that Eqs. (5) and (6) imply the Heisenberg–Kennard uncertianty relation, i.e., \(\sigma _{q_n}\sigma _{p_n}\ge \hbar /2\), \(n=1,\dots ,N\), where \(\sigma _{q_n}\) and \(\sigma _{p_n}\) are respectively the standard deviations of position and momentum of the nth spatial degree of freedom40. In this sense, the epistemic restriction constrains the allowed phase space distributions to satisfy the uncertainty relation. For this reason, we refer to \({\mathbb P}_{\{S,\rho \}}(p,q|\xi )\) defined in Eq. (7) as the “ER (epistemically restricted) phase space distribution” associated with a pair of real-valued functions \((S,\rho )\). We have recently shown in Refs.41,42 that the ER phase space distribution is not just a mathematical artefact, but can be indirectly reconstructed in experiment via the notion of weak momentum value defined as \(\frac{\langle q|\hat{p}|\psi \rangle }{\langle q|\psi \rangle }\)51,52,53. This, in turn, provides an interpretation of the complex weak momentum value in terms of the ER momentum fluctuation and quantum uncertainty. The classical limit of macroscopic physical regime is obtained when \(|\partial _qS |\gg |\frac{\xi }{2}\frac{\partial _q\rho }{\rho }|\), e.g., \(|\xi |\rightarrow 0\) or equivalently \(\hbar \rightarrow 0\) as per Eq. (5), and \(S \rightarrow S_{{\mathrm{C}}}\), so that the conditional distribution of momentum of Eq. (6) reduces back to that in conventional classical statistical mechanics given by Eq. (2), and the ER phase space distribution of Eq. (7) rolls back to the classical phase space distribution as \({\mathbb P}_{\{S,\rho \}}(p,q|\xi )\rightarrow {\mathbb P}_{\{S_{{\mathrm{C}}},\rho \}}(p,q)=\prod _{n=1}^N\delta (p_n-\partial _{q_n}S_{{\mathrm{C}}})\rho (q)\), lacking an epistemic restriction.

Assume further that the physical quantities in the above ER classical model are given by some real-valued functions of positions and momentum, O(q, p), having the same form as those in classical mechanics. We then obtain the following theorem expressing the equivalence between the conventional statistical average over the ER phase space distribution and the quantum expectation value.

Theorem 1

(Budiyono–Rohrlich40) Consider an ensemble of trajectories satisfying the epistemic restriction of Eqs. (6) and (5). The ensemble average of any classical physical quantity O(q, p) up to second order in momentum p over the ER phase space distribution of Eq. (7) is equal to the quantum expectation value as

$$\begin{aligned} \langle O\rangle _{\{S ,\rho \}}\doteq & {} \int \mathrm{d}^Nq~\mathrm{d}\xi ~\mathrm{d}^Np~O(q,p){\mathbb P}_{\{S ,\rho \}}(p,q|\xi )\chi (\xi )\nonumber \\= & {} \int \mathrm{d}^Nq~\mathrm{d}\xi ~\mathrm{d}^Np~O(q,p)\prod _{n=1}^N\delta \Big (p_n-\Big (\partial _{q_n}S +\frac{\xi }{2}\frac{\partial _{q_n}\rho }{\rho }\Big )\Big )\rho (q)\chi (\xi )\nonumber \\= & {} \langle \psi |\hat{O}|\psi \rangle , \end{aligned}$$
(8)

where \(\hat{O}\)is a Hermitian operator taking exactly the same form as that obtained by applying the Dirac canonical quantization procedure to O(q, p) with a specific ordering of operators, and the wave function \(\psi (q;t)\doteq \langle q|\psi (t)\rangle\)is defined as

$$\begin{aligned} \psi (q;t)\doteq \sqrt{\rho (q;t)}\exp (iS(q;t)/\hbar ). \end{aligned}$$
(9)

See the Methods section in Ref.40 for the proof. We note that when the physical quantity O(q, p) has cross terms between momentum of different degrees of freedom, i.e., \(p_np_m\), \(n\ne m\), the nonseparibility of \(\xi\) is indeed indispensable for obtaining Eq. (8).

A couple of important notes are in order here. Theorem 1 is developed in such a way that we can devise a classical statistical model so that the quantum expectation value is ‘reconstructed’ from the conventional statistical average over the ER phase space distribution of Eq. (7). We can however read Theorem 1 the other way around. Namely, the ER classical model can be seen as a ‘representation’ of quantum statistics in the ER classical phase space41. In this reading, first we are given a quantum wave function \(\psi (q)\) describing the preparation. We compute its amplitude \(\rho (q)\) and phase S(q) as

$$\begin{aligned} \rho (q)=\big |\psi (q)\big |^2~~~ \& ~~~S (q)=\frac{\hbar }{2i}\big (\log \psi (q)-\log \psi (q)^*\big ), \end{aligned}$$
(10)

complying with Eq. (9), based on which, we define the ER phase space distribution of Eq. (7). We then use the ER phase space distribution to transform the quantum expectation value into a conventional statistical phase space average as in Eq. (8)—by reading it from the right-hand side to the left.

Reading Eq. (8) as an epistemically restricted phase space representation of quantum expectation values, let us briefly discuss the apparent asymmetric treatment of q and p in the Theorem 1. Namely, while there is no restriction on the order of q, Theorem 1 only applies when p is at most second order. This asymmetry seems in particular peculiar from the perspective of quantum optics, wherein the quadrature phase space operators \((\hat{q},\hat{p})\) represent the field amplitudes oscillating out of phase with each other, so that they can switch positions. First, we emphasize that, as in quantum optics, the labelling of q and p in Theorem 1 is a matter of convention; namely, we can interchange the roles of the symbol q and p by working in the p (i.e., momentum) representation instead of in the q (i.e., position) representation. For a concrete example, consider the computation of the following expectation value \(\langle \psi |(\hat{q}^2+\hat{p}^n)|\psi \rangle\), where n is an integer larger than two. In this case, it is natural to work in the momentum representation and write the associated wave function in the polar form as \(\phi (p)\doteq \langle p|\psi \rangle =\sqrt{\rho (p)}e^{iS(p)/\hbar }\). One can then show that Eq. (8) in Theorem 1 still applies with the transformation \(q\leftrightarrow p\). However, in this case, \(\hat{q}\) must be at most second order. Hence, either way, within the ER phase space representation, one of the two quadrature phase space operators must be at most second order. Mathematically, the asymmetric treatment of q and p in Theorem 1 can be traced to their asymmetric roles in the epistemic restriction based on which the phase space representation is constructed. Below we shall take the convention that p plays the role of momentum which is limited to be at most second order.

As an important example of the application of Theorem 1, we have:

$$\begin{aligned} \langle (\eta _m-\langle \eta _m\rangle _{\{S,\rho \}})(\eta _n-\langle \eta _n\rangle _{\{S,\rho \}})\rangle _{\{S ,\rho \}}= & {} \langle \eta _m\eta _n\rangle _{\{S ,\rho \}}-\langle \eta _m\rangle _{\{S ,\rho \}}\langle \eta _n\rangle _{\{S ,\rho \}}\nonumber \\= & {} \langle \psi |\frac{1}{2}(\hat{\eta }_m\hat{\eta }_n+\hat{\eta }_n\hat{\eta }_m)|\psi \rangle -\langle \psi |\hat{\eta }_m|\psi \rangle \langle \psi |\hat{\eta }_n|\psi \rangle , \end{aligned}$$
(11)

where \(\eta =(p,q)\), \(m,n=1,\dots ,N\). Notice that the right-hand side is precisely the quantum covariance matrix of the phase space (quadrature) operators. Hence, reading Eq. (11) from the right-hand side to the left, the covariance matrix of position and momentum operators in quantum mechanics, which plays prominent roles in quantum optics and also in general CV quantum information processing, can be expressed as the ordinary classical statistical covariance matrix of the position and momentum over the ER phase space distribution.

It can be further shown that imposing the conservation of average energy to the ensemble of trajectories satisfying the epistemic restriction of Eqs. (5) and (6), i.e., \((\mathrm{d}/\mathrm{d}t)\langle H\rangle _{\{S,\rho \}}=0\), and also the conservation of trajectories (probability current), single out a unique dynamics for the wave function defined in Eq. (9) given by the unitary Schrödinger equation: \(i\hbar (\mathrm{d}/\mathrm{d}t)|\psi \rangle =\hat{H}|\psi \rangle\)40. \(\hat{H}\) is the quantum Hamiltonian having the same form as that obtained by applying the Dirac canonical quantization procedure to H(q, p) with a specific ordering of operators. In the macroscopic regime of classical limit so that Eq. (6) becomes Eq. (2), Eq. (8) gives back the conventional classical average of Eq. (3), and the Schrödinger equation reduces to the classical Hamilton–Jacobi equation of Eq. (1)40,42.

Computing the expectation values in a class of quantum circuits with an epistemically restricted ensemble of trajectories

We extend Theorem 1 to include a specific class of mappings of position and momentum variables to devise classical algorithms which efficiently compute the expectation values arising in a nontrivial class of CV quantum computational circuits or CV quantum processes on classical probabilistic computers.

Consider the following stochastic dynamics of ensemble of trajectories. Suppose that at an initial time we are given a pair of real-valued functions \((S(q),\rho (q))\), where \(\rho (q)\) is a normalized probability density of q. First, we draw a sample of the pair of variables \((q,\xi )\) from the joint probability density \(\rho (q)\chi (\xi )\), where \(\chi (\xi )\) satisfies Eq. (5). We then use \((q,\xi )\) to compute the value of p at this initial time as

$$\begin{aligned} p_n=\partial _{q_n}S(q)+\frac{\xi }{2}\frac{\partial _{q_n}\rho (q)}{\rho (q)}, \end{aligned}$$
(12)

\(n=1,\dots ,N\). Hence, in this way, we initially sample the phase space variables (p, q) from the joint probability distribution of Eq. (7), which is just the ER phase space distribution associated with a wave function \(\psi (q)=\sqrt{\rho (q)}\exp (iS(q)/\hbar )\), summarized in the previous section.

Next, assume that the phase space variables (p, q) evolves in time into a new phase space variable \((p',q')\) as

$$\begin{aligned} q\mapsto q'= & {} f(q,p),\nonumber \\ p\mapsto p'= & {} g(q,p), \end{aligned}$$
(13)

where \(f=(f_1,\dots ,f_N)^{{\mathrm{T}}}\) and \(g=(g_1,\dots ,g_N)^{{\mathrm{T}}}\) are vectors of functions on phase space, independent of \(\xi\). Here, we assume for simplicity that \(\xi\) is kept fixed during the mapping of phase space variables (i.e., during the time evolution). The above deterministic mapping induces the following transition probability over the phase space variables:

$$\begin{aligned} {\mathbb {P}}_{\{f,g\}}(p',q'|p,q)= & {} \prod _{n=1}^N\delta \big (q'_n-f_n(q,p)\big )\prod _{m=1}^N\delta \big (p'_m-g_m(q,p)\big ). \end{aligned}$$
(14)

At the end of the phase space transformation, we wish to compute the statistical phase space average of a physical quantity \(O(q',p')\). We therefore have to evaluate the following multidimensional integral:

$$\begin{aligned} \langle O\rangle _{\{(f,g);(S ,\rho )\}}\doteq & {} \int \mathrm{d}^Nq'~ \mathrm{d}^Np'~\mathrm{d}^Nq~\mathrm{d}\xi ~\mathrm{d}^Np~O(q',p'){\mathbb P}_{\{f,g\}}( p',q'|p,q){\mathbb P}_{\{S,\rho \}}(p,q|\xi )\chi (\xi )\nonumber \\= & {} \int \mathrm{d}^Nq'~\mathrm{d}^Np'~\mathrm{d}^Nq~\mathrm{d}\xi ~ \mathrm{d}^Np~O(q',p')\prod _{n=1}^N\delta \big (q'_n-f_n(q,p)\big ) \prod _{m=1}^N\delta \big (p'_m-g_m(q,p)\big )\nonumber \\&\times \prod _{k=1}^N\delta \Big (p_k-\Big (\partial _{q_k}S +\frac{\xi }{2}\frac{\partial _{q_k}\rho }{\rho }\Big )\Big )\rho (q)\chi (\xi ), \end{aligned}$$
(15)

where we have used Eqs. (7) and (14) to get the second equality. Note that the above computational scheme for average value can be seen as a classical stochastic process for 2N positions and momentum random variables, but with an initial phase space that is epistemically (statistically) restricted being sampled from the specific ER phase space distribution given by Eq. (7). Hence, assuming that the time evolution of the position and momentum variables described in Eq. (13) is efficiently tractable on a classical computer, and provided we can sample \(\xi\) from \(\chi (\xi )\), and (p, q) from \({\mathbb P}_{\{S,\rho \}}(p,q|\xi )\) of Eq. (7) efficiently, the above computational task can be run on a classical computer efficiently using Monte Carlo sampling method54. We show below that for a specific class of classical quantities O(q, p), and phase space mappings (f, g), Eq. (15) can be used to efficiently compute the expectation values arising in a wide important class of quantum circuits.

First, consider the case when f and g in Eq. (15) are linear in position and momentum variables; namely, we assume the following linear mapping of position and momentum variables:

$$\begin{aligned} q'= & {} f(q,p)=A\cdot q+B\cdot p+ I\cdot q_0,\nonumber \\ p'= & {} g(q,p)=C\cdot q+D\cdot p+ I\cdot p_0, \end{aligned}$$
(16)

where A, B, C and D are \(N\times N\)-matrices, I is identity matrix, \(q_0\) and \(p_0\) are N-elements column vectors, and \((A\cdot q)_m\doteq \sum _{n=1}^NA_{mn}q_n\), et cetera. Inserting Eq. (16) into Eq. (15), and evaluating the integration over \((p',q')\), we obtain

$$\begin{aligned} \langle O\rangle _{\{(f,g);(S ,\rho )\}}= & {} \int \mathrm{d}^Nq~\mathrm{d}\xi ~\mathrm{d}^Np~O'(q,p)\prod _{k=1}^N\delta \Big (p_k-\Big (\partial _{q_k}S +\frac{\xi }{2}\frac{\partial _{q_k}\rho }{\rho }\Big )\Big )\rho (q)\chi (\xi ). \end{aligned}$$
(17)

where \(O'(q,p)\doteq O\big (A\cdot q+B\cdot p+ I\cdot q_0,C\cdot q+D\cdot p+ I\cdot p_0\big )\).

Now, let us assume that \(O(q',p')\) in Eq. (15) is at most quadratic in the phase space variables \((p',q')\), i.e., it may contain terms like \(p'_mp'_n\), \(q'_mq'_n\), or \(p'_mq'_n\), but it does not contain cubic or higher order terms like \(q_mq_nq_l\), \(m,n,l=1,\dots ,N\). In this case, imposing the linear transformation of Eq. (16), the resulting transformed \(O'(q,p)\) in Eq. (17) must also be at most second order in (p, q). Importantly, \(O'(q,p)\) is therefore at most second order in p. This fact permits the application of Theorem 1 to show that the statistical average over the ER phase space distribution of Eq. (17) is equivalent to the quantum expectation value as

$$\begin{aligned} \langle O\rangle _{\{(f,g);(S ,\rho )\}}=\langle \psi |O'(\hat{q},\hat{p})|\psi \rangle . \end{aligned}$$
(18)

Here, the wave function \(\psi\) is defined as in Eq. (9), and the quantum observable is given by \(O'(\hat{q},\hat{p},)\doteq O\big (A\cdot \hat{q}+B\cdot \hat{p}+\hat{I}\cdot q_0,C\cdot \hat{q}+D\cdot \hat{p}+\hat{I}\cdot p_0\big )\). Hence, \(O'(\hat{q},\hat{p},)\) is obtained from \(O(\hat{q},\hat{p})\) via the following linear Affine transformation of the position and momentum operators:

$$\begin{aligned} \hat{q}'= & {} f(\hat{q},\hat{p})=A\cdot \hat{q}+B\cdot \hat{p}+\hat{I}\cdot q_0,\nonumber \\ \hat{p}'= & {} g(\hat{q},\hat{p})=C\cdot \hat{q}+D\cdot \hat{p}+\hat{I}\cdot p_0. \end{aligned}$$
(19)

Suppose further that the linear transformation of the position and momentum variables of Eq. (16) is symplectic55; namely, it preserves the canonical Poisson bracket relations: \([q'_n,p'_m]_{{\mathrm{PB}}}=[q_n,p_m]_{{\mathrm{PB}}}=\delta _{nm}\), \(n,m=1,\dots ,N\). Then, the corresponding Affine transformation of the position and momentum operators of Eq. (19) conserves the canonical commutation relations, i.e., \([\hat{q}_n',\hat{p}_m']=[\hat{q}_n,\hat{p}_m]=i\hbar \delta _{nm}\), \(n,m=1,\dots ,N\), so that the mapping \(O(\hat{q},\hat{p})\mapsto O'(\hat{q},\hat{p})\) can be implemented by a unitary transformation \(\hat{U}_{\{f,g\}}\) as55,56

$$\begin{aligned} O'(\hat{q},\hat{p})=\hat{U}^{\dagger }_{\{f,g\}}O(\hat{q},\hat{p})\hat{U}_{\{f,g\}}. \end{aligned}$$
(20)

Here \(\hat{U}_{\{f,g\}}\doteq \exp (-i\hat{H}\theta /\hbar )\), where \(\hat{H}\)—the quantum Hamiltonian—is the generator of the unitary transformation taking the form of a Hermitian operator at most quadratic in position and momentum operators, and \(\theta\) is a parameter with the dimension of time. Such unitaries are called Gaussian unitaries since they map Gaussian states onto Gaussian states: they only change the means and the covariances of the initial Gaussian states. The set of all Gaussian unitaries \(\hat{U}_{\{f,g\}}\) comprises a Clifford group and plays crucial roles in quantum optics and in general CV quantum information processing7,8,56,57.

Noting Eq. (20), Eq. (18) can thus be written as

$$\begin{aligned} \langle O\rangle _{\{(f,g);(S ,\rho )\}}=\langle \psi |\hat{U}_{\{f,g\}}^{\dagger }O(\hat{q},\hat{p})\hat{U}_{\{f,g\}}|\psi \rangle . \end{aligned}$$
(21)

Reading Eq. (21) from the right-hand side to the left, we therefore have the following theorem.

Theorem 2

The expectation values of quantum observables up to second order in position and momentum operators arising in any quantum circuits which start from arbitrary quantum states and then acted upon by unitary quantum gates generated by quadratic quantum Hamiltonians (Gaussian unitaries), can be computed efficiently on classical probabilistic computers using Monte Carlo sampling over the ER ensemble of trajectories.

We have shown that the classical algorithm above can compute the expectation value of quantum observables at most quadratic in position and momentum operators. It thus does not in general simulate the full quantum state. However, for the specific yet important class of quantum circuits which initiate from an arbitrary Gaussian state and evolve under an arbitrary Gaussian unitary inducing linear Affine transformation (in the Heisenberg picture) of the type of Eq. (19), our classical algorithm can be used to simulate the final quantum state. To see this, first recall that the initial Gaussian quantum states are transformed by the Gaussian unitaries into Gaussian states. Furthermore, noting that Gaussian states are completely determined by the means and covariances of the position and momentum operators56,57, to know the final Gaussian states, we only need to compute their means and covariances, i.e., \(\langle \psi _{{\mathrm{G}}}|\hat{U}_{{\mathrm{G}}}^{\dagger }\hat{\eta }_k\hat{U}_{{\mathrm{G}}}|\psi _{{\mathrm{G}}}\rangle\) and \(\langle \psi _{{\mathrm{G}}}|\hat{U}_{{\mathrm{G}}}^{\dagger }\frac{1}{2}(\hat{\eta }_k\hat{\eta }_l+\hat{\eta }_l\hat{\eta }_k)\hat{U}_{{\mathrm{G}}}|\psi _{{\mathrm{G}}}\rangle\), \(\eta =p,q\), and \(k,l=1,\dots ,N\), where \(|\psi _{{\mathrm{G}}}\rangle\) is the initial Gaussian state and \(\hat{U}_{{\mathrm{G}}}\) is the Gaussian unitary. The above quantum circuits for computing the quantum averages and quantum covariances of the phase space quadrature operators obviously fall into the scope of Theorem 2 so that they can be efficiently computed using our classical sampling algorithm. Moreover, the means and covariances of the phase space operators over the final Gaussian states can be further used to compute the probability of measurement outcome with respect to the Gaussian POVM.

We have thus the following corollary of Theorem 2.

Corollary 1

Quantum algorithms that start from the preparations of Gaussian states, followed by the operations of Gaussian unitaries and measurements over Gaussian POVM, can be simulated efficiently on classical probabilistic computers using Monte Carlo sampling over the ER ensemble of trajectories.

Notice that the classical simulatability of Gaussian quantum computations in Corollary 1 is just the CV version of the Gottesman-Knill theorem reported in Refs.7,8. Moreover, note that the CV Gottesman-Knill theorem does not apply for the computation of the expectation value when the initial quantum state is not Gaussian. In this regards, our classical algorithm summarized in Theorem 2 therefore extends the CV Gottesman-Knill theorem to include the computations of the expectation values of observables quadratic in position and momentum operators with nonGaussian initial quantum states.

Theorem 2 can also be extended to cover a class of nonlinear transformations of phase space operators induced by nonGaussian quantum operations as follows. We first note that to prove Theorem 2, we have used Theorem 1 which only requires that \(O(\hat{q},\hat{p})\) in Eq. (8) is at most second order in \(\hat{p}\); apart from that, it may contain terms with arbitrary degrees of \(\hat{q}\), e.g., \(\hat{q}^4\), \(\hat{p}\hat{q}^3\hat{p}\), et cetera. This means that we can devise an efficient classical algorithm which computes the quantum expectation values arising in quantum circuits starting from an arbitrary initial quantum state followed by the application of any quantum gate as long as the final quantum observable—after the Heisenberg evolution induced by the quantum circuits—does not have term third order or higher in momentum operators. For example, consider the following quantum circuit for computing the quantum expectation value \(\langle \psi |\hat{U}^{\dagger }_{\{f,g\}}O(\hat{q},\hat{p})\hat{U}_{\{f,g\}}|\psi \rangle\), where \(O(\hat{q},\hat{p})=\hat{p}^2+\hat{q}^4\) and \(\hat{U}_{\{f,g\}}\) is a nonGaussian quantum gate inducing nonlinear transformations \(\hat{p}\mapsto \hat{p}+\hat{q}^2\) and \(\hat{q}\mapsto \hat{q}\). One can see that under such transformation, the quantum observable after the Heisenberg evolution is maintained to be second order in \(\hat{p}\): \(O(\hat{q},\hat{p})\mapsto O'(\hat{q},\hat{p})=\hat{p}^2+\hat{p}\hat{q}^2+\hat{q}^2\hat{p}+2\hat{q}^4\). Hence, we can still proceed as before to efficiently compute the above quantum expectation value using the classical algorithm of Eq. (15) and evaluate the integral using Monte Carlo sampling over the associated ER phase space distribution. We therefore obtain the following theorem which is more general than Theorem 2.

Theorem 3

The quantum expectation values arising in any quantum circuits which start from arbitrary quantum states followed by the applications of unitary quantum gates inducing in general nonlinear transformations (in the Heisenberg picture) and yielding quantum observables at most quadratic in momentum operators, can be classically computed efficiently using Monte Carlo sampling over the ER ensemble of trajectories.

Let us summarize, before proceeding, the important ingredients for the classical algorithm based on the ER phase space representation of Eq. (15). First, given an initial pure quantum state \(|\psi \rangle\), we express it in the configuration representation to get the wave function \(\psi (q)=\langle q|\psi \rangle\), and compute its phase S(q) and amplitude \(\rho (q)\) as in Eq. (10). Hence, we assume that the wave function \(\psi (q)\), its phase S(q), and amplitude \(\rho (q)\), can be computed efficiently, and \(\rho (q)\) can be sampled efficiently. These requirements give a restriction on the form of the initial quantum states. For discrete variables case, such quantum states are called computationally trackable states in Ref.54. Moreover, we also assume that the computation of p in Eq. (12), and also the phase space transformation in Eq. (13), can be carried out efficiently on classical computers.

In the macroscopic classical limit, noting that the mapping of Eq. (13) is given by the classical Hamilton’s equation, and Eq. (7) reduces back to the phase space distribution of conventional classical statistical mechanics, then Eq. (15) gives back the computation of the average value in conventional classical statistical mechanics of Eq. (4). The epistemic restriction of Eq. (6) is thus playing a crucial role in the quantum circuits for computing expectation values. Namely, first, the epistemic restriction is a necessary ingredient for any form of quantum algorithm to compute the quantum expectation value without which we regain the genuine conventional classical algorithm. In other words, any such quantum algorithm (covertly) operates a specific epistemic restriction encoded in the Hilbert space in the form of noncommutative structure of the phase space operators. Moreover, for a specific class of quantum circuits computing the expectation values, we can disclose the epistemic restriction and exploit it to devise classical algorithms which efficiently compute the quantum expectation values.

Next, the transition probability in the classical algorithm of Eq. (15) can be sequentially concatenated, before computing the phase space average as follows:

$$\begin{aligned} \langle O\rangle _{\{(f_T,g_T);(f_{T-1},g_{T-1});\dots ;(f_1,g_1);(S ,\rho )\}}\doteq \int \mathrm{d}^Nq_T~\mathrm{d}^Np_T\dots \mathrm{d}^Nq_1~\mathrm{d}^Np_1~\mathrm{d}^Nq_0~\mathrm{d}\xi ~\mathrm{d}^Np_0~O(p_T,q_T)\nonumber \\ \times {\mathbb P}_{\{f_T,g_T\}}(p_T,q_T|p_{T-1},q_{T-1})\times \dots \times {\mathbb P}_{\{f_1,g_1\}}(p_1,q_1|p_0,q_0){\mathbb P}_{\{S,\rho \}}(p_0,q_0|\xi )\chi (\xi ).~~~ \end{aligned}$$
(22)

Assume that each transition probability takes the form of Eq. (14) in which all pairs of transformation \((f_t,g_t)\) has unitary implementation \(\hat{U}_t\), \(t=1,\dots ,T\). As long as the total transformation (in the Heisenberg picture) does not yield cubic or higher order term of momentum operator, one can proceed as before to show that the above ER classical stochastic process can be used to efficiently compute the quantum expectation value arising in the quantum circuit, i.e., we have

$$\begin{aligned} \langle O\rangle _{\{(f_T,g_T);(f_{T-1},g_{T-1});\dots ;(f_1,g_1);(S ,\rho )\}}=\langle \psi |\hat{U}_1^{\dagger }\dots \hat{U}_{T-1}^{\dagger }\hat{U}_T^{\dagger }\hat{O}\hat{U}_T\hat{U}_{T-1}\dots \hat{U}_1|\psi \rangle . \end{aligned}$$
(23)

We note an interesting and important point that this result still applies even if the initial quantum observable \(\hat{O}\) in Eq. (23) contains cubic or higher order terms of \(\hat{p}\), and/or some of the quantum gates \(\hat{U}_t\), \(t=1,\dots ,T\) generate cubic or higher order terms of \(\hat{p}\), as long as these terms are cancelled at the end of the processes so that the final observable \(\hat{O}'=\hat{U}_1^{\dagger }\dots \hat{U}_{T-1}^{\dagger }\hat{U}_T^{\dagger }\hat{O}\hat{U}_T\hat{U}_{T-1}\dots \hat{U}_1\) does not contain such term. This is comparable to the result obtained based on the \(s-\)ordered quasiprobability representation13, wherein any negativity is allowed in the quasiprobability distributions associated with the initial states, and/or generated during the intermediate quantum operations, as long as the final quantum states are nonnegatively represented.

Furthermore, bearing in mind that the above classical algorithm is developed in the form of classical stochastic processes, we can naturally generalize the algorithm to compute the quantum expectation values arising in quantum circuits which start from an arbitrary incoherent mixture (convex combination) of pure states, and evolves according to an arbitrary mixture of the (allowed class of) unitary quantum gates, as long as the associated classical mixing probabilities can be efficiently sampled. We only need to combine the sampling from the ER phase space distribution associated with the initial pure states given by Eqs. (7) and (9), with the sampling from their mixing probabilities. Moreover, the transition probability must now become a mixture of the transition probabilities generated by each unitary gate. Hence, incoherent mixing does not present significant fundamental difference.

Finally, let us discuss the convergence rate of the above classical algorithm for computing quantum expectation values using the Monte Carlo sampling from the ER phase space distribution. For this purpose, let us make more explicit the steps of the classical algorithm. For the sake of discussion, first, we ignore the unitary transformation and consider the computation of the quantum expectation value \(\langle \psi |\hat{O}|\psi \rangle\) of a Hermitian observable \(\hat{O}\) over a pure quantum state \(|\psi \rangle\). Our classical algorithm starts by sampling the phase space variables (p, q) from the ER phase space distribution \({\text{ P }}_{\{S,\rho \}}(p,q|\xi )\) associated with \(\psi (q)=\langle q|\psi \rangle =\sqrt{\rho (q)}e^{iS(q)/\hbar }\) given by Eq. (7). We then proceed to evaluate O(q, p) associated with \(\hat{O}\) for the above sampled value of (p, q). We independently repeat this protocol for sufficiently large \(K_{{\mathrm{C}}}\) number of times. Let us denote the value of O for the kth sample as \(\tilde{O}_k\), and compute the classical sample mean (average)

$$\begin{aligned} \langle O\rangle _{\{(S,\rho );K_{{\mathrm{C}}}\}}\doteq \frac{1}{K_{{\mathrm{C}}}}\sum _{k=1}^{K_{{\mathrm{C}}}}\tilde{O}_k. \end{aligned}$$
(24)

The law of large numbers and Eq. (8) then ensures that in the limit of infinite number of samples, the above classical sample mean approaches the quantum expectation value as

$$\begin{aligned} \lim _{K_{{\mathrm{C}}}\rightarrow \infty }\langle O\rangle _{\{(S,\rho );K_{{\mathrm{C}}}\}}=\int \mathrm{d}^Nq~\mathrm{d}\xi ~\mathrm{d}^Np~O(q,p){\mathbb P}_{\{S ,\rho \}}(p,q|\xi )\chi (\xi )=\langle \psi |\hat{O}|\psi \rangle . \end{aligned}$$
(25)

We note that \(\tilde{O}_k\), i.e., the value of the physical quantity O(p, q) computed for the kth sample, is in general not equal to one of the eigenvalues of the associated quantum observable \(\hat{O}\). O(p, q) is unbounded as can be seen from the computation of p in Eq. (12), and may take on any continuum real number even when the associated quantum observable \(\hat{O}\) only allows discrete spectrum of eigenvalues. Moreover, as discussed above, we can straightforwardly insert a quantum circuit generating a unitary transformation \(\hat{U}_{\{f,g\}}\) as long as the resulting transformed observable \(\hat{O}'=\hat{U}^{\dagger }_{\{f,g\}}\hat{O}\hat{U}_{\{f,g\}}\) is at most second order in \(\hat{p}\), and the associated mapping of the classical quantity \(O(p,q)\mapsto O'(p,q)\) can be obtained via a phase space mapping, \(q\mapsto f(p,q)\) and \(p\mapsto g(p,q)\), that is computationally trackable using classical computer.

How does the classical sample mean \(\langle O\rangle _{\{(S,\rho );K_{{\mathrm{C}}}\}}\) for a finite \(K_{{\mathrm{C}}}\) number of samples computed in Eq. (24) approaches the targeted quantum expectation value \(\langle \psi |\hat{O}|\psi \rangle\), or, how many samples \(K_{{\mathrm{C}}}\) are required so that the classical sample mean is within a tolerated small error from the quantum expectation value, with a sufficiently high probability? This important question on convergence rate of the estimation of the quantum expectation value by classical sample mean can be assessed by applying the Chebyshev’s inequality58. Namely, the probability that the classical sample mean is within an error bound \(\epsilon \ge 0\) from the true quantum expectation value, is given by

$$\begin{aligned} {\mathbb P}\Big (\big |\langle O\rangle _{\{(S,\rho );K_{{\mathrm{C}}}\}}-\langle \psi |\hat{O}|\psi \rangle \big |\le \epsilon \Big )\ge 1-\frac{\mathrm{Var}_{\{S,\rho \}}[O]}{K_{{\mathrm{C}}}\epsilon ^2}. \end{aligned}$$
(26)

Here, \(\mathrm{Var}_{\{S,\rho \}}[O]\) is the variance of the classical quantity O associated with the quantum observable \(\hat{O}\), over the ER phase space distribution associated with the quantum wave function \(\psi (q)=\sqrt{\rho (q)}e^{iS(q)/\hbar }\), defined as

$$\begin{aligned} \mathrm{Var}_{\{S,\rho \}}[O]\doteq & {} \langle (O-\langle O\rangle _{\{S,\rho \}})^2\rangle _{\{S,\rho \}}. \end{aligned}$$
(27)

We shall refer to this quantity as the ER classical variance, and assume that it is finite. Equation (26) shows that, the number of samples \(K_{{\mathrm{C}}}\) needed for the classical sample mean estimate \(\langle O\rangle _{\{(S,\rho );K_{{\mathrm{C}}}\}}\) to be within an error \(\epsilon\) from the true quantum expectation value \(\langle \psi |\hat{O}|\psi \rangle\), with a probability at least \(1-\delta\), \(\delta \ge 0\), is proportional to the ER classical variance as:

$$\begin{aligned} K_{{\mathrm{C}}}\sim \frac{\mathrm{Var}_{\{S,\rho \}}[O]}{\epsilon ^2\delta }. \end{aligned}$$
(28)

It is therefore instructive to study how the ER classical variance varies with the quantum states and the quantum observables, and how it is related to the epistemic restriction and the other signatures of nonclassicality. We leave this important and interesting problem for future study.

Comparison with the classical simulation based on quasiprobability representation

We note that our strategy to classically compute the quantum expectation values arising in a class of CV quantum circuits employing the ER phase space representation is close in spirit to the classical simulation of quantum algorithm using the quasiprobability phase space representations reported in Refs.9,10,11,12,13. Both methods transform the quantum circuits or algorithms into classical albeit unconventional stochastic processes and apply the classical sampling method to evaluate the resulting multidimensional integral on classical computers. Below we compare the two schemes.

First, we emphasize that the epistemic restriction in classical phase space directly reflects the quantum uncertainty relation40,41, and the ER phase space distribution defined in Eq. (7) is always non-negative. On the other hand, the construction of quasiprobability phase space representation relies heavily on the algebraic structure of the abstract space of linear Hermitian operators. The associated quasiprobability phase space distributions may yield negative value or highly irregular so that in general they cannot be seen as proper probabilities. Such negative quasiprobability, which is seen as the signature of nonclassicality, is difficult to grasp and its relation with the quantum uncertainty is not immediately clear. Negative quasiprobability is argued to be formally related with the fact that the transformation that maps the quantum states and operations to the quasiprobability distributions is (bi)linear in the state vector \(|\psi \rangle\)59,60,61. In contrast to this, as expressed in Eqs. (7) and (9), the ER phase space distribution is nonlinear in \(|\psi \rangle\).

Further, recall importantly that the classical algorithm based on the ER phase space representation can only efficiently compute the quantum expectation values arising in a class of quantum circuits with the final quantum observables at most second order in the momentum operators (Theorem 3). Namely, it cannot be used to compute the expectation values of arbitrary Hermitian operators, including the expectation values of arbitrary POVM giving the probability of measurement outcomes. This limitation is in a sense to be expected. If our classical algorithm would apply to all forms of quantum observables, then quantum speedup might not exist. Additionally, notice that Eqs. (8) or (15) are basically local hidden variable models for quantum expectation values, so that, according to Bell theorem for continuous variable62, they cannot be applied to compute the expectation values of arbitrary Hermitian operators. On the other hand, within the quasiprobability approach, one can in principle directly use the classical sampling algorithm to efficiently simulate the probability of measurement outcomes as long as the quasiprobability distributions associated with the initial quantum states, quantum operations, and POVM are all nonnegative9,10,11,12,13. In this sense, the negativity in the quasiprobability representation was then suggested as a nonclassical ingredient responsible for the quantum speedup9,11,28,29.

Noting this, below we confine our discussion to the computation of the quantum expectation values arising in a class of CV quantum circuits yielding observables at most second order in the momentum operators so that both the classical sampling algorithm based on the ER phase space representation and those based on the quasiprobability representation apply. In this case, as mentioned in Theorems 2 and 3, unlike the classical algorithms based on the quasiprobability approach using the Wigner function reported in Refs.10,11, the classical algorithm based on the ER phase space representation can be applied to a large class of such quantum circuits with the initial quantum states admitting negative Wigner function and with nonGaussian quantum operations inducing nonlinear transformations. Our results thus show that negativity of Wigner functions and nonlinearity of quantum transformation are not sufficient for quantum speedup.

Next, employing the general quasiprobability representation developed in Ref.63, Pashayan et al. in Ref.12 devised a classical sampling algorithm which mitigates the presence of negativity problem in quasiprobability approach. They showed that the classical simulation converges at a slower rate as the amount of the negativity in the quasiprobability distributions increases. Note however that, as for the general classical simulation approach based on quasiprobability representations, to run the algorithm, one needs to compute the associated quasiprobability distributions. For CV systems of large size, unless the initial quantum state and operations are factorizable into those of smaller systems, this computation of the quasiprobability distributions typically involves convoluted multidimensional integrations, which are in general computationally hard. In this sense, these quasiprobability distributions are in general difficult to sample. Moreover, to construct the classical algorithm, one needs to compute the total amount of negativity in the quasiprobability distributions which is also difficult to do for large nonfactoziable systems. The total amount of negativity is crucial in the convergence analysis of the Monte Carlo sampling technique.

By contrast, as can be seen in Eq. (7), to get the associated ER phase space distribution \({\mathbb P}_{\{S,\rho \}}(p,q|\xi )\) from the initial wave function \(\psi (q)=\sqrt{\rho (q)}e^{iS(q)/\hbar }\), we only need to perform spatial differentiations locally around each sample trajectory, which are in general easier to handle. The classical algorithm based on the ER phase space representation indeed requires that \(\psi (q)=\langle q|\psi \rangle\) can be computed efficiently, \(\rho (q)=|\psi (q)|^2\) can be sampled efficiently, and the mapping of Eq. (13) is trackable on classical computers. These requirements are arguably less restrictive than the requirement of product form initial quantum states and operations in the quasiprobability approach. Furthermore, the convergence rate is determined by the ER classical variance defined in Eq. (27) which can also be estimated using the Monte Carlo sampling technique.

Conclusion and discussion

Despite many continuous efforts and remarkable progresses in the last decades, the boundary and correspondence between classical and quantum computations is still not fully understood. This practically important basic problem in quantum computation is arguably deeply related to the longstanding foundational problem of quantum-classical divide. In this work, we employ a novel ER phase phase representation of quantum statistics40,41, to devise a classical algorithm that can efficiently compute the quantum expectation values arising in a class of CV quantum circuits which yield, after the Heisenberg time evolution, quantum observables at most second order in momentum. The classical algorithm is obtained by transforming the quantum circuits in the complex Hilbert space into classical stochastic processes in the classical phase space, but with the initial positions and momentum being sampled from a specific ER phase space distribution associated with the initial quantum state. We then use the Monte Carlo sampling technique to evaluate the resulting multidimensional integral. It remains an open challenging important problem to extend the method to discrete variable quantum computations.

The classical algorithm shows that the epistemic restriction, including the nonseparability of \(\xi\), is necessary for quantum computational circuits or algorithms. Namely, we can define quantum computations as those that benefits from a specific epistemic restriction in the classical phase space, encoded as quantum uncertainty relation or canonical commutation relation in the Hilbert space. In the context of computation, the epistemic restriction may thus be (paradoxically) liberating. Moreover, for a certain class of quantum circuits computing the expectation values, we can decode the quantum uncertainty relation back in the form of epistemic restriction in the classical phase space, and exploit it to devise an efficient classical algorithm to compute the quantum expectation values. The result suggests that the epistemic restriction may offer an intuitive conceptual tool to further study the boundary between quantum and classical computations. To this end, it is interesting to explore the general connection between the epistemic restriction and other signatures of nonclassicality suggested to bound the convergence rate of classical sampling algorithms, e.g., negativity in quasiprobability approach12, and interference in an approach based on sampling over Feynman-like paths64.

On the other hand, the epistemic restriction is not sufficient for all possible (universal) quantum computational algorithms. The classical algorithm of Eq. (15) does not simulate the dynamics of the quantum states (except for the Gaussian sector as mentioned in Corollary 1), and the random outcome in quantum measurement. As discussed in Ref.40, to recover the dynamics of the quantum states within the ER phase space representation, we need to follow the dynamics of the whole ensemble of trajectories which is required to respect the conservation of average energy. Moreover, to get a definite measurement outcome, we also need to keep track of each single random trajectory which is allowed to violate the conservation of energy. This observation suggests that the key behind the quantum speedup might be the ability of Nature to manage conserving the ensemble average energy in the presence of a global random variable, while allowing each trajectory in the ensemble to violate randomly the conservation of energy. Intuitively, to manage to do these dual tasks naively using only classical resource is computationally hard. By contrast, in classical mechanics, the average ensemble energy is automatically conserved since each single trajectory is deterministic conserving the energy. As an analogy, while each trajectory in the Bohmian mechanics65 violates the conservation of (local) energy (defined suitably to include the so-called quantum potential), the whole ensemble of Bohmian trajectories manage to mysteriously anticipate each other so that the ensemble average energy is conserved. In Bohmian mechanics, this is due to the presence of a physical wave function evolving under the Schrödinger equation, which co-orchestrates the whole ensemble of trajectories. We hope to further clarify this conjecture in the future.