Introduction

The extended Church-Turing thesis (ECT) states that every problem that can be efficiently computable with real physical devices are efficiently simulated with a Turing machine. It is expected that quantum computers would refute ECT by exploiting its inherent quantum supremacy. However, since scalable universal quantum computers that can perform actual quantum algorithms are not likely to be built in the foreseeable future, they are not “real physical devices” yet.

Boson sampling (BS)1 was introduced to defeat ECT with more feasible quantum devices, i.e., the linear optical network (LON) implementation. BS is considered a non-universal quantum computer with multi-photons in the multimode optical network. Aaronson and Arkhipov1 claimed that the transition amplitudes with no more than one photon per mode becomes hard to simulate with classical computers as the system scale increases.

The computational hardness of BS is from the hardness of matrix permanents. The transition amplitude from a pre-selected input state to a post-selected output state is determined by the permanent of a submatrix of a unitary matrix U in the LON. When no more than one photon is in both all input and output modes of the system, the amplitude can be classically simulated with Ryser’s formula2 (the best known algorithm for computing permanents). We recall the definition of the permanent of an M-dimensional square matrix A (Per[A]):

$${\rm{Per}}[A]=\sum _{\overrightarrow{\sigma }\in S}\,\prod _{i=1}^{M}\,{A}_{i,{\sigma }_{i}},$$
(1)

where A ij are the entries of A and the set S includes all the permutations of (1, 2, …, M), \(\{\overrightarrow{\sigma }\}\) (\(\overrightarrow{\sigma }=({\sigma }_{1},\ldots {\sigma }_{M})\) is an N-dimensional vector). The brute force computation of a matrix permanent in Eq. (1) requires N! terms in summation and each term is composed of the products of N elements of the matrix. Even though Ryser’s algorithm2 can perform the calculation in O(2N−1N2) arithmetic operations (it can be optimized further by Gray code as O(2N−1N) operations), the number of operations still increases exponentially with N (it was shown in Valiant3 and Aaronson4 that the computation of permanent is a #P-hard problem). Glynn5,6 derived a different algorithm that has the same order of computational cost with that of Ryser’s. Even though Jerrum et al.7 suggested a polynomial-time approximation algorithm for the permanents of matrices with non-negative elements, there exists no algorithms for arbitrary matrices that are more efficient than Ryser’s and Glynn’s yet. On the other hand, there have been some efforts in developing randomized algorithms for the permanents. Gurvits used Glynn’s formula to design a randomized algorithm8, and Aaronson and Hance9 generalized Gurvits’s sampling algorithm for matrices with either of repeated columns or repeated rows. A more generalized algorithm for matrices with repeated rows and columns, which can estimate the complexity of Fock state BS with multiple photons in both input and output modes, is introduced by Yung et al.10. When these algorithms are randomized, they estimate the matrix permanent with additive errors in polynomial runtimes. In this paper, we propose another generalized algorithm for matrices with repeated rows and columns. It is achieved by exploiting a series expansion of a product of variables regarding the linear combinations of variables11.

The classical minimal runtimes (\({{\mathscr{T}}}_{min}\)) of the algorithms mentioned above have interesting mathematical features, which render the algorithm to be related to a more general viewpoint of quantum complexity. The first observation is that the algorithm we derived here has the same \({{\mathscr{T}}}_{min}\) as that of the formula in Yung et al.10, such as Ryser’s and Glynn’s have the same \({{\mathscr{T}}}_{min}\). Considering the two algorithms arise from very different mathematical structures, we can regard that the obtained runtime is a rigorous criterion for the computational complexity of Fock state BS. The second obseration is that the functional form of \({{\mathscr{T}}}_{min}\) contains a summation of elementary symmetric polynomials. They have an intimate functional relation with the recently introduced coherence monotones, the coherence rank and generalized coherence concurrence12,13,14. This motivates us to define the generalized Fock state concurrence for a given state \(|\overrightarrow{n}\rangle \), which consists of the Fock state k-concurrence denoted by \({C}_{k}(\overrightarrow{n})\) with 0 ≤ k ≤ N, and the Fock state concurrence sum (the summation of \({C}_{k}(\overrightarrow{n})\) from k = 0 to k = N and denoted by \({C}_{S}(\overrightarrow{n})\)). The Fock state concurrence sum \({C}_{S}(\overrightarrow{n})\) is directly related to the amount of \({{\mathscr{T}}}_{min}\). We can state that the increase of \({C}_{S}(\overrightarrow{n})\) results in larger computational complexity, or \({C}_{S}(\overrightarrow{n})\) is a quantum resource for the complexity of the given system.

The concept of Fock state concurrence can also be compared with the Boltzmann entropy of the elementary quantum complexity \({S}_{B}^{q}\) introduced in Chin et al.15, which naturally emerges from the additive error bound for the approximated permanent estimator. By encompassing entropy and concurrence, our suggestion in this paper would provide the foundation for the quantum resource theory of linear optical quantum computing. In other words, by understanding the role of these quantum measures, we could find the origin of the quantum supremacy in quantum linear optics.

Results

First, the generalized Fock state concurrence and the concurrence sum are defined, and their physical relation with the generalized coherence concurrence14 is explained. Then an algorithm for computing the transition amplitudes of Fock state BS with multiple photons in input/output modes is proposed. By analyzing the minimal runtime of three algorithms (including ours) for computing the transition amplitudes of Fock state BS, we show that the Fock state concurrence sum is a quantum resource that determines the complexity of a given Fock state BS system.

The generalized Fock state concurrence family and the concurrence sum

Many theoretical analyses support the belief that quantum computers can perform some tasks faster than classical computers. Accordingly, it has been of particular interest to find the resources required for the quantum speedup. It is believed that entanglement is a critical resource for universal quantum computers16,17,18; however, the efficiency does not simply depend on the amount of entanglement19. It is also recently shown that the original Grover algorithm monotonically consumes coherence during the searching process13,20. There have been attempts to approach the problem in the reverse direction as well, i.e., to find conditions for a quantum system not to have any speedup. It is rigorously shown that nonnegative probability quasi-distributions (PQD) result in no quantum speedup21,22,23.

In the case of BS, the photon indistinguisability is considered the origin of the computational complexity in the Fock state BS1,24, and the degree of complexity is closely related to the majorization of the input-output photon distributions15. Whereas Rahimi-Keshari et al.25 approached this problem from the perspective of quasi-probability distributions (QPD), showing that the negativity of probability quasi-distribution (PQD) of linear optical networks is the necessary resource for the complexity.

In this section, we define a quantum measure from the multi-photon distribution patterns in multimode optical systems, the generalized Fock state concurrence and the Fock state concurrence sum. The generalized Fock state concurrence is a quantity analogous to the generalized entanglement concurrence26 and generalized coherence concurrence14. It will be shown in the later section that the Fock state concurrence sum becomes a resource that determines the complexity of Fock state BS.

Definitions

In the linear optical network of M optical modes into which N photons are injected, the Fock state vector is written as

$$|\overrightarrow{n}\rangle =|{n}_{1},{n}_{2},\ldots ,{n}_{M}\rangle ,\,\sum _{i}^{M}\,{n}_{i}=N$$
(2)

where n i represents the photon number for the ith mode (it is worth emphasizing that n i can be greater than 1 for our later discussion on the generalized Fock state BS). Then the coherence rank and k-concurrence of a Fock state \(|\overrightarrow{n}\rangle \) is defined as follows:

Definition 1. 

The Fock state coherence rank for a given Fock state \(\overrightarrow{n}\) is defined as the integer \({\alpha }_{\overrightarrow{n}}\), the number of nonzero elements for the particle distribution vector \(\overrightarrow{n}\).

Definition 2. 

The Fock state k-concurrence for a given Fock state \(\overrightarrow{n}\) is defined with the elementary symmetric polynomial as

$${C}_{k}(\overrightarrow{n})\equiv {(\frac{1}{(\genfrac{}{}{0ex}{}{N}{k})}{X}_{k}(\overrightarrow{n}))}^{\frac{1}{k}},$$
(3)

where \({X}_{k}(\overrightarrow{n})\) (kth elementary symmetric polynomial) is defined as

$${X}_{k}(\overrightarrow{n})=\sum _{{i}_{1} < {i}_{2} < \cdots < {i}_{k}=1}^{{\alpha }_{\overrightarrow{n}}}\,{n}_{{i}_{1}}{n}_{{i}_{2}}\cdots {n}_{{i}_{k}},\,(0\le k\le {\alpha }_{\overrightarrow{n}}\le M)$$
(4)

(we define \({X}_{0}(\overrightarrow{n})=1\), and \({X}_{1}(\overrightarrow{n})=N\) for any \(\overrightarrow{n}\)).

This is normalized so that \({C}_{k}(\overrightarrow{n})\) becomes 1 when \(\overrightarrow{n}\) is maximally coherent, i.e., \(\overrightarrow{n}=({\overrightarrow{1}}_{N},{\overrightarrow{0}}_{M-N})\). The Fock state k-concurrences \({C}_{k}(\overrightarrow{n})\) from 0 ≤ k ≤ N constitute the generalized Fock state concurrence family.

Definition 3. 

The Fock state concurrence sum for a given Fock state \(\overrightarrow{n}\) is defined as

$${C}_{S}(\overrightarrow{n})\equiv \frac{1}{{2}^{N}}\,\sum _{k=0}^{{\alpha }_{\overrightarrow{n}}}\,(\genfrac{}{}{0ex}{}{N}{k}){[{C}_{k}(\overrightarrow{n})]}^{k},$$
(5)

The factor 1/2N is multiplied for the normalization, i.e., \({C}_{S}(\overrightarrow{n})=1\) when \(\overrightarrow{n}\) is maximally coherent. Since X k are all Schur concave functions, which decrease as the Fock state vector \(\overrightarrow{n}\) is more majorized27, the k-concurrences are also Schur concave functions. Hence the concurrence sum \({C}_{S}(\overrightarrow{n})\) is also a Schur concave function.

To calculate the Fock state k-concurrence (and concurrence sum) of those states which are not expressed with a single photon distribution vector, we need to consider more comprehensive definitions than Definition 2. It can be achieved in a similar manner to the generalized concurrences of entanglement and coherence14,26, the situation is slightly different for our case though. The well-known convex-roof extention (see, e.g., Eltschka et al.28) is not exactly suitable here, for Definition 2 does not embrace the pure states that are superpositions of photon distribution vectors, i.e., when \(|\psi \rangle ={\sum }_{\overrightarrow{n}}\,{\psi }_{\overrightarrow{n}}|\overrightarrow{n}\rangle \) (\({\sum }_{\overrightarrow{n}}\,|{\psi }_{\overrightarrow{n}}{|}^{2}=1\)). Hence, we need two steps of extension for the generalized Fock state concurrence family:

Definition 4. 

The Fock state k-concurrence of a pure state \(|\psi \rangle ={\sum }_{\overrightarrow{n}}\,{\psi }_{\overrightarrow{n}}|\overrightarrow{n}\rangle \) is defined as

$${C}_{k}(|\psi \rangle )\equiv \sum _{\overrightarrow{n}}\,|{\psi }_{\overrightarrow{n}}{|}^{2}{C}_{k}(|\overrightarrow{n}\rangle ),$$
(6)

and the Fock state concurrence sum of |ψ〉 as

$${C}_{S}(|\psi \rangle )\equiv \sum _{\overrightarrow{n}}\,|{\psi }_{\overrightarrow{n}}{|}^{2}{C}_{S}(|\overrightarrow{n}\rangle \mathrm{).}$$
(7)

The Fock state k-concurrence of a mixed state ρ, which can be pure-state-decomposed as \(\rho ={\sum }_{a}\,{\rho }_{a}|{\psi }_{a}\rangle \langle {\psi }_{a}|\), is defined with the convex roof extension as

$${C}_{k}(\rho )\equiv \mathop{{\rm{\min }}}\limits_{\{{\rho }_{a},|{\psi }_{a}\rangle \}}\,{\rho }_{a}{C}_{k}(|{\psi }_{a}\rangle ),$$
(8)

and the Fock state concurrencce sum of ρ as

$${C}_{S}(\rho )\equiv \mathop{{\rm{\min }}}\limits_{\{{\rho }_{a},|{\psi }_{a}\rangle \}}\,{\rho }_{a}{C}_{S}(|{\psi }_{a}\rangle \mathrm{).}$$
(9)

Comparison of the Fock state concurrence with the single particle coherence concurrence

We can explain the intuitive relation between the Fock state coherence and the single photon coherence, which will clarify our concept of the Fock state concurrence.

The coherence as one of the fundamental non-classicalities is originated from the framework of superposition12,29, i.e., the partition of probability among several states for one quantum system. Since coherence is basis-dependent, we first need to fix a computational basis set. The quantification of coherence is possible under a given normalized basis set \(\{|i\rangle {\}}_{i=1}^{d}\), and we can state that a pure state is coherent in the basis set if and only if

$$|\psi \rangle =\sum _{i=1}^{k > 1}\,{\psi }_{i}|i\rangle \mathrm{.}$$
(10)

When k = 1, |ψ〉 is incoherent (the mixed state extension of coherence is straightforward. See Baumtratz et al.30). The degree of coherence is determined by the probability amplitude of the state, i.e.,

$$P(|\psi \rangle )=(|{\psi }_{1}{|}^{2},|{\psi }_{2}{|}^{2},\ldots ,|{\psi }_{d}{|}^{2})\,(\sum _{i=1}^{d}\,|{\psi }_{i}{|}^{2}=1).$$
(11)

The concept of majorization plays a crucial role here (for two nonincreasing real vectors \(\overrightarrow{x}\) and \(\overrightarrow{y}\) of dimension d, we state that \(\overrightarrow{x}\) is majorized by \(\overrightarrow{y}\) (or \(\overrightarrow{x}\prec \overrightarrow{y}\)) if and only if \({\sum }_{i=1}^{k}\,{x}_{i}\le {\sum }_{i=1}^{k}\,{y}_{i}\) for all k < d and \({\sum }_{i=1}^{d}\,{x}_{i}={\sum }_{i=1}^{d}\,{y}_{i}\)27). Indeed, for two pure states |ψ〉 and |ϕ〉, the following relation holds31:

$$\begin{array}{l}P(|\psi \rangle )\,{\rm{is}}\,{\rm{majorized}}\,{\rm{by}}\,P(|\varphi \rangle )\\ \begin{array}{rcl} & \iff & |\psi \rangle \,{\rm{is}}\,{\rm{more}}\,{\rm{coherent}}\,{\rm{than}}\,|\varphi \rangle ,\,{\rm{i}}{\rm{.e}}{\rm{}}\,.,\\ & & |\psi \rangle \,{\rm{can}}\,{\rm{be}}\,{\rm{transformed}}\,{\rm{to}}\,|\varphi \rangle \,{\rm{with}}\,{\rm{incoherent}}\,{\rm{operations}}\,({\rm{IO}}){\rm{.}}\end{array}\end{array}$$
(12)

Two extremal cases are when P(|ψ〉) = perm[(1, 0, …, 0)] (perm[\(\overrightarrow{v}\)] denotes any permutation vector of \(\overrightarrow{v}\)) and \(P(|\psi \rangle )=\frac{1}{\sqrt{d}}\mathrm{(1},1,\ldots ,\mathrm{1)}\). The state is incoherent for the former and maximally coherent for the latter. There are several coherence measures that satisfy (12) (see, e.g., Streltsov et al.32 for some examples).

One specific example is the d-slit experiment of a photon (Fig. 1(a)). When each slit is well-separated from the others, the photon state that passes through the slit is represented in the computational basis set {|ψ1〉, |ψ2〉, …, |ψ d 〉} (〈ψ i |ψ j 〉 = δ ij ), where |ψ i 〉 corresponds to the case when the photon passes the i-th slit. Then the photon state is given by

$$|{\rm{\Psi }}\rangle =\sum _{i=1}^{d}\,{c}_{i}|{\psi }_{i}\rangle \,(\sum _{i}\,|{c}_{i}{|}^{2}=1).$$
(13)

Therefore, we can state that the coherence of |ψ〉 is determined by

$$P(|{\rm{\Psi }}\rangle )=(|{c}_{1}{|}^{2},|{c}_{2}{|}^{2},\ldots ,|{c}_{d}{|}^{2}).$$
(14)
Figure 1
figure 1

(a) A single photon that passes through a well-separated multi-slit. (b) Multi-photons injected to a multi-mode linear optical network system U.

Full size image

When |c i |2 = 1 for some i, the state is incoherent and passes through the i-th slit deterministically. This state represents the particle-like property of the photon. On the other hand, when \(|{c}_{i}{|}^{2}=1/\sqrt{d}\) for all i, the state is maximally coherent and represents the wave-like property. The analysis can be generalized to the mixed state case by attaching a detector at the slit14,33. The coherence k-concurrence \({C}_{c}^{(k)}\) for a pure state |Ψ〉 is given by14

$${C}_{c}^{(k)}(|{\rm{\Psi }}\rangle )=d{(\frac{1}{(\genfrac{}{}{0ex}{}{d}{k})}\sum _{{i}_{1} < {i}_{2} < \cdots < {i}_{k}=1}^{d}{X}_{k}[P(|{\rm{\Psi }}\rangle )])}^{\frac{1}{k}},$$
(15)

where X k [P(|Ψ〉)] is the k-th elementary symmetric polynomial of P(|Ψ〉).

In the multi-photon case as to the Fock state BS, the probablity for each photon to be in a specific mode is 1. On the other hand, the probability distribution of the initial N-photon in M-modes are given by

$$P(|\overrightarrow{n}\rangle )=(\frac{{n}_{1}}{N},\frac{{n}_{2}}{N},\ldots ,\frac{{n}_{M}}{N})$$
(16)

for the photon distribution vector \(|\overrightarrow{n}\rangle =({n}_{1},{n}_{2},\ldots ,{n}_{M})\) (Fig. 1(b)). Therefore, Definition 2 of Fock state coherence k-concurrence is analogous to Eq. (15) with the replacement of (d, |Ψ〉) with (N, \(|\overrightarrow{n}\rangle \)).

The Fock state concurrence sum C S and the complexity of Fock state BS

In this section, we show that the Fock state concurrence sum C S defined in the former section plays a crucial role in the complexity of Fock state BS. To see the relation, we first derive a generalized algorithm for computing the transition amplitudes of Fock state BS with multiple photons in input/output modes. An intriguing fact about our algorithm is that the minimal runtime \({{\mathscr{T}}}_{{\rm{\min }}}\) for the algorithm is equal to that of another generalized algorithm presented in Yung et al.10. Furthermore, the functional form of \({{\mathscr{T}}}_{{\rm{\min }}}\) explicitly contains the Fock state concurrence sum C S . This implies that C S is a quantum measure that determines the computational complexity of a generalized Fock state BS system.

The derivation of a generalized algorithm for matrix permanents

In the linear optical network of M optical modes characterized by a unitary transformation \(\hat{U}\), the photon creation and annihilation operators \({\hat{a}}_{i}^{\dagger }\) and \({\hat{a}}_{i}\) in the i-th mode (i = 1, …, M) rotate under the acton of \(\hat{U}\) as

$$\hat{U}{\hat{a}}_{i}^{\dagger }{\hat{U}}^{\dagger }=\sum _{j=1}^{M}\,{U}_{ij}{\hat{a}}_{j}^{\dagger },\,\hat{U}{\hat{a}}_{i}{\hat{U}}^{\dagger }=\sum _{j=1}^{M}\,{U}_{ij}^{\ast }{\hat{a}}_{j}.$$
(17)

Scheel34 showed that the transition amplitude between the two Fock states \(|\overrightarrow{n}\rangle \) (= \(|{n}_{1},{n}_{2},\ldots ,{n}_{M}\rangle \), \({\sum }_{i}^{M}\,{n}_{i}=N\)) and \(|\overrightarrow{m}\rangle \) (= \(|{m}_{1},{m}_{2},\ldots ,{m}_{M}\rangle \), \({\sum }_{i}^{M}\,{m}_{i}=N\)) is proportional to the matrix permanent:

$$\langle \overrightarrow{m}|\hat{U}|\overrightarrow{n}\rangle =\frac{{\rm{Per}}([U{]}_{\overrightarrow{n},\overrightarrow{m}})}{\prod _{k=1}^{M}\,\sqrt{{n}_{k}!{m}_{k}!}},$$
(18)

where \({[U]}_{\overrightarrow{n},\overrightarrow{m}}\) is an N × N submatrix of U, which has n i of the i-th rows of U and m j of j-th row of U. \(N={\sum }_{i}^{M}\,{n}_{i}={\sum }_{i}^{M}\,{m}_{i}\) is the total number of photons.

The relation (18) holds for the arbitrary square complex matrix A35, i.e.,

$$\langle \overrightarrow{m}|\hat{A}|\overrightarrow{n}\rangle =\frac{{\rm{Per}}([A{]}_{\overrightarrow{n},\overrightarrow{m}})}{\prod _{k=1}^{M}\,\sqrt{{n}_{k}!{m}_{k}!}},$$
(19)

where \(\hat{A}\) can be expressed as \(\hat{A}=\exp \,[{\sum }_{i,j}\,{\hat{a}}_{i}^{\dagger }{(\mathrm{ln}{A}^{T})}_{ij}{\hat{a}}_{j}]\). This implies that a matrix permanent is obtained by calculating the corresponding transition amplitude between the given input-output Fock states. We will exploit this relation to obtain a generalized algorithm for matrix permanents.

The following lemma11 is useful for deriving our algorithm.

Lemma 1. 

For a vector \(\overrightarrow{x}=({x}_{1},{x}_{2},\ldots ,{x}_{M})\), the following identity holds:

$${x}_{1}^{{n}_{1}}\,\cdots \,{x}_{M}^{{n}_{M}}=\frac{1}{N{!2}^{N}}\sum _{{v}_{1}=0}^{{n}_{1}}\cdots \sum _{{v}_{M}=0}^{{n}_{M}}\,{(-1)}^{{N}_{v}}(\genfrac{}{}{0ex}{}{{n}_{1}}{{v}_{1}})\cdots (\genfrac{}{}{0ex}{}{{n}_{M}}{{v}_{M}}){[\sum _{i=1}^{M}({n}_{i}-2{v}_{i}){x}_{i}]}^{N},$$
(20)

where \({N}_{v}={\sum }_{i=1}^{M}\,{v}_{i}\) and \(N={\sum }_{i=1}^{M}\,{n}_{i}\).

A detailed proof is given in Kan11. The isomorphism between quantum states and multivariate polynomials from Theorem 3.6 of Aaronson and Arkhipov1 (see also Yung et al.10) connects the above lemma with the linear optical quantum system, which results in the following identity:

Theorem 1. 

There exists a generalized formula for the matrix permanent with repeated rows and columns that is expressed as

$${\rm{P}}{\rm{e}}{\rm{r}}({[A]}_{\overrightarrow{n},\overrightarrow{m}})=\frac{1}{{2}^{N}}\,\sum _{{v}_{1}=0}^{{n}_{1}}\cdots \sum _{{v}_{M}=0}^{{n}_{M}}\,{(-1)}^{{N}_{v}}(\genfrac{}{}{0ex}{}{{n}_{1}}{{v}_{M}})\cdots (\genfrac{}{}{0ex}{}{{n}_{M}}{{v}_{M}})\,\prod _{j=1}^{M}\,{[\sum _{i=1}^{M}({n}_{i}-2{v}_{i}){A}_{ij}]}^{{m}_{j}}.$$
(21)

Proof.

From Theorem 3.6 of Aaronson and Arkhipov1, the Fock state \(|\overrightarrow{n}\rangle ={\otimes }_{k=1}^{M}\,{({\hat{a}}_{k}^{\dagger })}^{{n}_{k}}/\sqrt{{n}_{k}!}\mathrm{|0}\rangle \) can be expanded using Lemma 1 as

$$|\overrightarrow{n}\rangle =\frac{1}{(\prod _{k=1}^{M}\,\sqrt{{n}_{k}!})N{!2}^{N}}\,\sum _{{v}_{1}=0}^{{n}_{1}}\cdots \sum _{{v}_{M}=0}^{{n}_{M}}\,{(-1)}^{{N}_{v}}(\genfrac{}{}{0ex}{}{{n}_{1}}{{v}_{1}})\cdots (\genfrac{}{}{0ex}{}{{n}_{M}}{{v}_{M}})\,{[\sum _{i=1}^{M}({n}_{i}-2{v}_{i}){\hat{a}}_{i}^{\dagger }]}^{N}|\overrightarrow{0}\rangle .$$
(22)

By substituting Eq. (22) into Eq. (19), we have

$${\rm{P}}{\rm{e}}{\rm{r}}({[A]}_{\overrightarrow{n},\overrightarrow{m}})=\frac{1}{N{!2}^{N}}\,\sum _{{v}_{1}=0}^{{n}_{1}}\cdots \sum _{{v}_{M}=0}^{{n}_{M}}\,{(-1)}^{{N}_{v}}(\genfrac{}{}{0ex}{}{{n}_{1}}{{v}_{1}})\cdots (\genfrac{}{}{0ex}{}{{n}_{M}}{{v}_{M}})\langle 0|{\hat{a}}_{1}^{{m}_{1}}\,\cdots \,{\hat{a}}_{M}^{{m}_{M}}{[\sum _{ij}({n}_{i}-2{v}_{i}){A}_{ij}{\hat{a}}_{j}^{\dagger }]}^{N}|0\rangle .$$
(23)

On the other hand, \({[{\sum }_{ij}({n}_{i}-2{v}_{i}){A}_{ij}{\hat{a}}_{j}^{\dagger }]}^{N}\) is expanded using the multinomial formula

$${(\sum _{i=1}^{M}{y}_{i})}^{N}=\sum _{{\sum }_{i=1}^{M}{s}_{i}=N}\,\frac{N!}{\prod _{i}\,{s}_{i}!}{y}_{1}^{{s}_{1}}\,\cdots \,{y}_{M}^{{s}_{M}},$$
(24)

as

$${[\sum _{ij}({n}_{i}-2{v}_{i}){A}_{ij}{\hat{a}}_{j}^{\dagger }]}^{N}=\sum _{{\sum }_{i=1}^{M}{s}_{i}=N}\,\frac{N!}{\prod _{i}\,{s}_{i}!}\prod _{j}\,{[\sum _{i}({n}_{i}-2{v}_{i}){A}_{ij}{\hat{a}}_{j}^{\dagger }]}^{{s}_{j}}.$$
(25)

Therefore, by substituting Eq. (25) into Eq. (23), we obtain Eq. (21) with the identity\(\langle \mathrm{0|}{\hat{a}}_{1}^{{m}_{1}}\,\cdots \,{\hat{a}}_{M}^{{m}_{M}}{({\hat{a}}_{1}^{\dagger })}^{{s}_{1}}\,\cdots \) \({({\hat{a}}_{M}^{\dagger })}^{{s}_{M}}\mathrm{|0}\rangle ={\prod }_{k=1}^{M}({m}_{k}!{\delta }_{{m}_{k}{s}_{k}})\). □

With A = U (unitary operator), the above formula is for the computation of transition amplitudes with multiple photons in both input and output modes.

By exploiting the symmetry in Eq. (20) 11, the number of terms can be reduced to about a half of that in Eq. (21). First, when at least one of {n k } is an odd number, n1 can be chosen to be an odd number without loss of generality. Then Eq. 21 is simplified as

$${\rm{P}}{\rm{e}}{\rm{r}}\,({[A]}_{\overrightarrow{n},\overrightarrow{m}})=\frac{1}{{2}^{N-1}}\sum _{{v}_{1}=0}^{({n}_{1}-1)/2}\cdots \sum _{{v}_{M}=0}^{{n}_{M}}\,{(-1)}^{{N}_{v}}(\genfrac{}{}{0ex}{}{{n}_{1}}{{v}_{1}})\cdots (\genfrac{}{}{0ex}{}{{n}_{M}}{{v}_{M}})\prod _{j=1}^{M}[\sum _{i}({n}_{i}-2{v}_{i}){A}_{ij}{]}^{{m}_{j}}$$
(26)

Second, when all {n k } are even numbers, one still can reduce the number of terms in the summation by dividing the first summation of v1 into a summation from 0 to n1 − 1 and v1 = n1. Since n1 − 1 is an odd number, the same symmetry that is used for Eq. (26) reduces the number of terms. \({\prod }_{k}\,({n}_{k}+\mathrm{1)/2}\) terms are required for the first case and \(({\prod }_{k}\,({n}_{k}+\mathrm{1)}-1)/2\) terms for the second case.

Now, we show that our formula is reduced to that of Glynn’s5,6 when \(\overrightarrow{n}=\overrightarrow{m}=\mathrm{(1},\ldots ,\,\mathrm{1)}\) and N = M. In this case (n i  − 2v i ) is either +1 or −1 for all i, and all the binomial coefficients become 1. Then Eq. (26) can be expressed as

$${\rm{P}}{\rm{e}}{\rm{r}}\,(A)=\frac{1}{{2}^{N-1}}\,\sum _{\overrightarrow{x}\in {\{-1,1\}}^{N}}\,(\prod _{i=1}^{N}\,{x}_{i})\prod _{j=1}^{N}\,(\sum _{k=1}^{N}\,{A}_{jk}{x}_{k})$$
(27)

which is the Glynn’s formula (see Method).

As pointed out by Gurvits8,9, Glynn’s form can be interpreted as an ensemble average over the random vector whose entries are ±1. Likewise, one can interpret Eq. (21) as a randomized algorithm in which v k is randomly generated among (0, 1, …, n k ) with probability \(p({v}_{k})=(\genfrac{}{}{0ex}{}{{n}_{k}}{{v}_{k}})/{2}^{{n}_{k}}\). Accordingly, Eq. 21 is rewritten as

$${\rm{Per}}\,({[A]}_{\overrightarrow{n},\overrightarrow{n}})=\sum _{\overrightarrow{v}=\overrightarrow{0}}^{\overrightarrow{n}}\,p(\overrightarrow{v})G(\overrightarrow{v}),$$
(28)

where \(p(\overrightarrow{v})\equiv p({v}_{1})\,\cdots \,p({v}_{M})\) and \({\sum }_{\overrightarrow{v}=\overrightarrow{0}}^{\overrightarrow{n}}\,p(\overrightarrow{v})=1\), and \(G(\overrightarrow{v})\equiv {(-\mathrm{1)}}^{{N}_{v}}\,{\prod }_{j=1}^{M}\,{[{\sum }_{i}({n}_{i}-2{v}_{i}){A}_{ij}]}^{{m}_{j}}\). \(G(\overrightarrow{v})\) is evaluated for each random instance of \(\overrightarrow{v}\) with the probability \(p(\overrightarrow{v})\), and then the matrix permanent is approximated as an average,

$${\rm{Per}}\,({[A]}_{\overrightarrow{n},\overrightarrow{m}})\simeq \frac{1}{{N}_{{\rm{Sample}}}}\,\sum _{i=1}^{{N}_{{\rm{Sample}}}}\,G({\overrightarrow{v}}^{(i)}),$$
(29)

where NSample is the number of samples.

Minimal classical runtime

The runtime for the classical simulation of Eq. (21) is obtained by identifying all the summations included in the algorithm as

$${\mathscr{T}}\,(\overrightarrow{n},\overrightarrow{m})={\mathscr{O}}\,(\prod _{i=1}^{M}\,({n}_{i}+\mathrm{1)}{\alpha }_{\overrightarrow{n}}{\alpha }_{\overrightarrow{m}}),$$
(30)

where \({\alpha }_{\overrightarrow{n}}\) and \({\alpha }_{\overrightarrow{m}}\) are the number of nonzero elements of \(\overrightarrow{n}\) and \(\overrightarrow{m}\) respectively. \({\prod }_{i=1}^{M}\,({n}_{i}+\mathrm{1)}\) comes from \({\sum }_{{v}_{1}=0}^{{n}_{1}}\,\cdots \,{\sum }_{{v}_{M}=0}^{{n}_{M}}\), and \({\alpha }_{\overrightarrow{n}}{\alpha }_{\overrightarrow{m}}\) comes from \({\prod }_{j=1}^{M}\,{[{\sum }_{i}({n}_{i}-2{v}_{i}){A}_{ij}]}^{{m}_{j}}\) in Eq. (21).

On the other hand, we can expand \(|\overrightarrow{m}\rangle \) instead of \(|\overrightarrow{n}\rangle \) with Kan’s series expansion as in Eq. (22). From this input-output symmetry, we obtain another runtime

$${\mathscr{T}}{}^{{\rm{^{\prime} }}}(\overrightarrow{n},\overrightarrow{m})={\mathscr{O}}\,(\prod _{i=1}^{M}\,({m}_{i}+1){\alpha }_{\overrightarrow{n}}{\alpha }_{\overrightarrow{m}}).$$
(31)

We can choose the shorter one between \({\mathscr{T}}\) and \({\mathscr{T}}\,^{\prime} \) for the optimal classical simulation. Therefore, the minimal running time for the algorithm, denoted by \({{\mathscr{T}}}_{min}(\overrightarrow{n},\overrightarrow{m})\), is given by

$${{\mathscr{T}}}_{min}(\overrightarrow{n},\overrightarrow{m})={\mathscr{O}}\,(min(\prod _{i=1}^{M}\,({n}_{i}+1),\prod _{j=1}^{M}\,({m}_{j}+1)){\alpha }_{\overrightarrow{n}}{\alpha }_{\overrightarrow{m}}).$$
(32)

A special case is when both n i and m i are not bigger than 1 for all i. Then \({\alpha }_{\overrightarrow{n}}={\alpha }_{\overrightarrow{m}}=N\) and \({\mathscr{T}}\,(\overrightarrow{n},\overrightarrow{m})={2}^{N}{N}^{2}\), which is the same runtime as that of Ryser’s formula.

The minimal runtime for our algorithm can be compared to that of another generalized algorithm suggested in Yung et al.10. Interestingly enough, the minimal runtime for the algorithm is exactly equal to that of ours (the same thing happens when we compare the runtime for Ryser’s and Glynn’s formula). A brief explanation is presented in Methods. This phenomena is intriguing since these two algorithms appear from very different mathematical backgrounds. While our algorithm is constructed from a series expansion of collective variables, the algorithm in Yung et al. is a direct generalization of Aaronson and Hance’s algorithm9. Two algorithms created from two totally different paths have the same classical runtime, from which we can surmise that the minimal runtime \({{\mathscr{T}}}_{min}\) is a credible criterion for the computational complexity of the generalized Fock state BS.

The minimal runtime and the Fock state concurrence sum

Now we are ready to see the functional relation of \({{\mathscr{T}}}_{{\rm{\min }}}\) with the Fock state concurrence sum C S . Actually, this relation is easily observed by reexpressing \({{\mathscr{T}}}_{{\rm{\min }}}(\overrightarrow{n},\overrightarrow{m})\) by expanding \({\prod }_{i=1}^{M}\,({n}_{i}+\mathrm{1)}\) along the order of n i as

$$\begin{array}{rcl}\prod _{i=1}^{M}\,({n}_{i}+\mathrm{1)} & = & 1+\sum _{i}\,{n}_{i}+\sum _{{i}_{1} < {i}_{2}}\,{n}_{{i}_{1}}{n}_{{i}_{2}}+\cdots \\ & & +\sum _{{i}_{1} < {i}_{2} < \,\cdots \, < {i}_{M}}\,{n}_{{i}_{1}}{n}_{{i}_{2}}\,\cdots \,{n}_{{i}_{M}}\mathrm{.}\end{array}$$
(33)

From the definition of the elementary symmetric polynomial (see Eq. (4)), we have

$$\prod _{i=1}^{M}\,({n}_{i}+\mathrm{1)}=\sum _{k=0}^{{\alpha }_{\overrightarrow{n}}}\,{X}_{k}(\overrightarrow{n}\mathrm{).}$$
(34)

Note that the summation is until \(k={\alpha }_{\overrightarrow{n}}\) because \({X}_{k}(\overrightarrow{n})=0\) for \(k > {\alpha }_{\overrightarrow{n}}\). As a result, \({{\mathscr{T}}}_{min}(\overrightarrow{n},\overrightarrow{m})\) is rewritten as

$${{\mathscr{T}}}_{min}(\overrightarrow{n},\overrightarrow{m})={\mathscr{O}}\,(min(\sum _{k=0}^{{\alpha }_{\overrightarrow{n}}}\,{X}_{k}(\overrightarrow{n}),\sum _{l=0}^{{\alpha }_{\overrightarrow{m}}}\,{X}_{l}(\overrightarrow{m})){\alpha }_{\overrightarrow{n}}{\alpha }_{\overrightarrow{m}}).$$
(35)

By using Definition 3, the minimal runtime is finally rewritten as

$${{\mathscr{T}}}_{min}(\overrightarrow{n},\overrightarrow{m})={\mathscr{O}}({2}^{N}(min\,({C}_{S}(\overrightarrow{n}),{C}_{S}(\overrightarrow{m}))){\alpha }_{\overrightarrow{n}}{\alpha }_{\overrightarrow{m}}),$$
(36)

which is a composition of the Fock state concurrence sum and Fock state coherence rank. This expression shows that in linear optics the Fock state concurrence sum is a critical resource that determines the computational complexity. We should emphasize that –in so far as we know– this is the first evidence that the summation of all the family members of concurrence can operate as an independent resource. Most works on the generalized concurrence in entanglement and coherence have focused on the role of some specific member as the resource for practical quantum processes (see, e.g., Sents et al.36, Girard et al.37, Chin13 and Chin14). On the other hand, as we have just seen, the generalized coherence concurrence of Fock state acts as a whole in the multimode linear optical system. In other words, not an individual member C k but the summation of the whole members C S becomes the deterministic resource for the process we are interested in.

As an example, when \(\overrightarrow{n}=(N\mathrm{,0},\ldots ,\mathrm{0)}\), we have \({C}_{k}(\overrightarrow{n})=0\) for k ≥ 2 and \({C}_{S}(\overrightarrow{n})=\frac{1+N}{{2}^{N}}\) (the minimal concurrence sum), which results in \({{\mathscr{T}}}_{{\rm{\min }}}(\overrightarrow{n},\overrightarrow{m})={\mathscr{O}}\mathrm{[(1}+N){\alpha }_{\overrightarrow{m}}]\le {\mathscr{O}}\mathrm{[(1}+N)N]\). For this case the runtime becomes polymonial. As another example, when \(\overrightarrow{n}=\overrightarrow{m}=\mathrm{(1},\ldots ,1,\mathrm{0,}\,\ldots ,\mathrm{0)}\), we have \({C}_{S}(\overrightarrow{n})={C}_{S}(\overrightarrow{m})=1\) (the maximal concurrence sum) and \({{\mathscr{T}}}_{{\rm{\min }}}(\overrightarrow{n},\overrightarrow{m})={\mathscr{O}}{\mathrm{[2}}^{N}{N}^{2}]\).

Eq. (36) also reveals an intriguing property of \({{\mathscr{T}}}_{min}\), which contrasts with that of the additive error bound \( {\mathcal E} \) for an approximated permanent estimator. In Chin et al.15, the Boltzmann and Shannon entropy of elementary quantum complexity is introduced to evaluate the quantum complexity of the given quantum particle distributions. And \( {\mathcal E} \) is explicitly expressed as the difference between the Boltzmann entropy and Shannon entropy of elementary quantum complexity. On the other hand, the relation between the entropies and \({{\mathscr{T}}}_{min}\) is implicit and only can be intuitively explained. Eq. (36) indicates that the generalized Fock state concurrence is another criterion for the computational complexity of linear optical systems. We can state that both entropy and concurrence are crucial measures (or resources) that directly determines the quantum complexity of linear optical computers.

Before closing this section, it is worth mentioning the role of the extended definitions for general states (Definition 4). With such definitions, we can calculate the concurrence sum for arbitary states, including coherent states and thermal states, etc. As a simple example, when ρ is a thermal state, i.e.,

$$\rho ={\rho }^{th}=\sum _{N=0}^{{\rm{\infty }}}\,\sum _{{\sum }_{i}{n}_{i}=N}^{{\rm{\infty }}}\,(\prod _{i=1}^{M}\,\frac{{\bar{n}}_{i}^{{n}_{i}}}{{({\bar{n}}_{i}+1)}^{{n}_{i}+1}})|\overrightarrow{n}\rangle \langle \overrightarrow{n}|,$$
(37)

where \({\bar{n}}_{i}\) represents the average photon number for each i, we have

$${C}_{k}({\rho }^{th})=\sum _{N=0}^{{\rm{\infty }}}\,\sum _{{\sum }_{i}{n}_{i}=N}^{{\rm{\infty }}}\,(\prod _{i=1}^{M}\,\frac{{\bar{n}}_{i}^{{n}_{i}}}{{({\bar{n}}_{i}+1)}^{{n}_{i}+1}}){C}_{k}(|\overrightarrow{n}\rangle ).$$
(38)

Since \({C}_{k}(|\overrightarrow{n}\rangle ) < 1\) except when \(|\overrightarrow{n}\rangle \) is maximally coherent, it is easy to see that\({C}_{S}({\rho }^{th}) < {C}_{S}(|\overrightarrow{n}\rangle =\) \(|{\overrightarrow{1}}_{N},{\overrightarrow{0}}_{M-N}\rangle =1\). We can surmise that the reason why the thermal state BS is simpler to simulate classically than the original BS38 is that the former has less quantum resource, i.e., concurrence sum, than the latter. This viewpoint would be compared more rigorously to that of Rahime-Keshari et al.38 in the future. We expect that it would reveal the role of concurrence sum in the complexity of various BS systems, such as Gaussian BS38,39,40 and Vibronic BS41,42.

Discussion

We expect our research to develop into two aspects, which are closely relevant to each other. First, the relation of the generalized concurrence \({C}_{k}(\overrightarrow{n})\) with the exact classical simulation of matrix permanents has many similarities with that of the Boltzmann entropy \({S}_{B}^{q}(\overrightarrow{n})\) with the randomized algorithm for approximated permanent computation in Chin et al.15. By delving into the role of \({C}_{k}(\overrightarrow{n})\) and \({S}_{B}^{q}(\overrightarrow{n})\) further, we could formulate the quantum resource theory of Fock state that is a useful tool for understanding the quantum computing power of linear optical computing. Second, our approach to the complexity problem of BS can be compared to that of Rahimi-Keshari et al.38, which investigated the role of single-mode nonclassicality in computational complexity. This viewpoint is different from Chin et al.15 that focused on the multimode quantum correlation. We expect that there exists a unified theory that embraces the partial interpretations of former works, and the Fock state resource theory including concurrence and entropy is a strong candidate for such a theory.

Methods

Glynn’s formula and its generalization

Here we briefly introduce Glynn’s formula for N × N matrix permanent computation and its generalization for the permanents of matrices that have repeated rows and columns.

Glynn’s formula with the random variable expectation for N × N matrix A is given by

$${\rm{P}}{\rm{e}}{\rm{r}}\,(A)=\frac{1}{{2}^{N-1}}\sum _{\overrightarrow{x}\in {\{-1,1\}}^{N}}\,(\prod _{i=1}^{N}\,{x}_{i})\prod _{j=1}^{N}\,(\sum _{k=1}^{N}\,{A}_{jk}{x}_{k}).$$
(39)

The summation of \(\overrightarrow{x}\) is over \(\overrightarrow{x}\in {\{-1,\mathrm{1\}}}^{N}\), or \(\overrightarrow{x}\in {\mathscr{X}}\equiv {\mathcal R} \mathrm{[2]}\times \cdots \times {\mathcal R} \mathrm{[2]}\), where \( {\mathcal R} [i]\) is a set that consists of the ith root of unity.

When A has repeated rows or columns, and the ith column (or row) is repeated n i -times, Eq. (39) is generalized to9

$${\rm{P}}{\rm{e}}{\rm{r}}\,(A)=\sum _{\overrightarrow{z}\in {\mathscr{X}}}\,{v}_{\overrightarrow{n}}^{2}(\prod _{i=1}^{N}\,{\bar{z}}_{i}^{{n}_{i}})\prod _{j=1}^{N}\,(\sum _{k=1}^{N}\,{A}_{jk}{z}_{k}),$$
(40)

where \({\mathscr{X}}\equiv {\mathcal R} [{n}_{1}+\mathrm{1]}\times \cdots \times {\mathcal R} [{n}_{N}+\mathrm{1]}\) and \({v}_{\overrightarrow{n}}\equiv \sqrt{{\prod }_{i=1}^{N}\,({n}_{i}!/{n}_{i}^{{n}_{i}})}\). The runtime for this algorithm is \(({\prod }_{k=1}^{N}\,({n}_{k}+\mathrm{1)}{\alpha }_{\overrightarrow{n}}N)={\mathscr{O}}\,({\sum }_{k=1}^{{\alpha }_{\overrightarrow{n}}}\,{X}_{k}(\overrightarrow{n})\,{\alpha }_{\overrightarrow{n}}N)\).

When A has repeated rows and columns, and the ith column is repeated n i -times and the jth column is repeated m j -times, the above equations are expressed more generally10 as

$${\rm{P}}{\rm{e}}{\rm{r}}\,(A)=\sum _{\overrightarrow{z}\in {\mathscr{X}}}\,{v}_{\overrightarrow{n}}^{2}(\prod _{i=1}^{N}\,{\bar{z}}_{i}^{{n}_{i}})\prod _{j=1}^{N}\,{(\sum _{k=1}^{N}{A}_{jk}{z}_{k})}^{{m}_{j}}.$$
(41)

Here, \({\mathscr{X}}\) is the same as that in Eq. (40). From the summation form of Eq. (41) and the symmetry between rows and columns, it is straightforward to see that the minimal runtime \({{\mathscr{T}}}_{min}(\overrightarrow{n},\overrightarrow{m})\) for Eq. (41) is equal to Eq. (32).