Graph isomorphism and Gaussian boson sampling

We introduce a connection between a near-term quantum computing device, specifically a Gaussian boson sampler, and the graph isomorphism problem. We propose a scheme where graphs are encoded into quantum states of light, whose properties are then probed with photon-number-resolving detectors. We prove that the probabilities of different photon-detection events in this setup can be combined to give a complete set of graph invariants. Two graphs are isomorphic if and only if their detection probabilities are equivalent. We present additional ways that the measurement probabilities can be combined or coarse-grained to make experimental tests more amenable. We benchmark these methods with numerical simulations on the Titan supercomputer for several graph families: pairs of isospectral nonisomorphic graphs, isospectral regular graphs, and strongly regular graphs.


I
The problem of graph isomorphism (GI) lies at an interesting point in the landscape of computational complexity theory.Though algorithms have been recently proposed which run in 'quasi-polynomial' time [1,2], it is still an open question in theoretical computer science whether there exists a polynomial-time algorithm that can determine whether two graphs are isomorphic; indeed, graph isomorphism is likely to belong to the class of NP-intermediate computational problems.Two other well-known problems which have similar status in the complexity landscape are integer factoring and the discrete logarithm problem.Famously, while no classically efficient algorithm for these two problems is known, they can be solved in polynomial time on quantum computers [3,4].Quantum algorithms with a superpolynomial runtime advantage have also been proposed for linear systems [5,6], semidefinite programming [7,8], knot invariants [9][10][11], and partition functions [12][13][14], among many others.Given these other success cases, it is natural to hypothesize that quantum hardware may also be useful for the graph isomorphism problem.
Over the last several years, several works have explored this problem, with quantum algorithms for tackling graph isomorphism proposed based on quantum annealing [15][16][17] and quantum graph states [18].However, the bulk of quantum algorithm proposals to distinguish non-isomorphic graphs have utilized the time-evolution of a quantum walker to calculate 'graph invariants' or 'graph certificates' which, ideally, produce identical results for two graphs if and only if they are isomorphic.Of the algorithms proposed, they differ mainly in the number of particles involved, the presence of interactions, localised perturbations, and construction of the GI certificate [19][20][21][22][23].It has subsequently been proven using this approach that conventional quantum walk algorithms, both discrete-time and continuous-time, with interactions and perturbations, cannot distinguish arbitrary non-isomorphic graphs [23][24][25].
To test the distinguishing ability of proposed quantum GI algorithms, a common benchmark has become their capacity to distinguish nonisomorphic strongly regular graphs (SRGs) with the same graph parameters.This provides an analytic approach to investigate the effectiveness of graph isomorphism proposals; if a particular certificate will always fail to distinguish two non-isomorphic SRGs, this can be shown to be because all elements of a certificate, as well as their multiplicities, are functions of SRG family parameters [26].
In this work, we present an approach to graph isomorphism which uses a near-term quantum computational device, namely a photonics-based Gaussian boson sampling apparatus [27,28].For this method, graphs are encoded into quantum-optical states of light -specifically Gaussian states -which are then subjected to photon-number-resolving measurements.Mathematically, we show that the resulting measurement outcome probabilities can be combined to give a complete set of graph invariants.Two graphs are isomorphic if and only if these graph invariants are equal.We also present several ways that these measurement probabilities can be combined and coarse-grained to obtain new quantities which can be used to distinguish nonisomorphic graphs.Finally, we perform classical numerical simulations of our proposed method on the Titan supercomputer.Using these results, we are able to distinguish 3852 out of 3854 nonisomorphic graphs using only a subset of measurement events.The remaining two graphs were distinguished by failing to satisfy a necessary condition introduced here as well.

M
Our main result is a necessary and sufficient condition to distinguish isospectral nonisomorphic graphs by virtue of comparing the probabilities of the measurement patterns of the graphs encoded in a Gaussian boson sampling (GBS) apparatus.We discovered the vital role played by a matrix function called the hafnian [29], applied to an adjacency matrix, for the GI problem.It leads to a complete set of graph invariants.The hafnian belongs to the family of matrix functions such as the determinant, permanent and pfaffian [30].It has been established that photon detection probabilities can be expressed in terms of the hafnians of a collection of graphs related to the original graph [31].Multiphoton detection probabilities are handled by introducing a new matrix product related to the Kronecker product and by showing how the output probabilities depend on the hafnian of the graph adjacency matrix as well.We further strengthen our graph invariant results by introducing the so-called symmetrized graphs invariants and showing that they correspond to coarse-grained measurement events in GBS.The measurement events are given by the stratification according to the total photon number and partitioned into the orbits of the permutation group.Their hafnian-based coarse-grained probabilities are again sufficient to distinguish isospectral nonisomorphic graphs.We extend these insights by deriving necessary conditions for isospectral graphs to be isomorphic by comparing the coarse-grained partition-averaged photon distribution from the Gaussian boson sampler.
Our method differs from previous quantum GI algorithm proposals.A great majority have utilized quantum walks, either using discrete-time quantum walks (DTQWs) [20,21] or continuous-time quantum walks (CTQWs) [23][24][25][26].Although the graph invariants constructed via quantum walk propagation on graph structures have shown success in distinguishing various families of SRGs, it has been proven that this distinguishing power is not universal -there will always exist graphs which (current) quantum walk-based algorithms cannot distinguish [32].In order to execute a quantum walk-based algorithm in a universal quantum photonic platform, it is necessary to implement a non-Gaussian operation as a vertex-dependent shift or via multiple interacting walkers.This is a major obstacle with the current and near-term technology our proposal does not suffer from.GBS is a Gaussian circuit followed by an array of photon-number-resolving detector (PNR) representing a non-Gaussian element.Unlike non-Gaussian unitary transformations, the PNRs are available in the state-of-the-art laboratories.
The most comprehensive simulations of quantum methods for GI were performed in [19] and [23].We successfully tested three types of isospectral graphs: pairs of isospectral nonisomorphic graphs (PINGs) as the first examples of such graphs [33], isospectral regular graphs [34] and mainly SRGs.There are numerous resources available detailing the SRG families containing more than one non-isomorphic graph [35,36]; as a result, SRGs have become a common benchmark in studying the distinguishing powers of the GI algorithms.Note that there may be other graph classes (such as k-equivalent graphs) which have been proven to be harder to distinguish than strongly regular graphs for particular quantum GI algorithms [32] -however, SRGs remain an ideal testing set, simply due to the large number of relatively small non-isomorphic graphs present in specific families [35,36].The largest tested and distinguished family by our approach was SRG (35,18,9,9) containing 3854 isospectral graphs.This family is supposedly tested in [23].However, the size of the family is mistakenly taken to be only 227 graphs (see Table I.).The same error appears in [19].Ironically, another SRG family considered there (SRG (35,16,16,8)), that happens to be complementary to SRG (35,18,9,9) and thus containing 3854 graphs as well, is counted properly and analyzed (see Table 1.).
In Section 3 we informally introduce the hardware setup (Gaussian boson sampler) where studied graphs are encoded.Section 4 contains all necessary definitions and previous results used in the paper including a detailed GBS description and a formal introduction of SRGs.Section 5 contains the main result and is split into four subsections: In 5.1 we gather several supporting results followed by the main results in Sections 5.2, 5.3 and 5.4.Section 6 contains the simulation results and Section 7 concludes with a scalability discussion and other open questions.

B
The basis for our graph isomorphism method is a near-term photonic quantum processor, specifically a GBS apparatus.This apparatus consists of three main components.First, squeezed states are generated in M quantum-optical modes.These states are then sent through an M -port linear-optical interferometer.Finally, a photon-number-resolving measurement is performed on each of the M output modes.The first two steps lead to the preparation of a zero-mean quantum-optical Gaussian state, which can be described efficiently using a covariance matrix.For a single mode, the covariance matrix has dimension 2 × 2, and encodes the covariances of the canonical quadrature operators (x, p) of that mode: with ξk ∈ {x k , pk }.For M modes, we have M pairs of quadrature operators and an 2M × 2M covariance matrix built from the set ξ ∈ {x 1 , p1 , . . ., xM , pM }.Multimode Gaussian states themselves are of limited interest in quantum computing.While they can be prepared with quantum hardware and exhibit entanglement, the covariance matrix scales linearly in the number of modes, so they can be efficiently simulated classically.However, when we introduce the photonnumber measurement, the story changes.A single photon-number measurement in mode k will return a nonnegative integer n k ∈ + , representing the number of photons which were detected.For measurement on M modes, we denote the collective photon-number output pattern by n = (n 1 , . . ., n M ) and call it a detection event.From [27], whenever n i = 1, ∀i the probability of this detection event is proportional to a function called the hafnian [29] (see Def. 1): where the matrix C is obtained from σ Q by basic matrix transformations (see Eq. ( 6)).
Unlike the simulation of Gaussian states, computing the hafnian is a #P-hard problem.In addition, approximating the GBS photon-number distribution is believed to be computationally hard [27].Thus, by combining Gaussian states with photon-number measurements -representing the wavelike and particle-like properties of light, respectively -we have a physical sampling apparatus whose behaviour is classically hard to replicate.The remainder of this paper will explore the question of how we can leverage this GBS device for the graph isomorphism problem, specifically, how we set the squeezing and interferometer parameters to represent the problem, and how to interpret the photon-number measurement outcomes to solve the problem.

N
In the following text the symbol k, denotes an all-ones rectangular matrix of size k × and k ≡ k,k .The following notation is extensively used: The symbol df = stands for 'defined' and a positive-definite matrix A will be denoted by A 0. Recall that any Gaussian n-dimensional real distribution with zero mean, denoted as G Σ , is given by Here, Σ is a positive definite matrix which is the covariance matrix of the Gaussian variables X 1 , . . ., X n .Since Σ is positive definite, there exists a unique positive definite matrix A such that Σ = A 2 .Let us change the variables y = A −1 x .That is, x = Ay.Hence the determinant of the Jacobian is det A. As det Σ = (det A) 2 we get that the density distribution for (Y 1 , . . ., Y n ) is the standard normal density distribution Therefore Y 1 , . . ., Y n are independent standard random variables.Assume that A = [a i j ] is a positive definite symmetric matrix.Then and Observe the well known fact that the odd moments [ Recall that two square matrices A, B are permutationally similar, if B = PAP , where P is a permutation matrix.In this case P −1 = P .Two Gaussian distributions corresponding to positive definite covariance matrices Σ, Σ ∈ n×n are called isomorphic, if Σ = P ΣP for some permutation P = P(σ).That is Denote by H N ⊃ H +,N ⊃ H ++,N the real space of N × N hermitian matrices, the closed cone of positive semidefinite hermitian matrices, and the open set of positive definite hermitian matrices.For is the hafnian of C [29].
Given a detection event n, its measurement probability was shown in [27] to be where which we view as a column vector (even though β is a complex conjugate of a complex β (entrywise), we consider β i as a new variable).We denote Then and σ = σ Q − 2M /2 is the covariance matrix in Eq. ( 1).Note that in order for σ to be a covariance matrix, C has to satisfy certain restrictions which will be elaborated on later.We call a GBS detection event corresponding to the measurement pattern n df = (n 1 , . . ., n M ) of an M -mode matrix C of size 2M pure if n i = n j , ∀i, j and mixed otherwise.For n = (1, . . ., 1) Eq. ( 4) reduces to (2) [27].

Isospectral graphs -Pairs of isospectral nonisomorphic graphs and strongly regular graphs.
In studying the graph isomorphism problem, it is convenient to consider a class of graphs known to be classically intractable to distinguish.An important tractable feature is the graph eigenvalues and the first examples of isospectral graphs were pairs of isospectral nonisomorphic graphs (PINGs) [33].The smallest connected example of a PING is on six vertices, see Fig. 1.PINGs may have some tractable features enabling F 1. PING on six vertices.
one to easily decide that they are not isomorphic.The situation gets complicated for graphs with symmetries such as regular graph or strongly regular graphs (SRGs), defined as follows [37].
Definition 2. Let G(V, E) be a a regular graph of degree k consisting of N vertices and adjacency matrix A, that is neither a complete graph (A = N − N ) nor a null graph.G is then said to be strongly regular with parameters SRG(N , k, λ, µ) if every pair of adjacent vertices have exactly λ common neighbours, and every pair of non-adjacent vertices have exactly µ common neighbours.
. Thus the SRG parameters are not independent.
Proof.Consider a vertex v in a graph with parameters SRG(N , k, λ, µ), and count in two ways the number of edges from vertices adjacent to v to vertices non-adjacent to v.
If multiple non-isomorphic strongly regular graphs share the same set of SRG parameters, we refer to this graph set as an SRG family, often simply denoted by the SRG parameters (N , k, λ, µ).Graphs within the same SRG family share various properties that are dependent only on the SRG parameters.The spectrum is rather special, consisting of just three eigenvalues with known multiplicities.SRG families contain nonisomorphic graphs which are isospectral, and difficult to distinguish using common classical measures [37].

G G
5.1.Multiphoton contributions in GBS -supporting results.We start with the exploration of how to interpret Eq. ( 4) for the detection events n where n i > 1 for some i.This corresponds to a multiphoton contribution of the output probability function.The multiphoton contributions play a vital role in our analysis.Example.Let n = (3, 2, 1, 4) and A an adjacency matrix of a simple weighted graph (without loops).Then .
The reason for introducing a new kind of structure is a compact expression for the probability of measurement of a mixed multiphoton event n as a hafnian function not unlike Eq. ( 2) for n = (1, . . ., 1).Definition 4. A 2M × 2M -dimensional real symmetric matrix R will be called GBS encodable if we can find a covariance matrix Ref. [31] introduced a necessary criterion for R to be GBS encodable.For some real symmetric R not satisfying the conditions a general procedure was created to produce a matrix related to R that is GBS- Even though this procedure is always guaranteed to succeed in creating a Gaussian covariance matrix, it is not a necessary condition.Here we strengthen this previous result by loosening the requirements on R.
Then the following conditions hold if and only if F 2 < 1 and F 22 0. Proof.
(1) Clearly we need to have that The assumption that . This is equivalent to λ N (F ) > −1.Therefore the claim follows.
(2) Use (1) to get the assumption that G 0 is equivalent to This inequality is equivalent to Observe next that Hence Corollary 3. Let R ∈ 2M ×2M be a nonzero real symmetric matrix with the following partition to M ×M blocks: Then there exists a Gaussian covariance matrix σ such that cR = X 2M [ 2M − (σ + 1 2 2M )] if and only if: Part (3) follows from the observation the singular values of R and Remark.The previous lemma was presented for the sake of completeness.In the rest of the paper we use the construction of GBS-encodable matrices introduced already in [31].
For future reference we recall the following straightforward result [31].

Lemma 4. Let M be even and C
where c > 0 and k ∈ .
is a symmetric and GBS-encodable matrix.Then the probability of sampling in the GBS event, Eq. (4), can be expressed as Proof.To prove this equality we assume that C = A ⊕ A where A is an arbitrary symmetric matrix of order M .Let p(n) be defined by the right-hand side of (4).We will show that p We now consider the quadratic form γ (C / ⊗ 2|n| )γ and focus on β (A/ ⊗ |n| )β (the 'second half' is treated in exactly the same way).We substitute in (4).We define its action for the quadratic form to be where we used Def. 3 by setting B i j = n i ,n j .We further set ∂ , ∀i (and similarly for the conjugated variables β i and the corresponding α in i ) and write It remains to show that p(n) = p(n).Indeed, the higher-order partial derivatives in (4) result in the same expression as the first-order ones in ( 13) whenever we set α = β = 0 at the end of the calculation.This is a consequence of the elementary properties of the differential operator, namely, and the chain rule for the n-th derivative given by Faà di Bruno's formula for . Now we put all the pieces together.The RHS of ( 14) is identified with (13) through (x k1 , . . ., x kn k ) → (α k1 , . . ., α kn k ) (or its conjugate) for a given 1 ≤ k ≤ M and so f ( But then, according to the LHS of ( 14) we may write (13) as The RHS of ( 16) is identified with the LHS of ( 15) by setting forms p(n) and from ( 13) we get . This follows from (15) and the definition of h(x).
Interestingly, many of the detection events have probability zero: Proof.Here we may assume c = 1 and so C = A ⊕ A. The probability expression (Eq.( 10)) contains haf [A/ ⊗ |n| ] and if n i > |n|/2 then A/ ⊗ |n| contains a zero matrix of size greater than M /2 (placed in the lower right corner of A/ ⊗ |n| ), see Definition 3. The hafnian of such a matrix is zero.This follows from the hafnian definition (Def. 1) where the hafnian of a 2M × 2M matrix is a sum of products of M entries a i j .Since i < j and none of the indices repeats for any summand then inevitably at least one of the a i j 's in every summand equals zero.
We prove several useful properties of the reduced Kronecker product / ⊗. where We then observe that dim [C / ⊗ 2|n| ] = 2|n| and the claim follows.

Remark.
For Lemma 8. Let A, B, C be matrices and assume that A/ ⊗C and B / ⊗C are defined.Then Proof.The ordinary Kronecker product satisfies . By removing columns and the corresponding rows from A, B on both sides of the expression we arrive at the claimed result.

Lemma 9. Let P π be a permutation matrix. Then the following diagram commutes for any matrix A and a detection event
where Pπ is another permutation matrix and m = π(n).
Proof.Following the lower route, the spanning basis β = (β 1 , . . ., for Ã.The new basis is then expanded by considering and so π(α) is a spanning basis of Ã/ ⊗ |m| .Going through the upper route, we observe that β → α by the action of Then, a permutation matrix exists transforming (18) into (17).Its construction is straightforward.The reordering (permutation) (α i1 , . . ., the overall transformation is an action of a permutation matrix we denoted by Pπ .
We use another result from [31] to prove the following lemma.
Proof.In order to encode an arbitrary graph we take two copies of a graph's adjacency matrix [31].Also, it is advantageous to rewrite (6) as The matrix A i commutes with X 2M [31] and so the eigenvalues of X 2M A i are products of eigenvalues of the constituents.Furthermore, using the second equality in (19) and σ Q,G i = σ A i + 2M /2 we conclude by a direct calculation that the eigenvalues of σ Q,G 1 and σ Q,G 2 coincide.The claim follows from the fact that the determinant is a product of eigenvalues.

GBS and a complete set of graph invariants.
Remark.We will use the transformation A → C = A⊕A for the application of GBS to the graph isomorphism problem.We recall C ∈ 2M ×2M .By 'doubling' A, one copy is conveniently spanned by (β 1 , . . ., β M ) whereas the second one by (β 1 , . . ., β M ).Moreover, to make the matrix C GBS-encodable (see Def. 4) we simply take The additional multiple of an identity on the diagonal does not affect the hafnian of A⊕ A as follows from Lemma 4 but it will become useful in the next sections.

GBS and Moments of Multivariate Gaussians.
The moments µ n 1 ,...,n 2M (Σ) of a (zero-mean) 2M -dimensional multivariate real normal distribution N (0, Σ) are given by the following formula: where Σ is the covariance matrix.This follows from the fact that exp [ 1  2 x Σx ] is the moment-generating function of the multivariate normal.
Let R = c(A ⊕ A + k 2M ) be GBS encodable for k sufficiently high so that R 0. Then is exactly the moment µ n (R) of the 2M -dimensional (zero-mean) multivariate normal distribution N (0, R) if we ignore the prefactor (n! det σ Q ) −1 .For clarity, we have changed variables so that x j = β j and x M + j = β j for j = 0, . . ., M .From the above equations, it is clear that the different photon-counting probabilities of a GBS setup are directly related to various moments of a multivariate normal distribution.Importantly, however, they do not give us all moments (µ n 1 ,...,n m ,n m+1 ,...,n 2M (R)), but rather the smaller set (µ n 1 ,...,n M ,n 1 ,...,n M (R)).This is something we need to be careful of.
The moment-generating function factorizes: where we have used the notation x (M ) df = (x 1 , . . ., x M ) and x (2M ) df = (x M +1 , . . ., x 2M ).We set c = 1 (since we omit the determinant prefactor where it otherwise plays a role) and also k = 0.This step will cost us the positive-definiteness of A but at the moment this is just a formality to properly define the moment generating function.We would have set k = 0 afterwards anyway to recover the correct probability expression.The moments of this factorized distribution are then Rewriting, we find where µ n 1 ,...,n M (A) and µ n M +1 ,...,n 2M (A) are moments of the M -dimensional normal distributions N (0, A).
Connecting back to photon-counting probabilities, we recover c and conclude that, for the considered case of block-diagonal R, Finally, we note that the moments are exactly the hafnian of some appropriate matrix, so where the second equality also follows from Lemma 4 and 8.
Proposition 11.Suppose we have two isospectral graphs G 1 and G 2 .Assume we can encode the adjacency matrices A i of either graph into a Gaussian boson sampling setup.Then these graphs are isomorphic iff the hafnians are related by a permutation, haf , for all n.Furthermore, the permutation π must be the same for all n.
Proof of ⇒: Suppose G 1 and G 2 are isomorphic.Equivalently, their adjacency matrices are related by a permutation where P π(i)i = δ i,π(i) for some permutation π.If we encode these adjacency matrices directly into the covariance matrices of two Gaussian states, then graph isomorphism is equivalent to the multivariate normal distributions corresponding to these two Gaussian states being related by a permutation of coordinates: All moments of these 2M -dimensional distributions must correspondingly be related by the permutation π ⊕ π, where R i = A i ⊕ A i are made into GBS encodable matrices (we keep on omitting the c factors).Looking back to Eq. ( 24), we have These moments must be equal for any choices p = (n 1 , ..., n M ) and q = (n M +1 , ..., n 2M ).Thus, we conclude that In particular, for p = q = n, where n is arbitrary, we get We now use the fact that adjacency matrices A contain only 0s or 1s, so haf [A/ ⊗ |n| ] ≥ 0 for any possible A or n.This leads to which proves the statement.Proof of ⇐: Eq. ( 32) immediately implies Eq. ( 31), even when the hafnians are not positive.Furthermore, Eqns.( 27)-( 31) are all equivalent.Hence, the graphs having adjacency matrices A 1 and A 2 are isomorphic.are the same for the two graphs for all possible n.

GBS and Symmetrized
To avoid a notational clash in this section we use n σ instead of σ(n) used in the previous sections.We start with the following result.
Theorem 13.The following statements are equivalent for two Gaussian distributions with zero mean and positive definite covariance matrices Σ, Σ ∈ n×n : (1) The two Gaussian distributions are isomorphic.
(3) The matrices Σ −1 and (Σ ) −1 are permutationally similar.( 4) For each homogeneous symmetric polynomial p(x ) of even degree the expected value of p(x ) is the same for the two Gaussian distributions.(5) The symmetrized moments of the two Gaussian distributions are the same.
Proof.For a given Gaussian distribution with zero mean and the covariance matrix Σ let us denote by H Σ the distribution with the following density By G Σ we denote the density function of the normal distribution exp Clearly if Σ and Σ are permutationally similar then h Σ (x ) = h Σ (x ) for each x .So we need to prove the other direction.Let p(x ) be a monomial n of even degree.We denote m = (m 1 , . . ., m n ).Consider the moments Thus if two Gaussian distributions are isomorphic, then H Σ = H Σ and the corresponding density functions are the same.In particular, the moments of H Σ and H Σ are the same.
Our first main result is the claim that if the moments of H Σ and H Σ are the same then H Σ = H Σ .It is not true that the equality of the moments yield that the distribution are the same.However, it is true for all distributions that have the form of H Σ .Indeed a sufficient condition is M (u) = [exp [〈u, X〉]] is a well-defined vector for all u ∈ n [38].This condition is satisfied for H Σ , where Σ is a positive definite matrix.
It is left to show that if H Σ = H Σ then Σ and Σ are permutationally similar.We first analyze the behavior of the n! numbers for a fixed but arbitrary vector x ∈ n and σ ∈ S n .
Let orb Σ be all pairwise distinct matrices of the form P (σ)ΣP(σ) for σ ∈ S n .Recall that | orb Σ| divides n! and n!/| orb Σ| is the cardinality of the automorphism group, i.e., all σ ∈ S n such that P (σ)ΣP(σ) = Σ.Let A 1 , A 2 be two real symmetric matrices of order n.Define two corresponding quadratic forms Hence for generic, (randomly selected x ) we have that f 1 (x ) = 2 (x ).Similarly: let A 1 , . . ., A k be k pairwise distinct symmetric matrices.Set Since any two pairs of matrices in orb Σ are pairwise distinct it follows that for a generic x we will have exactly | orb Σ| distinct values in (33) and each value is repeated n!/| orb Σ| times.
Assume now that for each x ∈ n in (33) we have the equality h Σ (t x ) = h Σ (t x ): for some fixed t ≥ 0. Fix x in general position.Then for t = 1 we get that the number of distinct values in (33) is | orb Σ| for Σ and | orb Σ | for Σ , respectively.Let for σ ∈ S n .In the equality (35) set t = k for k = 0, 1, . . ., n!.Thus we have the equalities: for k = 1, . . ., n!.These equalities yield that the two multisets {a(σ), σ ∈ S n } and {a (σ), σ ∈ S n } are the same.Hence the n! moments of discrete distributions equally distributed on n! points given in (34) for Σ and Σ are the same.Hence these two multisets are the same.First it yields that | orb Σ| = | orb Σ |.Moreover there exists P(σ) such that x P (σ)Σ −1 P(σ)x = x (Σ ) −1 x .Moreover, for each Σ i = P i ΣP i in the orbit of Σ (under the action of the group of permutations) we have a permutation x .Now if we change x to y we still have the same equality y P (σ)Σ −1 P(σ)y = y (Σ ) −1 y.This finally shows that P (σ)Σ −1 P(σ) = Σ .So indeed the covariance matrices are permutationally similar.
Let e i = (δ i,1 , . . ., δ i,M ), where δ i, j is the Kronecker delta function.Assume that n = (n 1 , . . ., n M ) ∈ M + .Let B = [b i, j ] be a real symmetric matrix and recall the "n-th moment corresponding to B" from the beginning of this section as ), the n-th moment of the Gaussian distribution given by the covariance B, is equal to µ(n, B) up to a multiplicative constant.) To proceed, we also recall the generalization of the classical Leibniz's formula of the derivative of the product of m functions in one variable: . Then Leibniz's formula yields the multilinear Leibniz's formula: where n a 1 , a 2 , . . ., a m = M j=1 n j a 1, j , a 2, j , . . ., a m, j (40) and a i = (a i,1 , . . ., a i,M ) and i ∈ [m].We now apply this formula for f (x ) = 1 2 x Bx and m = |n|/2.Then in (39) we need to consider only the case where We have two kinds of a i .Namely, either a i = 2e p or a i = e p + e q , where 1 ≤ p < q ≤ n.Let us discuss briefly all possibilities for the decomposition of n as n = |n|/2 i=1 a i .We claim that the set {a 1 , . . ., a |n|/2 } corresponds to the following multigraph G = G(a 1 , . . ., a |n|/2 ) with multiple edges and self loops.Each a i = e p + e q , where 1 ≤ p < q ≤ n corresponds to an edge {p, q}.Each a i = 2e p corresponds to a self loop on vertex p (the degree of a self loop is 2).So A(G) = [c pq (G)], the adjacency matrix of G, is a symmetric matrix whose entries are nonnegative integers with the following properties.Each diagonal entry c pp (G) is an even integer.c pp (G)/2 is the number of a i of the form 2e p .For 1 ≤ p < q ≤ M the integer c pq (G) is the number of a i of the form e p + e q .
We denote by S M ,0 all M × M symmetric matrices with zero diagonal.Assume that A ∈ S M ,0 .Note that for f = 1/2(x Ax ) we get that ∂ 2 i f = 0 for each i.
For each ( j 1 , . . ., j k ) ∈ F k (n) and i ∈ [M ] denote m i ( j 1 , . . ., j k ) the number of j l that are equal to i.So M i=1 m i ( j 1 , . . ., j k ) = k.Lemma 15.Let A ∈ S M ,0 and t ∈ are given.Assume that |n| is even.Then Here µ(0, A) = 1.
Proof.We consider the formula (41 Let k be the number of terms in {a 1 , . . ., a |n|/2 } of the form 2e i .The contribution of all terms {a 1 , . . ., a |n|/2 } for which k = 0 is µ(n, A).Let us assume that k ∈ [|n|/2].Without loss of generality we can assume that n i !m i ( j 1 ,..., j k )!2 m i ( j 1 ,..., j k ) (n i −2m i ( j 1 ,..., j k ))! in front of µ(n − 2 k l=1 e j l , A).This comes from the equality We need to see this equality on the level of the derivative with respect to the variable x i .Let m i = m i ( j 1 , . . ., j k ).If m i = 0, then (n − 2 k j=1 e j l ) i = n i for the coordinate i we have obvious equality.Assume now that m i ≥ 1.Then on the left-hand side of the above equality we the factor On the right-hand side we have the factors: We are now ready to give the proof for the main result of this section.
Example.Let us illustrate Lemma 15 and Theorem 12 on an example of a graph whose adjacency matrix is A⊕ A of size 2M = 6 for orbit of n = (2, 3, 3).Since |n| = 8 it is clear that only the fourth power (|n|/2 = 4) of (x (M ) ) A(t)x (M ) survives an encounter with the partial derivatives.So, following Lemma 15, we write Each t k coefficient corresponds to a polynomial of the matrix entries in the exponential of some µ(n − 2 k l=1 e j l , A) ≡ µ(m, A).In accordance with Eq. ( 43) we get, for example for t 2 , m = (2, 1, 1) since 1 16 = a 12 a 13 .
As the final step, we symmetrize the orbit represented by n (in this case the orbit size equals 3) which causes a permutation of indices in (45b).
It is advantageous to stratify the measurement events of an M -mode interferometer according to the total photon number |n| ≥ 0. Once M and |n| are fixed, all possible detection events can be split into the orbits O i (equivalence classes under permutation) that partition the set of all events for a fixed M and |n|.We choose the class representative to be a detection event n = (n j ) M j=1 such that n i ≤ n j , ∀i, j and denote by G n its stabilizer.Clearly G n ⊂ G = S M and the orbits are generated by the left action of the coset G/G n .In order to find the orbits with a great number of detection events (presumably the most likely ones) we count the orbit size according to , where k j are the multiplicities of the j-th photon events satisfying j=0 jk j = |n| and ≤ M .The probability of measurement of a given pattern (n 1 , . . ., n M ) is given by p(n) in Eq. ( 10) (or, more precisely, by its doubled version, Eq. ( 26), see Lemma 4 and the remark on page 10), where How does the number of orbits increase with |n|?This is equivalent to the question of integer |n| partition, that is, in how many ways one can write where the order of the sum plays no role and the number of parts is 1 ≤ m ≤ M .We naturally order the parts such that λ i ≤ λ i+1 , ∀i.Suppose M ≥ |n| first.No closed formula is known but the generating function for integer partition provides the number of orbits for a given |n|.Also, very precise estimates have been uncovered and the growth of the number of orbits is exponential in |n|.For M < |n|, not all number partitions can be realized and the counting is given by the generating function for the number of integer partitions into exactly M parts.Note that we only partition even numbers in this paper, since GBS assigns zero probability for odd |n|.
Curiously, if we try to coarse-grain the probability distribution further and introduce the partition probability where the first sum on the RHS is over the partitions of |n| and the second sum over the orbit elements.We find for all |n| whenever the graphs G 1 , G 2 are isospectral.
Proof.Any undirected graph G on 2M vertices can be encoded as a pure covariance matrix whose circuit decomposition consists of an array of M single-mode squeezing transformation S(r k ) (0 ≤ k ≤ M ) followed by an M -mode linear interferometer U [31].For each S(r k ) we find and so where |(n, M )〉 carry all completely symmetric representations of su(M ) (each representation labeled by n/2).Given λ k (A) and 0 < c < 1/ A 2 for G's adjacency matrix A ( [31], see also Lemma 2) we can write cλ k = tanh r k .Therefore α in (r 1 , . . ., r M ) = α in (λ 1 , . . ., λ M ).Since the interferometer U preserves the number of particles |n| the partition probability p G (|n|) is unaffected by it.Then However, the RHS is independent on the graph (depends only on λ k common for isospectral graphs) and the claim follows.
Remark.Note that the similar argument does not hold for the less coarse-grained probability p G (O n ) in Eq. ( 46) since the interferometer 'mixes' the orbits.
Even though the coarse-grained probability, Eq. ( 48), cannot be used to distinguish nonisomorphic graphs, not all hope is lost.Possible strategies and the closely related problem of scalability is discussed in Sec. 7.

Modifying the results for C = A ⊕ A and beyond.
Given A of even order consider p(n, C) in (4) and µ(n, C).If for the two graphs A and B we have the equalities for the symmetrized sums what can we say?If instead of considering just the matrix A we will consider the matrix A ⊕ A then we can conclude that our bigger graph is a disjoint union of two isomorphic graphs.So if the union of two isomorphic graphs is isomorphic to the union of another two isomorphic graphs, then the two graphs are also isomorphic.If we consider the functions µ(m, A ⊕ A).Note that m = (m 1 , . . ., m 2M ) = (n, n ) where n = (n 1 , . . ., n M ), n = (n M +1 , . . ., n 2M ).It now follows that Therefore, if we have equalities for the symmetrized sums, we get a whole hierarchy of necessary conditions by considering A ⊕k for k > 2.
5.4.Partition-averaged photon distribution as a necessary condition for graph isomorphism.The main result of this section will be a simpler necessary condition for isospectral graphs to be isomorphic.

Definition 5.
Let A be the adjacency matrix of an M -vertex graph G and 1 ≤ k ≤ M .Then we call the partition-averaged photon distribution of the k-th detector the function and its coarse-grained version reads Theorem 18.The partition-averaged photon distributions introduced in Definition 5 of two isomorphic graphs are identical up to a permutation of output modes which can be verified in polynomial time.

Lemma 19.
Let G A and G Ã be isomorphic graphs.Then the output probability distribution from GBS with encoded graphs is related by a permutation.
We prove the same statement by using Eq. ( 4) where C = c(A ⊕ A), C = c( Ã ⊕ Ã) and we can ignore c here by setting c = 1.We introduce P df = P π ⊕ P π and write But that implies that the probability of a pure event remains the same since Pβ by definition merely relabels the output modes and the partial derivatives do not care: For mixed detection events the situation is different.If one of the n i 's in is different from the rest, the corresponding partial derivative breaks the symmetry and unlike the pure case one concludes that However, if we permute the derivative variables (symbolically written as ∂ β i ,β i → ∂ (P π β i ),(P π β i ) ), we find the desired equality Next, using map (18), we rewrite the both sides of the last equation as where δ df = (α, α) and P df = Pπ ⊕ Pπ was introduced in Lemma 9. We used Eq. ( 10), Lemma 8 and Lemma 9 (the upper route in the commutative diagram to go from the LHS of (56) to the LHS of (57)).But since δ and P is a permutation, the hafnian is preserved and the output probability distribution is merely permuted.
To conclude the proof we notice that the overall detection probability is a sum of invariant (for pure events) or permuted (for the mixed ones) probability distributions where the permutation is the same for all mixed n's.
To simplify the notation in the rest of the section we write haf 2 G (n) ≡ haf 2 [A/ ⊗ |n| ] in Eq. ( 52).Given the stratification into orbits, it is advantageous to collect n k together with the factorial coefficients and the (squared) hafnians of a graph G to N and haf G , respectively, and rewrite (52) as where n G is M -tuple of numbers.
Example.Let us illustrate (58) for a graph G on M = 4 vertices and for |n| = 2.There are two orbits represented by (0, 0, 0, 2) and (0, 0, 1, 1).Since the graph is doubled, we have M = 4 detectors and then We can clearly identify the sums on the RHS of (52).Note that due to Corollary 16 the hafnians of the (0, 0, 0, 2) orbit are zero and so are the corresponding contributions to n G .
The following proof is best viewed together with the above example.
Proof of Theorem 18. Lemma 19 shows that, if graphs G 1 and G 2 are isomorphic, then the ordered set of hafnians for one graph is a permutation of the same ordered set for the other graph.This statement translates into a permutation of haf G introduced earlier: haf G 2 = π(haf G 1 ).Note that the pure orbit elements are fixed points of π.Now we observe that the i-th row of N by construction coincides with the sequence assembled from the i-th elements of all n (see (59)).So instead of swapping the rows of N we correspondingly swap these sequences in the argument of all haf 2 G forming haf G .But this transformation is a permutation of the set of all n's for a fixed |n| since it preserves the photon number.Hence where P π swaps the rows of N.
and thanks to Lemma 10 we can rewrite the equality as But the LHS is n G 2 and the action of a permutation P π on the RHS implies that it is equal to π(n G 1 ).Therefore which is the main result.
To conclude the proof we observe that it takes only a polynomial number of steps to uncover how the partition-averaged photon distribution is permuted.We order n G 1 and n G 2 in an increasing order and if the two ordered sets differ the graphs cannot be isomorphic One could be tempted to argue that the opposite is true (that is, if the partition-averaged photon distributions the same then the graphs are isomorphic).The following counterexample shows that there is no hope for the converse of Theorem 18.

Example
Remark.Note that since the hafnian sets differ in the previous example we know that the graphs are not isomorphic.It just can't be concluded from comparing the partition-averaged photon distributions for |n| = 4 and it can't even be concluded from (46) since p G 1 (O n ) = p G 2 (O n ) for all orbits for |n| = 4 (including n = (1, 1, 1, 1) again!).The first differences both in 〈n k 〉 G and p G (O n ) appear for some orbits of |n| = 8.Interestingly, 〈n k 〉 G is always uniform for SRGs and when it differs for two nonisomorphic SRGs, it differs in a magnitude.

S
In the following section, we present the results of the GBS quantum GI algorithm applied to various SRG families and other isospectral graphs.The algorithm itself is presented in the Appendix.Among other graphs, we examine the SRG (35,18,9,9) family, and show that, using various detection patterns, the GBS fully distinguishes all 3854 graphs in this family.Due to the large number of photon event permutations required to calculate the probability of detection, and the classically intractable graph hafnian calculation, the results were computed in parallel using the Python Hafnian library [39] and the Titan supercomputer 1 .
Example.Let us start with the smallest connected PING in Fig. 1.The hafnians of the adjacency matrices coincide so we have to look to all possible GBS-measurable submatrices.These correspond to all measurement patterns with at most one photon per mode (in this example we won't study the multiphoton contributions coming from A/ ⊗ n ).Hence, we can measure in total 6 4 = 15 graphs on 4 vertices as well as 2 vertices (the subgraphs with an odd number of vertices have zero perfect matchings and therefore zero hafnian).The hafnians of the latter (let's call them 2-hafnians) do not differ but the 4-hafnian sets do differ: 4-haf A G 1 = (0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2) ≡ (0 7 , 1 7 , 2 1 ), (63a) Example.Consider another example of a PING [33] in Fig. 2, this time on nine vertices.Their hafnians are F 2. PING on nine vertices.
zero since the number of vertices is odd but also all (8,6,4,2)-hafnians are identical (that is, all subgraphs accessible to GBS have the same number of perfect matchings).We thus have to change the strategy and systematically investigate multiphoton detection events by using the stratification according to the overall photon number and analyze all detection events, both pure and mixed.This is the way we will proceed in the upcoming examples.Here it turns out that the first differences between the two graphs happen for |n| = 6.Table 1 on page 22 summarizes the result.The leftmost column contains all partitions of |n| = 6 (see (47)) where each partition is represented by a naturally ordered orbit representative.The orbit size is in the second column in accordance with the discussion preceding Eq. ( 46).In experimental terms, an orbit consists of a measurement pattern and all of its permutations.The two rightmost columns contain the hafnians haf [A/ ⊗ |n| ].We notice a difference in three orbits (greyed): (1, 1, 1, 1, 2), (1, 1, 2, 2) and (1, 1, 1, 3) (the zeros omitted).Also note that the last four rows corresponding to events which do not occur as predicted by Lemma 6.For another graph G 3 , isomorphic to G 2 , we get the same hafnians for all orbits as an additional check.
Another interesting piece of information is the actual partition-averaged photon distribution for a given orbit provided by (51) or (52) (we omit the determinants in this example).For non-SRGs the partitionaveraged photon distribution is typically non-flat.To illustrate (51) let's choose an orbit where no difference was found: O n for n = (1, 2, 3).The first two plots in Fig. 3 are clearly different (that is, non-permutationally invariant).In accordance with the result of Section 5.4, namely Theorem 18, this is enough to decide that the two graphs are not isomorphic.For a graph G 3 isomorphic to G 2 we notice a mere permutation of bars in the rightmost panel of Fig. 3 again in accordance with Theorem 18.
A similar conclusion follows from the analysis of (52) and the situation is depicted in Fig. 4 for |n| = 8.
Example.Regular isospectral nonisomorphic graphs appear to be somewhere between PINGs and SRGs in terms of the difficulty to distinguish them.We analyzed a pair of graphs on ten vertices introduced in [34, page 110, (a) and (b)].Neither the coarse-grained photon distribution, Eq. ( 52), nor its probability 4. Coarse-grained partition-averaged photon distribution, Eq. 52, for G 1 , G 2 and G 3 G 2 for all orbits contributing to |n| = 8.The x axis labels the detectors.Note that since we omitted the determinantal prefactor and set c = 1 the distribution is not normalized.
Example (SRG(29,14,6,7)).Fig. 5 summarizes the path to distinguish all 41 isospectral graphs.As the starting point we took orbit O n for n = (1, 1, 1, 1).This orbit has no distinguishing power.This is indicated by a single square box containing all graphs.One could start with a measurement event containing more single photons but the problem is, as the number of vertices increases, the orbit size grows rapidly making the simulations rather resource-expensive.Also, it is desirable to find the orbit with the smallest possible |n| distinguishing all graphs to (i) heuristically assess the performance of our algorithm and (ii) make sure that the physical resources needed to run the algorithm are not excessive.This is because the smaller |n| is the less squeezing we need in an actual experiment.Sometimes, however, a smaller |n| does not guarantee a faster simulation.In the current example the orbit of n = (2, 2, 2, 2, 2, 2) where |n| = 12 is much more computationally feasible than n = (1, 2, 3, 4) where |n| = 10.This is because |O n | = 475020 for the former SRG (29,4,6,17) 1-41 Example (SRG (35,18,9,9)).The analysis of the biggest family of SRGs we studied is summarized in Fig. 6.Starting from the top, the orbit representatives on the left indicate how successful they are in distinguishing the graphs.The right box contains the number of newly distinguished graphs whereas the left box contains the number of remaining graphs.The effect is cumulative so for example the orbit of n = (4, 4, 4, 4, 4) together with n = (4, 4, 4, 4) distinguishes 2443 graphs.Our computational resources were not enough to distinguish the two remaining graphs using Theorem 12.The necessary condition developed in Sec.5.3 was used instead.
Example (SRG(16,6,2,2)).The two graphs can be easily distinguished by our method but in this case we illustrate the probability function for all partitions and their orbits up to |n| = 14.In Fig. 7 we plot Eq. 46 for orbits whose probability is nonzero (so their number is less than given by partitioning |n|).The x axis labels these orbits and in the plot we indicate the actual partitioning by white and gray background.Even for a fixed |n| some orbits are more likely than others.We observed that the probability is correlated to the size of the orbit.This confirms our intuition from the paragraph before Eq.(46).
Q1.The main result of this paper is a necessary and sufficient condition for two isospectral graphs to be isomorphic.We found a complete set of graph invariants.However, this is only half satisfactory because we don't know where the difference between two graphs 'kicks in'.Without this knowledge we can use the iff condition only in one direction.The ideal situation would be to have a deterministic or probabilistic criterion for the existence of such a threshold orbit as a polynomial function of the graph size.The numerical experiments are favorable as far as the polynomial growth goes -there is no indication that the threshold value grows fast.As the SRG example at the end of Section 5.4 shows, this is actually a more subtle problem: the set of hafnians may be different which is one sufficient condition as shown in Section 5.2 but their sum of squares is forming p G i (O n ) (the coarse-grained probability as another sufficient condition) is the same.The latter often comes 'later', that is, for higher orbits than the former.7. Probabilities of various orbits.Each point is a probability of an orbit O n given by Eq. ( 46) for one of the graphs from the SRG(16,6,2,2) family.The prefactor in ( 46) is 1/ det σ Q,G = ((−1 + 4c 2 ) 15 (−1 + 36c 2 )) 1/2 where we set c=1/6.9.Orbits whose probability is zero are omitted.again different.We call this a branching effect.It seems to be a probabilistic effect but from our experiments it holds overwhelmingly.This results in a dramatic increase of orbits with different hafnian sets as we increase |m| and greatly helps in the practical use of our algorithm.The intuition behind this behavior is that when the hafnian sets differ for an orbit of n, then another orbit of m > n (obtained by adding two new rows/columns or copies thereof to the adjacency matrix corresponding to n) contains as its subgraphs all the graphs that already had different hafnian sets.The question to answer is how likely this effect is to occur.

Q3.
A problem closely related to the previous question is how long it takes to approach a true probability distribution within a given precision for a chosen orbit of interest.This is typically answered by the methods that are standard in random graph and probability theory.
7.2.Scalability discussion.As discussed in the previous subsection, we don't know at what point two non-isomorphic graphs start to differ when sampled using a GBS device.Presumably this critical |n| grows with the graph size and the numerical evidence suggests that the growth is not exponential or even fast in general in the graph size.Nonetheless, even a moderate rate of growth of the threshold orbit n with the graph size M could be fatal for large graphs.This is because the physical interpretation of 0 < c < A 2 that directs the amount of squeezing -a quantity directly related to |n|.As we see from Fig. 7, the class with a (nearly) maximal probability (a single orbit in fact) is n = (0) M and this is a generic case.We can tune the c parameter by increasing squeezing such that the probability of a different (desired) orbit O n increases.This typically makes the vacuum contribution smaller but still dominating the probability landscape.What happens is that the desired orbit's probability p G A (O n ) increases with respect to the vacuum and other lower contributions but the probability flattens and therefore its values are inevitably impractically low to gather enough statistics reasonably fast.To put it differently, there is no significant peak (or a concentration effect) for the desired orbit O n .The situation is a bit alleviated by a heuristic observation that beyond the threshold orbit O n , where the difference in the hafnian set is first observed, the orbits O m where n < m overwhelmingly add to distinguishability.But it is tempting to avoid this 'probability dilution' altogether.
The rise of the number of orbits for a fixed |n| is equivalent to the integer partition problem briefly discussed towards the end of Section 5.2.Let us look closer at this problem.Any graph G A can be encoded as a pure covariance matrix whose circuit decomposition consists of an array of M squeezers S followed by an interferometer U on M modes.Since the eigenvalues (more precisely the singular values) of A are related to the squeezing parameters, two isospectral graphs, G A , G Ã, have the same S -the fact already used in Lemma 17.The input to the interferometer is given by (50).Suppose that we know or suspect that a given |n| contains orbits that are different if two graphs are not isomorphic.Since an interferometer is a passive unitary operator we know the input state responsible for it -it is the state whose coefficients are α in (r 1 , . . ., r M ) from (50).So the task becomes to prepare such a state.This could be a computationally hard task.Even though ) i=1 α in (r 1 , . . ., r M ) |(n, M )〉 i living in the completely symmetric subspaces are, in general, entangled.There are two problems, though.First, generation of such states could potentially require a circuit of a great depth.Second, even if n scales favorably with M , the number of coefficients α in (r 1 , . . ., r M ) can be overwhelming to work with for M 0 and so we may run into the issue of tractability to describe the necessary unitary operation.If these issues were resolved we would gained a complete control over the output distribution behavior.But we would also switch from GBS to a generalized (multiphoton) boson sampling (BS) [40].The orbit representative is a state with a given photon number per input mode.The difference compared to BS is that we do not require the input/output state to be in the 0,1 subspace per mode.Hence we arrived (by a detour) to the output probability function whose form is most likely governed by some permanent function [41] -a form which will most likely be a variation on the reduced Kronecker product we have introduced in Sec.5.1.
The use of a fixed input photon distribution will dramatically change the odds of detecting the difference in the probability distribution.We can take a look at Fig. 7 for, say, |n| = 14 and p G A (O n ) for all orbits in this partition will be considerably higher with their mutual ratio preserved.We leave for a future exploration the question if the probability distribution is always skewed (or even concentrated) such that some events (preferably the ones where there is a difference) are overwhelmingly likely than the others.What can we do if this is not the case?The thing we cannot do is to coarse-grain the probability more than in Eq. (46) due to lemma 17.Then, one option would be to coarse-grain the probability more than in (46) but less than in (48).The reason for this effort is to have a favorable scaling.Recall that the number of partitions of |n| grows exponentially with |n|.If we partially coarse-grain a given orbit we may keep the scaling polynomial and still detect a difference for non-isomorphic graphs.The question is how to split a given orbit.At this point we can offer only certain heuristic rules based on our simulations.

A -A
Algorithm 1 GBS graph isomorphism algorithm: returns the GI invariant of adjacency matrix A with respect to orbit o, considering hafnians' to the nth power.Generate the reduced Kronecker product of matrix A 7: Ap ← kron_reduced(A, p) 8: Append the hafnian to the nth power Calculate the sum of hafnians n i=1 X m i ], where (m 1 , . . ., m n ) ∈ n + and n i=1 m i is odd, are zero.A polynomial p(x ), x = (x 1 , . . ., x n ) ∈ n is called symmetric if for each permutation σ : [n] → [n] = (1, . . ., n) we have the equality p(x ) = p(x σ ), where x σ = (x σ(1) , . . ., x σ(n) ).Denote by S n the symmetric group of bijections σ : [n] → [n].Denote by P n ⊂ n×n the group of n × n permutation matrices.So P(σ)x = x σ .

Q2.
When a threshold orbit n is reached we made the following observation: upon examining orbits m for |m| > |n| satisfying m > n (m i ≥ n i , ∀i and at least one n i = m i ) the set of hafnians is

17 : 2
Calculate the mean photon distribution for the orbit16: phDist ← array[len(o)] for i ∈ len(perms) do 18: for j ∈ len(o) do 19: phDist[j] ← phDist[j] + perms[i, j]*result[i] Function to return the reduced Kronecker product of matrix A, given a sequence of integers n indicating the multi-mode photon detection event.
2and in total m 1 × m 2 matrices B = [B i j ] where each B i j is an n i × n j matrix.Then the reduced Kronecker product C = A/ ⊗B will denote a block matrix C partitioned as B and the blocks areC = [a i j B i j ].Remark.The dimension of A/ ⊗B is the dimension of the matrix B. Then A/ ⊗B is a submatrix of the Kronecker tensor product of two matrices A ⊗ B = [a i j B].The matrix B in this paper will always be assembled from B i j = n i ,n j ∈ n i ×n j where n is a measurement pattern.In this case we will write B = |n| ∈ |n|×|n| where ⊗ n 1 .Also note that if dim B i j = 1, ∀i, j the reduced Kronecker product becomes the Hadamard (Schur) product.