1. Introduction

The classical simulation of quantum mechanics is a critically important field for science in general. Not only does this pursuit allow us to push the boundaries of what is possible to simulate, it also enables new insights into notions of complexity and information theory. The field of quantum computing itself has its origins in classical simulation, with early observations by Feynman and Manin that quantum mechanics, seemingly, is incompatible with classical mechanics, in a particular computational sense. Indeed, the advent of quantum complexity theory and several pioneering works in the 1990s suggested that the class of problems efficiently computable on a (randomized) classical computer (BPP) is a strict subset of those efficiently computable on a quantum computer (BQP), up to standard complexity theoretic assumptions [1].

Quantum simulation has itself been driven by experimental developments in the field of quantum computing, for as sizes in the lab increase, increasingly better classical simulation tools are required for verification [2–4]. As a result, a plethora of techniques have been developed for such purposes in finite-dimensional systems, including state vector simulation [5], stabilizer state simulation [6–8], tensor network methods [9–12], and many other optimized routines within these broad classes [13–16].

Infinite-dimensional systems are typically treated quite distinctly from finite-dimensional ones in terms of simulation, although there are certain similarities. For example, a generalization of the Clifford group to infinite dimensions, via the Weyl relations, allows for a stabilizer formalism to be employed [17], which shows that Gaussian computations can be simulated efficiently, in a similar sense to the Gottesman–Knill theorem in finite dimensions [18,19]. Conversely, the phase space formalism, originally developed for continuous variable (CV) systems [20], can also be defined in finite dimensions, which connects the negativity of the discrete Wigner function to the ability to attain a “quantum speedup” [21–23]. In addition, and more recently, tensor network methods for Bosonic systems have also been explored [24–26]. Improvements to non-Gaussian simulation continue to increase the threshold for a quantum advantage [27–29] and our understanding of the intrinsic complexity of simulating these systems [30–32]. We review other relevant simulation techniques in Section 2.

1.1 Contributions

In this work, we outline an approach to the simulation of infinite-dimensional systems, inspired by stabilizer rank simulations in finite-dimensional systems with “magic” [8]. First we show (Theorem 1) that any state in the Fock basis can be decomposed, to arbitrary accuracy, as a finite superposition of coherent states. Then, with such a decomposition, we demonstrate that the cost of classically simulating linear optics (Theorem 2), and Fock basis measurements (Section 3.4), scale linearly with respect to the number of terms in the decomposition. This results in a simulation time complexity for computing Fock state amplitudes in the setting of Boson sampling, scaling as $O(m^2 2^n)$ for $n$ photons in $m$ modes, scale similarly in $n$ with techniques relying upon permanent calculations [though the memory overhead in our approach is worse, $O(m2^n)$ versus $O(m)$] [30,31]. In addition, we note that this approach can simulate bunched $n$ photon systems $|n\rangle |0\rangle ^{\otimes m-1}$ in a time that scales linearly in the number of photons, $O(nm^2)$ (note, permanent computations also reduce to similar complexity in this case). This allows up to a logarithmic number of N00N states $(|N0\rangle + |0N\rangle )$—fundamental for quantum metrology and sensing [33]—to be simulated efficiently, in a linear optical setting. Moreover, the framework is flexible in the sense that it can, in principle, simulate both linear and nonlinear optics, finite photon number states, coherent states, and squeezed states in the same manner (though these do not all have the same simulation complexity, i.e., some are less efficient).

In Section 5, we also discuss trade-offs that can be employed, interchanging accuracy (fidelity) for simulation cost (memory and time). Further, we outline an approach that can be applied to circuits with structure, partitioning the system according to the number of entangling operations, resulting in a time and memory saving for certain simulations.

This work provides a well-defined notion of classical complexity for the simulation of infinite-dimensional quantum systems, related to the efficiency of a decomposition to coherent states. Indeed, coherent states are often considered the “most classical” of quantum states [34], and the simulation cost in this framework is linear in the number of these classical states. As such, we construct a quantum resource theory for coherent state basis simulations (Section 4), which further relates this “non-classicality” to: (i) the number of photon additions and (ii) the amount of squeezing required to construct a state. We compare this resource to the non-Gaussianity measures, the Wigner negativity, and the stellar rank.

2. Background

2.1 Hilbert Space, Operators, and States

Here we introduce necessary prerequisites and notation for what is to follow in Section 3. The starting point is a separable infinite-dimensional Hilbert space, $\mathcal {H}_\infty$, such as the Hilbert space for a particle in one dimension (the space of square integrable functions over the real line). By definition, such a Hilbert space admits a countable orthonormal basis, $\mathcal {H}_\infty =\mathrm {span}\{|n\rangle \}_{n=0}^\infty$, which we will call the Fock basis. Each basis element, $|n\rangle$, will be called a single mode Fock state (or simply Fock state). Although we have not specified the system per se, for ease of language, we will refer to $|n\rangle$ as the Fock state of $n$ photons. Also note that while in this section we will focus on the single-mode case, the discussion generalizes straightforwardly to the multi-mode case, $\mathcal {H} = \mathcal {H}_\infty ^{\otimes m} = \mathrm {span}\{|n_1, \dots, n_m\rangle \}_{n_i=0}^\infty$, where $m$ is the number of modes.

With the basis in hand, we can construct relevant operators for quantum optics. In particular, we define the annihilation and creation operators,

These operators have other remarkable properties. Namely, it is easy to show they admit a non-trivial commutation relation $[\hat {q}, \hat {p}]=i/2$ (where the right-hand side has an implicit identity operator). It is further possible to show these operators have a continuous and real spectra [36], i.e.,

This brings our attention to the displacement operator, $\hat {D}(\alpha )$. This can be used to construct an uncountably infinite family of quantum states, known as coherent states. In particular, by starting at the so-called vacuum state $|0\rangle$ (the Fock state of 0 photons) and applying the displacement operator, one arrives at the coherent state

Coherent states themselves have several useful and interesting properties. First, they form a basis over the Hilbert space

Coherent states are often referred to as the “most classical” states of a CV system [34]. This is because they minimize the uncertainty relations for both quadratures, $\hat {q},\hat {p}$, with equal variance in either observable. Squeezed states, however, also minimize the uncertainty relation, while not equalizing the variance of each observable. For example, one could minimize the position uncertainty at the cost of increasing the momentum uncertainty. Define the “squeeze” operator as [38]

Before delving into more details, it is worth listing some useful identities involving displacement and squeeze operators, and their canonical transformations. The displacement operators are unitary and form a one-parameter family,

Similarly, the squeeze operator is unitary,

2.2 QPD Simulation

Central to the discussion of optics simulation is the Wigner function, which provides a phase space formalism for infinite-dimensional quantum systems. Every state $\rho$ of $m$ modes has an associated Wigner function, $W:\mathbb {C}^{m}\rightarrow \mathbb {R}$, which is informationally equivalent to the state. For a single-mode state $\rho$, the Wigner function is defined as [39]

The Wigner function for a coherent state $|\beta \rangle$ is a Gaussian:

Gaussian operations (including measurements) are defined as those which map Gaussian states to Gaussian states, which can be implemented efficiently, in time $O(poly(m))$ for $m$ modes [40]. This is achieved by representing the state by its Wigner function, which for pure Gaussian states, is a Gaussian function, thus depending only on $O(m^2)$ numbers to specify the first two moments of the distribution (means and variances of the quadratures). Gaussian operations include displacement operators, beam splitters, phase shifts, as well as measurements of the position or momentum operators (homodyne) [41], or measuring in the coherent state basis (heterodyne).

We can see an example where QPD positivity does not hold. The simplest case, perhaps, is a cat state, which is a normalized state of the form $|\psi _\pm \rangle = c(|\alpha \rangle \pm |-\alpha \rangle )$, where the $+$ ($-$) sign is often called an even (odd) cat state, due to the fact that in the Fock basis, it only contains even (odd) photon numbers. Due to the non-orthogonality, the normalization condition $\langle \psi _\pm |\psi _\pm \rangle =1$ is non-trivial, from computing overlaps via Eq. (4): $2|c|^2 = 1/(1 \pm e^{-2|\alpha |^2})$. In Supplement 1B we show that the Wigner function for an off-diagonal term $|\alpha \rangle \langle \beta |$ is

 

Fig. 1. Wigner function for even cat state. The Wigner function evaluated at $\alpha = \alpha _r + i\alpha _i$ for a normalized even cat state $|2+2j\rangle +|-2-2j\rangle$.

Download Full Size | PDF

We also briefly mention for completeness that the theory of efficient simulation can be generalized beyond the Wigner function, for the Wigner function is just one QPD in a continuum [43]. In all of these cases, the efficient QPD simulation is only possible when the relevant QPDs (states, operators, measurements) represent a genuine probability distribution, otherwise a sign problem of sorts is present, rendering the efficient classical simulation by these means, in general, infeasible. Indeed, QPD negativity is often considered a quantum resource and, as discussed in more detail in Supplement 1A, the classical simulation complexity can be shown to scale exponentially in the number $n$ of non-Gaussian initial states, $\mathcal {N}^{2n}$, where $\mathcal {N}=\int d\alpha |W(\alpha )|\ge 1$ measures the amount of negativity ($\mathcal {N}=1$ if and only if the Wigner function is non-negative).

2.2.1 Linear Combination of Gaussians

References [42,44] generalize the discussion above to allow simulations where the state is represented as a linear combination of Gaussian Wigner functions, $W_\rho (\alpha ) = \sum _{g} c_g W_g(\alpha )$. In particular, it is demonstrated in Ref. [42] how to represent many states of interest in this framework, including GKP and Fock states, the latter of which can describe $|n\rangle$ with $n+1$ Gaussian Wigner functions (to arbitrary accuracy). Update rules for Gaussian operations are given in [44], and non-Gaussian measurements can be performed by representing the measurement outcome as a linear combination of Gaussians, meaning that (Gaussian) Boson sampling can be simulated with time complexity scaling as [$O(2^n)$] $O(4^n)$, for $n$ photons in the output.

2.3 Stellar Rank

States in the Hilbert space can be constructed using the stellar representation [45,46]. In particular, every single-mode state $|\psi \rangle = \sum _{n=0}^\infty \psi _n |n\rangle$ has an associated stellar function

An additional and equally remarkable result is that the classical simulation complexity is intimately related to the stellar rank [48]. In particular, a Fock state with stellar rank $r$ can be simulated under Gaussian operations/measurements with complexity $O(r^3 2^r + poly(m))$, for $m$ modes. Note for an initial state of $n$ photons (e.g., $|1\rangle ^{\otimes n}$), the stellar rank is $r=n$, as the stellar rank is additive under tensor products for pure states.

If one additionally wants to perform a measurement in a non-Gaussian basis (such as the Fock basis), this can be done at an additional cost related to the stellar rank of the output measurement, $r = r_\psi + \sum _{k=1}^m r_k$, where $r_k$ is the stellar rank for a particular output result in the $k$th mode, and $r_\psi$ the stellar rank of the initial state [49]. For an $n$ photon simulation (such as Boson sampling), $r = 2n$, as the input and output each have rank $n$.

2.4 Direct Fock State Simulation of Linear Optics

To simulate systems of strictly finite photon number, there exists another, more “direct” approach, by simulating in the Fock basis. Simulating Fock states is notoriously tricky, in fact, this forms the basis for the classical hardness of Boson sampling [30]. The general setup of interest in the present work is that of linear optics, where $n$ Fock states are injected over $m$ spatial modes (which are physical, e.g., fiber optics or waveguides). That is, the initial state is of the form $|n_1, n_2, \dots, n_n\rangle |0\rangle ^{\otimes (m-n)}$. The typical use case is for single-photon states, i.e., $n_i=1$. Linear optical (LO) operations—beam splitters and phase shifts—are then applied between and on these spatial modes, after which Fock basis measurements are performed. The best classical algorithms known for computing single amplitudes [30,50] and outputting measurement results [31] run in time (space) $O(n 2^n)$ [$O(m)$] for $n$ single-photon inputs, based on calculations of the permanent (or slightly better in some cases [32]). Full Fock state simulation of LO is however much more expensive in general, due to the size of the Hilbert space; for $n$ photons distributed among $m$ modes, the dimension is

A general LO transformation is generated by a “bilinear” Hamiltonian, of the form

Since such transformations preserve photon number, the evolution of a state can be described by evolving its creation operators:

Of course, in systems with symmetry, not every mode will be populated, and the actual cost can be significantly less (the simulation cost is linear in the number of Fock states with non-zero amplitude). One can review Ref. [53] for a more detailed exposé on the current state of the art for Fock simulation.

In the next section, we will outline a simulation approach using coherent states, where the simulation cost can be much better than the full Fock dimension, even if all Fock states are populated. In particular, while the simulation cost, in the worse case, is exponential in the number of photons, it is only quadratic in the number of modes.

3. Simulation by Coherent State Rank

Here we outline a general strategy for simulating infinite-dimensional systems using a finite rank decomposition over coherent states, i.e., states of the form $\sum _{i=1}^k c_i|\alpha _i^{(1)}, \dots, \alpha _i^{(m)}\rangle$ for $m$ modes, where each $|\alpha _i^{(j)}\rangle$ is a coherent state, and $k$ is the “coherent rank.” In general, we will allow for such a decomposition to be approximate, in the following sense.

Definition 1. The “approximate coherent rank” of a state $|\psi \rangle$ is the smallest integer $k$ such that for any $\epsilon >0$, there exists a coherent rank $k$ state $|\tilde {\psi }\rangle$ where $|\langle \psi |\tilde {\psi }\rangle |^2 > 1 - \epsilon$.

With such a representation, the memory requirements are simply from storing $k(m+1)$ complex numbers, i.e., each of the $c_i$, and the $\alpha _i^{(j)}$. Overall, we follow a similar approach as in finite rank stabilizer simulations, see, e.g., [8]. In particular, we first show how to decompose Fock states into an approximate finite rank coherent state representation. We then discuss operations that are “free” within this approach, i.e., those that do not increase the rank, and can be applied in time $O(k)$. Following this, we consider “resourceful” operations (that increase the rank), and then various methods to perform Fock basis measurements.

We will see, in many situations, this will outperform the direct Fock evolution [Eq. (21)], potentially exponentially so, as the time complexity for computing amplitudes scales similarly as methods relying upon the permanent.

3.1 Fock Basis Decomposition

To simulate a Fock state using coherent states, we must first outline a strategy for converting a (possibly multi-mode) Fock state to a coherent state superposition. This relies upon the following theorem.

Theorem 1. A single-mode state in the Fock basis containing at most $N$ photons, $|\psi \rangle =\sum _{n=0}^Na_n |n\rangle$, has approximate coherent rank at most $N+1$. One explicit construction for this is

Proof. We show this by construction. Consider the target state $|\psi \rangle =\sum _{n=0}^N a_n|n\rangle$, where $\sum _n |a_n|^2=1$ (the protocol can also be applied to unnormalized “states”—one can always first normalize it and later multiply by the inverse—but for ease of exposition, we assume it is normalized here). Using Eq. (1), we wish to satisfy the relations

This leads us to the following corollary.

Corollary 1. The state

Proof. This follows immediately from direct application of Eq. (22). □

We can observe that for a single-photon state, $N=1$, this is an odd (normalized) cat state of the form $|\epsilon \rangle - |-\epsilon \rangle$, where one can see by inspection that this indeed tends to $|1\rangle$ as $\epsilon \rightarrow 0$ (all even Fock terms are 0 by construction, and all weight will concentrate on the $|1\rangle$ state for small enough $\epsilon$). For higher photon number Fock states, this idea generalizes; we see that $|\tilde {N}\rangle$ only has non-zero amplitude on $|N\rangle, |2N+1\rangle, |3N+2\rangle,\dots$.

After submitting the first version of this manuscript, we became aware that similar results as Corollary 1 have previously been observed. For example, in Ref. [54], a projector onto select Fock states is constructed as a sum of $N+1$ phase shifts, $\sum _{k=0}^N e^{(\hat {n}-N) 2 \pi i k /(N+1)}$, which when acting on a coherent state $|\epsilon \rangle$ gives Eq. (28), upon normalization. In addition, decompositions of this type were also used in Refs. [55,56].

This shows that to approximate a Fock state of the form $|n_1, n_1, \dots, n_m\rangle$ results in a superposition of $\prod _{i=1}^m (n_i+1)$ coherent states. For example, the initial state for Boson sampling, $|1\rangle ^{\otimes n}|0\rangle ^{\otimes (m-n)}$, results in a superposition of $2^n$ coherent states. We will see that under certain “free operations” that do not increase the rank (Section 3.2), this can be much more efficient than storing each individual Fock amplitude, as was discussed in Section 2.4. Moreover, a simulation of a system with initial state $|n\rangle |0\rangle ^{\otimes (m-1)}$ is only a superposition of $n+1$ coherent states. Thus, under free operations, such a system is efficient to simulate, whereas in the Fock basis, this could result in $d_{n,m}$ [Eq. (17)] amplitudes to store.

We also mention here that Theorem 1 is technically only an upper bound on the approximate coherent rank, since we do not rule out more efficient representations (i.e., the optimality of Theorem 1 is an open question). It is however easy to verify for the Fock state $|1\rangle$, the approximate coherent rank is exactly 2.

While the above construction was presented for finite Fock state representations (i.e., up to $N$ photons), it can, in certain cases, be used to approximate states with support over the entire Hilbert space. As a pertinent example, we apply this to squeezed vacuum states $|\zeta \rangle$ (as used in Gaussian Boson sampling [57]), which by Eq. (5) is

 

Fig. 2. Squeezed vacuum state approximation. Infidelity of approximate squeezed vacuum state $|\tilde {\zeta }\rangle$ with respect to reference target $|\zeta \rangle$ of Eq. (29) (note, here we set $\phi =0$ in $\zeta =re^{i\phi }$). The approximation $|\tilde {\zeta }\rangle$ is a superposition of $N$ coherent states as described in Supplement 1C.

Download Full Size | PDF

Lastly, we show we can rigorously upper bound the coherent state rank, in particular, we have the following proposition.

Proposition 1. A single-mode state $|\psi \rangle$ with finite stellar rank$r$,

Proof. Using Eq. (7), we have

It is important to note that in the stellar rank formalism, the state $| G_{\psi } \rangle$ is a general Gaussian state and therefore, it may have nontrivial squeezing $\zeta \neq \,0$, which can result in an arbitrarily large coherent rank. However, for $\zeta =O(1)$, by keeping terms up to $n$th order in the Taylor expansion of $\hat {S}(\zeta )$, we can approximate $| \zeta \rangle$ with approximate coherent rank $2n+1$.

3.2 Free Operations

We consider the class of operations that do not increase the rank, $k$, for states of the form $\sum _{i=1}^k c_i|\alpha _i^{(1)}, \dots, \alpha _i^{(m)}\rangle$. Generally speaking, such operations are those that map a tensor product of $m$ coherent states into another tensor product of $m$ coherent states. In particular, such operations generate no mode entanglement when acting on coherent states, which contains the set of linear optical operations. In line with the resource theory literature, we call such operations “free operations.”

The first obvious example is a single-mode displacement operator, which to update the state, takes time $O(k)$, simply by updating the relevant single-mode coherent state for each term in the sum, as well as the amplitudes: $\hat {D}(\beta )|\alpha \rangle = e^{i \mathrm {Im}( \bar {\alpha }\beta )}|\alpha + \beta \rangle$.

Another obvious and relevant example is a single-mode phase shift, as used in linear optics, $\hat {P}(\phi ) = e^{i\phi \hat {n}}$. It is easy to see such an operation acts as $\hat {P}(\phi )|\alpha \rangle = |e^{i\phi }\alpha \rangle$.

In fact, the result can be generalized to that of all LO operations.

Theorem 2. All linear optical unitary transformations $\hat {U}$, with respective transfer matrix $u$ [Eq. (19)], are free. The time complexity for updating a rank $k$ state is $O(k\ell ^2)$, where $\hat {U}$($u$) acts non-trivially on $\ell \le m$ modes.

Proof. Consider the action of such a unitary on $|\alpha _1, \dots, \alpha _m\rangle$:

Using that $\hat {U}e^{\hat {X}}\hat {U}^\dagger = e^{\hat {U}\hat {X}\hat {U}^\dagger }$ by unitarity, we can write

For LO transformations therefore, the update can be performed by simple matrix-vector multiplication,

The above shows the coherent rank of a state is unchanged by any LO operation. Further, if $\hat {U}$ only acts non-trivially on $\ell \le m$ modes, then the relevant submatrix of $u$ has size $\ell \times \ell$, resulting in computational cost to update a rank $k$ state, $O(k \ell ^2)$. □

We can use Eq. (33) to show for a beam splitter described by a unitary $\hat {B}$ with transfer matrix $u$,

Lastly, while not a unitary operation, it is worth noting that, essentially by definition, the annihilation operator is also free (in fact, it is a free operation as opposed to a free unitary, since it can consume resource, see Supplement 1F), as it just modifies each complex amplitude $c_i\rightarrow c_i\alpha _{i}^{(j)}$, when applying $\hat {a}_j$ (therefore taking time $O(k)$ to update the state). In addition, any function $f(\hat {a})$ can be similarly implemented (assuming it can be computed classically in polynomial time). In general, the creation operator is not free, as we will see below. However, noting that $\langle \phi | \hat {a}|\psi \rangle = \overline {\langle \psi |\hat {a}^\dagger |\phi \rangle }$, the simulation cost for computing output amplitudes/probabilities in circuits composed of free unitaries and only creation operators is equivalent to the time-reversed circuit containing entirely free (unitary and annihilation) operators. See Fig. 3 for an example.

 

Fig. 3. Example of a free circuit and its inverse. The application of only linear optical (LO) unitaries, displacements, and annihilation operators are free, meaning the coherent rank of $|\psi \rangle$ is unchanged under such a circuit (left). Similarly, if one has a circuit of only LO unitaries, displacements and creation operators (right), the output probabilities can be computed with the same cost as the time-reversed (free) circuit. Circuits with both creation and annihilation operators however incur additional cost for implementing the non-free operations.

Download Full Size | PDF

3.3 Resourceful Operations

Any operation that is not classed as free is “resourceful.” Such an operation, when acting on a multi-mode coherent state will generate a non-trivial superposition. One example is a single-mode squeezing operator, Eq. (5), which when acting on the (coherent) vacuum generates a squeezed vacuum state, Eq. (29), and hence squeezing is not a free operation (see Supplement 1C for one decomposition). Another example (though not unitary) is the creation operator, $\hat {a}^\dagger$, which is clearly not free. In fact, as we will see below in Corollary 2, $\hat {a}^\dagger$ causes the rank of a state to double. We can first make a definition analogous to Definition 1 at the operator level.

Definition 2. The (approximate) “operator coherent rank” for an $m$-mode operator $\hat {A}$ is $\ell$ if and only if $\hat {A}|\vec {\alpha }\rangle$ has (approximate) coherent state rank at most $\ell$, for all $|\vec {\alpha }\rangle$.

In other words, the action of an operator with coherent rank $\ell$ on a coherent state results in a state of coherent rank at most $\ell$ (and this is achieved for at least one coherent state). One can readily observe that an operator with (approximate) operator coherent rank $\ell$ can be implemented on a state with (approximate) coherent rank $k$, resulting in a state of at most (approximate) coherent rank $k\ell$.

An operator composed as a sum of displacements, $\hat {U} = \sum _{i=1}^{\ell } a_i \hat {D}_i(\alpha _i)$, has operator coherent rank at most $\ell$, and similarly for an operator represented in the coherent state basis, $\hat {U}=\sum _{i,j=1}^{\ell } a_{i,j}|\alpha _i\rangle \langle \beta _i |$.

One potential construction for this is to use a similar Fourier decomposition as in Theorem 1, but at the level of an operator. Consider an operator acting on at most $N$ photons: $\hat {A}=\sum _{n,m=0}^N A_{n,m} |n\rangle \langle m |$. This could result, for example, from the truncation of an operator to the maximum $N$ photon subspace. Then one can use the following proposition.

Proposition 2. The operator $\hat {A}=\sum _{n,m=0}^N A_{n,m} |n\rangle \langle m |$ in the Fock basis up to $N$ photons has approximate operator coherent rank at most $N+1$. One explicit representation is

Proof. This follows almost exactly that of Theorem 1. □

Such a decomposition can also be applied to multi-mode operators. One clear issue with such a method is that if $N$ is large, the state rank $k$ will increase very quickly $k\rightarrow k(N+1)$, meaning very few operators in this form can be applied. Nevertheless, this allows one, in principle, to achieve universality with displacements and SNAP operations $\exp (i\theta _n|n\rangle \langle n |)=|n\rangle \langle n |(e^{i\theta _n}-1) + \mathbb {I}$ [58], which has approximate operator coherent rank at most $n+2$.

A second, perhaps more practical approach, is to recompute each coherent state in the superposition to a desired accuracy in the Fock basis, and then reexpress the state via Theorem 1. For example, in the single-mode case, $\hat {U}|\psi \rangle = \sum _{i=1}^k c_i \hat {U}|\alpha _i\rangle$, and one can truncate each $\hat {U}|\alpha _i\rangle = \sum _{n\ge 0} a_n |n\rangle \approx \sum _{n=0}^{N_i} a_n |n\rangle$ to the $N_i$ photon subspace (to desired accuracy), which can be (approximately) written as $N_i+1$ coherent states. In this example, the output state has coherent rank at most $k + \sum _{i=1}^k N_i$.

Lastly, we show some results pertaining to the implementation of creation operators. We will see a creation operator can be implemented doubling the rank, or in general, an operator composed of $n$ such operators will increase the rank by a factor of $n+1$.

Theorem 3. Any single-mode operator composed of only creation and annihilation operators, of the form $\hat {A}:=\sum _{i=1}^p b_i \prod _{j=1}^{l_i} (\hat {a}^{\dagger })^{n_{i,j}} (\hat {a})^{m_{i,j}}$, where there are at most $n$ creation operators in total for each term in the sum (i.e., $n \ge \sum _{j=1}^{l_i} n_{i,j}, \forall i$), has approximate operator coherent rank at most $n+1$.

Proof. Such an operator acts on a coherent state as

This leads immediately to the following observations (we also give an alternative proof of this in Supplement 1D using post-selection).

Corollary 2. The operator $\hat {a}^\dagger$ has approximate operator coherent rank $2$ and an operator $\sum _{l=0}^n b_l (\hat {a}^{\dagger })^l$ has approximate operator coherent rank at most $n+1$. The classical simulation of photon-added coherent states therefore has a cost linear in the number of photon additions.

Theorem 3 allows one to implement operators generated by the creation/annihilation operators. For example, squeezing operations [Eq. (5)] can be implemented to arbitrary accuracy by expanding $\hat {S}(\zeta )$ to order $n$ [terms up to $(\bar {\zeta }\hat {a}^2 - \zeta \hat {a}^{\dagger 2})^n$], which will increase the states coherent rank by a factor of at most $2n+1$.

In general, applying multiple resourceful operations in a circuit causes an exponential overhead for the classical simulation complexity. For example, $n$ independent applications of the creation operator (potentially on different modes) will generically result in an $O(2^n)$ increase in the coherent state rank, by Corollary 2. This is similar to the case of stabilizer/magic simulations, where only a logarithmic number of $T$ gates can be tolerated [8]. Of course, in some situations, this can be improved upon as demonstrated above, where (for example) applying $n$ creation operators sequentially on a single mode only increases the simulation cost linearly (instead of exponential) in $n$, meaning a polynomial (instead of logarithmic) number can be applied.

3.4 Measurements

The last ingredient in the coherent state rank simulation paradigm is that of measurements. Here we primarily consider Fock basis measurements, but will also mention properties of the coherent state basis measurement. There are several ways in which to perform a measurement in this approach, each of which has different pros and cons.

For consistency, we define an $m$ mode state of rank $k$ as before: $|\psi \rangle = \sum _{i=1}^k c_i|\alpha _i^{(1)}, \dots, \alpha _i^{(m)}\rangle$. This could correspond to the state after evolving a Fock state through a linear optical network, for example. We will see that computing the probability for a particular Fock/coherent basis measurement can be achieved in a time linear in the rank $k$, but for actually outputting measurement results according to the desired distribution, there are additional complications.

3.4.1 Computing Individual Amplitudes

Before discussing probabilistic measurement (sampling), we first observe that the computation of individual amplitudes/probabilities can be performed in time $O(mk)$ by evaluating the overlaps

We further note that certain output amplitudes can be computed much more efficiently in the setting of Boson sampling by simulating the time-reversed circuit. For example, since $\langle n, 0, \dots, 0| \hat {U}|1,\dots, 1\rangle =\overline {\langle 1, \dots, 1| \hat {U}^\dagger |n, 0, \dots, 0\rangle }$, and the state $|n, 0, \dots, 0\rangle$ has rank $n+1$, the simulation time complexity to compute this output amplitude is $O(m^2 n)$, instead of $O(m^2 2^n)$ for the forward circuit. In general, the cost to compute some particular transition amplitude is determined by which of the input or output has the most efficient coherent rank representation:

To calculate the equivalent conditional coherent basis normalization of $\langle \beta _1, \dots, \beta _p|\psi \rangle$ via Eq. (38), one replaces the $\hat {a}_j^{\dagger n_j}, \hat {a}_\ell ^{n_\ell }$ with displacement operators $\hat {D}(\beta _j), \hat {D}(-\beta _\ell )$, which as they are free, the cost of computing this norm is linear in the rank, $O(m^2 k)$.

3.4.2 Metropolis Sampling Algorithm

For sake of an example, consider the case of interest being Boson sampling-like simulations, with $n$ single-photon inputs in $m$ modes. The coherent rank of such a state is $k=2^n$ and the Fock dimension is $d_{n,m}$ [Eq. (17)]. In the most “naive” approach to outputting measurement samples from the desired distribution, one can simply compute measurement probabilities $p_{\vec {n}}$ until some (random) threshold $r<1$ is reached, $\sum _{\vec {n}} p_{\vec {n}}>r$. Although computing each individual measurement probability is linear in $k$, $O(m k)$, this procedure may require the evaluation of $O(d_{n,m})$ probabilities, i.e., resulting in overall time complexity $O(m 2^n d_{n,m})$. This can therefore be more costly (in time) than doing the full Fock basis simulation, which scales as $O(n d_{n,m}$) for outputting measurement samples, as discussed in Section 2.4 [the memory overhead will still, in general, be better in the coherent representation however, $O(m2^n)$ versus $O(d_{n,m})$].

We can improve upon this “naive” sampling by a heuristic approach, using a Markov chain. This could also be applied to the coherent state basis, though as it is a continuous distribution, we focus on Fock measurements here for simplicity. Similar methods have been proposed and used for Boson sampling [59] and in stabilizer rank simulation [8]. We will outline a very high level version of this type of algorithm as an example for the reader, but mention that this can be implemented in many different ways, each likely with different benefits and drawbacks.

The main idea is to set up a Markov chain, where each state in the chain is a measurement outcome. One can propose new outcomes and accept them with the appropriate probability. For example, we have the following.

While this method certainly can result in an efficient sampling for practical cases of interest, it has the usual issues and suffers from the same constraints as in most applications of Metropolis sampling. We already briefly discussed the issue of the thermalization time above. Another issue worth addressing, as also described in Ref. [8], is the Markov chain must be irreducible; given two configurations (measurement results) with non-zero probability, there needs to be a path between them, where the path is defined by the type of updates in step 3 above. While in random circuits (such as for Boson sampling) it is likely to have this property (since all measurement results will likely have a non-zero amplitude), in circuits with structure, many amplitudes may be identically 0, leaving the possibility of disconnected regions. Though this can be mitigated to some extent through numerical techniques (dynamically increasing the search area, random restarts, etc.), the Metropolis sampling approach may not always be a feasible one.

3.4.3 Conditional Probability Sampling

Here we describe an approach to measurement using conditional probabilities. As before, in principle, this could be applied to the coherent state basis also, though the continuous nature adds additional technical details. We will describe the procedure in words.

Starting with mode 1, output the measurement result $n_1$ with probability $P^{(1)}_{n_1}=\|\langle n_1|\psi \rangle \|^2$ (where $\|\phi \|^2=\| |\phi \rangle \|^2 = \langle \phi |\phi \rangle$). To implement this, one can pick a random number $r\in (0,1)$, and output $n_1$ such that $\sum _{n=0}^{n_1}P^{(1)}_n>r$. Then, with the $m-1$ mode (normalized) state $|\psi _{2,\dots, m}\rangle =\frac {1}{{\|\langle n_1|\psi \rangle \|}} \langle n_1|\psi \rangle$, one applies the same algorithm to the next mode, outputting $n_2$ with probability $\| \langle n _2|\psi _{2, \dots m} \rangle \|^2$. Continue until all $m$ modes have been measured.

The cost of this procedure to output result $(n_1, \dots, n_m)$ is $O(\sum _{i=1}^m n_i\mathcal {N}_i)$, where $\mathcal {N}_i$ is the cost of computing the norm of the projected state $\langle n _1, \dots, n_{i}|\psi \rangle$. As discussed in Section 3.4.1, the worst case scaling for such a computation is $O(mk^2)$, meaning the cost to output a single measurement sample is $O(Nmk^2)$, where $N=\sum _i n_i$ is the total number of photons. The best or average cases however may be significantly better than this, see Eq. (38). In general though, this is a huge computational overhead, e.g., for single-photon inputs, $k=2^n$. (Note, in the case of coherent basis measurements—heterodyne—the norm can be computed in a time linear in $k$.)

There are however two possible ways to improve upon this. First, since the overlap between coherent states decreases exponentially in their separation, $|\langle \alpha |\beta \rangle |^2 = e^{-|\alpha - \beta |^2}$, one may be able to pick a representation where all states are approximately orthogonal, in which case one can output a single measurement sample of $N$ photons in time $O(Nk)$ [if the states in the decomposition are approximately orthogonal, the norm can be computed in time $O(k)$]. One example as to how this can be achieved follows the Fourier decomposition in Eq. (22). Consider approximating the Fock state $|n\rangle$. Note one can increase the error parameter $\epsilon$ (which decreases the overlap between states in the decomposition) and the expense of increasing the rank of the state (e.g. take $N>n$ terms, where the amplitudes $a_i=\delta _{i,n}$ for $i=0, \dots, N-1$). The parameter $\epsilon$ does not, per se, have to be small, since the error goes as $N! \epsilon ^{2N+1}/(2N+1)!$ and so can be compensated by increasing $N$. This can allow one to construct a decomposition with a greater number of terms, but where the coherent states in the superposition are approximately orthonormal. Since unitarity preserves the overlap of the initial states in the decomposition, one only needs to guarantee the orthogonality in the initial states construction.

A second approach is to use a method of estimating the norm, as demonstrated in Ref. [8], for a stabilizer rank decomposition. This relies upon the following proposition.

Proposition 3. For sufficiently large $L$, the expected overlap of coherent state $|\beta _1, \dots, \beta _m\rangle$, where $|\beta _i|<L$, with an $m$ mode (generally unnormalized) “state” $|\psi \rangle =\sum _{i=1}^k c_i|\alpha _i^{(1)}, \dots, \alpha _i^{(m)}\rangle$ is

Proof. We perform the calculation for a single mode (the result generalizes straightforwardly). We wish to compute

In the $m$ mode case, we just get a product of integrals, and since the trace distributes over tensor products $\mathrm {Tr}[A\otimes B] = \mathrm {Tr}[A]\mathrm {Tr}[B]$, the result is the same, but with normalization $L^{2m}$. □

Proposition 3 implies yet another method for computing the norm linearly in the rank, in principle. In particular, since each overlap squared $|\langle \beta _1, \dots, \beta _m|\psi \rangle |^2$ costs $O(mk)$ steps to compute (as previously discussed), the cost for estimating the norm is $O(Tmk)$, where $T$ is the number of random samples required to estimate the norm to desired accuracy. Combining this with the above description for conditional sampling, the cost for outputting a single measurement sample of $N$ photons is $O(NTm k)$.

Of course, there are some potential practical issues to address. First, the choice of $L$ is critical; if it is too small, the estimate will not converge to the correct value, and if it is too large, the number of samples $T$ required will also be large. To pick $L$, one could start with a small value and increase it until convergence, but of course this costs additional resources. In Fig. 4(A), we show an example of estimating the norm for a single-mode state using Proposition 3. We see in this case, $L\approx 3$ suffices to get within 1% accuracy of the true value.

 

Fig. 4. Estimating the norm via Proposition 3. (A) We estimate the norm squared for a rank three unnormalized state $\frac {1}{2} |2\rangle$ (using Corollary 1 with $\epsilon =0.35$), taking 750,000 samples per data point. (B) Numerical computation of the standard deviation from 750,000 samples, for states of the form $\frac {1}{2}|n\rangle ^{\otimes m}$, for $n=1,2$ in the legend (i.e., states with norm squared 1/4). We use Corollary 1 to write $|n\rangle ^{\otimes m}$ as a superposition of coherent states with rank $(n+1)^m$, where we use error parameter $\epsilon =0.35$. We also vary the sampling radius $L$. The straight lines are least squares fits to the data in the form $\sigma = A L^{c m}$, with $c$ shown in the legend.

Download Full Size | PDF

Another obvious issue is the number of samples required. First, let us define the random variable $X = L^{2m} |\langle \vec {\beta } |\psi \rangle |^2$, where $|\vec {\beta }\rangle = |\beta _1, \dots, \beta _m\rangle$. By Proposition 3, we know $\langle X\rangle _\beta \rightarrow \|\psi \|^2$ for a large enough sampling radius $L$. The variance, for sufficiently large $L$, is given by

It is interesting to note that the finite-dimensional (qubit) analog of this does not suffer such a sampling overhead, since the classical states—stabilizer states—form a 2-design and thus sampling fluctuations (standard deviation) are suppressed by the Hilbert space dimension [7,8,60]. In contrast, it is known that any set of CV states can at best form a 1-design (indeed, coherent states form a 1-design, as they resolve the identity) [61,62].

3.4.4 Summary of Measurement Results

We sum up the results of this section in Table 1 for Fock basis measurements, but direct the reader to the relevant sections above for a more nuanced picture. Note, the first row does not correspond to outputting statistics according to the measurement distribution, it is the cost for simply computing a single probability (the other rows do however correspond to outputting results according to the measurement distribution).

Table 1. Measurement Procedure Cost for $m$-Mode Rank-$k$ State

View Table | View all tables in this article

3.5 Comment on Errors

Since the framework described allows for a certain degree of imprecision in a decomposition (e.g., Fock states are only approximately represented, as in Theorem 1, albeit to arbitrary precision), a natural question is the propagation of errors. If the initial approximate state $|\tilde {\psi }\rangle$ has target fidelity $|\langle \tilde {\psi }|\psi \rangle |^2 = 1-\epsilon$ with the ideal state $|\psi \rangle$, under free unitary operations, the overlap remains unchanged. For an initial product state $|\tilde {\psi }\rangle ^{\otimes n}$, the fidelity is $(1-\epsilon )^n$. To achieve a total fidelity $1 - \delta$, one can take $\epsilon = 1 - (1-\delta )^{1/n} \approx \delta /n$ for small $\delta$ (i.e., it only requires linearly decreasing precision per state). Of course, if one implements resourceful operations inexactly (as in Corollary 2), errors can be introduced during the simulation; however, these can be controlled easily by picking the error parameters appropriately.

In Fig. 5, we show the output measurement probabilities from a Boson sampling like simulation, for sizes up to 10 photons in 10 modes, comparing the results from an exact simulation in the Fock basis (as outlined in Section 2.4) with simulation using the coherent rank decomposition. We see the results match very well, for all output states, where the initial fidelity is set to over 99.9% for each single-photon state.

 

Fig. 5. Boson sampling simulation. Starting with initial state $|1\rangle ^{\otimes n}$, where $n=6,8$, and $10$, we simulate the output from a Haar random unitary transfer matrix. In the Fock basis, this results in storing all 462, 6435, and 92378 complex amplitudes, respectively (i.e., the dimension of the Hilbert space for the $n$ photons in $n$ modes). The coherent state decomposition uses Corollary 1 to represent $|1\rangle$ as an odd cat state, with error parameter $\epsilon =0.2$, which has fidelity $|\langle 1|\tilde {1}\rangle |^2 > 0.999$. This results in $(n+1)2^n= 448, 2304,$ and $11264$ complex numbers stored in memory, respectively. We compute the Fock probabilities via Eq. (37), which are enumerated on the $x$-axis, sorted by the probability.

Download Full Size | PDF

Lastly we discuss the issue of numerical stability. We take the canonical example of starting in a state with $n$ single photons $|1\rangle ^{\otimes n}$ (though any other $n$ photon configuration would work equally well). In the coherent framework, this is represented as a product of odd cat states, which for small error $\epsilon$ can be written as

When computing a transition amplitude (again with a slight abuse of notation) $\langle \vec {n}|\hat {U}|1\rangle ^{\otimes n}$, we see the $\epsilon$ dependence can be almost dropped entirely:

In more general scenarios, this kind of “bookkeeping” may not be desired/feasible, in which case there could, in principle, be numerical precision issues at large sizes. This of course depends strongly on the desired accuracy. In Fig. 5, we see $\epsilon = 0.2$ gives high accuracy (above $99{\% }$) to all output probabilities. In Eq. (40), this yields a prefactor of $2.5^n$, which would likely pose no numerical problems at the accessible system sizes (i.e., those that fit in memory). Overall, the value $\epsilon$ should be chosen carefully to avoid numerical issues at large sizes. Fortunately, the fidelity of the approximate state converges quickly as $O(\epsilon ^{2(n+1)}/(n+1)!)$ for an $n$-photon single-mode state, as shown in Theorem 1. Unless very large sizes and extremely high precision are required, one can likely avoid any such issues.

4. Resource Theory of Coherent State Rank

The framework introduced in this work can be phrased as a quantum resource theory [63]. As a brief reminder, a resource theory is a triple, $(\mathcal {F}, \mathcal {O}, \mathcal {R})$, where $\mathcal {F}$ is the set of “free states”, a subset of quantum states $\rho \in \mathcal {S}(\mathcal {H})$ that are devoid of the resource of interest [$\mathcal {S}(\mathcal {H})$ is the state space, i.e., the set of all density matrices]; A$\mathcal {O}$ is the set of “free operations”, a subset of quantum channels, $\mathcal {E}: \mathcal {S}(\mathcal {H}) \rightarrow \mathcal {S}(\mathcal {H})$ that do not generate any resource when applied to states in $\mathcal {F}$; and $\mathcal {R}:\mathcal {S}(\mathcal {H}) \rightarrow \mathbb {R}_{\ge 0 }$ is a functional quantifying the amount of resource in quantum states, with certain properties such as: (i) faithfulness, i.e., vanishing on the resource-free states only, $\mathcal {R}(\rho ) = 0 \iff \rho \in \mathcal {F}$ and (ii) monotonicity under free operations, i.e., free operations can at most consume the resource but never increase it, $\mathcal {R}(\mathcal {E}(\rho )) \leq \mathcal {R}(\rho ) \,\forall \mathcal {E} \in \mathcal {O}$.

In our formulation, free states are multi-mode coherent states of the form $|\alpha _1, \dots, \alpha _m\rangle$, where $|\alpha _i\rangle$ is a coherent state, Eq. (1). Free operations are those described in Section 3.2, which map (multi-mode) coherent states to coherent states (up to the amplitude), such as beam splitters, phase shifts, displacement operators. Then, akin to the definition of “magic" or “non-stabilizerness" in stabilizer rank simulations [64–66], a resource measure for an arbitrary $m$ mode state $|\psi \rangle$ is the logarithm of the rank of its minimal decomposition into free states. That is, $\mathcal {R}(|\psi \rangle )=\log k$, where $k$ is the smallest positive integer such that $|\psi \rangle =\sum _{i=1}^k c_i |\alpha _i^{(1)}, \dots, \alpha _i^{(m)}\rangle$ (logarithm base is unimportant, but we will use base 2 here). Properties (i) and (ii) are immediately satisfied by this definition. In addition, it is easy to check this satisfies subadditivity, $\mathcal {R}(|\psi \rangle \otimes |\phi \rangle ) \le \mathcal {R}(|\psi \rangle ) + \mathcal {R}(|\phi \rangle )$, and in particular, $\mathcal {R}(|\psi \rangle \otimes |\phi \rangle ) = \mathcal {R}(|\psi \rangle )$ for $|\phi \rangle \in \mathcal {F}$.

We emphasize that our focus in this work is primarily on resource quantification and not on state conversion. This allows us to restrict to pure states and free unitaries, while acknowledging the fact that free unitaries cannot consume any resource, see Supplement 1F. A similar notion related to our coherent rank resource framework is the so called “degree of non-classicality,” which has previously been studied in, e.g., Refs. [55,67].

Moreover, the coherent state rank is closely related to the notion of coherence rank in the resource theory of coherence [68]. For example, if the coherent states in the expansion of a quantum state are nearly orthogonal, then the two are equivalent. However, generically, coherent states are not orthogonal and therefore a direct connection to resource theory of coherence is slightly more involved. The quantification of coherence in infinite-dimensional systems is not a new problem. Two approaches are quite common here. The first assumes Fock states as the “incoherent basis" and quantifies the coherence of coherent states with respect to this, e.g., Refs. [69,70]. However, this has the disadvantage of being at odds with the traditional notion of non-classicality in quantum optics, where coherent states are the most classical, while Fock states are admittedly non-classical [34]. The second approach defines coherent states as the preferred incoherent basis. This includes the resource theory of superposition [71] and subsequent works [72], including an operational framework [73]. Our resource theory framework falls under the so-called passive linear optics framework and the operational characterizations introduced in the references above are broadly applicable to our work.

In the present work, we consider the coherent state rank as the quantity of interest, since this directly determines the classical simulation complexity. We add one comment here however, that since the decomposition of an arbitrary state $|\psi \rangle$ to a coherent state basis is typically not exact for finite $k$ [e.g., see Theorem 1 and Eq. (3)], we will allow the above resource theory to hold approximately, in the sense of Definition 1.

We can contrast the coherent state rank resource theory to another related resource theory, namely one using the monotone of the Wigner logarithmic negativity [74], defined as $\mathcal {W}(|\psi \rangle ) = \log \int d\vec {\alpha } |W_{|\psi \rangle } (\vec {\alpha })|$, where $\vec {\alpha } = (\alpha _1, \dots, \alpha _m)$ for an $m$-mode pure state $|\psi \rangle$. In this framework, Gaussian states (including coherent states) are free, as the Wigner function is everywhere non-negative. First note that this implies $\mathcal {W}(|\psi \rangle )=0=\mathcal {R}(|\psi \rangle )$ for $|\psi \rangle \in \mathcal {F}$. Moreover, in Fig. 6, we plot both resource measures for Fock states, which shows the two definitions are consistent (in fact, are related approximately by a constant factor). Here we use the results of Theorem 1 that shows one can approximate $|n\rangle$ using $n+1$ coherent states. However, it is also clear the two definitions will not generally align with one another, in the sense one can find states $|\psi \rangle, |\phi \rangle$ such that $\mathcal {W}(|\psi \rangle )<\mathcal {W}(|\phi \rangle )$, but $\mathcal {R}(|\psi \rangle )>\mathcal {R}(|\phi \rangle )$, or vice versa (this is easy to show in the case where we consider exact decompositions, but harder to see when we allow approximate states, due to the difficulty of proving the optimality of the approximate rank $k$ of a particular state decomposition, though we expect the statement to remain true).

 

Fig. 6. Resource monotones for Fock states. We plot the measures of resource for the Fock states $|n\rangle$, using the coherent state rank definition, $\mathcal {R}$, and the Wigner negativity resource (with logarithm base 2). We also perform a curve fit (dash line), which shows the Wigner resource is approximated by $\mathcal {W}(n)\approx 0.48 \mathcal {R}(n)$. In Supplement 1E, we demonstrate by exact calculation that indeed the result is only approximate and likely the similar growth is due to the fact $W_n(\alpha )=W_n(|\alpha |)$ has $n$ roots in $|\alpha |$.

Download Full Size | PDF

Similarly, while the coherent rank decomposition is related to the stellar rank via Proposition 1, we wish to make it clear these are also distinct theories. In the context of a resource theory, the stellar rank $r$ is a genuine resource monotone where Gaussian states and operations are free ($r=0$), and simulation complexity scales as $O(2^r)$ (Section 2.3). By Proposition 1, one can see the approximate coherent rank resource measure

In terms of simulation complexity, the classical requirements for Fock input state $|1\rangle ^{\otimes n}|0\rangle ^{\otimes (m-n)}$ under linear optics and with coherent basis measurements scales as $O(2^n)$ in both the stellar and coherent rank frameworks, though the former scales as $O(4^n)$ for Fock basis measurements (and similarly for linear combinations of Gaussians in Section 2.2.1). Indeed, it is easy to find differences. As another example, the state $|n\rangle |0\rangle ^{m-1}$ can be simulated in a time linear in $n$ within the coherent state framework, whereas this is a state with stellar rank $n$ and simulation complexity therefore scaling as $2^n$ in the stellar rank paradigm (Section 2.3). Another key example is that the cat state $|ix\rangle -|-ix\rangle$ has coherent rank exactly 2, whereas the stellar rank is infinite [47]. However, Gaussian Boson sampling is typically much more efficient in the stellar framework (though as discussed below, for small enough squeezing or low enough target fidelity, the two can be equivalent, but when the squeezing is large, the stellar formalism will be most efficient).

Nevertheless, in all three cases, Wigner logarithmic negativity, coherent state rank, and the stellar rank, the resource measure, $\mathcal {W}$, $\mathcal {R}$ and $r$ relate to the simulation cost. For coherent state rank and stellar rank, the cost is exponential, scaling as $2^{\mathcal {R}, r}$, for free operations. For the Wigner negativity, as discussed in Supplement 1A, the simulation complexity scales as $4^{\mathcal {W}}$, although this has additional difficulties/errors as it is a sampling approach.

5. Discussion

We have developed a paradigm for the simulation of quantum optics, by decomposition to coherent states, in a manner akin to stabilizer state decompositions in finite-dimensional systems. In this framework, linear optical operations are free, as are computing amplitudes in the Fock or coherent state bases. For Boson sampling type simulations, this can lead to exponential improvements over naive Fock space simulation, with a simulation cost comparable to methods relying upon the permanent.

In Table 2 we summarize some of the complexity results for a few examples using the finite coherent state rank decomposition discussed in this work. Starting with an initial state of $n$ photons, we document the cost for storing this in memory (space), and the update time for a $m$-mode LO operation (e.g., random Boson sampling circuit). Here, the measurement column is the time cost to compute the probability of a particular measurement outcome, in the Fock or coherent state basis. All results in the table should be interpreted in the “big O” style, where, e.g., lower-order terms and constant factors related to implementation details are hidden. In the fourth row down, $|0.882\rangle$ refers to a squeezed vacuum state (e.g., Gaussian Boson sampling [GBS]) with squeezing parameter $\zeta =0.882$, corresponding to mean photon number $\langle \hat {n}\rangle \approx 1$. We find such a state can be approximated to fidelity above 0.99 using eight coherent states, which is used here (for lower/higher target fidelity, fewer/more states are needed in general). Also note, for larger (smaller) squeezing, typically more (fewer) coherent states will be required in the decomposition.

Table 2. Complexity Overview for $n$ Photon Simulations

View Table | View all tables in this article

In all cases in Table 2, the simulation requirements here can be much less than the equivalent full Fock space simulation (Section 2.4), which scales as the full Fock dimension in the worst case (such as in Boson sampling). Notably, the simulation cost here is only quadratic in the number of modes, even if all modes of an $n$-photon system are populated. This allows us to represent such systems significantly more efficiently, based on the decomposition used.

The method introduced is trivially parallelizable, since different compute nodes can store a subset of the terms in a coherent state decomposition and update them independently under linear optical operations. Moreover, this method potentially allows one to construct a trade-off between simulation complexity and accuracy, e.g., by reducing the number of terms in a coherent state decomposition, at the expense of fidelity. In the above table for example, we use eight coherent states to approximate a particular squeezed state to fidelity 0.99, however, fewer terms are required in general for a lower target fidelity; for this case, $\zeta =0.882$, only two coherent states are required to achieve fidelity above 0.9, yielding a cost to simulate GBS scaling as $2^n$. This allows a heuristic approach to simulating Boson sampling up to a target fidelity.

In addition, this construction allows for a tensor-network-like approach, potentially resulting in memory (and time) saving, in circuits with structure [9,11]. For example, in systems that have only a few entangling operations across a particular partition, we can split the state as a sum of tensor product states across this partition. Starting with a product state, one can always write $|\psi \rangle = |\psi _L\rangle \otimes |\psi _R\rangle$. Let us assume for simplicity equal partitions, so each partition size (number of modes) is half of the full one (i.e., dimension $\sqrt {d}$ instead of $d$). Instead of storing one state over the entire Hilbert space $O(d)$, we use the structure to store two states, but each only over their respective half partition, $O(\sqrt {d})$. In the coherent rank picture for say Boson sampling with $n$ photons, the cost to store this in memory is only $O(2^{n/2})$, instead of $O(2^n)$ for a state over the full Hilbert space. Then, any operation acting only on the left or right half can be applied trivially. For entangling operations (e.g., two-mode operations acting across the partition), if they can be decomposed as $\hat {U} = \sum _{i=1}^kc_i \hat {u}_i \otimes \hat {v}_{i}$, where $\hat {u}_i, \hat {v}_i$ are free (e.g., displacements) acting on each partition, then the state becomes a superposition of size $k$, but preserving the partition structure in each term. In the coherent rank decomposition picture, writing $|\psi _L\rangle = \sum _{m=1}^{k_L} a_m |\vec {\alpha }_m\rangle$, $|\psi _R\rangle = \sum _{n=1}^{k_R} b_n |\vec {\beta }_n\rangle$, we initially just store $k_L + k_R$ states in memory [here, $|\vec {\alpha }_m\rangle$ ($|\vec {\beta }_n\rangle$) is only defined over the modes in the left (right) partition]. Upon application of $\hat {U}$ across the partition, we get