Asymptotic evaluation of bosonic probability amplitudes in linear unitary networks in the case of large number of bosons

I Introduction

Linear optical networks play an important role in the quantum manifestations of light. Indeed, the well-known Hong-Ou-Mandel (HOM) dip HOM (see also Refs. CST ; LM ) is a direct manifestation of the quantum indistinguishability of photons. Recently LB a generalization of the HOM effect and difference in behavior of bosons and fermions was theoretically studied in the general setting of Bell multiport beam splitters (see also, for instance, Refs. Mult1 ; Mult2 ). Many important results on quantum interference phenomena in multiport devices for bosons and fermions were recently discovered. There is a zero transmission law ZTL in Bell multiports due to symmetry of the network matrix. Moreover, a generalization of the suppression laws and many-particle interferences beyond the boson and fermion statistics in unitary linear networks are found MPI . Recently the experimental advances have allowed to verify the HOM effect and the zero transmission law with three photons on a tritter HOM2zero . There are other important experimental advances in quantum interference experiments with indistinguishable particles in linear multiport devices, for instance, the recent three-photon quantum interference experiment on an integrated eight-mode optical device MPH . Moreover, the experimental efforts are now also directed at building the so-called boson sampler E1 ; E2 ; E3 ; E4 (see also below). Finally, a linear bosonic network is the central part of an ingenious proposal of quantum computation based on linear optics KLM .

It is known that the outcome probability amplitudes in unitary bosonic networks (e.g. in optical multiports) are expressed through matrix permanents, somewhat similar to the fermonic amplitudes which are given by matrix (a.k.a. Slater) determinants. Considering $N$ non-interacting bosons in an unitary network of $M$ input and output modes, one has to compute the matrix permanents of the $N\times N$-dimensional complex matrices composed of repeated rows and columns of the network matrix to find the bosonic transition amplitudes in the network. The permanent of a matrix Minc can be obtained by taking the well-known Laplace expansion formula for matrix determinant and setting to “$+$” the signatures of all permutations in the summation.

The relation of bosonic amplitudes to matrix permanents was explored a long time ago C in connection with the quantum fields theory. Recently it has attracted a lot of renewed attention. One reason is that the similarity between fermions and bosons does not go along if one tries to compute matrix permanent: unlike it is for matrix determinant, computation of the permanent of an arbitrary matrix is $\#P$-complete, which is the result of a classic paper in the computational complexity theory Valiant . This means that no algorithm polynomial in matrix size can compute matrix permanents. The fastest known algorithm for computing the permanent of an arbitrary (complex) $n\times n$-dimensional matrix is based on Ryser’s formula Ryser and requires $\mathcal{O}(n^{2}2^{n})$ flop operations. An interesting “physical” proof of the matrix permanent complexity result, based on the linear optical computing proposal KLM , was recently discovered A1 . Moreover, even the problem of approximating the permanent of an arbitrary complex matrix to a polynomial multiplicative error was also shown A1 to be a $\#P$-hard problem.There is an important exception of permanents of matrices with positive elements, which can be effectively approximated JSV , but such permanents are not connected to the quantum transition amplitudes in linear networks. Complexity of the matrix permanent was analyzed in the context of quantum mechanics in Ref. TT , where a method based on quantum measurement was proposed to directly measure the matrix permanent in a bosonic quantum multiport device. It was shown that the permanent of an arbitrary matrix can be expressed as a quantum observable. However, the catch of the method lies in an exponential number of necessary measurements, since the variance of the observable giving matrix permanent is exponential in matrix size TT (for comparison it was shown that matrix determinant can be found in just a single measurement).

Deep connection between the complexity of bosonic networks and that of matrix permanents was used in the recent proposal of a new model of quantum computer with noninteracting identical bosons, which, though not being an universal quantum computer, nevertheless can perform computations considered to be hard on a classical computer A1 ; AA . Such new quantum computer was compared to classical Galton’s board where, instead of classical balls and a single entry point, identical bosons are launched into different modes of a linear network. A very crucial difference, however, is that in the quantum network case the bosonic probabilities themselves are not known beforehand, since they cannot be effectively computed on a classical computer (for a sufficiently large network). One has to run the actual sampling experiments to find them. In the technical part, the proposal depends on the hardness to approximate the matrix permanent of a large $N\times N$-dimensional submatrix of an arbitrary unitary $M\times M$-dimensional matrix and some numerically tested conjectures. This proposal has generated experimental efforts to build the necessary bosonic network E1 ; E2 ; E3 ; E4 (see also the related experiments of Refs. HOM2zero ; MPH ).

It is believed that the hardness of computing the permanent of an arbitrary matrix is related to the matrix rank. Computation of the matrix permanent of an arbitrary (complex) $n\times n$-dimensional matrix of rank $R_{n}$ would apparently require at least on the order of $n^{R_{n}}$ flop operations on a classical computer (see, for instance, Ref. Gurvits ; moreover, this estimate for a matrix with repeated columns or rows is derived in Appendix D). In terms of bosonic networks, different limits are possible in this respect, which can be roughly divided by the relation between the number of bosons $N$ (i.e. matrix size) and the number of modes $M$ in the network (i.e. the maximum rank). For instance, it is shown AA that even an approximate computation of the matrix permanents describing linear bosonic networks is $\#P$-hard in the limit of $N\to\infty$ and $M\gg N^{2}$. The latter inequality is related to the “boson birthday bound”, i.e. the probability of two bosons to land into the same mode is negligible in this limit (and there is an experimental confirmation BosBirthExp ).

But even for Bell multiports there is no analytical approach capable to give probabilities of individual events when the latter are greater than zero (however, probabilities of averaged output events, such as the probability to find a specific number of bosons in one output port or a specific number of occupied output ports was in fact approximated by an analytical formula based on the classical consideration ZTL ). Such an analytical approach is especially lacking in the limit of a large matrix size, when the computational complexity makes numerical evaluation practically hard. One would not expect the existence of an asymptotic analytical approximation in the limit $N\to\infty$ and $M\geq N$, since it would contradict the $\#P$-hardness of numerical approximation of the permanent in this limit (indeed, one could then run a numerical scheme mimicking the analytical approach). An asymptotic analytical approximation could only exist for $M\ll N$. Such an approach is developed below by reducing the summation of the Laplace expansion of matrix permanent to a multidimensional integral of the saddle-point type in the limit $N\to\infty$ and $M$ finite. The very existence of an asymptotic approximation is important and stems from the fact that it is not dependent on the matrix size $N$ (i.e. the number of bosons in the network), but only on the solution of some $2M-1$ bilinear equations (the matrix scaling problem) giving the saddle points of the integral.

Practical usage of the proposed approach crucially depends on effective solution of two other problems of different type. One is the unitary matrix scaling problem, i.e. for a given unitary matrix $U$ to find all diagonal scaling matrices $X$ and $Y$ with the complex-valued elements such that the product matrix $XUY$ has given row and column sums (these matrices give the saddle points of the integral). The other problem is to derive a formula for multidimensional saddle point method with the multiple, i.e. coalescing, saddle points. Only the case of simple saddle points has a general solution so far FB and a special case of two coalescing saddle points on the real axis is treated FP . However, the coalescing saddle points represent exceptional cases rather than the general case. For simple saddle points, which is the general case, an explicit formula for bosonic probability amplitudes is derived below. As we use the well-known Stirling formula for the factorial in one of the steps, the approach is restricted to the case when there is no input or output mode without at least one boson in it (though this latter limitation can be lifted, for simplicity sake it is not pursued in this work). The asymptotic formula for bosonic transition amplitudes is tested on two-mode network (i.e. the beam splitter) where it shows a very good accuracy within the above described limitations. Three-mode network (i.e. the tritter) is also considered.

The rest of the text is organized as follows. In section II an integral approximation of bosonic probability amplitudes in unitary networks is derived. The integral is evaluated by the saddle point method and an explicit formula is derived when the saddle points are simple. Some details of the calculations are relegated to B and C. For completeness, a short derivation of bosonic probability amplitudes as matrix permanents is given in A. Comparison with the classical particles on a network is also discussed. In section III the general formula is tested on two-mode network (the beam splitter) where the sources of error are identified and discussed. Moreover, three-mode network (the tritter) is also considered. Section IV contains summary of the results. The computational complexity of the permanent of a matrix with repeated rows and/or columns, i.e. for which the asymptotic approximation is developed, is considered in detail in Appendix D.

II The saddle point method for probability amplitudes in linear bosonic networks

We consider the transition probability amplitudes of $N$ bosons between the input ($|f_{1}\rangle,\ldots,|f_{M}\rangle$) and output ($|g_{1}\rangle,\ldots,|g_{M}\rangle$) modes of a $M$-mode unitary network, as in Fig. 1, which are given by (see, for instance, Refs. TT ; Scheel )

Here matrix $U[n_{1},\ldots,n_{M}|m_{1},\ldots,m_{M}]$ consists of repeated rows and columns of the network matrix $U_{kl}=\langle g_{l}|f_{k}\rangle$ (the order being insignificant), with $n_{k}$ duplicates of the $k$th row and $m_{l}$ of the $l$th column, satisfying $\sum_{k=1}^{M}n_{k}=\sum_{k=1}^{M}m_{k}=N$ (for more details, see A). Recall that the permanent of a $N\times N$-dimensional matrix $A$ is given by summation over all possible permutations $\tau$ in the product of $N$ different matrix elements Minc , i.e.

To derive an asymptotic approximation for the bosonic amplitude given in Eq. (1) we regroup of the summation in Eq. (2) in such a way that the nature of the matrix $U[n_{1},\ldots,n_{M}|m_{1},\ldots,m_{M}]$ as a matrix with repeated rows and columns is used. This regrouping leads also to an interesting interpretation of the permanent of $U[n_{1},\ldots,n_{M}|m_{1},\ldots,m_{M}]$ as an average over a lattice of contingency tables, see Fig. 2, with probabilities given by a bi-multivariate generalization of the hypergeometric distribution, known as the Fisher-Yates distribution in mathematical statistics Fisher ; Yates (see also a number of reviews Refs. Good ; DE ; DG and the references therein). This representation is discussed below in detail.

By successive application of the Laplace expansion for permanents Minc one can express the boson probability amplitude (1) as an average over the lattice of $M\times M$-dimensional matrices $S_{kl}$ ($S_{kl}\in\{0,1,2,3\ldots\}$) with given row and column sums, equal here to the distribution of bosons in the input and output modes: $\sum_{l=1}^{M}S_{kl}=n_{k}$ and $\sum_{k=1}^{M}S_{kl}=m_{l}$. Indeed, the first application of the Laplace expansion consists of dividing the matrix $U[n_{1},\ldots,n_{M}|m_{1},\ldots,m_{M}]$ into two parts, and the permanent into a sum over the products of partial permanents, one involving $n_{1}$ first rows (all equal) and the other involving the rest of matrix rows. The summation runs over all possible partitions of the column indices between the two permanents, such that one permanent contains $n_{1}$ column indices ($w_{1},\ldots,w_{n_{1}}$, with $S_{11}$ of them being equal to $1$, $S_{12}$ being equal to $2$, etc), while the other permanent contains $N-n_{1}$ column indices ($w_{n_{1}+1},\ldots,w_{N}$). Thus, application of the Laplace expansion as above described gives (with $(w_{1},\ldots,w_{N})$ being a permutation of $(\underbrace{1,\ldots,1}_{m_{1}},\underbrace{2,\dots,2}_{m_{2}},\ldots,\underbrace{M,\ldots,M}_{m_{M}})$)

where we have used that permanent of the first submatrix is equal to $n_{1}!\prod_{l=1}^{M}U_{1l}^{S_{1l}}$ and, due to the permutational invariance of matrix permanent, the r.h.s. contains functions of numbers of repeated columns and not the column indices themselves. Repeating the above procedure, by taking out successively each set of repeated rows, we obtainEq. (4) could also be deduced from Eq. (1) and the formula for bosonic probability amplitude derived in Ref. Urias by application of Wick’s theorem; note that using the matrix permanent and the Laplace expansion is an impressive shortcut to that involved derivation.

Eq. (4) has an interesting statistical interpretation.Besides physical interpretation as the sum over quantum probability amplitudes of all possible transitions through the network, i.e. a variant of R. Feynman’s path integral formula. Indeed, the ratio of factorials which appears on the r.h.s. in Eq. (4), divided by $N!$, is know as the Fisher-Yates distribution Fisher ; Yates (see also the reviews Refs. Good ; DE ; DG ). It appears in applied mathematical statistics, namely in Fisher’s exact test of independence of two properties, and uses the so-called contingency tables (here $S_{kl}$). The Fisher-Yates distribution gives the conditional probability of getting a matrix $S_{kl}$ (the contingency table) of the joint frequencies of two statistically independent properties, given the row and column sums are equal to their marginal frequencies (the margins, in our case $n_{1},\ldots,n_{M}$ and $m_{1},\ldots,m_{M}$, respectively). Since the two properties are independent, a simple exercise in combinatorics (see also footnote “d” ) leads to the following probability formula of the Fisher-Yates distribution (using a shortcut notation $\{v_{i}\}$ for a set of indexed variables: $v_{1},v_{2},\ldots$)

Note that the delta symbols in Eq. (4) restrict the summation to matrices $S_{kl}$ with given margins and precisely under these constraints the probabilities of Eq. (5) sum to 1. The matrix permanent of Eq. (4) is thus multiplied by $N!$ the value of the characteristic function $\chi(\{\lambda_{kl}\}|\{n_{k},m_{l}\})\equiv\langle\exp\{\sum_{k,l=1}^{M}\lambda_{kl}S_{kl}\}\rangle$ of the Fisher-Yates distribution at the parameters $\lambda_{kl}$ equal to logarithms of the elements of $U$:

Eq. (6) has the following physical interpretation. The bosonic transition amplitude in a quantum network, between the input $|n_{1},\ldots,n_{M}\rangle_{f}$ and output $|m_{1},\ldots,m_{M}\rangle_{g}$ Fock states, is an average of products of amplitudes of “elementary processes”, as depicted in Fig. 2, corresponding to contingency table $S_{kl}$ (i.e. a matrix satisfying the constraints $\sum_{l=1}^{M}S_{kl}=n_{k}$ and $\sum_{k=1}^{M}S_{kl}=m_{l}$), assuming the mutual statistical independence of distribution of bosons in the input and output modes.There are $C_{in}=\frac{N!}{\prod_{k=1}^{M}n_{k}!}$ ways to distribute bosons over the input modes, $C_{out}=\frac{N!}{\prod_{l=1}^{M}m_{l}!}$ ways to distribute bosons over the output modes, and, as the two distributions are independent, in total $C_{in}C_{out}$ ways to distribute bosons over the input and output modes. For a table $S_{kl}$ we select $C_{S}=\frac{N!}{\prod_{k,l=1}^{M}S_{kl}!}$ combinations from these distributions, thus $P(\{S_{kl}\}|\{n_{k},m_{l}\})=\frac{C_{S}}{C_{in}C_{out}}$, i.e. Eq. (5).

Eq. (6), however, does not seem to be of any help for numerical evaluation of matrix permanents, since the number of contingency tables scales exponentially with $N$ (precisely, their number scales exponentially with the margins for a fixed table size and, in fact, the problem of counting the contingency tables is $\#P$-hard, see, for instance Refs. DG ; Bend ; Barv ; GM ; Barv2 ).

On the other hand, computation of the permanent of a matrix with repeated rows and/or columns can be effectively carried out by using the available in this case reductions in Ryser’s algorithm (see, for instance, Ref. Thesis ). Indeed, it is shown in Appendix D that a modified Ryser algorithm requiring just $\mathcal{O}(N^{M+1})$ flops is available in this case (this algorithm was used for obtaining Fig. 7 of section III.2 below).

II.2 Approximating the bosonic probability amplitude by a multidimensional integral

On the other hand, the fact that the number of contingency tables in Eq. (6) scales exponentially with $N$ is an indication on possibility of an asymptotic approach. Indeed, since the lattice of contingency tables $S_{kl}$ is exponential in $N$, if we divide $S_{kl}$ by $N$, the resulting matrix $p$, $p_{kl}\equiv S_{kl}/N$, will belong as $N\to\infty$ to a dense lattice in the continuous convex set of matrices with real-valued entries and given row and column sums. Therefore, a sum over such a lattice can be replaced by a multidimensional integral with $p_{kl}$ as the integration variables. Moreover, a large parameter $N$ would appear in the exponent of the integrand, thus allowing for an asymptotic evaluation of the permanent. This is the approach pursued in the following.

First of all, we need to approximate the Fisher-Yates distribution (5) by a manageable smooth function of the integration variables $p_{kl}$. Using an approximate formula for the multinomial coefficient, given by Eq. (59) of B, we obtain for $n_{k},m_{k},S_{kl}\geq 1$:

where we have denoted by $\mathcal{I}$ the mutual information function, namely

with the Shannon entropy function denoted by $\mathcal{H}$.

The sum in Eq. (4) can be replaced by an integral as $N\to\infty$, since the difference between the elements of two neighboring $p$-matrices (the $N^{-1}$-scaled contingency tables $S_{kl}$) is of order $1/N$. The Kronecker delta symbols must be replaced by the $N^{-1}$-scaled Dirac delta functions:

There are only $2M-1$ independent constrains in the product of $2M$ Kronecker deltas in Eq. (4), since the contingency table elements $S_{kl}$ sum to $N$, giving both the sum of the row sums and of the column sums. Hence, the integration domain is $(M-1)^{2}$-dimensional. Using these observations and some elementary algebra we obtain from Eqs. (4), (7), and (9):

Here we have introduced an integration measure $d\mu(\{p_{kl}\})$ over an $(M-1)^{2}$-dimensional subspace in the convex set of all matrices with positive elements constrained only by the $2M-1$ Dirac delta functions from Eq. (9).Note that an arbitrary subset of $2M-1$ delta functions can be used, see also C. It reads

The error of the approximation in Eq. (10) is estimated to have a multiplicative order $\sim 1/N$, since this is the order of our approximation of the Fisher-Yates distribution by a smooth function in Eq. (7), whereas replacing a finite sum by an integral brings also an error on the order of difference between the values of two nearest lattice points, i.e. $\Delta p_{kl}\sim 1/N$. The multidimensional integral in Eq. (10) is in the standard form used for asymptotic expansion in powers of $1/N$ by the saddle point method (called also the steepest descent method). However, since the integral representation in Eq. (10) has already an error of order $1/N$, only the leading term of the resulting asymptotic expansion is meaningful.

II.3 The matrix scaling problem giving the saddle points

At this stage, let us recall the general formula for the leading term of an $n$-dimensional integral, given by the saddle point approximation, when the saddle points are simple FB :

where the summation is over contributing saddle points $z_{j}$. The determinant in the denominator of Eq. (12) is of the Hessian matrix, i.e. the matrix composed of the second-order derivatives of $\phi(z)$, taken at the respective saddle point. The saddles $z_{j}$ are found by a deformation, as allowed by analyticity of $\phi(z)$, of the integration domain in the extended complex-valued space of $z$, such that the deformed domain is contained in the steepest descent regions of the integrand. The next term in the asymptotic expansion, as compared to the leading term of Eq. (12), has the relative order of $1/N$. Let us now apply the result (12) to our specific case and derive the saddle point approximation of the permanent.

The saddle points (matrices $p$, in our case) are found as extremals of the multivariate function in the exponent of the integrand. In our case the function reads

with, however, only $(M-1)^{2}$ independent variables out of the total $M^{2}$ matrix elements $p_{kl}$. Using the Lagrange multipliers $\lambda_{k}$ and $\mu_{l}$ one can equivalently look for extremals of the augmented function

Equating the differential of $\mathcal{F}$ to zero we obtain that the saddle points have the following general form

II.3.2 Calculation of the Hessian. The main result

In our case, the Hessian matrix is with respect to some independent $(M-1)^{2}$ variables from the $M^{2}$ elements of matrix $p$. Therefore, an extension of the saddle point method to the constrained integration is needed, which runs as follows. We rewrite the constrains on variables $p_{kl}$ in Eq. (11) as a set of linear equations by introducing a matrix $C_{j,kl}$, where the enumeration order of the double index $(k,l)$ is as follows $(k,l)=\{(1,1),\ldots,(1,M),(2,1),\ldots(2,M),\ldots,(M,1),\ldots,(M,M)\}$, i.e. index $k$ runs slower than index $l$. The constraints can be rewritten as followsThe specific subset of $2M-1$ constraints is in accord with the selected measure in Eq. (11).

Matrix $C$ in Eq. (19) has rank equal to $2M-1$. It can be partitioned into a $(2M-1)\times(2M-1)$-dimensional nonsingular submatrix $C^{(I)}$ and a submatrix $C^{(II)}$. These two matrices induce a similar partition of the elements $p_{kl}$, treated as a vector with double index $(k,l)$. Using a vector notation $\overline{p}$, we can cast the system of constraints given by Eq. (19) as follows

Eq. (20) allows to extract independent variables from the $M^{2}$ elements of $p$ and calculate the needed Hessian. First, by introducing two vectors, $\overline{\xi}$ consisting of $2M-1$ dependent integration variables and $\overline{\eta}$ of $(M-1)^{2}$ independent ones, as follows $\overline{\xi}=C^{(I)}\overline{p}^{(I)}$ and $\overline{\eta}=\overline{p}^{(II)}$, we satisfy the constraints by fixing the value of $\overline{\xi}$ according to Eq. (20) (i.e. by integrating over $\overline{\xi}$ using the Dirac delta functions in Eq. (11)) and obtain the rest of the measure $d\mu$ as follows

Second, due to linearity of constraints (20), the determinant of the matrix of second derivatives of $\phi(\{p_{kl}\})$ with respect to $(M-1)^{2}$ independent variables can be evaluated from the full matrix of second derivatives with respect to all variables $p_{kl}$ by using Eq. (20). We get the following result

Here $[\widetilde{B},I]$ stands for the block matrix constructed from the transposed $(2M-1)\times(M-1)^{2}$-dimensional matrix $B$ and the $(M-1)^{2}\times(M-1)^{2}$-dimensional matrix unit $I$. Furthermore, the determinant on the r.h.s. of Eq. (25) can be further simplified by using an identity which generalizes Sylvester’s identity for determinant of a block matrix (see C for details) valid for a nonsingular matrix $A$:

where matrix $B$ is as in Eq. (25). In our case the Hessian matrix $A$ in the second differential of $\phi=\mathcal{I}(\{p_{kl}\})-\sum_{k,l=1}^{M}p_{kl}\ln U_{kl}$ is diagonal, i.e.

thus Eq. (26) applies for $p_{kl}\neq 0$. Note that the exponent with a large parameter $N$ in the integral on the r.h.s. of Eq. (10) evaluated at a saddle point $p$ such that some $p_{kl}=0$ would be infinite (thus this case is ruled out). Indeed, using Eq. (16), we get at a saddle point $p_{kl}=x_{k}U_{kl}y_{l}$:

Now, using Eqs. (12), and (21)-(28) into Eq. (10) and noticing that

we obtain a formula for the leading term approximation to the matrix permanent in the case of simple saddle points (our main result)

Here the sum over all contributing saddle points $p^{(s)}_{kl}=x^{(s)}_{k}U_{kl}y^{(s)}_{l}$ is implied and matrix $D^{\prime}$ in the denominator is as follows

The above discussion implies that the symmetries of bosonic probability amplitudes are preservedThis is also manifested by exact cancellation of probability amplitudes in the generalized HOM effect, Figs. 3 and 6 of section III. by the saddle point approximation. For instance, the inversion symmetry: ${}_{g}\langle m_{1},\ldots,m_{M}|n_{1},\ldots,n_{M}\rangle_{f}={}_{f}\langle n_{1},\ldots,n_{M}|m_{1},\ldots,m_{M}\rangle_{g}^{*}$.

II.4 Comparison with classical identical particles on a network

Let us compare the transition probabilities of indistinguishable bosons with the transition probabilities of classical particles (which we consider identical). In the classical case, the elementary process of Fig. 2 means redistribution of $n_{k}$ identical classical particles from the $k$th input mode into the $M$ output modes with the output distribution given by the same contingency table $S_{kl}$. The difference is that the probabilities are multiplied and summed up, thus instead of the amplitude of an elementary quantum process, given by the product $\prod_{l=1}^{M}U_{kl}^{S_{kl}}$ in Fig. 2, we have the probability of an elementary classical process, given by $\prod_{l=1}^{M}|U_{kl}|^{2S_{kl}}$. As the particles are identical (i.e. the paths of the individual particles through the network are not traced), the total probability of such an elementary process is given by the latter product multiplied by the number of redistributions of $n_{k}$ identical particles from the $k$th input mode into $M$ output modes, i.e. by the factor $\frac{n_{k}!}{\prod_{l=1}^{M}S_{kl}!}$. Summing up over all such probabilities (i.e. over the elementary processes from different input modes) and identifying the Fisher-Yates distribution in the summation, we obtain the transition probability of $N$ classical particles through $M$-mode network (see also Ref. MPI )

where the input and output distributions are $\{n_{1},\ldots,n_{M}\}$ and $\{m_{1},\ldots,m_{M}\}$, respectively (cf. with the quantum probability amplitude given by Eqs. (1) and (6)). For instance, for Bell multiports $|U_{kl}|^{2}=\frac{1}{{M}}$ and we obtain the resulting probability as follows (see also Ref. ZTL )

One can apply the saddle point approximation also to the classical probability given by Eq. (33). Using this simple observation, we can compare complexity of the saddle point approximation in the classical and quantum cases. A very important difference is spotted immediately: the classical analog of the matrix scaling problem is formulated for a matrix of positive elements $A_{kl}\equiv|U_{kl}|^{2}$ (note that matrix $A$ is doubly stochastic, i.e. its row and column sums are equal: $\sum_{k=1}^{M}A_{kl}=\sum_{l=1}^{M}A_{kl}=1$). It is known that the matrix scaling problem for a positive matrix has a unique positive solution, since it is equivalent to a minimization problem of a convex function MatScal ; MatScal2 ; MScalUnique . Note that the corresponding saddle point belongs to the integration domain over the contingency tables $p_{kl}=S_{kl}/N$, i.e. $x_{k}A_{kl}y_{l}<1$, since both $A_{kl}\geq 0$ and $x_{k},y_{k}>0$ (the positive solution is constrained by the margins). Therefore, it is the only contributing saddle point in the classical case. This is remarkably different from the quantum case, where, as is discussed below, generally there are more than one complex-valued contributing saddle points. They loose interpretation of the dominating “real processes” of Fig. 2 (since the corresponding contingency table $S_{kl}$ is complex). However, the saddle points describe in a simpler way the quantum interferences between exponentially many of such real processes in Eq. (6).

By using the explicit form of the saddle point, substituting Eq. (35) into Eq. (29) and the resulting expression into Eq. (33) one recovers the exact result given by Eq. (34) for a classical analog of Bell multiports.

III Testing accuracy of the saddle point approximation

To apply the saddle point approximation (29) we first have to solve the matrix scaling problem (16) which is a bilinear system of equations in $x$ and $y$. One can reduce the number of variables by half by resolving one of the equations in Eq. (16), for instance, $y_{l}=\sum_{k=1}^{M}U^{*}_{kl}({n_{k}}/{Nx_{k}})$. By introducing a vector $R$ containing all $M-1$ independent variablesMultiplication of all $x$-variables by a complex number $\lambda$, $x_{k}\to\lambda x_{k}$, does not change the saddle points, since it induces the inverse scaling of the $y$-variables: $y_{k}\to y_{k}/\lambda$. and a set of $M-1$ vectors $Z^{(l)}$, defined as follows:

we obtain from Eq. (29) a reduced system in the following form

Eq. (37) is hard to solve analytically for more than two modes (despite considerable efforts, the solution has not been found even for the simple case of Bell multiports). Moreover, all solutions of Eq. (37) are needed, since all saddle points contribute to the approximation in general. On the other hand, it is relatively easy to find solutions to Eq. (37) numerically (for instance, by the all-purpose nonlinear equations solver available in MATLAB with a random initial guess to find all possible solutions). It was found that the total number of solutions is dependent on $\{n_{1}/N,...,n_{M}/N\}$ and $\{m_{1}/M,...,m_{M}/N\}$ and that there can be symmetries leading to degeneracies. The absence of a formula for the total number of solutions prevents analysis of the computational complexity of Eq. (37).

Below we consider accuracy of the saddle point approximation in two cases: the beam splitter and the tritter, where the beam splitter allows an analytical solution, while already the tritter case requires a numerical solution.

The case of beam-splitter, $M=2$, is analytically solvable. As is known Berns , two and three dimensional networks are uniquely defined by moduli $|U_{kl}|$ of the unitary matrix elements, whereas the $2M-1$ phases are scaled out by changing the unimportant phases of the input and output states. Let us consider the symmetric beam splitter, which is given by the following matrix

Then Eq. (37) leads to a quadratic equation for $R_{1}=\frac{\sqrt{n_{2}}}{\sqrt{n_{1}}}\frac{x_{1}}{x_{2}}$:

There are two saddle points $p$ whose $x$ component in Eq. (15) readsTo obtain $x_{1,2}$ from $R_{1}$ we have taken into account the scale invariance $x\to\lambda x$ and $y\to y/\lambda$.

The $y$ component is given by $y_{1}=(-\frac{n_{1}}{Nx_{1}}+\frac{n_{2}}{Nx_{2}})/\sqrt{2}$ and $y_{2}=(\frac{n_{1}}{Nx_{1}}+\frac{n_{2}}{Nx_{2}})/\sqrt{2}$. In Eq. (40) we have introduced a phase $\phi$, however, it is a real value only under the condition that $\gamma^{2}<1$, i.e. when

with $\Delta n=n_{2}-n_{1}$ and $\Delta m=m_{2}-m_{1}$. Under condition (41) the $y$ component can be given as follows

Since $4m_{1}m_{2}(1-\sigma^{2})=4n_{1}n_{2}(1-\gamma^{2})=N^{2}-(\Delta n)^{2}-(\Delta m)^{2}$, the same threshold condition (41) applies to phases $\psi$ and $\delta$. When condition (41) is violated, the phases $\phi$, $\psi$, and $\delta$ become complex-valued ($\phi$ and $\psi$ become imaginary). This corresponds to a phase transition in the bosonic probability amplitudes (see below).

The saddle points are given by Eq. (15). We get:

valid in the whole domain $|\Delta n|\leq N$, $|\Delta m|\leq N$. When condition (41) is satisfied, i.e. when $|\gamma|\leq 1$, both saddle points are complex-valued and contribute to the integral in Eq. (29). When it is violated, the saddle points become real valued. However, one of them does not contribute to the integral, since it escapes from the integration domain $0\leq p_{kl}\leq 1$. Which one of the saddle points contributes depends on the sign of $\gamma$, i.e. the sign of $\Delta m$, and the ratio $n_{1}/n_{2}$ (see also Fig. 3(c) below).

Finally, after a simple algebra, the determinant given by Eq. (32) becomes

where we have used that $16n_{1}n_{2}m_{1}m_{2}=\left(N^{2}-(\Delta n)^{2}\right)\left(N^{2}-(\Delta m)^{2}\right)$. Substituting Eqs. (40), (42), and (47) into Eq. (29) and the resulting approximation into Eq. (1) we obtain the saddle-point approximation to bosonic probability amplitudes of the symmetric beam-splitter (38).

III.1.2 Comparison with the exact result

To compare with the exact result the following expression for bosonic probability amplitude for network matrix of Eq. (38) will be used (see also Ref. LM )

The correspondence of the exact result (48) with the saddle point approximation can be divided into three regions. In the first region condition (41) is satisfied. This region contains the generalized HOM effect CST ; LM . In this case the amplitudes of two contributing terms in Eq. (29), corresponding to two saddle points, have the same moduli and the only difference lies in their relative phase. The corresponding domain in the two-dimensional plane with coordinates $\Delta n$ and $\Delta m$ is the inside of the circle $(\Delta n)^{2}+(\Delta m)^{2}\lesssim N^{2}$. There is cancellation of the probability amplitudes for $n_{1}=n_{2}=N/2$ and odd values of $m_{1}$, Fig. 3(a) (see also Refs. CST ; LM ). In the saddle point approach, the cancellation is due to symmetry of the two saddle-point contributions with the only difference being their relative phase given by $(-1)^{m_{1}}$ (note: the cancellation is captured exactly by the saddle-point approximation).

On the other hand, there is another regime: the exponential decay of the probability amplitude as $m_{1}$ approaches either or $N$, see Fig. 3(c). This regime has not been studied previously (for instance, the approximation of Ref. LM only captures the oscillating regime). It appears when condition (41) is violated. In this case, there is just one contributing saddle point, the one which has the smallest moduli contribution to the permanent (this is similar to what occurs in the saddle-point approximation to the Airy function).

The third region is about the circle $(\Delta n)^{2}+(\Delta m)^{2}\approx N^{2}$. This region contains two coalescing saddle points and cannot be approximated by Eq. (29) valid for the simple saddle points only (their contributions diverge on this circle, which is due to the determinant (47) approaching zero). Outside this region, which is restricted to narrow neighborhoods of the points $m_{1}=10$ and $m_{1}=50$ in Fig. 3(c), the saddle point approximation has a very good accuracy as is shown in Fig. 3(b), where the accuracy of the approximation (29) is compared to that for the binomial coefficient, given by Eq. (59) of B and used to build the approximation of the Fisher-Yates distribution (7). It is seen that the relative error is approximately twice as that in the approximation of the binomial coefficient.

The transition from the two contributing saddle points, with an oscillating probability amplitude as function of $m_{1}$, to a single contributing saddle point, with an exponentially decaying probability amplitude, is similar to the Airy function behavior, thus the related integral in Eq. (29) can be, in principle, expressed through a linear combination of the Airy function and its first derivative. Such results are available for the real-valued coalescing saddle points (for instance, in Refs. FB ; FP ). However, the method needs to be generalized to the complex-valued case before it could be applied to the integrals approximating the bosonic probability amplitudes.

Finally, it is interesting to observe that the saddle point approximation can give correct results down to the very small number of bosons, for instance, it correctly predicts the HOM effect HOM (though, obviously, small $N$ violate the assumption $N\gg 1$). Several results for small number of bosons are collected in Fig. 4, where we have $2\leq N\leq 9$.

III.1.3 Scaling of the relative error of the saddle point approximation

Let us verify that the relative error scales as $1/N$ for fixed $n_{1}$ and $m_{1}$. The scaling of the relative error of the saddle point approximation is given in Figs. 5 and 6. Comparing Fig. 5 with Fig. 6 one can notice the oscillations of the relative error around the law of the inverse proportionality in the latter case. The origin of these oscillations is unclear. For instance, they are not due to approaching the boundary circle $(\Delta n/N)^{2}+(\Delta m/N)^{2}=1$, since all data points from the same figure represent one and the same point in the $(\Delta n/N,\Delta m/N)$-square. The only explanation is a very complicated general dependence of bosonic probability amplitude on $N$ for a fixed set of distributions $\{n_{1}/N,...,n_{M}/N\}$ and $\{m_{1}/M,...,m_{M}/N\}$ due to the fact that the phases of the individual saddle point contributions are multiplied by $N$.

III.2 The three-mode (tritter) case

Let us now consider three mode network (the tritter). The canonical symmetric tritter, used, for instance, in the recent experiment HOM2zero , has the following network matrix

Nonlinear system in Eq. (37) with $U$ from Eq. (49) seems to be unsolvable analytically, however, numerical solution contains at most six different saddle points. Moreover, numerical simulations with random three-mode unitary matrices $U$ has shown that six is the maximal number of saddle points for any three-mode network, where in most cases all saddle points contribute to the approximation and the respective vector parameters $x$ and $y$ have the following “most probable” form (after fixing one of the $x$-vector elements due to the scale invariance $x\to\lambda x$ of Eq. (37))

where $\phi_{k}$ and $\psi_{k}$ are real values (phases). This “most probable” form is very similar to the general solution in the beam splitter case of section III.1. However, the corresponding contributions to the approximation from such saddle points are not always of the same moduli as distinct from two-mode network. The approximation is not even qualitatively correct for small number of bosons. Moreover, it is found that the relative error always has oscillations reminiscent of those in Fig. 6 (what can explain the poor performance of the approximation for small $N$ in this case). Behavior of the relative error is illustrated in Fig. 7. Computation of the matrix permanent is carried out by a modified Ryser’s algorithm (similar as in Ref. Thesis ). Such algorithm has only polynomial in $N$ complexity as is shown in Appendix D.

IV Conclusion

We have shown that an asymptotic evaluation of matrix permanents giving bosonic probability amplitudes in unitary linear networks is possible for large number of bosons $N$ and fixed network size $M$, such that $N\gg M$. The asymptotic approximation reduces the problem of evaluation of permanents of $N\times N$-dimensional matrices with repeated rows and columns to a solution of a matrix scaling problem for $M\times M$-dimensional network matrix, which is a system of bilinear equations in $2M-1$ variables. For simple saddle points, an explicit formula for bosonic probability amplitudes is derived.

The asymptotic approximation has been compared with the exact analytical result available for two-mode network, i.e. the beam-splitter, and has been found to have good accuracy correlated with accuracy of the approximation of multinomial coefficient (used as a building block of the saddle-point approximation). The approximation error is studied also for three-mode network, i.e. the tritter, where the saddle points were found numerically. Interestingly, in the beam splitter case, the approximation correctly reproduces behavior of probability amplitudes even for small number of bosons, for instance, it reproduces the original HOM effect. The relative error of the approximation is found to scale inversely with the number of bosons, however, the scaling is plagued by oscillations about the inverse scaling law (the origin of which is unclear). These oscillations degrade the approximation for small number of bosons in the case of the tritter, where, in contrast to the beam splitter case, the approximation performs poorly for small number of bosons.

There are various regimes of behavior of bosonic probability amplitudes in unitary networks, which are dependent on the number of contributing saddle points. For instance, in the beam-splitter case there are two regimes: (i) the oscillating regime, when two saddle point contribute to bosonic probability amplitude and the generalized Hong-Ou-Mandel effects take place, and (ii) the regime of exponential decay of bosonic probability amplitudes, when only one saddle point contributes.

Practical application of the method is conditioned on solution of the matrix scaling problem giving the saddle points (where the whole set of solutions is generally required). Another problem of different type has to be solved before practical application of the method is attempted. One has to derive a general formula for the saddle point method applicable when the saddle points coalesce. There is also an important topological problem of identifying the contributing saddle points, which is a highly nontrivial one in general FB , but may have general solution for the type of integrals appearing in the saddle point approximation of matrix permanents.

It is unlikely that there is analytical solution to the matrix scaling problem for $M$-mode network with $M>2$, since the total number of different solutions grows rapidly with $M$. For $M=2$ there are at most two saddle points. Numerical simulations with random network matrices indicate that for $M=3$ there are at most six different saddle points, while for $M=4$ there are at most twenty different saddle points. Therefore, the matrix scaling problem may have an exponential in $M$ computational complexity (note that this complexity is an attribute of a quantum network: for classical particles on a network there is just one contributing saddle point).

Acknowledgements

The author would like to thank the referee for many valuable comments that resulted in substantial improvement of the presentation. This work was supported by the CNPq and FAPESP of Brazil.

Appendix A Bosonic amplitudes expressed as matrix permanents

Consider a $M$-dimensional quantum unitary network where $N$ bosons are launched into the input modes $|f_{1}\rangle,\ldots,|f_{M}\rangle$ and are detected in the output modes $|g_{1}\rangle,\ldots,|g_{M}\rangle$. The network is given by a $M\times M$-dimensional unitary matrix $U$, $U^{\dagger}U=I$, which transforms the input modes into the output modes,

i.e. between two orthogonal bases of a $M$-dimensional single-particle Hilbert space $H$. The goal is to express the bosonic transition amplitude between two Fock states $|n_{1},\ldots,n_{M}\rangle_{f}$ and $|m_{1},\ldots,m_{M}\rangle_{g}$, giving the input and output states of $N$ bosons:

where $\sum_{k=1}^{M}n_{k}=\sum_{k=1}^{M}m_{k}=N$ and it is implied that the two sets $\{i_{1},\ldots,i_{N}\}$ and $\{j_{1},\ldots,j_{N}\}$ are composed of repeated mode indices, e.g.

On the r.h.s.’s of Eqs. (52) and (53) there are unnormalized symmetric states of $N$ particles in the tensor product of $N$ single-particle Hilbert spaces $H\otimes H\otimes\ldots\otimes H$. Such a state is given by a sum over all permutations $\tau$ of $N$ indices of the single-particle states, divided by the number of all permutations, e.g.

The bosonic transition amplitude between the Fock states of Eqs. (52) and (53) is given by a double sum over the two sets of permutations of indices in the inner product of the $f$ and $g$ states, i.e. the indices of elements of the network matrix $U$. This double sum is converted to a single one over all permutations of the column indices by transferring one of the two permutations to the co-product indices, i.e.

where $\sigma^{\prime}\equiv\sigma\cdot\tau^{-1}$ runs over all permutations of $N$ elements. One immediately recognizes a matrix permanent on the r.h.s. of Eq. (56), where the matrix consist of repeated rows and columns of the network matrix $U$ (where the order is insignificant due to permutational invariance of matrix permanent). Therefore, we introduce the notation $U[n_{1},\ldots,n_{M}|m_{1},\ldots,m_{M}]$ for such a matrix and obtain the resulting bosonic amplitude proportional to permanent of this matrix, as in Eq. (1) of section II (see also Refs. C ; TT ; Scheel ).

Appendix B Approximating the multinomial coefficient

An approximation of the multinomial coefficient, and hence of the Fisher-Yates distribution, can be based on an exact formula for the factorial for $n\geq 1$Mortici :

where $\theta_{n}$ is tightly bounded: $1/6<\theta_{n}<0.177$. Interestingly, Eq. (57) can be extended to all $n\geq 0$ by carefully defining the limit $0^{0}=1$ and redefining the lower bound to $\theta_{0}=\frac{1}{2\pi}<\frac{1}{6}$. After some algebraic manipulations, the multinomial coefficient becomes

where $\mathcal{H}$ is the expected Shannon entropy function $\mathcal{H}(\{\frac{n_{k}}{N}\})\equiv-\sum_{k=1}^{M}\frac{n_{k}}{N}\ln\left(\frac{n_{k}}{N}\right)$. Since $\theta$ in Eq. (58) is always divided by $N\gg 1$, the simplest approximation is to drop it, thereby making an error of order $\mathcal{O}(N^{-1})$ and restricting ourselves to $n_{k}\geq 1$. Assuming the latter, we obtain

Finally, an even better approximation of the multinomial coefficient (for all $n\geq 0$) could be obtained if an uniform nonzero $\theta_{n}$ is selected, for instance $\theta_{n}=1/6$, i.e. similar as in Gosper’s approximation of the factorial Gosper . Such a formula, though being more complicated, would then improve the approximation of the Fisher-Yates distribution of section II. However, for the sake of simplicity, we do not pursue this approach in the present work.

Appendix C Proofs of the determinant identities of section II

Let us first proof a generalization of Sylvester’s determinant identity, i.e. Eq. (26) of section II. To this goal, one can use the following auxiliary Gaussian integral

where the real part of the quadratic form in the exponent is positive definite. The Gaussian integral in Eq. (62) can be easily evaluated and we obtain

On the other hand, one can also use the Fourier representation of the Dirac delta functions in the integrand of Eq. (60) and, by interchanging the integration order, evaluate the integral $J$ on the whole $x$ space and then take the inverse Fourier transform. As all integrals are Gaussian they are easily evaluated. We obtain

Comparison of Eqs. (63) and (64) gives the determinant identity of Eq. (26) for nonsingular matrices $A$ with a positive definite real part. The validity can be extended to arbitrary nonsingular matrices $A$ by uniqueness of the analytic continuation in the complex plane, by noticing that the r.h.s.’s of Eqs. (63) and (64) are analytic functions of the elements of $A$.

Finally, let us simplify the expression for a $(2M-1)\times(2M-1)$-dimensional principal minor of $D$ (31). To this goal the following determinant identity valid for $2\times 2$-block matrices can be used

which is a generalization of the formula for determinant of $2\times 2$-dimensional matrices. For the $(2M-1)\times(2M-1)$-dimensional principal minors of $D$ we obtain

Appendix D On the computational complexity of the permanent of a matrix with repeated rows and/or columns

One can reduce the number of operations in Ryser’s formula Ryser giving the matrix permanent when the matrix has repeated rows or columns, (see also Appendix B in Ref. Thesis ). Indeed, Ryser’s formula uses the inclusion and exclusion principle of Sylvester, it can be cast as

where $S_{R}\subset\{1,\ldots,N\}$ and has $N-R$ elements (hence, the first term corresponds to $S_{0}$). In Eq. (70) the $R$th term is a sum over the products of row sums of the matrices extracted from $A$ by crossing out $R$ columns. Now let us consider the $N\times N$-dimensional matrix $U[n_{1},\ldots,n_{M}|m_{1},\ldots,m_{M}]$ of section II. Introducing the notation $U_{k_{i}l_{j}}$ for its $(i,j)$th element (where $i$ and $j$ run from 1 to $N$, while $1\leq k_{i},l_{j}\leq M$) we get from Eq. (70)

with the summation over $\{S_{R}\}$ being reduced to that over $\{r_{1},\ldots,r_{M}\}$ with $r_{1}+\ldots+r_{M}=R$ ($r_{l}$ is the number of the $l$th column duplicates crossed out from $U[n_{1},\ldots,n_{M}|m_{1},\ldots,m_{M}]$). One can now easily estimate the number of floating point operations (flops) necessary for computing the permanent of $U[n_{1},\ldots,n_{M}|m_{1},\ldots,m_{M}]$. Indeed, computation of the $R$th term by Eq. (71) involves $T_{R}$ summations over $\{r_{1},\ldots,r_{M}\}$ and $N$ multiplications of a sum comprised of between 1 and $M$ elements $U_{kl}$. The worst case for the last two operations is thus $MN$ flops. On the other hand, the number $T_{R}$ can expressed as

where $P_{N}(z)=\prod_{k=1}^{M}(1+z+\ldots+z^{m_{k}})$. Thus $T_{R}$ is the $R$th term in the Taylor expansion of $P_{N}(z)$ about $z=0$ and with $\Delta z=1$. While $T_{R}$ seems to be given by a complicated dependence on $R$ and $\{m_{l}\}$, their sum, i.e. the total number of flops in the summations over all sets $\{S_{R}\}$, has a simple expression. Indeed, using the fact that the Taylor expansion for $P_{N}(z)$ has only $N+1$ terms, we get

Using this result and the previous estimates on the number of flops in the product and summation inside each sum over $\{r_{1},\ldots,r_{M}\}$, as in Eq. (71), we get that the number of necessary flops $\mathcal{F}$ in Ryser’s formula can be reduced to a value satisfying

(note that by setting $M=N$ and $m_{k}=1$ we get the well-known upper estimate on the number of flops necessary for computing the permanent of an arbitrary matrix by Ryser’s formula: $\mathcal{N}=\mathcal{O}(N^{2}2^{N})$). Let us now analyze the worst case which is obtained by uniformly distributing the bosons over the modes, i.e. when $m_{k}={N}/{M}$ (this maximizes the product in Eq. (74)). We obtain the upper estimate on the number of flops (assuming $N\gg M$, and $M$ fixed)

This result conforms with the general guess, stated in the Introduction, on the number of necessary flops for computing the permanent of a $N\times N$-dimensional matrix of rank $M$.

Finally we note that complete characterization of a network for a given input $\{n_{1},\ldots,n_{M}\}$ is given by probabilities of all possible distributions $\{m_{1},\ldots,m_{M}\}$ of particles in the output modes. Hence, to characterize $N$ bosons on a $M$-mode network with $N\gg M$ and $M$ fixed one has to compute

permanents. Thus the total number of necessary flops to characterize $N$ bosons on a $M$-mode network with $N\gg M$ and $M$ fixed is $\mathcal{F}_{Total}=\mathcal{O}(N^{2M})$.

References