Asymptotic Gaussian law for noninteracting indistinguishable particles in random networks.

For N indistinguishable bosons or fermions impinged on a M-port Haar-random unitary network the average probability to count n 1, … n r particles in a small number r ≪ N of binned-together output ports takes a Gaussian form as N ≫ 1. The discovered Gaussian asymptotic law is the well-known asymptotic law for distinguishable particles, governed by a multinomial distribution, modified by the quantum statistics with stronger effect for greater particle density N/M. Furthermore, it is shown that the same Gaussian law is the asymptotic form of the probability to count particles at the output bins of a fixed multiport with the averaging performed over all possible configurations of the particles in the input ports. In the limit N → ∞, the average counting probability for indistinguishable bosons, fermions, and distinguishable particles differs only at a non-vanishing particle density N/M and only for a singular binning K/M → 1, where K output ports belong to a single bin.

Introduction

Indistinguishable identical particles show correlated behavior due to their quantum statistics even in the absence of interactions: indistinguishable bosons show bunching, e.g., leave the balanced beam splitter in the same port, as demonstrated in the famous experiment with single photons[1] and recently also with massive bosons[2], while indistinguishable fermions show anti-bunching[3]. Fermionic (bosonic) behavior can also be emulated with bosons (fermions) by using entangled particles[4, 5]. Scaling up from two indistinguishable particles on a beam splitter to many particles in large-size multiports significantly increases complexity of behavior, as shown theoretically[6-11] and demonstrated in a series of spectacular experiments[11-16]. Calculating the probabilities or reproducing the statistics of output configurations of noninteracting indistinguishable bosons in a large optical multiport requires exponential resources in the number of particles, which is the essence of the Boson Sampling idea[17] (and its generalizations[18, 19]), where with a few dozens of indistinguishable photons one could have the computational supremacy over the best of current digital computers[17]. The proof-of-principle experiments of several groups open a way to build such a device[20-27].

The behavior of identical particles in multiports is an interplay between the quantum statistics and interference. On the one hand, indistinguishable particles in multi-port networks show correlated behavior beyond the quantum statistics, e.g., in the symmetric (Bell type) multiports complex multi-particle interference results in common forbidden output configurations both for bosons and fermions[10], confirmed recently with photons[16]. On the other hand, the quantum statistics shows up as the two-particle correlations at the Anderson localization in propagation through a disordered media[28-30], induces the “photon clouding” in a unitary optical multiport[24], defines the moments of the output distribution in the scattering of identical particles in chaotic cavities[31], and is responsible for the generalized boson bunching and fermion anti-bunching[32].

A natural question arises: Is there a universal statistics-dependent law in behavior of noninteracting indistinguishable identical bosons (fermions) in multi-port networks? It is shown below that indeed the probability of counting indistinguishable particles in binned-together output ports of a unitary M-port, averaged over the Haar-random unitary matrix representing the multiport, has a statistics-dependent asymptotic Gaussian form as N ≫ 1, where the quantum statistics enters through the particle density N/M. (Our main interest is the bosonic case, since indistinguishable bosons and distinguishable particles can share the same port in a multiport, whereas we also consider fermions to identify the contribution of the quantum statistics). This quantum asymptotic law is reminiscent of the well-known asymptotic law for a multinomial distribution (the de Moivre-Lagrange-Laplace theorem[33] below), which governs an analogous setup with distinguishable (or classical) particles. In the thermodynamic limit N → ∞, the average probabilities of counting indistinguishable bosons, fermions, and distinguishable particles differ only for a finite particle density α = N/M and a singular binning, i.e., K/M → 1 with K output ports belonging to a single bin.

Results

Asymptotic Gaussian law for distinguishable particles in a random network

Consider a unitary M-port network U Fig. 1, whose output ports are partitioned into r bins having K ≡ (K 1, …, K ) ports, with N noninteracting identical particles at the input. Let us begin with the case of distinguishable particles impinging on a multiport. Here by the term “distinguishable particles” we mean the identical particles in different states with respect to the degrees of freedom not affected by a multiport[1, 8, 34], such as the arrival time in the case of photons (e.g., particles sent one at a time through the multiport). The probability for a single particle from input port k to land into bin i reads , where p(l|k) = |U |2 is the probability of the transition k → l. Below we will be interested in the average probability in a random unitary multiport (except where stated otherwise, here and below the term “average” and the notation 〈…〉 means the average over the Haar-random unitary matrix U). A random unitary optical multiport can be experimentally realized with a very high fidelity[25] and without explicit matrix calculations[35]. We have 〈p 〉 = q  ≡ K /M, since 〈|U |2〉 = 1/M (see Supplementary Information). For identical particles sent one at a time through a random multiport, the average probability to count n ≡ (n 1, …, n ) particles in the output bins becomes a multinomial distribution

Figure 1

Unitary network with the output ports grouped into bins. A quantum network, having a unitary matrix U, with N indistinguishable identical particles at its input and binned-together output ports (three bins in this case). We are interested in the probability of counting n = (n 1, n 2, n 3) particles in the output bins. (Two or more bosons as well as classical particles may share the same input port).

Eq. (1) turns out to be also a good approximation to the average counting probability for N distinguishable bosons or fermions at the input (as the numerical results presented below show), e.g., as in propagation of distinguishable particles in the disordered media[28-30] or chaotic cavities[31]. However, Eq. (1) cannot give the exact average probability for these cases, since for N simultaneous distinguishable particles at the input there is an extra factor (see Supplementary Information) due to the correlations between the matrix elements |U |2.

The classic result in the probability theory, below referred to as the de Moivre-Lagrange-Laplace theorem[33], states that given q 1, …, q , such that q  ≠ 0 and ∑ q  = 1, the probability of Eq. (1) takes asymptotically the following multivariate Gaussian form as N → ∞ (more details in Methods)

In other words, for N ≫ 1 and r ≪ N, the numbers n 1, …, n decompose as follows , where y = (y 1, …, y ) is a vector of random variables, constrained by , with the joint probability density ρ(y) of a Gaussian form[36] (see also Methods):

One important note. For a fixed unitary multiport U, averaging the transition probability over the input port gives (the averaging over the input port k and over the Haar-random unitary matrix U give the same result). Hence Eqs (1) and (2) apply also for any unitary multiport with the averaging performed over all possible configurations of distinguishable particles in the input ports (i.e., over k = (k 1, …, k ), 1 ≤ k  ≤ M).

An intriguing question is: What is the quantum version of the above asymptotic result, i.e., what form the analogous asymptotic average probability has if N indistinguishable identical particles are at the input? Since the averaging over the random multiport matrix is performed, one would expect that the multiport-specific quantum interference effects do not show up in the result, whereas the difference in counting the classical (distinguishable) and quantum (indistinguishable) particles in the output bins is solely due to the quantum statistics. It turns out that in the quantum case the asymptotic form of the counting probability is the Gaussian law of Eq. (2) modified by the quantum statistics through the particle density α ≡ N/M.

Asymptotic Gaussian law for indistinguishable particles in a random network

Let us consider bosons and fermions simultaneously (for fermions only up to one particle can be found per network port, as in Fig. 1). For indistinguishable bosons (fermions), the unitary invariance of the Haar measure makes the average probability uniform over the input/output configurations (see Supplementary Information), i.e., for an input k = (k 1, …, k ) and an output l = (l 1, …, l ) configurations the average probability of the transition k → l is just the inverse of the number of Fock states of N bosons (fermions) in M ports: . Summing up the number of the transitions corresponding to a given counting result n of N bosons (fermions) in r bins, we get the following quantum equivalent of Eq. (1) (with the upper signs for bosons and the lower ones for fermions)where the classical 〈P ((n)〉 and quantum factors are separated. It turns out that the quantum factor in Eq. (4) can be approximated as N → ∞ in the same way as the classical probability in Eq. (1), i.e., there is also an asymptotic Gaussian law for Eq. (4). Namely, for any partition of the output ports into r bins, the probability 〈P ((n)〉 has the following asymptotic form as N → ∞ (see Methods for the derivation)with α = N/M. For the Gaussian law (5) to be a good approximation to the probability (4) for N ≫ 1 the number of bins must satisfy r ≪ N, min(K ) (see Methods).

It follows from Eq. (5) that for N ≫ 1 the ith bin occupation number n admits the following decomposition , where (y 1, …, y ) is the same set of Gaussian random variables as in Eq. (3). Therefore, as N → ∞ we get x  ≡ n /N → q , i = 1, …, r for indistinguishable bosons (for a finite density α < ∞), fermions, and distinguishable particles. For bosons there is also the asymptotically infinite density case, N → ∞ for a fixed M ≫ 1, resulting in a qualitatively different (non-deterministic) limit . It is known[37] that in this case the output probabilities p ((l|k; U) and p ((l|k; U) of the transition k → l for distinguishable particles and indistinguishable bosons, respectively, differ even asymptotically as N → ∞ (except for the trivial cases, such as U being a permutation matrix).

The asymptotic Gaussian laws (2) and (5) apply for an arbitrary assignment of output ports to bins, while keeping the total number of ports in each bin fixed. Since the average probability 〈P ((n)〉 (4) depends only on the number of Fock states corresponding to the particle counts n, the asymptotic Gaussian law (5) is valid also for an arbitrary assignment of output configurations l = (l 1, …, l ), allowed by statistics, to bins of configurations enumerated by a vector index n = (n 1, …, n ) (replacing the vector of particle counts in the bins of the output ports), if the total number of configurations l in bin n converges as N → ∞ to the number of configurations corresponding to counting n particles in r bins of the output ports.

Similarly as in the classical case, the probability formula (4) and the asymptotic Gaussian law (5) apply also to any given unitary multiport with the averaging performed over the input configurations allowed by the quantum statistics. Indeed, the quantum transition probability, k → l, satisfies the time-inversion relation p ((l|k; U) = p ((k|l; U †), thus averaging it over the allowed input configurations k produces a result uniform in l and equal to that obtained by averaging over the Haar-random unitary U (see Methods).

Comparison of the asymptotic Gaussian laws with the corresponding probabilities

To compare the asymptotic Gaussian law with the corresponding probability of particle counting let us first give the formula for the respective probability P(n; U) in an arbitrary multiport U. Suppose for the moment that the particles are distinguishable (one can trace the path of each particle through a network by the distinct values of some internal degree(s) of freedom of the particles). Introduce for each output bin a N-dimensional positive semi-definite Hermitian matrix , with k = (k 1, …, k ) being the vector of the input ports of the particles, and assume that a 1, …, a is a partition of (1, …, N) into r subsets, such that vector a contains n indices identifying the particles which end up in the ith bin. Then the probability of counting n = (n 1, …, n ) particles in the output bins becomes

In the quantum case, the fact that it is impossible to trace the individual paths of indistinguishable particles through a network results in the cross-particle interference terms, weighted by a function J(σ) acting on the group of permutations σ of N objects[32, 34] (see Supplementary Information): J ((σ) = 1 for indistinguishable bosons and J ((σ) = sgn(σ) for fermions. For an input configuration k with s particles (0 ≤ s  ≤ N for bosons and 0 ≤ s  ≤ 1 for fermions) in input port j, j = 1, …, M, we obtain (see Methods)

Let us first test the convergence of the average probability 〈P(n)〉 to the corresponding asymptotic Gaussian law as N → ∞. Taking the three bin partition, see Fig. 1, as an illustrative example we numerically compare Eq. (1) to Eq. (2) and Eq. (4) to Eq. (5) (separately for bosons and fermions) by computing the maximum of the absolute value of the difference as a function of n = (n 1, n 2, N − n 1 − n 2) divided by the maximum value of 〈P(n)〉 (both the maximum of the Gaussian approximation in Eqs (2) and (5) and the maximum of the average probability 〈P(n)〉 in Eqs (1) and (4) decrease with N for fixed α and q 1, …, q ). The results are presented in Table 1, where we start with N = 12 and M = 36 (i.e., α = 1/3), increase N by factor 2 and compute the respective M and the partition (K 1, K 2, K 3) for fixed q 1 = q 3 = 1/4 and q 2 = 1/2 (these values are chosen to have integer partitions for all M).

Table 1

The numerically computed for the distinguishable particles, bosons, and fermions for the three bin partition with q 1 = q 3 = 1/4 and q 2 = 1/2.

M	K 1	K 2	N	max∣G^(D)(n)-⟨P^(D)(n)⟩∣max⟨P^(D)(n)⟩	max∣G^(B)(n)-⟨P^(B)(n)⟩∣max⟨P^(B)(n)⟩	max∣G^(F)(n)-⟨P^(F)(n)⟩∣max⟨P^(F)(n)⟩
36	9	18	12	0.096443	0.148050	0.075394
72	18	36	24	0.065787	0.102118	0.037101
144	36	72	48	0.045826	0.068456	0.025218
288	72	144	96	0.031497	0.047765	0.015706
576	144	288	192	0.022142	0.032836	0.010365
1152	288	576	384	0.015580	0.023029	0.006918

The maximum is taken over all distributions n = (n 1, n 2, n 3) of N particles over the output bins.

To see how well the average probabilities (1) and (4) and the asymptotic Gaussian laws (2) and (5) compare with the corresponding particle counting probabilities for a randomly selected multiport matrix U, we numerically estimate the standard deviation employing the following formulawith a large number of randomly selected networks (the error controlled by comparing the results for independent data sets). Figure 2(a,b) gives the comparison of the Gaussian laws with the corresponding average counting probabilities of N = 8 particles in two equal-size bins q 1 = q 2 = 1/2 of the Haar-random multiports with M = 24. The filled regions in Fig. 2(b) give the standard deviation Δ numerically estimated using Eq. (8). Both the standard deviation Δ and also its relative value Δ/〈P(n)〉 decrease with N at the peak n  = q N of the Gaussian law (5), see Fig. 2(c,d), moreover, the standard deviation for bosons rapidly settles on that of distinguishable particles with increase of N.

Figure 2

The Gaussian laws vs. the probabilities of particle counting in two output bins. Panel (a) shows the average probability (4) (filled circles) vs. the Gaussian law (5) (open circles), where bosons (fermions) correspond to the solid (dashed) line, and the classical average probability (1) (filled squares) vs. the Gaussian law (2) (open squares). Panel (b) shows the average probabilities (filled markers) and the Gaussian laws (open markers) together with the standard deviations; for a fixed n, the height of each filled region is twice the respective Δ of Eq. (8). Panel (c) shows the standard deviation Δ (8) at n = N/2 for bosons (the solid line), fermions (the dashed line), and distinguishable particles (the dash-dot line). Panel (d) shows the relative standard deviation Δ/〈P〉 with 〈P〉 ≡ 〈P(N/2, N/2)〉. Here q 1,2 = 1/2 and α = 1/3. In panels (a) and (b) N = 8, M = 24, and 0 ≤ n ≤ 8 corresponds to one of the bins.

Symmetries in the multiport matrix and a large deviation from the average probability

A multiport matrix U with symmetries may result in the counting probability P(n; U) having a large deviation (as compared with the standard deviation) from the respective average probability for particular input configurations. Let us consider as an example the Fourier multiport [6, 38]. In this case most of the output configurations l are forbidden for a cyclicly symmetric input configuration k [10, 16]. For instance, with the help of numerical simulations it is found that for M = N 2 and an input satisfying the cyclic symmetry k  = k  + N (mod M), the equipartition of output ports into two bins (q = 1/2) produces the actual quantum distribution satisfying , i.e., equal to that of Eq. (1). Nevertheless, averaging over all allowed input configurations reproduces 〈P ((n)〉, as it should be. In general, there are unitary matrices having symmetries and possessing many continuous parameters, e.g., the block-Hadamard unitarywhere V is an arbitrary unitary matrix. The unitary in Eq. (9) has forbidden laws for any V. For instance, for N = M/2 particles at the input and the equipartition of the output ports into two bins by the order of their indices, the N-dimensional Hermitian matrices which define the counting probability (7) are drawn from the following 2N-dimensional Hadamard onesdepending on the indices of the occupied input ports. As the result, some of the bin configurations (n 1, n 2) are forbidden for certain input configurations (e.g., for an even number of particles N, only even number of bosons n 1,2 can be detected in either bin, whereas fermions are always divided into two halves n 1 = n 2 = N/2).

The Haar measure for the M-port unitary matrices has M 2 independent real parameters, where 2M − 1 of them are due to the invariance by multiplication with a diagonal unitary matrix[39]. The example in Eq. (10) shows that there exist continuous subsets of the unitary matrices which result in large deviations of the counting probability P(n; U) from the respective average. In general, for a fixed U we have in total (using Eqs (19) and (33) from Methods)values of P(n; U), whereas at most of the values could be independent (for fermions, when K  ≥ N is not satisfied for all i = 1, …, r, some of the values of P ((n; U) are zero by the Pauli exclusion principle). We get the upper bound as N → ∞, with r and α being fixed. Hence, the relative dimension of the subset of the unitary matrices with prescribed independent values of P(n; U) decreases with N at least for r ≥ 4 bins.

Quantum to classical transition in the thermodynamic limit N → ∞

Quantum Statistical Mechanics predicts that, in the thermodynamic limit, for a system of weakly interacting identical particles the classical behavior appears at the vanishing density[40] (α → 0 in our notations). Is there an analog of this rule for noninteracting indistinguishable particles in a random multiport? Let us call the boundary scaling M vs. N the largest scaling leading to the average probabilities in Eqs (1) and (4) being different in the limit N → ∞. Since the average transition probability 〈p(l|k)〉 decreases exponentially with N [17] (see Supplementary Information), to have a finite average probability 〈P(n)〉 one has to sum an exponential in N number of such transition probabilities. This observation together the fact that Eq. (5) applies to an arbitrary assignment of output configurations to the bins of configurations seem to lead to the conclusion that as N → ∞ the difference between the quantum (4) and classical (1) average probabilities vanishes (for a bounded particle density α). However, Eq. (5) and the above conclusion do not hold for a singular binning, the simplest of such is obtained by assigning K < M ports to bin 1 and the rest to bin 2 (or any number of bins), with K/M → 1 as M → ∞ (since q  ≠ 1 is the applicability condition). In this case, setting for the occupied input ports k 1, …, k , we obtain from Eq. (7) the probability to detect all N input particles in K output ports as follows[32] where per(…) stands for the matrix permanent[41]. The probability in Eq. (12) satisfies[32] P ( ≤ P ( ≤ P ( for any multiport and any K < M (this order of the probabilities is reflected also in Fig. 2(a,b) at the end-points n = 0 and n = N). The average values of the probabilities in Eq. (12) were reported before[32]. We have (we use Eq. (1) for distinguishable particles, which is a good approximation, Fig. 2)where the upper (lower) sign stands for bosons (fermions). Let us assign K output ports to a single bin, such that m ≡ M − K is fixed as K, M → ∞ and assume that the input configuration contains up to one particle per input port. We get in the limit N → ∞ for N/M → α ≠ 0 (see Methods):

The result of Eq. (14) implies different survival probability of non-interacting identical particles in some lossy networks: even as N → ∞ the probability to detect N indistinguishable bosons (fermions) and distinguishable particles at the output of a “crowded” M-port network with M = O(N) can be different. Indeed, an arbitrary lossy linear network A (AA † ≠ I) can be equivalently represented by three consecutive multiports according to the singular-value decomposition A = UDV [42], where U and V are two unitary multiports, whereas corresponds to M beam splitters having the transmission coefficients η 1, …, η connecting the multiports U and V, see Fig. 3(a). The commutation rules for the creation and annihilation operators require that 0 ≤ η  ≤ 1 for a lossy multiport (see the Supplement to ref. 32). For a network with m strongly lossy channels η 1 = … = η  = 0 and K = M − m transparent ones (η  = … = η  = 1), by averaging over the Haar-random U, i.e., the path-dependence of the loss of particles (whereas V has no effect on the survival probability) we arrive at the survival probability given by Eq. (13) with the asymptotic results in Eq. (14). The conclusion extends to the networks having some of the loss rates satisfying η ≪ 1, while the others 1 − η ≪ 1, Fig. 3(b).

Figure 3

The survival probability in a lossy network. Panel (a): an equivalent representation of a lossy M-port linear network A = UDV with a unitary 2M-port network, where U and V are some unitary M-ports and the diagonal matrix corresponds to M beamsplitters (with the transmission coefficients η 1, …, η ) placed between them. Panel (b): the average survival probability of N indistinguishable bosons (above the middle dash-dot line) and fermions (below the middle dash-dot line) in a random 4N-port network with m = 3 lossy channels having the loss rates η 1,2,3 = 0.1 (dots on the solid lines) and η 1,2,3 = 0.3 (open circles on the dashed lines). The dash-dot lines give the analytical average survival probability of Eq. (13) for m = 3 completely lossy channels η 1,2,3 = 0, where from top the to bottom we have bosons, distinguishable particles, and fermions, respectively. The averages are computed over 9000 random unitary matrices U.

Discussion

The existence of forbidden events common for bosons and fermions in the symmetric Bell multiports[10] could be viewed as a manifestation of the primary role of the quantum interference over the quantum statistics. Now this view is complemented by the existence of the asymptotic statistics-dependent Gaussian law for indistinguishable particles in multiports for the counting probability in the binned-together output ports averaged either over the Haar-random multiports, or, for a fixed multiport, over the allowed input configurations of the particles (a quantum analog of the de Moivre-Lagrange-Laplace theorem[33] for the multinomial distribution describing classical particles at the input). The asymptotic Gaussian law can have direct applications for noninteracting indistinguishable bosons or fermions in the setups where randomness plays a key role, such as the multiphoton propagation through disordered media[28-30] and the multi-particle scattering in chaotic cavities[31].

The rapid advance of the quantum technology has scaled up the number of identical particles and the size of employed networks in current experiments[11-16], where a wealth of complex behavior is reported. A lot of attention is devoted to indistinguishable photons in multiports[19-27] due to the Boson Sampling idea[17] as a near-future feasible “sampling computer”, an alternative to the universal quantum computer[43] for the demonstration of the quantum supremacy over the digital computers. Recently it was suggested[44] that there might be also some decision problems easily solvable on a Boson Sampling device but not on the digital computers, where binning of the output configurations in a random multiport was proposed for devising such classically hard decision problems. While the complexity of behavior in the linear bosonic networks asymptotically challenges the digital computers, the search for simple asymptotic laws, such as the asymptotic Gaussian law of the present work, brings the dividends of understanding such a complex behavior.

Methods

Derivation of the quantum asymptotic Gaussian law

Derivation of the de Moivre-Lagrange-Laplace theorem can be based on the Stirling’s formula , where for n ≥ 1[45] and θ 0 = 1/(2π). It allows one to approximate the multinomial distribution of Eq. (1) as follows (here for n  ≠ 0)with x  ≡ n /N. Employing a meticulous error estimate[36] as N → ∞ one can then derive the asymptotic result of Eq. (2). The main steps are the following. Observe that the Kullback-Leibler divergence in the exponent in Eq. (15) has a simple expansion for x close to q (its mean value)

For N ≫ 1 and for small |x  − q | the first term on the right-hand side in Eq. (16) is a sufficient approximation for D (x||q) when the former is multiplied by a large negative factor (−N) and placed in the exponent in Eq. (15). In its turn, the factor x in the denominator on the right-hand-side of Eq. (15) can be replaced by q within the same approximation. The result is the stated Gaussian law

The above steps, besides clarifying a basic idea behind the proof of the classical asymptotic result, allow one to justify the analogous quantum asymptotic result by application of exactly the same approximations to the quantum average probability (4). Indeed, in the quantum case, additionally to the classical multinomial distribution of Eq. (15), there is the quantum factor

To estimate the latter, let us use the following resultvalid for all n ≥ 0 for the positive sign and for 0 ≤ n < m for the negative sign. Eq. (19) follows from the second-order Euler-Maclaurin summation formula[46] where and {x} is the fractional part of x. Setting f(x) = ln(1 ± x/m) and using Eq. (20) one can derive Eq. (19). With the use of Eq. (19) the quantum factor is approximated uniformly over n as follows (for fermions N/M ≤ δ < 1)

Using simple algebraic manipulations one can expose in the logarithm of the right-hand-side of Eq. (21) a Kullback-Leibler divergence multiplied by a large negative factor, similar as in Eq. (15). Setting α = N/M, introducing , and observing that 0 ≤ X  ≤ 1 and (an analog of x in Eqs (15)–(16)) we have

Expanding the Kullback-Leibler divergence (the third term) on the right-hand-side of Eq. (22) as in Eq. (16), using that and replacing X by q in the second term (for similar reasons as in the classical case) one arrives at the following Gaussian approximation for the quantum factor

Multiplying the two factors in Eqs (17) and (23) we arrive at the asymptotic Gaussian of Eq. (5).

The counting probability in the output bins of a given multiport

Let us recall that the probability to detect N identical particles in the output ports l = (l 1, …, l ), corresponding to occupations m = (m 1, …, m ), for an input k = (k 1, …, k ), corresponding to occupations s = (s 1, …, s ), reads[32, 34] where a unitary matrix U expands the input ports of a multiport over the output ones and . In Eq. (24) the function J(σ) of the permutation σ of N objects, describing the state of partial distinguishability of the particles, is given as followswhere is an operator representation of σ in the Hilbert space of the internal states of N particles and ρ ( is their internal state[32, 34] (more details in Supplementary Information). Consider now a partition of the output ports of a M-port unitary network U into r bins with K = (K 1, …, K ) ports. Summing up the probabilities of Eq. (24) with the occupations m = (m (1), …, m () corresponding to counting of n = (n 1, …, n ) particles in r output bins giveswhere we have introduced for output bin i a N-dimensional positive semi-definite Hermitian matrixand a partition a 1, …, a of the indices (1, …, N) into r subsets, such that vector a contains n indices. The partition represents the permutations of port indices complementary to the group G of symmetries of the input configuration s (with s! elements) in the group of permutations of N objects. In the case of distinguishable particles we have[34] which allows to further simplify the expression in Eq. (26) by computing the sum over σ. We get the factor s! (which cancels the same factor in the denominator) resulting in

For bosons, in computation of the probability by Eq. (26) one cannot avoid an exponential in N number of the floating point operations. The numerical computation can be carried out using Ryser-Glynn’s formula for the matrix permanent[47], by first computing the sum over σ for a fixed partition a (i.e., a product of the matrix permanents[41]) and then the sum over all possible partitions.

Averaging the probability over the input configurations allowed by the statistics

Unitarity of the multiport matrix U results in equivalence of the averaging the probability p(l|k; U) of the transition k → l over the Haar-random unitary U to the averaging over all the input configurations k, allowed by the statistics, for an arbitrary (fixed) unitary U. Indeed, denoting the respective occupations by s and m, from Eq. (24) we have (see also ref. 48)where U(k|l) is a N-dimensional submatrix of U on the rows k and columns l, per(…) stands for the matrix permanent[41] and the N-dimensional matrix U[s|m] is built by selecting rows and columns of U with the repetitions s and m, respectively. In both cases, the invariance properties of the matrix determinant and permanent imply the following time-inversion symmetry p ((l|k; U) = p ((k|l; U †). On the other hand due to unitarity we must havewhere in the case of bosons we sum over all m satisfying and in the case of fermions over the sets of distinct indices l = (l 1, …, l ). Using the time-invariance symmetry together with Eq. (31) in the definition of the average over all the (allowed) input configurations, we conclude, for instance, thatwhere we have counted the number of Fock states of N bosons in M ports. Since (see Supplementary Information) we obtain the announced result. Similar for fermions.

The asymptotic average counting probability for a singular binning

Let us prove the asymptotic results in Eq. (14). For N → ∞ and N/M → α we obtain (by expanding the logarithm of the left-hand side)

On the other hand, by applying Eq. (19) to the product in the quantum average in Eq. (13) we get

The first factor on the right-hand side of Eq. (34) can be estimated similar as in Eq. (33)

Using Eqs (33)–(35) into Eq. (13) we arrive at the asymptotic result stated in Eq. (14).