# Quantum-inspired permanent identities

Ulysse Chabaud  Institute for Quantum Information and Matter, California
Institute of Technology, Pasadena, CA 91125, USA<EMAIL_ADDRESS>Abhinav
Deshpande  Institute for Quantum Information and Matter, California Institute
of Technology, Pasadena, CA 91125, USA  Saeed Mehraban  Computer Science,
Tufts University, Medford, MA 02155, USA<EMAIL_ADDRESS>
###### Abstract

The permanent is pivotal to both complexity theory and combinatorics. In
quantum computing, the permanent appears in the expression of output
amplitudes of linear optical computations, such as in the Boson Sampling
model. Taking advantage of this connection, we give quantum-inspired proofs of
many existing as well as new remarkable permanent identities. Most notably, we
give a quantum-inspired proof of the MacMahon master theorem as well as proofs
for new generalizations of this theorem. Previous proofs of this theorem used
completely different ideas. Beyond their purely combinatorial applications,
our results demonstrate the classical hardness of exact and approximate
sampling of linear optical quantum computations with input cat states.

###### Contents

  1. 1 Introduction
  2. 2 Background
    1. 2.1 Notations and preliminary material
    2. 2.2 Boson Sampling
    3. 2.3 The MacMahon master theorem
  3. 3 Main results
  4. 4 Discussion
  5. 5 Quantum-inspired proofs of permanent identities
    1. 5.1 Glynn’s formula
    2. 5.2 The Glynn–Kan formula
    3. 5.3 The Cauchy–Binet theorem
    4. 5.4 Generating functions
  6. 6 New quantum-inspired permanent identities
    1. 6.1 Generalizations of the MacMahon master theorem
    2. 6.2 New generating functions
  7. 7 Boson Sampling with input cat states

##  1 Introduction

The permanent of an $m\times m$ matrix $A=(a_{ij})_{1\leq i,j\leq m}$ is a
combinatorial function defined as [1, 2]:

$\mathrm{Per}(A)=\sum_{\sigma\in\mathcal{S}_{m}}\prod_{k=1}^{m}a_{k\sigma(k)},$
(1)

where $\mathcal{S}_{m}$ is the symmetric group over $m$ symbols. While the
closely related determinant can be computed efficiently using many methods
such as Gaussian elimination, Valiant famously proved the $\\#$P-hardness of
computing the permanent [3].

Interestingly, the permanent appears in the expression of the output
amplitudes of linear optical quantum computations with noninteracting bosons
[4, 5], as in the Boson Sampling model of quantum computation [6]. This
connection has lead to several linear optical proofs of existing and new
classical complexity results: computation of the permanent is $\\#$P-hard [7],
(inverse polynomial) multiplicative estimation of the permanent of positive
semidefinite matrices is in $\textsf{BPP}^{\textsf{NP}}$ [8], multiplicative
estimation of the permanent of orthogonal matrices is $\\#$P-hard [9], and
computation of a class of multidimensional integrals is $\\#$P-hard [10]. It
has also lead to the introduction of a quantum-inspired classical algorithm
for additive estimation of the permanent of positive semidefinite matrices
[11]. Moreover, an approach inspired by quantum tomography recently showed
that (subexponential) multiplicative estimation of the permanent of positive
semidefinite matrices is NP-hard [12].

Beyond its importance for complexity theory, the permanent has numerous
applications for solving combinatorial problems [1, 2] and identities for the
permanent have been instrumental in these applications. For example, the
MacMahon master theorem [13], which relates the permanent to a coefficient of
the Taylor series of a determinant, is an invaluable tool for proving
combinatorial identities [14, 15, 16]. Similarly, Ryser’s formula [17] and
Glynn’s formula [18, 19, 20] are routinely used to compute the permanent more
efficiently than the naive brute-force approach.

While linear optics has been used as a tool to explore the classical
complexity of the permanent, previous work suggests that it can also be a
useful way to obtain simple proofs of theorems about the permanent: for
instance, [5] shows that the permanent of a unitary matrix $U$ lies in the
(closed) complex unit disk by expressing $|\mathrm{Per}(U)|^{2}$ as an output
probability of a linear optical sampling computation, while [6] derives simple
permanent identities using various representations of the same linear optical
sampling computation. This begs the following questions:

  1. 1.

Can we use linear optics to prove existing remarkable permanent identities,
such as the MacMahon master theorem?

  2. 2.

Can we use linear optics to derive new remarkable permanent identities?

In this work, we show that the answer to both questions is yes.

We give quantum-inspired proofs of several important permanent identities in
section 5. In particular, we show that the MacMahon master theorem can be
understood as two different ways of computing an inner product between two
Gaussian quantum states.

We also derive new quantum-inspired identities for the permanent. Our main
results are summarized in section 3 and proven in section 6. These include
generalizations of the MacMahon master theorem (Theorems 1 and 2), new
generating functions for the permanent (Theorem 3) and a formula for the sum
of two permanents (Theorem 4).

As a bonus, our findings also have consequences for the classical complexity
of exact and approximate sampling of linear optical sampling computations with
input cat states (Theorem 5), which we discuss in section 7. We rigorously
prove that the corresponding quantum probability distributions are as hard to
sample as the original Boson Sampling distribution [6], for all cat state
amplitudes in the exact case and for small enough amplitudes in the
approximate case. Until now, a formal proof was available only in the exact
case [21].

##  2 Background

###  2.1 Notations and preliminary material

$\bm{0}=(0,\dots,0)$ | $\bm{p}\leq\bm{q}\Leftrightarrow\forall k\in\\{1,\dots,m\\},\;p_{k}\leq q_{k}$
---|---
$\bm{1}=(1,\dots,1)$ | $-\bm{z}=(-z_{1},\dots,-z_{m})$
$|\bm{p}|=p_{1}+\dots+p_{m}$ | $\bm{z}^{*}=(z_{1}^{*},\dots,z_{m}^{*})$
$\bm{p}!=p_{1}!\dots p_{m}!$ | $\|\bm{z}\|^{2}=|z_{1}|^{2}+\dots+|z_{m}|^{2}$
$\ket{\bm{p}}=\ket{p_{1}}\otimes\dots\otimes\ket{p_{m}}$ | $\bm{z}^{\bm{p}}=z_{1}^{p_{1}}\dots z_{m}^{p_{m}}$
$\bm{p}+\bm{q}=(p_{1}+q_{1},\dots,p_{m}+q_{m})$ | $\partial_{\bm{z}}^{\bm{p}}=\partial_{z_{1}}^{p_{1}}\dots\partial_{z_{m}}^{p_{m}}$
$\bm{p}\oplus\bm{q}=(p_{1},\dots,p_{m},q_{1},\dots,q_{m})$ |  $d^{m}\bm{z}d^{m}\bm{z}^{*}=d\Re{z_{1}}d\Im{z_{1}}\dots d\Re{z_{m}}d\Im{z_{m}}$
Table 1: Multi-index notations used in this paper, for $m\in\mathbb{N}^{*}$,
$\bm{p}=(p_{1},\dots,p_{m})\in\mathbb{N}^{m}$,
$\bm{q}=(q_{1},\dots,q_{m})\in\mathbb{N}^{m}$ and
$\bm{z}=(z_{1},\dots,z_{m})\in\mathbb{C}^{m}$.

We use bold math for multi-index expressions (see Table 1 above). We denote by
$\mathbb{T}=\\{z\in\mathbb{C},|z|=1\\}$ the complex unit circle. For all
$m\in\mathbb{N}^{*}$, all $\bm{z}=(z_{1},\dots,z_{m})\in\mathbb{C}^{m}$ and
all $\bm{p}=(p_{1},\dots,p_{m})\in\mathbb{N}^{m}$, we use the notation
$[\bm{z}^{\bm{p}}]$ to denote the coefficient of
$\bm{z}^{\bm{p}}=z_{1}^{p_{1}}\dots z_{m}^{p_{m}}$ in an analytic expression.
For all $m\in\mathbb{N}^{*}$, all
$\bm{p}=(p_{1},\dots,p_{m})\in\mathbb{N}^{m}$ and all
$\bm{q}=(q_{1},\dots,q_{m})\in\mathbb{N}^{m}$, we denote by
$A_{\bm{p},\bm{q}}$ the matrix obtained from an $m\times m$ matrix $A$ by
first repeating its $i^{th}$ row $p_{i}$ times (deleting the row if $p_{i}=0$)
for all $i\in\\{1,\dots,m\\}$ and then repeating its $j^{th}$ column $q_{j}$
times (deleting the column if $q_{j}=0$) for all $j\in\\{1,\dots,m\\}$. By
convention, we set the permanent of non-square matrices to $0$ and
$\mathrm{Per}(A_{\bm{0},\bm{0}})=1$.

Let $m\in\mathbb{N}^{*}$ denote the number of modes. We denote (unnormalized)
quantum states using the Dirac ket notation. These are vectors in an infinite-
dimensional Hilbert space spanned by the orthonormal Fock basis

$\\{\ket{\bm{p}}=\ket{p_{1}}\otimes\dots\otimes\ket{p_{m}}\\}_{\bm{p}=(p_{1},\dots,p_{m})\in\mathbb{N}^{m}}.$
(2)

Hereafter, we label Fock states using $p$, $q$, $\bm{p}$, $\bm{q}$, and
$\bm{k}$, coherent states using $\alpha$, $\beta$, $\bm{\alpha}$, and
$\bm{\beta}$, and two-mode squeezed states using $\lambda$, $\mu$,
$\bm{\lambda}$, and $\bm{\mu}$. In particular, coherent states are defined as

$\ket{\alpha}=e^{-\frac{1}{2}|\alpha|^{2}}\sum_{p\geq
0}\frac{\alpha^{p}}{\sqrt{p!}}\ket{p},$ (3)

for all $\alpha\in\mathbb{C}$, and (unnormalized) two-mode squeezed states as

$\ket{\lambda}=\sum_{p\geq 0}\lambda^{p}\ket{pp},$ (4)

for all $\lambda\in\mathbb{C}$ with $|\lambda|<1$. Moreover, cat states are
defined as

$\ket{\mathrm{cat}_{\alpha}}=\frac{e^{\frac{1}{2}|\alpha|^{2}}}{2\sqrt{\sinh(|\alpha|^{2})}}(\ket{\alpha}-\ket{-\alpha}),$
(5)

for all $\alpha\in\mathbb{C}$.

We make use of the following inner products involving these states: for all
$p\in\mathbb{N}$, $\alpha,\beta\in\mathbb{C}$, and $\lambda\in\mathbb{C}$ with
$|\lambda|<1$,

$\innerproduct{p}{\alpha}=e^{-\frac{1}{2}|\alpha|^{2}}\frac{\alpha^{p}}{\sqrt{p!}}$
(6)
$\innerproduct{\alpha}{\beta}=e^{-\frac{1}{2}|\alpha|^{2}-\frac{1}{2}|\beta|^{2}+\alpha^{*}\beta},$
(7) $\innerproduct{pp}{\lambda}=\lambda^{p},$ (8)
$\bra{\lambda}(\ket{\alpha}\otimes\ket{\beta})=e^{-\frac{1}{2}|\alpha|^{2}-\frac{1}{2}|\beta|^{2}+\lambda^{*}\alpha\beta},$
(9)

which can be readily checked in the Fock basis. Moreover, for all
$\bm{p},\bm{q}\in\mathbb{N}^{m}$ [5, 6],

$\bra{\bm{p}}\hat{U}\ket{\bm{q}}=\frac{\mathrm{Per}(U_{\bm{p},\bm{q}})}{\sqrt{\bm{p}!\bm{q}!}},$
(10)

where $\hat{U}$ is a passive linear operation over $m$ modes whose action on
the creation operators of the modes is described by the unitary matrix
$U=(u_{ij})_{1\leq i,j\leq m}$ as

$\hat{U}\hat{a}^{\dagger}_{j}\hat{U}^{\dagger}=\sum_{i=1}^{m}u_{ij}\hat{a}^{\dagger}_{i},$
(11)

for all $j\in\\{1,\dots,m\\}$. As their name indicates, passive linear
operations do not change the total number of photons [22]: for all
$n\in\mathbb{N}$,

$\hat{U}\Pi_{n}=\Pi_{n}\hat{U},$ (12)

where $\Pi_{n}:=\sum_{|\bm{p}|=n}\ket{\bm{p}}\\!\bra{\bm{p}}$ is the $m$-mode
projector onto states with total photon number equal to $n$. Passive linear
operations map tensor products of coherent states to tensor products of
coherent states [22]: for all $\bm{\alpha}\in\mathbb{C}^{m}$,

$\hat{U}\ket{\bm{\alpha}}=\ket{U\bm{\alpha}}.$ (13)

Recall that coherent states form an overcomplete basis:

$\int_{\bm{\alpha}\in\mathbb{C}^{m}}\ket{\bm{\alpha}}\\!\bra{\bm{\alpha}}\frac{d^{m}\bm{\alpha}d^{m}\bm{\alpha}^{*}}{\pi^{m}}=\hat{\mathbb{I}},$
(14)

where $\hat{\mathbb{I}}$ is the identity operator over $m$ modes. We will also
make use of the following Gaussian inner product: for all
$\bm{\lambda},\bm{\mu}\in\mathbb{C}^{m}$ with $|\lambda_{k}|<1$ and
$|\mu_{k}|<1$ for all $k\in\\{1,\dots,m\\}$, and for all passive linear
operation $\hat{U}$ over $2m$ modes,

$\displaystyle\bra{\bm{\lambda}^{*}}\hat{U}\ket{\bm{\mu}}$
$\displaystyle=\int_{\bm{\beta}\in\mathbb{C}^{2m}}\bra{\bm{\lambda}^{*}}\ket{\bm{\beta}^{*}}\\!\bra{\bm{\beta}^{*}}\hat{U}\ket{\bm{\mu}}\frac{d^{2m}\bm{\beta}d^{2m}\bm{\beta}^{*}}{\pi^{2m}}$
(15)
$\displaystyle=\int_{\bm{\beta}\in\mathbb{C}^{2m}}\bra{\bm{\lambda}^{*}}\ket{\bm{\beta}^{*}}\left(\bra{\bm{\mu}}{U^{\dagger}\bm{\beta}^{*}}\rangle\right)^{*}\frac{d^{2m}\bm{\beta}d^{2m}\bm{\beta}^{*}}{\pi^{2m}}$
$\displaystyle=\int_{\bm{\beta}\in\mathbb{C}^{2m}}e^{\sum_{k=1}^{m}\lambda_{k}\beta_{k}^{*}\beta_{m+k}^{*}}e^{\sum_{k=1}^{m}\mu_{k}(U^{T}\bm{\beta})_{k}(U^{T}\bm{\beta})_{m+k}}e^{-\sum_{j=1}^{2m}\beta_{j}^{*}\beta_{j}}\frac{d^{2m}\bm{\beta}d^{2m}\bm{\beta}^{*}}{\pi^{2m}}$
$\displaystyle=\int_{\bm{\beta}\in\mathbb{C}^{2m}}\exp\left[-\frac{1}{2}\begin{pmatrix}\bm{\beta}\\\
\bm{\beta}^{*}\end{pmatrix}^{T}\\!\\!\\!V_{U}(\bm{\lambda},\bm{\mu})\begin{pmatrix}\bm{\beta}\\\
\bm{\beta}^{*}\end{pmatrix}\right]\frac{d^{2m}\bm{\beta}d^{2m}\bm{\beta}^{*}}{\pi^{2m}}$
$\displaystyle=\frac{1}{\sqrt{\mathrm{Det}(V_{U}(\bm{\lambda},\bm{\mu}))}},$

where we have used the overcompleteness of coherent states (14) over $2m$
modes in the first line, the action of passive linear operations on coherent
states (13) in the second line, the overlap between two-mode squeezed states
and coherent states (9) in the third line, and where we have introduced in the
fourth line the $(4m)\times(4m)$ symmetric matrix

$V_{U}(\bm{\lambda},\bm{\mu}):=\begin{pmatrix}-UV_{\bm{\mu}}U^{T}&I_{2m}\\\
I_{2m}&-V_{\bm{\lambda}}\end{pmatrix},$ (16)

where for all $\bm{w}\in\mathbb{C}^{n}$,

$V_{\bm{w}}:=\begin{pmatrix}0_{m}&\mathrm{Diag}(\bm{w})\\\
\mathrm{Diag}(\bm{w})&0_{m}\end{pmatrix}.$ (17)

Note that we associate mode $1$ with mode $m+1$, mode $2$ with mode $m+2$, and
so on, i.e. $\ket{\bm{\lambda}}=\bigotimes_{k=1}^{m}\ket{\lambda_{k}}$, where
$\ket{\lambda_{k}}$ is an unormalized two-mode squeezed state over modes $k$
and $m+k$.

Finally, we note that any $m\times m$ matrix $B$ with $\|B\|\leq 1$, where
$\|\cdot\|$ is the spectral norm, may be embedded as the top-left submatrix of
a $(2m)\times(2m)$ unitary matrix $U$ (see [6, Lemma 29] for an explicit
construction). We will use this fact multiple times throughout the paper to
extend identities proven for unitary matrices to the case of generic matrices.

###  2.2 Boson Sampling

Boson Sampling is a sub-universal model of quantum computation introduced by
Aaronson and Arkhipov (AA) in [6], which takes as input a Fock state
$\ket{\bm{1}\oplus\bm{0}}=\ket{1}^{\otimes n}\otimes\ket{0}^{\otimes m-n}$,
evolves it according to a passive linear operation $\hat{U}$ over $m$ modes
with unitary matrix $U$, and measures the photon number of each output mode.

With Eq. (10), the probability of detecting
$\bm{p}=(p_{1},\dots,p_{m})\in\mathbb{N}^{m}$ output photons is given by:

$\displaystyle P_{\text{BS}}(\bm{p}|n):=$ $\displaystyle|\langle p_{1}\dots
p_{m}|\hat{U}(\ket{1}^{\otimes n}\otimes\ket{0}^{\otimes m-n})|^{2}$ (18)
$\displaystyle=$
$\displaystyle\frac{\left|\mathrm{Per}(U_{\bm{p},\bm{1}\oplus\bm{0}})\right|^{2}}{\bm{p}!}.$

This model of quantum computation, while not believed to be universal, is
already capable of outperforming its classical counterparts: AA showed that
the output probability distribution $P_{\text{BS}}$ is hard to sample exactly
classically for $m\geq 2n$, or the polynomial hierarchy of complexity classes
collapses to its third level [6].

Moreover, under additional plausible conjectures, AA showed that this collapse
holds for $m=\Theta(n^{5}\log^{2}n)$ even if only an efficient classical
algorithm for approximate sampling exists (i.e. a classical algorithm which
samples efficiently from a probability distribution that has a small total
variation distance with the ideal Boson Sampling output probability
distribution $P_{\text{BS}}$). They further conjectured that this approximate
hardness should hold for $m=\Theta(n^{2})$.

The approximate hardness proof for Boson Sampling is based on a matrix-hiding
argument which requires a specific regime of parameters $m$ and $n$: given
$\delta>0$, the relation between $m$ and $n$ should be such that the
distribution $\mathcal{H}^{n,m}$ of $n\times n$ submatrices of $m\times m$
Haar-random unitary matrices multiplied by $\sqrt{m}$ is $O(\delta)$-close in
total variation distance to the distribution $\mathcal{G}^{n\times n}$ of
$n\times n$ matrices of i.i.d. Gaussians (see Theorem 35 in [6]). In
particular, AA showed that, for any $\delta>0$,
$\|\mathcal{H}^{n,m}-\mathcal{G}^{n\times n}\|_{\text{TV}}=O(\delta)$ when
$m\geq\frac{n^{5}}{\delta}\log^{2}\frac{n}{\delta}$.

This result was later refined in [23] with the bound
$\|\mathcal{H}^{n,m}-\mathcal{G}^{n\times
n}\|_{\text{TV}}\leq\frac{n^{3}}{2(m-n)}$, which gives a total variation
distance $O(\delta)$ for $m\geq\frac{n^{3}}{\delta}$.

Finally, in the case of orthogonal matrices rather than unitary matrices, it
was shown in [24] that (see the equation following Eq. (2.31) in [24], with
$n\rightarrow m$ and $p=q\rightarrow n$):

$\displaystyle D_{\text{KL}}(\mathcal{H}^{n,m},\mathcal{G}^{n\times n})$
$\displaystyle\leq\frac{n^{3}}{4(m-n)}-\frac{(m-n)-(n+1)}{2}\left[\frac{n^{3}+O(n^{2})}{2(m-n)^{2}}+\frac{n^{4}+O(n^{3})}{3(m-n)^{3}}\right]$
(19) $\displaystyle\leq
O\\!\left(\frac{n^{2}}{m}+\frac{n^{4}}{m^{2}}+\frac{n^{5}}{m^{3}}\right),$

where $D_{\text{KL}}$ is the Kullback–Leibler distance and where we obtain the
second line under the assumption $n=o(m)$. In particular, for $\delta>0$,
choosing $m\geq\frac{n^{2}}{\delta^{2}}$ implies
$D_{\text{KL}}(\mathcal{H}^{n,m},\mathcal{G}^{n\times n})\leq O(\delta^{2})$,
and thus $\|\mathcal{H}^{n,m}-\mathcal{G}^{n\times n}\|_{\text{TV}}\leq
O(\delta)$ by Pinsker’s inequality (see Eq. (1.3) in [24]). It was also argued
in [24] that a similar result should hold for unitary matrices.

Summarising these results, approximate hardness of Boson Sampling is proven
for passive linear operations described by orthogonal matrices in the regime
$m=\Theta(n^{2})$ [24] and for unitary matrices in the regime
$m=\Theta(n^{3})$ [23].

###  2.3 The MacMahon master theorem

The MacMahon master theorem is an important result in combinatorics which
relates the permanent to the determinant:

######  Theorem (MacMahon master theorem [13]).

Let $\bm{z}=(z_{1},\dots,z_{m})$ be formal variables. For any $m\times m$
matrix $A$,

$\sum_{\bm{p}\in\mathbb{N}^{m}}\frac{\bm{z}^{\bm{p}}}{\bm{p}!}\mathrm{Per}(A_{\bm{p},\bm{p}})=\frac{1}{\mathrm{Det}(I-ZA)},$
(20)

where $Z=\mathrm{Diag}(\bm{z})$.

This theorem is particularly useful to derive short proofs of combinatorial
identities: it expresses the permanent of an $m\times m$ matrix $A$ with rows
and columns repeated in the same way as the coefficient

$\mathrm{Per}(A_{\bm{p},\bm{p}})=\bm{p}![\bm{z}^{\bm{p}}]\left(\frac{1}{\mathrm{Det}(I-ZA)}\right),$
(21)

where $Z=\mathrm{Diag}(\bm{z})$ and $\bm{p}\in\mathbb{N}^{m}$, while the same
permanent may also be expressed as the coefficient

$\mathrm{Per}(A_{\bm{p},\bm{p}})=\bm{p}![\bm{z}^{\bm{p}}](A\bm{z})^{\bm{p}}.$
(22)

Picking a specific matrix $A$ and a pattern $\bm{p}$ and computing the above
expressions yields combinatorial identities, a famous example being the short
proof of Dixon’s identity [25, 14],
$\sum_{k=0}^{2n}(-1)^{k}\binom{2n}{k}^{3}=(-1)^{n}\binom{3n}{n,n,n}$, by
taking

$A=\begin{pmatrix}0&1&-1\\\ -1&0&1\\\ 1&-1&0\end{pmatrix},$ (23)

and $\bm{p}=(2n,2n,2n)$ for $n\in\mathbb{N}^{*}$.

Various generalizations of the MacMahon master theorem have been introduced
over the years [26, 27, 28, 29, 30]. In physics, this theorem plays an
important role in the quantum theory of angular momentum [31] and is also
oftentimes interpreted as an instance of the boson-fermion correspondence
[27].

##  3 Main results

In this section, we summarize our main findings, which we prove in section 6.
We obtain the following generalization of the MacMahon master theorem:

######  Theorem 1.

Let $\bm{x}=(x_{1},\dots,x_{m})$ and $\bm{y}=(y_{1},\dots,y_{m})$ be formal
variables. For all $m\times m$ matrices $A$ and $B$,

$\sum_{\bm{p},\bm{q}\in\mathbb{N}^{m}}\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}})\mathrm{Per}(B_{\bm{q},\bm{p}})=\frac{1}{\mathrm{Det}(I-XAYB)},$
(24)

where $X=\mathrm{Diag}(\bm{x})$ and $Y=\mathrm{Diag}(\bm{y})$.

We further show that this generalization extends to $N$ matrices:

######  Theorem 2.

Let $N\geq 2$. For all $k\in\\{1,\dots,N\\}$, let
$\bm{z}_{k}=(z_{k1},\dots,z_{km})$ be formal variables. For all $m\times m$
matrices $A^{(1)},\dots,A^{(N)}$,

$\displaystyle\sum_{\bm{p}_{1},\dots,\bm{p}_{N}\in\mathbb{N}^{m}}\prod_{k=1}^{N}\frac{\bm{z}_{k}^{\bm{p}_{k}}}{\bm{p}_{k}!}\mathrm{Per}(A^{(1)}_{\bm{p}_{1},\bm{p}_{2}})\mathrm{Per}(A^{(2)}_{\bm{p}_{2},\bm{p}_{3}})\dots\mathrm{Per}(A^{(N)}_{\bm{p}_{N},\bm{p}_{1}})$
(25)
$\displaystyle\quad\quad\quad\quad\quad\quad=\frac{1}{\mathrm{Det}(I-Z_{1}A^{(1)}\dots
Z_{N}A^{(N)})},$

where $Z_{k}=\mathrm{Diag}(\bm{z}_{k})$ for all $k\in\\{1,\dots,N\\}$.

As a corollary of Theorem 1, we obtain:

######  Corollary 1.

For any $m\times m$ matrix $A$ and all $\bm{p},\bm{q}\in\mathbb{N}^{m}$ with
$|\bm{p}|=|\bm{q}|=n\in\mathbb{N}$,

$\mathrm{Per}(A_{\bm{p},\bm{q}})=\frac{\bm{p}!\bm{q}!}{n!}[\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}]\left(\bm{x}^{T}A\bm{y}\right)^{n}.$
(26)

As a consequence of Corollary 1, we obtain the following family of generating
functions for the permanent:

######  Theorem 3.

Let $f(z)=\sum_{n=0}^{+\infty}f_{n}z^{n}$ be a series. Let
$\bm{x}=(x_{1},\dots,x_{m})$ and $\bm{y}=(y_{1},\dots,y_{m})$ be formal
variables. For any $m\times m$ matrix $A$,

$f(\bm{x}^{T}A\bm{y})=\sum_{\begin{subarray}{c}n\in\mathbb{N},\bm{p},\bm{q}\in\mathbb{N}^{m}\\\
|\bm{p}|=|\bm{q}|=n\end{subarray}}f_{n}n!\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}}),$
(27)

Equivalently,

$\partial^{\bm{p}}_{\bm{x}}\partial^{\bm{q}}_{\bm{y}}f(\bm{x}^{T}A\bm{y})\big{|}_{\bm{x}=\bm{y}=\bm{0}}=\mathrm{Per}(A_{\bm{p},\bm{q}})\,\partial^{n}_{z}f(z)\big{|}_{z=0},$
(28)

for all $n\in\mathbb{N}$ and all $\bm{p},\bm{q}\in\mathbb{N}^{m}$ such that
$|\bm{p}|=|\bm{q}|=n$. As a result, when $f_{n}\neq 0$,

$\mathrm{Per}(A_{\bm{p},\bm{q}})=\underset{\bm{x},\bm{y}\in\mathbb{T}^{m}}{\mathbb{E}}\left[\frac{\bm{p}!\bm{q}!}{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}\frac{f(\bm{x}^{T}A\bm{y})}{\partial^{n}_{z}f(z)\big{|}_{z=0}}\right],$
(29)

where the average is over random vectors with complex coefficients of modulus
$1$.

From this theorem we derive various permanent identities, the most notable one
being a formula for the sum of two permanents:

######  Theorem 4.

For all $m\times m$ matrices $A$ and $B$, all $n\in\mathbb{N}$ and all
$\bm{p},\bm{q}\in\mathbb{N}^{m}$ such that $|\bm{p}|=|\bm{q}|=n$,

$\displaystyle\mathrm{Per}(A_{\bm{p},\bm{q}})+\mathrm{Per}(B_{\bm{p},\bm{q}})$
(30)
$\displaystyle\quad=\sum_{k=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{k}}{\binom{n-1}{k}}\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{a}+\bm{b}+\bm{c}=\bm{p}\\\
\bm{a}^{\prime}+\bm{b}^{\prime}+\bm{c}^{\prime}=\bm{q}\\\
|\bm{a}|=|\bm{b}|=|\bm{a}^{\prime}|=|\bm{b}^{\prime}|=k\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\frac{\bm{p}!\bm{q}!}{\bm{a}!\bm{b}!\bm{c}!\bm{a}^{\prime}!\bm{b}^{\prime}!\bm{c}^{\prime}!}\mathrm{Per}(A_{\bm{a},\bm{a}^{\prime}})\mathrm{Per}(B_{\bm{b},\bm{b}^{\prime}})\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}}).$

In particular, when $\bm{p}=\bm{q}=\bm{1}$,

$\mathrm{Per}(A)+\mathrm{Per}(B)=\sum_{k=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{k}}{\binom{n-1}{k}}\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{a}+\bm{b}+\bm{c}=\bm{1}\\\
\bm{a}^{\prime}+\bm{b}^{\prime}+\bm{c}^{\prime}=\bm{1}\\\
|\bm{a}|=|\bm{b}|=|\bm{a}^{\prime}|=|\bm{b}^{\prime}|=k\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\mathrm{Per}(A_{\bm{a},\bm{a}^{\prime}})\mathrm{Per}(B_{\bm{b},\bm{b}^{\prime}})\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}}).$
(31)

Section 5 is primarily devoted to quantum-inspired proofs of existing
permanent identities, including a new proof of the MacMahon master theorem.
Along the way, we obtain the following inner product formula:

######  Lemma 1.

For all $\alpha\in\mathbb{C}$, all $\bm{p}\in\mathbb{N}^{m}$, all $n\leq m$,
and any passive linear operation $\hat{U}$ over $m$ modes with unitary matrix
$U$,

$\displaystyle\langle\bm{p}|\hat{U}(\ket{\mathrm{cat}_{\alpha}}^{\otimes
n}\otimes\ket{0}^{\otimes m-n})$
$\displaystyle=\frac{\alpha^{n}}{\sqrt{\sinh^{n}(|\alpha|^{2})}}\langle\bm{p}|\hat{U}(\ket{1}^{\otimes
n}\otimes\ket{0}^{\otimes m-n})$ (32)
$\displaystyle=\frac{\alpha^{n}}{\sqrt{\sinh^{n}(|\alpha|^{2})}}\frac{\mathrm{Per}(U_{\bm{p},\bm{1}\oplus\bm{0}})}{\sqrt{\bm{p}!}}.$

Lemma 1 has the following direct implications for the hardness of Boson
Sampling with input cat states [21], which we discuss in section 7:

######  Theorem 5.

Let $m=\text{poly }n\geq 2n$ and $\alpha\neq 0$. Boson Sampling with input
$\ket{\mathrm{cat}_{\alpha}}^{\otimes n}\otimes\ket{0}^{\otimes m-n}$ with
$\alpha\in\mathbb{C}$ is hard to sample exactly classically unless the
polynomial hierarchy of complexity classes collapses to its third level.

Moreover, assuming $|\alpha|=\mathcal{O}(n^{-1/4}\log^{1/4}m)$ Boson Sampling
with input cat states is as hard to sample approximately as Boson Sampling
with input single-photons, i.e. Boson Sampling with input cat states is hard
to sample approximately classically unless the polynomial hierarchy of
complexity classes collapses to its third level, in the same regime as Boson
Sampling, modulo the complexity conjectures introduced by AA in [6].

This result extends the arguments of [21]—which gave formal proof of classical
hardness in the exact sampling case—to the approximate sampling case.

##  4 Discussion

In the next sections, we introduce quantum-inspired proofs of permanent
identities. This approach allows us to give quantum-operational
interpretations of seminal results, such as the MacMahon master theorem [13].
In particular, we show that this theorem can be seen as two facets of the same
bosonic Gaussian amplitude.

This quantum-inspired approach also yields a breadth of new permanent
identities. We give some examples of combinatorial applications of these
identities in section 6 and we anticipate that they have many more. Beyond
these purely combinatorial applications, it would be interesting to
investigate whether our new identities may be used to obtain more efficient
classical algorithms for computing or estimating the permanent. In particular,
our Theorem 3 provides new estimators for the permanent which could be of
interest, by minimizing the variance of these estimators over the choice of
the analytic function $f$.

We use the formalism of linear optics with noninteracting bosons, but our
approach can be applied more generally to linear optics with other types of
particles. For instance, we expect our proof techniques to lead to determinant
identities in the fermionic case [32], immanant identities in the case of
partially distinguishable particles [33, 34], and additional permanent
identities in the case of generalized bosons [35]. Moreover, graphical
languages are currently being developed for linear optical quantum
computations [36, 37], which could lead to graphical proofs of remarkable
identities in combination with our approach.

Finally, our Theorem 5 gives solid complexity-theoretic foundations for the
hardness of Boson Sampling with input cat states. Generation of such states
has progressed tremendously in the recent years, thanks to circuit QED in
particular [38, 39, 40]. We hope that these foundations will motivate an
experimental demonstration of quantum speedup based on Schrödinger’s cat
states.

##  5 Quantum-inspired proofs of permanent identities

In this section, we derive quantum-inspired proofs of existing permanent
identities. Along the way, we obtain a generalization of Glynn’s formula [20]
for the permanent of matrices with repeated rows and columns in Eq. (43), as
well as a generalization of the Glynn–Kan formula [41] for the permanent of
matrices with repeated rows and columns in Eq. (54).

###  5.1 Glynn’s formula

Glynn’s formula for the permanent of an $m\times m$ matrix $A=(a_{ij})_{1\leq
i,j\leq m}$ is [20]:

$\mathrm{Per}(A)=\frac{1}{2^{m-1}}\sum_{\begin{subarray}{c}\bm{x}\in\\{-1,1\\}^{m}\\\
x_{1}=1\end{subarray}}x_{1}\dots
x_{m}\prod_{i=1}^{m}\left(\sum_{j=1}^{m}a_{ij}x_{j}\right).$ (33)

By symmetry, it is equivalent to the identity

$\mathrm{Per}(A)=\frac{1}{2^{m}}\sum_{\bm{x}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}\prod_{i=1}^{m}\left(\sum_{j=1}^{m}a_{ij}x_{j}\right).$ (34)

###### Proof of Glynn’s formula.

To prove this identity using quantum-mechanical tools, let us introduce the
unnormalized cat state
$\ket{\tilde{\mathrm{cat}}_{\alpha}}:=\frac{1}{2\alpha}(\ket{\alpha}-\ket{-\alpha})$,
for $\alpha\in\mathbb{C}$. Using the Fock basis expansion of coherent states
(6), we have $\lim_{\alpha\rightarrow
0}\ket{\tilde{\mathrm{cat}}_{\alpha}}=\ket{1}$ in trace distance. Hence, using
the Fock basis expansion of $\hat{U}$ (10),

$\mathrm{Per}(U)=\lim_{\alpha\rightarrow 0}\langle 1\dots
1|\hat{U}|\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}\rangle,$
(35)

where $\hat{U}$ is a passive linear operation over $m$ modes with unitary
matrix $U=(u_{ij})_{1\leq i,j\leq m}$. We compute the right hand side of this
equation:

$\displaystyle\langle 1\dots
1|\hat{U}|\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}\rangle$
$\displaystyle=\frac{1}{(2\alpha)^{m}}\sum_{\bm{x}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}\langle 1\dots 1|\hat{U}|x_{1}\alpha\dots x_{m}\alpha\rangle$
$\displaystyle=\frac{1}{(2\alpha)^{m}}\sum_{\bm{x}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}\langle 1\dots 1|(\alpha U\bm{x})_{1}\dots(\alpha U\bm{x})_{m}\rangle$
(36)
$\displaystyle=\frac{1}{(2\alpha)^{m}}\sum_{\bm{x}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}e^{-\frac{1}{2}|\alpha|^{2}\|U\bm{x}\|^{2}}\prod_{i=1}^{m}(\alpha
U\bm{x})_{i}$
$\displaystyle=\frac{e^{-\frac{m}{2}|\alpha|^{2}}}{2^{m}}\sum_{\bm{x}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}\prod_{i=1}^{m}\left(\sum_{j=1}^{m}u_{ij}x_{j}\right),$

where we used the action of $\hat{U}$ on coherent states (13) in the second
line, the Fock basis expansion of coherent states (6) in the third line, and
the fact that $U$ is unitary in the last line. Taking the limit when
$\alpha\rightarrow 0$ yields Glynn’s formula for unitary matrices. ∎

With the same calculations, we may obtain a more general version of Eq. (5.1):
for all $\alpha\in\mathbb{C}$ and $\bm{p}\in\mathbb{N}^{m}$,

$\displaystyle\langle p_{1}\dots
p_{m}|\hat{U}|\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}\rangle$
$\displaystyle=\frac{1}{(2\alpha)^{m}}\sum_{\bm{x}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}e^{-\frac{1}{2}|\alpha|^{2}\|U\bm{x}\|^{2}}\prod_{i=1}^{m}\frac{\left[(\alpha
U\bm{x})_{i}\right]^{p_{i}}}{\sqrt{p_{i}!}}$ (37)
$\displaystyle=\frac{\alpha^{|\bm{p}|-m}e^{-\frac{m}{2}|\alpha|^{2}}}{2^{m}\sqrt{\bm{p}!}}\sum_{\bm{x}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}\prod_{i=1}^{m}\left(\sum_{j=1}^{m}u_{ij}x_{j}\right)^{\\!p_{i}}\\!\\!\\!\\!\\!.$

Note that when $|\bm{p}|<m$ all products in the sum have at least one $x_{i}$
missing. Hence, by symmetry, $\langle p_{1}\dots
p_{m}|\hat{U}|\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}\rangle=0$
when $|\bm{p}|<m$. With the Fock basis expansion of $\hat{U}$ (10), taking the
limit when $\alpha\rightarrow 0$ yields a version of Glynn’s formula for
unitary matrices with repeated rows:

$\mathrm{Per}(U_{\bm{p},\bm{1}})=\frac{\delta_{|\bm{p}|,m}}{2^{m}}\sum_{\bm{x}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}\prod_{i=1}^{m}\left(\sum_{j=1}^{m}u_{ij}x_{j}\right)^{\\!p_{i}}\\!\\!\\!\\!\\!,$
(38)

where $\delta$ is the Kronecker symbol.

So far, all identities are derived for unitary matrices $U$. In order to
retrieve the same identity for a generic matrix $A=(a_{ij})_{1\leq i,j\leq n}$
of size $n$, we can embed $\frac{1}{\|A\|}A$ (or $A$ directly if $A=0$) as a
submatrix of a unitary matrix $U$ of size $m=2n$ [6, Lemma 29] and compute:

$\frac{\mathrm{Per}(A_{\bm{q},\bm{1}})}{\sqrt{\bm{q}!}}=\|A\|^{|\bm{q}|}\lim_{\alpha\rightarrow
0}(\bra{q_{1}\dots q_{n}}\otimes\bra{0}^{\otimes
n})\hat{U}(\ket{\tilde{\mathrm{cat}}_{\alpha}}^{\otimes
n}\otimes\ket{0}^{\otimes n}).$ (39)

To do so, we compute a slightly more general inner product: for
$\alpha\in\mathbb{C}$, $\bm{p}\in\mathbb{N}^{m}$ and $n\leq m$,

$\displaystyle\langle p_{1}\dots
p_{m}|\hat{U}(\ket{\tilde{\mathrm{cat}}_{\alpha}}^{\otimes
n}\otimes\ket{0}^{\otimes m-n})$ (40)
$\displaystyle\quad\quad=\frac{1}{(2\alpha)^{n}}\sum_{\bm{x}\in\\{-1,1\\}^{n}}x_{1}\dots
x_{n}\langle p_{1}\dots p_{m}|\hat{U}|x_{1}\alpha\dots x_{n}\alpha\,0\dots
0\rangle$
$\displaystyle\quad\quad=\frac{1}{(2\alpha)^{n}}\sum_{\bm{x}\in\\{-1,1\\}^{n}}x_{1}\dots
x_{n}e^{-\frac{1}{2}|\alpha|^{2}\|U(\bm{x}\oplus\bm{0})\|^{2}}\prod_{i=1}^{m}\frac{\left[(\alpha
U(\bm{x}\oplus 0))_{i}\right]^{p_{i}}}{\sqrt{p_{i}!}}$
$\displaystyle\quad\quad=\frac{\alpha^{|\bm{p}|-n}e^{-\frac{n}{2}|\alpha|^{2}}}{2^{n}\sqrt{\bm{p}!}}\sum_{\bm{x}\in\\{-1,1\\}^{n}}x_{1}\dots
x_{n}\prod_{i=1}^{m}\left(\sum_{j=1}^{n}u_{ij}x_{j}\right)^{\\!p_{i}}\\!\\!\\!\\!\\!.$

In particular, setting $m=2n$ and $\bm{p}=\bm{q}\oplus\bm{0}$ for
$\bm{q}\in\mathbb{N}^{n}$,

$(\bra{q_{1}\dots q_{n}}\otimes\bra{0}^{\otimes
n})\hat{U}(\ket{\tilde{\mathrm{cat}}_{\alpha}}^{\otimes
n}\otimes\ket{0}^{\otimes
n})=\frac{\alpha^{|\bm{q}|-n}e^{-\frac{n}{2}|\alpha|^{2}}}{2^{n}\sqrt{\bm{q}!}}\\!\\!\\!\\!\sum_{\bm{x}\in\\{-1,1\\}^{n}}\\!\\!\\!\\!x_{1}\dots
x_{n}\prod_{i=1}^{n}\left(\sum_{j=1}^{n}u_{ij}x_{j}\right)^{\\!q_{i}}\\!\\!\\!\\!\\!.$
(41)

Choosing $u_{ij}=\frac{1}{\|A\|}a_{ij}$ for $1\leq i,j\leq n$ and letting
$\alpha$ go to $0$ proves the claim: the left hand side converges to
$\frac{1}{\sqrt{\bm{q}!}}\mathrm{Per}(U_{\bm{q}\oplus\bm{0},\bm{1}\oplus\bm{0}})=\frac{1}{\|A\|^{|\bm{q}|}\sqrt{\bm{q}!}}\mathrm{Per}(A_{\bm{q},\bm{1}})$,
while the right hand side converges to
$\frac{\delta_{|\bm{q}|,n}}{\|A\|^{|\bm{q}|}2^{n}\sqrt{\bm{q}!}}\sum_{\bm{x}\in\\{-1,1\\}^{n}}x_{1}\dots
x_{n}\prod_{i=1}^{n}\left(\sum_{j=1}^{n}a_{ij}x_{j}\right)^{\\!q_{i}}$. Hence,

$\mathrm{Per}(A_{\bm{q},\bm{1}})=\frac{\delta_{|\bm{q}|,n}}{2^{n}}\sum_{\bm{x}\in\\{-1,1\\}^{n}}x_{1}\dots
x_{n}\prod_{i=1}^{n}\left(\sum_{j=1}^{n}a_{ij}x_{j}\right)^{\\!q_{i}}\\!\\!\\!\\!\\!.$
(42)

A similar generalization of Glynn’s formula for matrices with repeated rows
(or columns) based on roots of unity has previously appeared in [42].

Moreover, for any $m\times m$ matrix $A=(a_{i,j})_{1\leq i,j\leq m}$ and all
$\bm{p},\bm{q}\in\mathbb{N}^{m}$ with $|\bm{p}|=|\bm{q}|=n$, we have the
formula

$\mathrm{Per}(A_{\bm{p},\bm{q}})=\frac{\bm{q}!}{n^{m}}\sum_{\bm{x}\in\\{1,e^{\frac{2i\pi}{n}},\dots,e^{\frac{2i(n-1)\pi}{n}}\\}^{m}}\bm{x}^{-\bm{q}}(A\bm{x})^{\bm{p}},$
(43)

which is a version of Glynn’s formula for matrices with repeated rows and
columns. This formula, which has previously appeared in [43, 44, 45], can be
easily proven by expanding the product
$(A\bm{x})^{\bm{p}}=\prod_{i=1}^{m}\left(\sum_{j=1}^{m}a_{ij}x_{j}\right)^{p_{i}}$
and using properties of roots of unity. We give a quantum-inspired proof in
what follows.

For all $n\geq 1$, all $q\leq n$, and all $\alpha\in\mathbb{C}$, we define the
following unnormalised superposition of coherent states:

$\ket{\sigma^{q,n}_{\alpha}}:=\frac{\sqrt{q!}}{n\alpha^{q}}\sum_{k=0}^{n-1}e^{-\frac{2ikq\pi}{n}}\ket{e^{\frac{2ik\pi}{n}}\alpha}.$
(44)

Using the definition of coherent states (3), we obtain

$\displaystyle\ket{\sigma^{q,n}_{\alpha}}$
$\displaystyle=\frac{\sqrt{q!}e^{-\frac{1}{2}|\alpha|^{2}}}{n\alpha^{q}}\sum_{p\geq
0}\frac{\alpha^{p}}{\sqrt{p!}}\left(\sum_{k=0}^{n-1}e^{\frac{2ik(p-q)\pi}{n}}\right)\ket{p}$
(45) $\displaystyle=\sqrt{q!}e^{-\frac{1}{2}|\alpha|^{2}}\sum_{l\geq
0}\frac{\alpha^{ln}}{\sqrt{(ln+q)!}}\ket{ln+q}.$

where we used $\sum_{k=0}^{n-1}e^{\frac{2ik(p-q)\pi}{n}}=n$ if $p=q\mod n$ and
$0$ otherwise in the second line. Hence, $\lim_{\alpha\rightarrow
0}\ket{\sigma^{q,n}_{\alpha}}=\ket{q}$ in trace distance. Using the Fock basis
expansion of $\hat{U}$ (10), we thus have, for all
$\bm{p},\bm{q}\in\mathbb{N}^{m}$ with $|\bm{p}|=|\bm{q}|=n\geq 1$,

$\frac{\mathrm{Per}(U_{\bm{p},\bm{q}})}{\sqrt{\bm{p}!\bm{q}!}}=\lim_{\alpha\rightarrow
0}\langle p_{1}\dots
p_{m}|\hat{U}\bigotimes_{j=1}^{m}\ket{\sigma^{q_{j},n}_{\alpha}},$ (46)

where $\hat{U}$ is a passive linear operation over $m$ modes with unitary
matrix $U=(u_{ij})_{1\leq i,j\leq m}$. We compute the right hand side of this
equation:

$\displaystyle\langle p_{1}\dots
p_{m}|\hat{U}\bigotimes_{j=1}^{m}\ket{\sigma^{q_{j},n}_{\alpha}}$
$\displaystyle=\frac{\sqrt{\bm{q}!}}{n^{m}\alpha^{n}}\sum_{k_{1},\dots,k_{m}=0}^{n-1}\prod_{j=1}^{m}e^{-\frac{2ik_{j}q_{j}\pi}{n}}\langle
p_{1}\dots
p_{m}|\hat{U}\bigotimes_{j=1}^{m}\ket{e^{\frac{2ik_{j}\pi}{n}}\alpha}$ (47)
$\displaystyle=\frac{\sqrt{\bm{q}!}}{n^{m}\alpha^{n}}\sum_{\bm{x}\in\\{1,e^{\frac{2i\pi}{n}},\dots,e^{\frac{2i(n-1)\pi}{n}}\\}^{m}}\\!\\!\bm{x}^{-\bm{q}}\langle
p_{1}\dots p_{m}|\hat{U}\ket{\alpha x_{1}\dots\alpha x_{m}}$
$\displaystyle=\frac{\sqrt{\bm{q}!}e^{-\frac{m}{2}|\alpha|^{2}}}{\sqrt{\bm{p}!}n^{m}}\sum_{\bm{x}\in\\{1,e^{\frac{2i\pi}{n}},\dots,e^{\frac{2i(n-1)\pi}{n}}\\}^{m}}\\!\\!\bm{x}^{-\bm{q}}(U\bm{x})^{\bm{p}}.$

where we used the definition of $\ket{\sigma_{\alpha}^{q,n}}$ (44) in the
first line and where the derivation for the last line is identical to that of
Eqs. (5.1-37). Letting $\alpha$ go to $0$ proves Eq. (43) for unitary
matrices. Finally, the case of a generic nonzero $m\times m$ matrix $A$ is
obtained as in Eq. (40) by embedding $\frac{1}{\|A\|}A$ as a submatrix of a
unitary matrix $U$ of size $2m$ and computing

$\left(\langle p_{1}\dots p_{m}|\otimes\bra{0}^{\otimes
m}\right)\hat{U}\left(\bigotimes_{j=1}^{m}\ket{\sigma^{q_{j},n}_{\alpha}}\otimes\ket{0}^{\otimes
m}\right)\\!.$ (48)

###  5.2 The Glynn–Kan formula

A symmetrized version of Glynn’s formula for the permanent of an $m\times m$
matrix $A$ has been recently derived [41] under the name Glynn–Kan formula:

$\mathrm{Per}(A)=\frac{1}{4^{m}m!}\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}y_{1}\dots y_{m}(\bm{x}^{T}A\bm{y})^{m}.$ (49)

While the Glynn–Kan formula gives a slower classical algorithm than the Glynn
formula for computing the permanent, it is motivated by a connection with
quantum algorithms for estimating the permanent [41].

###### Proof of the Glynn–Kan formula.

To prove this identity using quantum-mechanical tools, we again use the
unnormalized cat state
$\ket{\tilde{\mathrm{cat}}_{\alpha}}:=\frac{1}{2\alpha}(\ket{\alpha}-\ket{-\alpha})$,
as in the previous section. Using the Fock basis expansion of coherent states
(6), we have $\lim_{\alpha\rightarrow
0}\ket{\tilde{\mathrm{cat}}_{\alpha}}=\ket{1}$ in trace distance. Hence, using
the Fock basis expansion of $\hat{U}$ (10),

$\mathrm{Per}(U)=\lim_{\alpha\rightarrow
0}\langle\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}|\hat{U}|\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}\rangle,$
(50)

where $\hat{U}$ is a passive linear operation over $m$ modes with unitary
matrix $U$. We compute the right hand side of this equation:

$\displaystyle\langle\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}|\hat{U}|\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}\rangle$
$\displaystyle\quad\quad=\frac{1}{4^{m}|\alpha|^{2m}}\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}y_{1}\dots y_{m}\langle x_{1}\alpha\dots
x_{m}\alpha|\hat{U}|y_{1}\alpha\dots y_{m}\alpha\rangle$
$\displaystyle\quad\quad=\frac{1}{4^{m}|\alpha|^{2m}}\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}y_{1}\dots y_{m}\langle x_{1}\alpha\dots x_{m}\alpha|(\alpha
U\bm{y})_{1}\dots(\alpha U\bm{y})_{m}\rangle$
$\displaystyle\quad\quad=\frac{1}{4^{m}|\alpha|^{2m}}\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}y_{1}\dots
y_{m}e^{-\frac{1}{2}|\alpha|^{2}\|\bm{x}\|^{2}}e^{-\frac{1}{2}|\alpha|^{2}\|U\bm{y}\|^{2}}\prod_{i=1}^{m}e^{x_{i}\alpha^{*}(\alpha
U\bm{y})_{i}}$ (51)
$\displaystyle\quad\quad=\frac{e^{-m|\alpha|^{2}}}{4^{m}|\alpha|^{2m}}\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}y_{1}\dots y_{m}e^{|\alpha|^{2}\bm{x}^{T}U\bm{y}}$
$\displaystyle\quad\quad=\frac{e^{-m|\alpha|^{2}}}{4^{m}}\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}y_{1}\dots
y_{m}\sum_{k=0}^{+\infty}\frac{|\alpha|^{2k-2m}}{k!}(\bm{x}^{T}U\bm{y})^{k},$

where we used the action of $\hat{U}$ on coherent states (13) in the third
line, the overlap between coherent states (7) in the fourth line, and the fact
that $U$ is unitary in the fifth line. For $k<m$, all products in the
expansion of $(\bm{x}^{T}U\bm{y})^{k}$ have at least one $x_{i}$ missing.
Hence, by symmetry, the terms for $k<m$ vanish and we obtain

$\langle\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}|\hat{U}|\tilde{\mathrm{cat}}_{\alpha}\dots\tilde{\mathrm{cat}}_{\alpha}\rangle=\frac{e^{-m|\alpha|^{2}}}{4^{m}}\\!\\!\\!\\!\\!\\!\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}\\!\\!\\!\\!x_{1}\dots
x_{m}y_{1}\dots
y_{m}\sum_{l=0}^{\infty}\frac{|\alpha|^{2l}}{(m+l)!}(\bm{x}^{T}U\bm{y})^{m+l}.$
(52)

Finally, taking the limit when $\alpha$ goes to $0$ yields

$\mathrm{Per}(U)=\frac{1}{4^{m}m!}\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}y_{1}\dots y_{m}(\bm{x}^{T}U\bm{y})^{m}.$ (53)

This proves the Glynn–Kan formula for a unitary matrix $U$. Once again, the
proof extends straightforwardly to any matrix $A$ of size $n$ by embedding
$\frac{1}{\|A\|}A$ as a submatrix of a unitary matrix $U$ of size $2n$ and
computing $\lim_{\alpha\rightarrow
0}(\bra{\tilde{\mathrm{cat}}_{\alpha}}^{\otimes n}\otimes\bra{0}^{\otimes
n})\hat{U}(\ket{\tilde{\mathrm{cat}}_{\alpha}}^{\otimes
n}\otimes\ket{0}^{\otimes n})$. ∎

Similar to Eq. (43), we obtain a version of the Glynn–Kan formula for matrices
with repeated rows and columns using roots of unity as

$\mathrm{Per}(A_{\bm{p},\bm{q}})=\frac{\bm{p}!\bm{q}!}{n^{2m}n!}\sum_{\bm{x},\bm{y}\in\\{1,e^{\frac{2i\pi}{n}},\dots,e^{\frac{2i(n-1)\pi}{n}}\\}^{m}}\bm{x}^{-\bm{p}}\bm{y}^{-\bm{q}}(\bm{x}^{T}A\bm{y})^{n},$
(54)

for any $m\times m$ matrix $A=(a_{ij})_{1\leq i,j\leq m}$ and all
$\bm{p},\bm{q}\in\mathbb{N}^{m}$ with $|\bm{p}|=|\bm{q}|=n$. This formula,
which appears to be new, can be easily proven using our Corollary 1 by
expanding the product $(\bm{x}^{T}A\bm{y})^{n}$ and using properties of roots
of unity. We give a quantum-inspired proof in what follows.

For all $n\geq 1$, all $q\leq n$, and all $\alpha\in\mathbb{C}$, we again use
the unnormalised superposition of coherent states as in Eq. (44):
$\ket{\sigma^{q,n}_{\alpha}}:=\frac{\sqrt{q!}}{n\alpha^{q}}\sum_{k=0}^{n-1}e^{-\frac{2ikq\pi}{n}}\ket{e^{\frac{2ik\pi}{n}}\alpha}$.
Recall that $\lim_{\alpha\rightarrow 0}\ket{\sigma^{q,n}_{\alpha}}=\ket{q}$ in
trace distance. Using the Fock basis expansion of $\hat{U}$ (10), we thus
have, for all $\bm{p},\bm{q}\in\mathbb{N}^{m}$ with $|\bm{p}|=|\bm{q}|=n\geq
1$,

$\frac{\mathrm{Per}(U_{\bm{p},\bm{q}})}{\sqrt{\bm{p}!\bm{q}!}}=\lim_{\alpha\rightarrow
0}\bigotimes_{j=1}^{m}\bra{\sigma^{p_{j},n}_{\alpha}}\hat{U}\bigotimes_{j=1}^{m}\ket{\sigma^{q_{j},n}_{\alpha}},$
(55)

where $\hat{U}$ is a passive linear operation over $m$ modes with unitary
matrix $U=(u_{ij})_{1\leq i,j\leq m}$. We compute the right hand side of this
equation:

$\displaystyle\bigotimes_{j=1}^{m}\bra{\sigma^{p_{j},n}_{\alpha}}\hat{U}\bigotimes_{j=1}^{m}\ket{\sigma^{q_{j},n}_{\alpha}}$
$\displaystyle=\frac{\sqrt{\bm{p}!\bm{q}!}}{n^{2m}|\alpha|^{2n}}\sum_{\begin{subarray}{c}k_{1},\dots,k_{m}=0\\\
l_{1},\dots
l_{m}=0\end{subarray}}^{n-1}\prod_{j=1}^{m}e^{-\frac{2ik_{j}p_{j}\pi}{n}}e^{-\frac{2il_{j}q_{j}\pi}{n}}\bigotimes_{j=1}^{m}\bra{e^{\frac{2ik_{j}\pi}{n}}\alpha}\hat{U}\bigotimes_{j=1}^{m}\ket{e^{\frac{2il_{j}\pi}{n}}\alpha}$
(56)
$\displaystyle=\frac{\sqrt{\bm{p}!\bm{q}!}}{n^{2m}|\alpha|^{2n}}\sum_{\bm{x},\bm{y}\in\\{1,e^{\frac{2i\pi}{n}},\dots,e^{\frac{2i(n-1)\pi}{n}}\\}^{m}}\\!\\!\bm{x}^{-\bm{p}}\bm{y}^{-\bm{q}}\langle
x_{1}\alpha\dots x_{m}\alpha|\hat{U}|y_{1}\alpha\dots y_{m}\alpha\rangle$
$\displaystyle=\frac{\sqrt{\bm{p}!\bm{q}!}e^{-m|\alpha|^{2}}}{n^{2m}}\sum_{\bm{x}\in\\{1,e^{\frac{2i\pi}{n}},\dots,e^{\frac{2i(n-1)\pi}{n}}\\}}\\!\\!\bm{x}^{-\bm{p}}\bm{y}^{-\bm{q}}\sum_{l=0}^{\infty}\frac{|\alpha|^{2l}}{(n+l)!}(\bm{x}^{T}U\bm{y})^{n+l}.$

where we used the definition of $\ket{\sigma_{\alpha}^{q,n}}$ (44) in the
first line and where the derivation for the last line is identical to that of
Eqs. (5.2-52). Letting $\alpha$ go to $0$ proves Eq. (54) for unitary
matrices. Once again, the proof extends straightforwardly to any nonzero
matrix $A$ of size $m$ by embedding $\frac{1}{\|A\|}A$ as a submatrix of a
unitary matrix $U$ of size $2m$ and computing

$\left(\bigotimes_{j=1}^{m}\bra{\sigma^{p_{j},n}_{\alpha}}\otimes\bra{0}^{\otimes
m}\right)\hat{U}\left(\bigotimes_{j=1}^{m}\ket{\sigma^{q_{j},n}_{\alpha}}\otimes\ket{0}^{\otimes
m}\right)\\!.$ (57)

###  5.3 The Cauchy–Binet theorem

The Cauchy–Binet theorem for the permanent expresses the permanent of the
product of two $m\times m$ matrices $A$ and $B$ as a sum of products involving
permanents of submatrices of these two matrices [1, 2]: for all
$\bm{p},\bm{q}\in\mathbb{N}^{m}$,

$\mathrm{Per}((AB)_{\bm{p},\bm{q}})=\sum_{\bm{k}\in\mathbb{N}^{m}}\frac{1}{\bm{k}!}\mathrm{Per}(A_{\bm{p},\bm{k}})\mathrm{Per}(B_{\bm{k},\bm{q}}).$
(58)

###### Proof of the Cauchy–Binet theorem.

Using the Fock basis expansion of passive linear operations (10) gives a quick
quantum-inspired proof of this identity: for $U$ and $V$ two $m\times m$
unitary matrices and for all $\bm{p},\bm{q}\in\mathbb{N}^{m}$,

$\displaystyle\mathrm{Per}((UV)_{\bm{p},\bm{q}})$
$\displaystyle=\sqrt{\bm{p}!\bm{q}!}\bra{\bm{p}}\widehat{UV}\ket{\bm{q}}$
$\displaystyle=\sqrt{\bm{p}!\bm{q}!}\bra{\bm{p}}\hat{U}\hat{V}\ket{\bm{q}}$
$\displaystyle=\sqrt{\bm{p}!\bm{q}!}\bra{\bm{p}}\hat{U}\left(\sum_{\bm{k}\in\mathbb{N}^{m}}\ket{\bm{k}}\\!\bra{\bm{k}}\right)\hat{V}\ket{\bm{q}}$
(59)
$\displaystyle=\sum_{\bm{k}\in\mathbb{N}^{m}}\sqrt{\bm{p}!\bm{q}!}\bra{\bm{p}}\hat{U}\ket{\bm{k}}\\!\bra{\bm{k}}\hat{V}\ket{\bm{q}}$
$\displaystyle=\sum_{\bm{k}\in\mathbb{N}^{m}}\frac{1}{\bm{k}!}\mathrm{Per}(U_{\bm{p},\bm{k}})\mathrm{Per}(V_{\bm{k},\bm{q}}),$

where we used the fact that Fock states form a basis in the third line and
where we used Eq. (10) once in the first line and twice in the last line.

In order to retrieve the formula for generic matrices $A$ and $B$, we can
embed $\frac{1}{\|A\|}A$ as a submatrix of a unitary matrix $U$ and
$\frac{1}{\|B\|}B$ as a submatrix of a unitary matrix $V$ and compute
$\mathrm{Per}((UV)_{\bm{p}\oplus\bm{0},\bm{q}\oplus\bm{0}})$. ∎

###  5.4 Generating functions

The permanent may be seen as a monomial coefficient in the Taylor expansion of
various functions [1, 2]. For example, the permanent of an $m\times m$ matrix
$A=(a_{ij})_{1\leq i,j\leq m}$ with rows repeated according to
$\bm{p}\in\mathbb{N}^{m}$ and columns repeated according to
$\bm{q}\in\mathbb{N}^{m}$ is given by the coefficient of
$\bm{z}^{\bm{q}}=z_{1}^{q_{1}}\dots z_{m}^{q_{m}}$ in

$\bm{q}!(A\bm{z})^{\bm{p}}=\prod_{i=1}^{m}q_{i}!\left(\sum_{j=1}^{m}a_{ij}z_{j}\right)^{\\!p_{i}}\\!\\!\\!\\!\\!.$
(60)

This may be thought of as the ‘monomial version’ of Glynn’s formula in Eq.
(42). Formally:

$\sum_{\bm{q}\in\mathbb{N}^{m}}\frac{\bm{z}^{\bm{q}}}{\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}})=(A\bm{z})^{\bm{p}},$
(61)

for formal variables $\bm{z}=(z_{1},\dots,z_{m})$.

###### Proof of ‘Glynn’s monomial formula’.

To prove this relation with quantum mechanical tools, we fix
$\bm{p}\in\mathbb{N}^{m}$, $\bm{\alpha}\in\mathbb{C}^{m}$, a unitary matrix
$U$ of size $m$, and we compute:

$\displaystyle\sum_{\bm{q}\in\mathbb{N}^{m}}\frac{\bm{\alpha}^{\bm{q}}}{\bm{q}!}\mathrm{Per}(U_{\bm{p},\bm{q}})$
$\displaystyle=\sqrt{\bm{p}!}e^{\frac{1}{2}\|\bm{\alpha}\|^{2}}\bra{\bm{p}}\hat{U}\left(\sum_{\bm{q}\in\mathbb{N}^{m}}\ket{\bm{q}}\\!\bra{\bm{q}}\right)\ket{\bm{\alpha}}$
$\displaystyle=\sqrt{\bm{p}!}e^{\frac{1}{2}\|\bm{\alpha}\|^{2}}\bra{\bm{p}}\hat{U}\ket{\bm{\alpha}}$
$\displaystyle=\sqrt{\bm{p}!}e^{\frac{1}{2}\|\bm{\alpha}\|^{2}}\bra{\bm{p}}(U\bm{\alpha})\rangle$
(62)
$\displaystyle=e^{\frac{1}{2}\|\bm{\alpha}\|^{2}-\frac{1}{2}\|U\bm{\alpha}\|^{2}}(U\bm{\alpha})^{\bm{p}}$
$\displaystyle=(U\bm{\alpha})^{\bm{p}},$

where we used the Fock basis expansions of coherent states (6) and of
$\hat{U}$ (10) in the first line, the fact that Fock states form a basis in
the second line, the action of $\hat{U}$ on coherent states (13) in the third
line, Eq. (6) again in the fourth line, and the fact that $U$ is unitary in
the last line. Once again, the relation for a generic nonzero matrix $A$ is
obtained by embedding $\frac{1}{\|A\|}A$ as a submatrix of a unitary matrix
$U$. ∎

Another generating function for the permanent is due to Jackson [46]: for
formal variables $\bm{x}=(x_{1},\dots,x_{m}),\bm{y}=(y_{1},\dots,y_{m})$,

$\sum_{\bm{p},\bm{q}\in\mathbb{N}^{m}}\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}})=e^{\bm{x}^{T}A\bm{y}}.$
(63)

Note that this relation implies Corollary 1, by expanding the Taylor series of
the exponential and considering the $\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}$
coefficient when $|\bm{p}|=|\bm{q}|=n\in\mathbb{N}^{*}$ (an alternative proof
of this result is given in the next section):

$\mathrm{Per}(A_{\bm{p},\bm{q}})=\frac{\bm{p}!\bm{q}!}{n!}[\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}](\bm{x}^{T}A\bm{y})^{n}.$
(64)

###### Proof of Jackson’s formula.

Jackson’s formula may be derived using quantum mechanical tools in a similar
way: for $\bm{\alpha},\bm{\beta}\in\mathbb{C}^{m}$ and a unitary matrix $U$,

$\displaystyle\sum_{\bm{p},\bm{q}\in\mathbb{N}^{m}}\frac{\bm{\alpha}^{\bm{p}}\bm{\beta}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(U_{\bm{p},\bm{q}})$
$\displaystyle=e^{\frac{1}{2}\|\bm{\alpha}\|^{2}+\frac{1}{2}\|\bm{\beta}\|^{2}}\bra{\bm{\alpha}^{*}}\left(\sum_{\bm{p}\in\mathbb{N}^{m}}\ket{\bm{p}}\\!\bra{\bm{p}}\right)\hat{U}\left(\sum_{\bm{q}\in\mathbb{N}^{m}}\ket{\bm{q}}\\!\bra{\bm{q}}\right)\ket{\bm{\beta}}$
(65)
$\displaystyle=e^{\frac{1}{2}\|\bm{\alpha}\|^{2}+\frac{1}{2}\|\bm{\beta}\|^{2}}\bra{\bm{\alpha}^{*}}\hat{U}\ket{\bm{\beta}}$
$\displaystyle=e^{\frac{1}{2}\|\bm{\alpha}\|^{2}+\frac{1}{2}\|\bm{\beta}\|^{2}}\bra{\bm{\alpha}^{*}}U\bm{\beta}\rangle$
$\displaystyle=e^{\frac{1}{2}\|\bm{\beta}\|^{2}-\frac{1}{2}\|U\bm{\beta}\|^{2}}e^{\bm{\alpha}^{T}U\bm{\beta}}$
$\displaystyle=e^{\bm{\alpha}^{T}U\bm{\beta}},$

where we used the Fock basis expansions of coherent states (6) and of
$\hat{U}$ (10) in the first line, the fact that Fock states form a basis in
the second line, the action of $\hat{U}$ on coherent states (13) in the third
line, the overlap between coherent states (7) in the fourth line, and the fact
that $U$ is unitary in the last line. Once again, the relation for a generic
nonzero matrix $A$ is obtained by embedding $\frac{1}{\|A\|}A$ as a submatrix
of a unitary matrix $U$. ∎

We conclude this section with arguably one of the most remarkable permanent
identities, the MacMahon master theorem [13], which relates the permanent and
the determinant through a generating function (see section 2.3):

$\sum_{\bm{p}\in\mathbb{N}^{m}}\frac{\bm{z}^{\bm{p}}}{\bm{p}!}\mathrm{Per}(A_{\bm{p},\bm{p}})=\frac{1}{\mathrm{Det}(I-ZA)},$
(66)

where $Z=\mathrm{Diag}(\bm{z})$, with $\bm{z}=(z_{1},\dots,z_{m})$ formal
variables.

###### Proof of the MacMahon master theorem.

To prove this relation with quantum mechanical tools, we show that the
MacMahon master theorem describes two different ways of computing an inner
product between two Gaussian states.

For $\lambda\in\mathbb{C}$ and $|\lambda|<1$, we make use of (unnormalized)
two-mode squeezed states of the form $\ket{\lambda}=\sum_{p\geq
0}\lambda^{p}\ket{pp}$. We write $\ket{\bm{\lambda}}$ with
$\bm{\lambda}=(\lambda_{1},\dots,\lambda_{m})\in\mathbb{C}^{m}$ a tensor
product of $m$ two-mode squeezed states (note that we are associating mode $1$
with mode $m+1$, mode $2$ with mode $m+2$, and so on). For $U$ a unitary
matrix and $\bm{\lambda},\bm{\mu}\in\mathbb{C}^{m}$, with $|\lambda_{k}|<1$
and $|\mu_{k}|<1$ for all $k\in\\{1,\dots,m\\}$, we compute, using the Fock
basis expansion of $\hat{U}$ (10):

$\sum_{\bm{p}\in\mathbb{N}^{m}}\frac{\bm{\lambda}^{\bm{p}}\bm{\mu}^{\bm{p}}}{\bm{p}!}\mathrm{Per}(U_{\bm{p},\bm{p}})=\sum_{\bm{p}\in\mathbb{N}^{m}}\bm{\lambda}^{\bm{p}}\bm{\mu}^{\bm{p}}\bra{\bm{p}}\hat{U}\ket{\bm{p}}.$
(67)

Using the Fock basis expansion of two-mode squeezed states (8), we have

$\bm{\lambda}^{\bm{p}}=\bra{\bm{\lambda}^{*}}\bm{p}\bm{p}\rangle.$ (68)

Moreover, a quick computation in Fock basis shows that

$\bm{\mu}^{\bm{p}}\bra{\bm{p}}\hat{U}\ket{\bm{p}}=\bra{\bm{p}\bm{p}}\hat{I}\otimes\hat{U}\ket{\bm{\mu}}.$
(69)

Hence, Eq. (67) rewrites

$\displaystyle\sum_{\bm{p}\in\mathbb{N}^{m}}\frac{\bm{\lambda}^{\bm{p}}\bm{\mu}^{\bm{p}}}{\bm{p}!}\mathrm{Per}(U_{\bm{p},\bm{p}})$
$\displaystyle=\sum_{\bm{p}\in\mathbb{N}^{m}}\bra{\bm{\lambda}^{*}}\bm{p}\bm{p}\rangle\\!\bra{\bm{p}\bm{p}}\hat{I}\otimes\hat{U}\ket{\bm{\mu}}$
(70)
$\displaystyle=\bra{\bm{\lambda}^{*}}\left(\sum_{\bm{p}\in\mathbb{N}^{m}}\ket{\bm{p}\bm{p}}\\!\bra{\bm{p}\bm{p}}\right)\hat{I}\otimes\hat{U}\ket{\bm{\mu}}$
$\displaystyle=\bra{\bm{\lambda}^{*}}\hat{I}\otimes\hat{U}\ket{\bm{\mu}},$

where in the last line we used
$\bra{\bm{\lambda}^{*}}\left(\sum_{\bm{p}\in\mathbb{N}^{m}}\ket{\bm{p}\bm{p}}\\!\bra{\bm{p}\bm{p}}\right)=\bra{\bm{\lambda}^{*}}$,
which can be checked in Fock basis.

We now compute the Gaussian inner product
$\bra{\bm{\lambda}^{*}}\hat{I}\otimes\hat{U}\ket{\bm{\mu}}$. From the action
of passive linear operations on creation operators (11), the unitary matrix
associated to the passive linear operation $\hat{I}\otimes\hat{U}$ is $I\oplus
U$. With the Gaussian overlap derived in Eqs. (15-17), we thus obtain:

$\bra{\bm{\lambda}^{*}}\hat{I}\otimes\hat{U}\ket{\bm{\mu}}=\frac{1}{\sqrt{\mathrm{Det}(V_{I\oplus
U}(\bm{\lambda},\bm{\mu}))}},$ (71)

where

$V_{I\oplus U}(\bm{\lambda},\bm{\mu})=\begin{pmatrix}-(I\oplus
U)V_{\bm{\mu}}(I\oplus U)^{T}&I_{2m}\\\
I_{2m}&-V_{\bm{\lambda}}\end{pmatrix},$ (72)

where for all $\bm{w}\in\mathbb{C}^{m}$,

$V_{\bm{w}}=\begin{pmatrix}0_{m}&\mathrm{Diag}(\bm{w})\\\
\mathrm{Diag}(\bm{w})&0_{m}\end{pmatrix}.$ (73)

We have

$\displaystyle\\!\mathrm{Det}(V_{I\oplus U}(\bm{\lambda},\bm{\mu}))$
$\displaystyle=\mathrm{Det}(I-(I\oplus U)V_{\bm{\mu}}(I\oplus
U)^{T}V_{\bm{\lambda}})$ $\displaystyle=\mathrm{Det}\\!\left(I-(I\oplus
U)\\!\begin{pmatrix}0_{m}&\mathrm{Diag}(\bm{\mu})\\\
\mathrm{Diag}(\bm{\mu})&0_{m}\end{pmatrix}\\!(I\oplus
U)^{T}\\!\\!\begin{pmatrix}0_{m}&\mathrm{Diag}(\bm{\lambda})\\\
\mathrm{Diag}(\bm{\lambda})&0_{m}\end{pmatrix}\\!\right)\\!\\!\\!$
$\displaystyle=\mathrm{Det}\\!\begin{pmatrix}I-\mathrm{Diag}(\bm{\mu})U^{T}\mathrm{Diag}(\bm{\lambda})&0_{m}\\\
0_{m}&I-U\mathrm{Diag}(\bm{\mu})\mathrm{Diag}(\bm{\lambda})\end{pmatrix}$ (74)
$\displaystyle=\mathrm{Det}(I-\mathrm{Diag}(\bm{\mu})U^{T}\mathrm{Diag}(\bm{\lambda}))\,\mathrm{Det}(I-U\mathrm{Diag}(\bm{\mu})\mathrm{Diag}(\bm{\lambda}))$
$\displaystyle=\mathrm{Det}(I-\mathrm{Diag}(\bm{\lambda})\mathrm{Diag}(\bm{\mu})U)^{2},$

where the last line is obtained by using Sylvester’s determinant theorem
$\mathrm{Det}(I+MN)=\mathrm{Det}(I+NM)$, together with the transpose for the
first determinant and
$\mathrm{Diag}(\bm{\mu})\mathrm{Diag}(\bm{\lambda})=\mathrm{Diag}(\bm{\lambda})\mathrm{Diag}(\bm{\mu})$
for the second determinant. Combining Eqs. (70), (71) and (5.4), we obtain,
$\bm{\lambda},\bm{\mu}\in\mathbb{C}^{m}$, with $|\lambda_{k}|<1$ and
$|\mu_{k}|<1$ for all $k\in\\{1,\dots,m\\}$,

$\sum_{\bm{p}\in\mathbb{N}^{m}}\frac{\bm{\lambda}^{\bm{p}}\bm{\mu}^{\bm{p}}}{\bm{p}!}\mathrm{Per}(U_{\bm{p},\bm{p}})=\bra{\bm{\lambda}^{*}}\hat{I}\otimes\hat{U}\ket{\bm{\mu}}=\frac{1}{\mathrm{Det}(I-\mathrm{Diag}(\bm{\lambda})\mathrm{Diag}(\bm{\mu})U)},$
(75)

where we used the value at $\bm{\lambda}=\bm{\mu}=0$ and the fact that the
left hand side is a continuous function of $\bm{\lambda}$ and $\bm{\mu}$ to
determine the sign of the square root. Replacing
$(\lambda_{1}\mu_{1},\dots,\lambda_{m}\mu_{m})$ by the formal variables
$\bm{z}=(z_{1},\dots,z_{m})$ concludes the proof. The relation for a generic
nonzero matrix $A$ of size $n$ is obtained by embedding $\frac{1}{\|A\|}A$ as
a submatrix of a unitary matrix $U$ of size $2n$ and taking
$\bm{z}=(z_{1},\dots,z_{n},0,\dots,0)$. ∎

##  6 New quantum-inspired permanent identities

In this section, we derive new quantum-inspired identities for the permanent
and we give some combinatorial applications of these identities.

###  6.1 Generalizations of the MacMahon master theorem

In this section, we introduce new generalizations of the MacMahon master
theorem (see section 2.3).

We first derive new quantum-inspired identities involving the permanent of
even-sized matrices: for $M$ a $(2m)\times(2m)$ matrix,

$\mathrm{Per}(M)=[z^{m}]\left(\frac{1}{4^{m}}\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}\frac{x_{1}\dots
x_{m}y_{1}\dots
y_{m}}{\sqrt{\mathrm{Det}(I-zV_{\bm{x}}MV_{\bm{y}}M^{T})}}\right),$ (76)

and, for formal variables
$\bm{x}=(x_{1},\dots,x_{m}),\bm{y}=(y_{1},\dots,y_{m})$,

$\sum_{\bm{p},\bm{q}\in\mathbb{N}^{m}}\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(M_{\bm{p}\oplus\bm{p},\bm{q}\oplus\bm{q}})=\frac{1}{\sqrt{\mathrm{Det}(I-V_{\bm{x}}MV_{\bm{y}}M^{T})}},$
(77)

where for all $\bm{w}=(w_{1},\dots,w_{m})$,

$V_{\bm{w}}=\begin{pmatrix}0_{m}&\mathrm{Diag}(\bm{w})\\\
\mathrm{Diag}(\bm{w})&0_{m}\end{pmatrix}.$ (78)

As a direct consequence, when $M=A\oplus B$, with $A$ and $B$ two $m\times m$
matrices, we obtain Theorem 1:

$\sum_{\bm{p},\bm{q}\in\mathbb{N}^{m}}\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}})\mathrm{Per}(B_{\bm{p},\bm{q}})=\frac{1}{\mathrm{Det}(I-XAYB^{T})},$
(79)

where $X=\mathrm{Diag}(x_{1},\dots,x_{m})$ and
$Y=\mathrm{Diag}(y_{1},\dots,y_{m})$. In particular,
$\mathrm{Per}(A)\mathrm{Per}(B)$ is given by the $x_{1}\dots x_{m}y_{1}\dots
y_{m}$ coefficient of $1/\mathrm{Det}(I-XAYB^{T})$. As a consequence, we
obtain Corollary 1: for all $n\in\mathbb{N}^{*}$ and all
$\bm{p},\bm{q}\in\mathbb{N}^{m}$ with $|\bm{p}|=|\bm{q}|=n$,

$\mathrm{Per}(A_{\bm{p},\bm{q}})=\frac{\bm{p}!\bm{q}!}{n!}[\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}]\left(\bm{x}^{T}A\bm{y}\right)^{n}.$
(80)

We also generalize Eq. (79) to the case of $N$ matrices to obtain Theorem 2:
for all $n\times n$ matrices $A^{(1)},\dots,A^{(N)}$,

$\displaystyle\sum_{\bm{p}_{1},\dots,\bm{p}_{N}\in\mathbb{N}^{m}}\prod_{k=1}^{N}\frac{\bm{z}_{k}^{\bm{p}_{k}}}{\bm{p}_{k}!}\mathrm{Per}(A^{(1)}_{\bm{p}_{1},\bm{p}_{2}})\mathrm{Per}(A^{(2)}_{\bm{p}_{2},\bm{p}_{3}})\dots\mathrm{Per}(A^{(N)}_{\bm{p}_{N},\bm{p}_{1}})$
(81)
$\displaystyle\quad\quad\quad\quad\quad\quad=\frac{1}{\mathrm{Det}(I-Z_{1}A^{(1)}\dots
Z_{N}A^{(N)})},$

where $\bm{z}_{k}=(z_{k1},\dots,z_{km})$ are formal variables and
$Z_{k}=\mathrm{Diag}(\bm{z}_{k})$ for all $k\in\\{1,\dots,N\\}$.

###### Proofs of Theorem 1, Corollary 1, and Theorem 2.

In section 5.1, we have obtained a linear optical proof of Glynn’s formula for
the permanent using the fact that a cat state $\ket{\mathrm{cat}_{\alpha}}$ of
small amplitude $\alpha$ approximates a single-photon Fock state $\ket{1}$.
Hereafter, we consider another approximation of Fock states using
superpositions of two-mode squeezed states. Such superpositions were recently
studied in [47] in the context of quantum sensing.

To prove the identity in Eq. (76), we introduce
$\ket{\lambda^{-}}:=\frac{1}{2\lambda}\left(\ket{\lambda}-\ket{-\lambda}\right)$,
for $|\lambda|<1$, where $\ket{\lambda}=\sum_{p\geq 0}\lambda^{p}\ket{pp}$ is
an unnormalized two-mode squeezed state. We have $\lim_{\lambda\rightarrow
0}\ket{\lambda^{-}}=\ket{11}$ in trace distance. Hence, for $U$ a
$(2m)\times(2m)$ unitary matrix we have with the Fock basis expansion of
$\hat{U}$ (10):

$\mathrm{Per}(U)=\lim_{\lambda\rightarrow 0}\prescript{\otimes
m}{}{\bra{\lambda^{-}}}\hat{U}\ket{\lambda^{-}}^{\otimes m}.$ (82)

For all $\lambda\in\mathbb{R}$, with $|\lambda|<1$, let us compute

$\prescript{\otimes m}{}{\bra{\lambda^{-}}}\hat{U}\ket{\lambda^{-}}^{\otimes
m}=\frac{1}{4^{m}\lambda^{2m}}\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}y_{1}\dots y_{m}g_{U}(\lambda\bm{x},\lambda\bm{y}),$ (83)

where we have defined for all $\bm{\lambda},\bm{\mu}\in\mathbb{C}^{m}$ (note
that we associate mode $1$ with mode $m+1$, mode $2$ with mode $m+2$, and so
on)

$g_{U}(\bm{\lambda},\bm{\mu}):=\bra{\bm{\lambda}^{*}}\hat{U}\ket{\bm{\mu}}.$
(84)

With the Gaussian overlap from Eq. (15), we have

$g_{U}(\bm{\lambda},\bm{\mu})=\frac{1}{\sqrt{\mathrm{Det}(V_{U}(\bm{\lambda},\bm{\mu}))}},$
(85)

where

$V_{U}(\bm{\lambda},\bm{\mu})=\begin{pmatrix}-UV_{\bm{\mu}}U^{T}&I_{2m}\\\
I_{2m}&-V_{\bm{\lambda}}\end{pmatrix},$ (86)

where for all $\bm{w}\in\mathbb{C}^{m}$,

$V_{\bm{w}}=\begin{pmatrix}0_{m}&\mathrm{Diag}(\bm{w})\\\
\mathrm{Diag}(\bm{w})&0_{m}\end{pmatrix}.$ (87)

Moreover,

$\mathrm{Det}(V_{U}(\bm{\lambda},\bm{\mu}))=\mathrm{Det}(I-V_{\bm{\lambda}}UV_{\bm{\mu}}U^{T}),$
(88)

so with Eq. (85)

$g_{U}(\bm{\lambda},\bm{\mu})=\frac{1}{\sqrt{\mathrm{Det}(I-V_{\bm{\lambda}}UV_{\bm{\mu}}U^{T})}}.$
(89)

Setting $\bm{\lambda}=\bm{\mu}=\lambda\bm{1}$ for $\lambda\in\mathbb{R}$ and
letting $\lambda$ go to $0$ in Eq. (83), we obtain that $\mathrm{Per}(U)$ is
given by

$\frac{1}{4^{m}}\\!\sum_{\bm{x},\bm{y}\in\\{-1,1\\}^{m}}x_{1}\dots
x_{m}y_{1}\dots y_{m}[\lambda^{2m}]g_{U}(\lambda\bm{x},\lambda\bm{y}),$ (90)

where $[\lambda^{2m}]g_{U}(\lambda\bm{x},\lambda\bm{y})$ is the $\lambda^{2m}$
coefficient in the Taylor expansion (in $\lambda$) of

$g_{U}(\lambda\bm{x},\lambda\bm{y})=\frac{1}{\sqrt{\mathrm{Det}(I-\lambda^{2}V_{\bm{x}}UV_{\bm{y}}U^{T})}}.$
(91)

This concludes the proof of Eq. (76). Once again, the relation for a generic
nonzero matrix $A$ is obtained by embedding $\frac{1}{\|A\|}A$ as a submatrix
of a unitary matrix $U$.

Now with Eq. (84), for all $\bm{\lambda},\bm{\mu}\in\mathbb{C}^{n}$,

$\displaystyle g_{U}(\bm{\lambda},\bm{\mu})$
$\displaystyle=\sum_{\bm{p},\bm{q}\in\mathbb{N}^{m}}\bm{\lambda}^{\bm{p}}\bm{\mu}^{\bm{q}}\bra{\bm{p}\bm{p}}\hat{U}\ket{\bm{q}\bm{q}}$
(92)
$\displaystyle=\sum_{\bm{p},\bm{q}\in\mathbb{N}^{m}}\frac{\bm{\lambda}^{\bm{p}}\bm{\mu}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(U_{\bm{p}\oplus\bm{p},\bm{q}\oplus\bm{q}}),$

where we used the Fock basis expansion of $\hat{U}$(10) in the second line.
Combining Eqs. (84), (89) and (92) we obtain

$\sum_{\bm{p},\bm{q}\in\mathbb{N}^{m}}\frac{\bm{\lambda}^{\bm{p}}\bm{\mu}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(U_{\bm{p}\oplus\bm{p},\bm{q}\oplus\bm{q}})=\bra{\bm{\lambda}^{*}}\hat{U}\ket{\bm{\mu}}=\frac{1}{\sqrt{\mathrm{Det}(I-V_{\bm{\lambda}}UV_{\bm{\mu}}U^{T})}},$
(93)

where for all $\bm{w}\in\mathbb{C}^{n}$,

$V_{\bm{w}}=\begin{pmatrix}0_{m}&\mathrm{Diag}(\bm{w})\\\
\mathrm{Diag}(\bm{w})&0_{m}\end{pmatrix}.$ (94)

This concludes the proof of Eq. (77). Once again, the relation for a generic
nonzero matrix $A$ is obtained by embedding $\frac{1}{\|A\|}A$ as a submatrix
of a unitary matrix $U$.

When $U=A\oplus B$, with $A$ and $B$ two $m\times m$ matrices, we have
$U_{\bm{p}\oplus\bm{p},\bm{q}\oplus\bm{q}}=A_{\bm{p},\bm{q}}\oplus
B_{\bm{p},\bm{q}}$ and the permanent of a block-diagonal matrix is the product
of the permanents of the blocks, so Eq. (93) gives

$\sum_{\bm{p},\bm{q}\in\mathbb{N}^{m}}\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}})\mathrm{Per}(B_{\bm{p},\bm{q}})=\frac{1}{\sqrt{\mathrm{Det}(I-V_{\bm{x}}(A\oplus
B)V_{\bm{y}}(A^{T}\oplus B^{T}))}}.$ (95)

Writing $X=\mathrm{Diag}(x_{1},\dots,x_{m})$ and
$Y=\mathrm{Diag}(y_{1},\dots,y_{m})$, we have

$\displaystyle\\!\\!\\!\\!\\!\\!\mathrm{Det}(I-V_{\bm{x}}(A\oplus
B)V_{\bm{y}}(A^{T}\oplus B^{T}))$
$\displaystyle=\mathrm{Det}\\!\left(I-\begin{pmatrix}0_{m}&X\\\
X&0_{m}\end{pmatrix}\\!\begin{pmatrix}A&0_{m}\\\
0_{m}&B\end{pmatrix}\\!\begin{pmatrix}0_{m}&Y\\\
Y&0_{m}\end{pmatrix}\\!\begin{pmatrix}A^{T}&0_{m}\\\
0_{m}&B^{T}\end{pmatrix}\right)\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(96) $\displaystyle=\mathrm{Det}\\!\begin{pmatrix}I-XBYA^{T}&0_{m}\\\
0_{m}&I-XAYB^{T}\end{pmatrix}$
$\displaystyle=\mathrm{Det}(I-XBYA^{T})\,\mathrm{Det}(I-XAYB^{T})$
$\displaystyle=\mathrm{Det}(I-XAYB^{T})^{2},$

where the last line is obtained by using Sylvester’s determinant theorem
$\mathrm{Det}(I+MN)=\mathrm{Det}(I+NM)$ together with the transpose for the
first determinant. Combining Eqs. (95) and (96) concludes the proof of Theorem
1.

This theorem reduces to the MacMahon master theorem when $A=I$ or $B=I$, since
$\mathrm{Per}(I_{\bm{p},\bm{q}})=\bm{p}!\delta_{\bm{p},\bm{q}}$. Another
particular case of interest is when $B=J_{m}$, where $J_{m}$ is the all-$1$
matrix of size $m\times m$, which satisfies $\mathrm{Per}(J_{m})=m!$. In this
case, we obtain, for all $\bm{x},\bm{y}\in\mathbb{C}^{m}$,

$\displaystyle\sum_{\begin{subarray}{c}n\in\mathbb{N},\bm{p},\bm{q}\in\mathbb{N}^{m}\\\
|\bm{p}|=|\bm{q}|=n\end{subarray}}\frac{n!\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}})$
$\displaystyle=\frac{1}{\mathrm{Det}(I-AYJ_{m}X)}$ (97)
$\displaystyle=\frac{1}{\mathrm{Det}(I-A\bm{y}\bm{x}^{T})}$
$\displaystyle=\frac{1}{1-\bm{x}^{T}A\bm{y}},$

where we used the fact that $A\bm{y}\bm{x}^{T}$ is a rank-one matrix to
compute the determinant. This implies Corollary 1: for all
$n\in\mathbb{N}^{*}$ and all $\bm{p},\bm{q}\in\mathbb{N}^{m}$ with
$|\bm{p}|=|\bm{q}|=n$,

$\displaystyle\mathrm{Per}(A_{\bm{p},\bm{q}})$
$\displaystyle=\frac{\bm{p}!\bm{q}!}{n!}[\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}]\left(\frac{1}{1-\bm{x}^{T}A\bm{y}}\right)$
(98)
$\displaystyle=\frac{\bm{p}!\bm{q}!}{n!}[\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}]\left(\bm{x}^{T}A\bm{y}\right)^{n},$

where we used the Taylor series $\frac{1}{1-z}=\sum_{k=0}^{+\infty}z^{k}$.

We note that Theorem 1—which is a generalization of the MacMahon master
theorem to two matrices—may also be obtained by combining the MacMahon master
theorem in Eq. (66) with the Cauchy–Binet theorem in Eq. (58), applied to the
matrix $AYB^{T}$. As it turns out, we may apply the same proof technique
inductively in order to generalize the MacMahon master theorem to $N$ matrices
$A^{(1)},\dots,A^{(N)}$, for $N\geq 1$ and obtain Theorem 2: assuming that Eq.
(81) holds for some $N\geq 2$, we have, for all $m\times m$ matrices
$B^{(1)},\dots,B^{(N)}$,

$\displaystyle\sum_{\bm{p}_{1},\dots,\bm{p}_{N}\in\mathbb{N}^{m}}\prod_{k=1}^{N}\frac{\bm{z}_{k}^{\bm{p}_{k}}}{\bm{p}_{k}!}\mathrm{Per}(B^{(1)}_{\bm{p}_{1},\bm{p}_{2}})\mathrm{Per}(B^{(2)}_{\bm{p}_{2},\bm{p}_{3}})\dots\mathrm{Per}(B^{(N)}_{\bm{p}_{N},\bm{p}_{1}})$
(99)
$\displaystyle\quad\quad\quad\quad\quad\quad=\frac{1}{\mathrm{Det}(I-Z_{1}B^{(1)}\dots
Z_{N}B^{(N)})},$

where $\bm{z}_{k}=(z_{k1},\dots,z_{km})$ are formal variables and
$Z_{k}=\mathrm{Diag}(\bm{z}_{k})$ for all $k\in\\{1,\dots,N\\}$.

Let $A^{(1)},\dots,A^{(N+1)}$ be $m\times m$ matrices,
$\bm{z}_{N+1}=(z_{N+1,1},\dots,z_{N+1,m})$ be formal variables and
$Z_{N+1}=\mathrm{Diag}(\bm{z}_{N+1})$. Setting
$B^{(N)}=A^{(N)}Z_{N+1}A^{(N+1)}$ and $B^{(k)}=A^{(k)}$ for all
$k\in\\{1,\dots,N\\}$, we obtain with Eq. (99):

$\displaystyle\frac{1}{\mathrm{Det}(I-Z_{1}A^{(1)}\dots
Z_{N}A^{(N)}Z_{N+1}A^{(N+1)})}$ (100)
$\displaystyle\quad=\sum_{\bm{p}_{1},\dots,\bm{p}_{N}\in\mathbb{N}^{m}}\prod_{k=1}^{N}\frac{\bm{z}_{k}^{\bm{p}_{k}}}{\bm{p}_{k}!}\mathrm{Per}(A^{(1)}_{\bm{p}_{1},\bm{p}_{2}})\dots\mathrm{Per}(A^{(N-1)}_{\bm{p}_{N-1},\bm{p}_{N}})\mathrm{Per}((A^{(N)}Z_{N+1}A^{(N+1)})_{\bm{p}_{N},\bm{p}_{1}}).$

On the other hand, with the Cauchy–Binet theorem from Eq. (58),

$\displaystyle\mathrm{Per}((A^{(N)}Z_{N+1}A^{(N+1)})_{\bm{p}_{N},\bm{p}_{1}})$
$\displaystyle=\\!\\!\\!\\!\\!\sum_{\bm{p}_{N+1}\in\mathbb{N}^{m}}\\!\frac{1}{\bm{p}_{N+1}!}\mathrm{Per}(A^{(N)}_{\bm{p}_{N},\bm{p}_{N+1}})\mathrm{Per}((Z_{N+1}A^{(N+1)})_{\bm{p}_{N+1},\bm{p}_{1}})$
(101)
$\displaystyle=\\!\\!\\!\\!\\!\sum_{\bm{p}_{N+1}\in\mathbb{N}^{m}}\\!\frac{\bm{z}_{N+1}^{\bm{p}_{N+1}}}{\bm{p}_{N+1}!}\mathrm{Per}(A^{(N)}_{\bm{p}_{N},\bm{p}_{N+1}})\mathrm{Per}(A^{(N+1)}_{\bm{p}_{N+1},\bm{p}_{1}}),$

where we used the fact that multiplying any single row of $M$ by a variable
$z$ changes $\mathrm{Per}(M)$ to $z\mathrm{Per}(M)$. Combining Eqs. (100) and
(101) completes the induction step and the proof of Theorem 2. ∎

Applications. Recall that the MacMahon master theorem is a useful tool for
deriving combinatorial identities (see section 2.3): it expresses the
permanent of an $m\times m$ matrix $A$ with rows and columns repeated in the
same way as the coefficient

$\mathrm{Per}(A_{\bm{p},\bm{p}})=\bm{p}![\bm{z}^{\bm{p}}]\left(\frac{1}{\mathrm{Det}(I-ZA)}\right),$
(102)

where $Z=\mathrm{Diag}(\bm{z})$ and $\bm{p}\in\mathbb{N}^{m}$, while the same
permanent may also be expressed as the coefficient

$\mathrm{Per}(A_{\bm{p},\bm{p}})=\bm{p}![\bm{z}^{\bm{p}}](A\bm{z})^{\bm{p}}.$
(103)

Our generalizations of the MacMahon master theorem may be used in a similar
fashion to obtain simple proofs of combinatorial identities. Let us illustrate
this with an example: for $n\in\mathbb{N}$ and $a,b\in\mathbb{C}$, setting

$S_{n}(a,b):=\sum_{k=0}^{n}\binom{n}{k}^{2}a^{k}b^{n-k},$ (104)

we aim to prove the relation

$S_{n}(a,b)^{2}=\sum_{l=0}^{n}\binom{2l}{l}\binom{n+l}{2l}(-1)^{n-l}(a-b)^{2n-2l}S_{l}(a^{2},b^{2}),$
(105)

which for $a=b=1$ directly implies the well-known
$\sum_{k=0}^{n}\binom{n}{k}^{2}=\binom{2n}{n}$.

###### Proof.

By Theorem 1 we have

$\mathrm{Per}(A_{\bm{p},\bm{q}})\mathrm{Per}(B_{\bm{p},\bm{q}})=\bm{p}!\bm{q}![\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}]\left(\frac{1}{\mathrm{Det}(I-XAYB^{T})}\right),$
(106)

where $X=\mathrm{Diag}(x_{1},\dots,x_{m})$ and
$Y=\mathrm{Diag}(y_{1},\dots,y_{m})$. Setting

$A=B=\begin{pmatrix}1&a\\\ 1&b\end{pmatrix},$ (107)

for $a,b\in\mathbb{C}$, $n\in\mathbb{N}^{*}$, and $\bm{p}=\bm{q}=(n,n)$, we
obtain with Eq. (103):

$\displaystyle\mathrm{Per}(A_{\bm{p},\bm{q}})$
$\displaystyle=(n!)^{2}[z_{1}^{n}z_{2}^{n}](z_{1}+az_{2})^{n}(z_{1}+bz_{2})^{n}$
(108) $\displaystyle=(n!)^{2}\sum_{k+l=n}\binom{n}{k}\binom{n}{l}a^{k}b^{l}$
$\displaystyle=(n!)^{2}S_{n}(a,b).$

On the other hand, with Eq. (106) we have

$\displaystyle\mathrm{Per}(A_{\bm{p},\bm{q}})^{2}$
$\displaystyle=(n!)^{4}[x_{1}^{n}x_{2}^{n}y_{1}^{n}y_{2}^{n}]\left(\frac{1}{\mathrm{Det}(I-\mathrm{Diag}(x_{1},x_{2})A\mathrm{Diag}(y_{1},y_{2})A^{T}}\right)$
$\displaystyle=(n!)^{4}[x_{1}^{n}x_{2}^{n}y_{1}^{n}y_{2}^{n}]\left(\frac{1}{1-x_{1}y_{1}-a^{2}x_{1}y_{2}-x_{2}y_{1}-b^{2}x_{2}y_{2}+(a-b)^{2}x_{1}x_{2}y_{1}y_{2}}\right)$
$\displaystyle=(n!)^{4}\sum_{k=n}^{2n}[x_{1}^{n}x_{2}^{n}y_{1}^{n}y_{2}^{n}]\left(x_{1}y_{1}+a^{2}x_{1}y_{2}+x_{2}y_{1}+b^{2}x_{2}y_{2}-(a-b)^{2}x_{1}x_{2}y_{1}y_{2}\right)^{k}$
(109)
$\displaystyle=(n!)^{4}\sum_{k=n}^{2n}\sum_{j=0}^{k}\binom{k}{j}(-1)^{k-j}(a-b)^{2k-2j}\\!\\!\\!\\!\\!\\!\\!\\!\sum_{n_{11}+n_{12}+n_{21}+n_{22}=j}\frac{j!a^{2n_{12}}b^{2n_{22}}}{n_{11}!n_{12}!n_{21}!n_{22}!}\\!\\!\\!\\!\\!\\!$
$\displaystyle\quad\quad\times[x_{1}^{n}x_{2}^{n}y_{1}^{n}y_{2}^{n}]\left((x_{1}x_{2}y_{1}y_{2})^{k-j}x_{1}^{n_{11}+n_{12}}x_{2}^{n_{21}+n_{22}}y_{1}^{n_{11}+n_{21}}y_{2}^{n_{12}+n_{22}}\right),\\!\\!\\!\\!\\!\\!$

where we used the Taylor expansion of $\frac{1}{1-z}$ in the third line and
the multinomial theorem in the last line. The indices in the above expression
must satisfy $n_{11}=n_{22}$, $n_{12}=n_{21}$, which implies
$j=2n_{11}+2n_{12}$, i.e. $j$ is even. Moreover, $n=k-j+n_{11}+n_{12}=k-j/2$.
Relabeling $j=2l$ and $n_{11}=p$ we have $k=l+n$ and we obtain

$\displaystyle\mathrm{Per}(A_{\bm{p},\bm{q}})^{2}$
$\displaystyle=(n!)^{4}\sum_{l=0}^{n}\binom{n+l}{2l}(-1)^{n-l}(a-b)^{2n-2l}\sum_{p=0}^{l}\frac{(2l)!a^{2l-2p}b^{2p}}{p!^{2}(l-p)!^{2}}$
(110)
$\displaystyle=(n!)^{4}\sum_{l=0}^{n}\binom{2l}{l}\binom{n+l}{2l}(-1)^{n-l}(a-b)^{2n-2l}S_{l}(a^{2},b^{2}).$

With Eq. (108) this proves

$S_{n}(a,b)^{2}=\sum_{l=0}^{n}\binom{2l}{l}\binom{n+l}{2l}(-1)^{n-l}(a-b)^{2n-2l}S_{l}(a^{2},b^{2}).$
(111)

∎

###  6.2 New generating functions

We have encountered several generating functions for the permanent of an
$m\times m$ matrix $A$ with differently repeated rows and columns. Jackson’s
formula gives [46]:

$e^{\bm{x}^{T}A\bm{y}}=\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}n\in\mathbb{N},\bm{p},\bm{q}\in\mathbb{N}^{m}\\\
|\bm{p}|=|\bm{q}|=n\end{subarray}}\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}}),$
(112)

for formal variables $\bm{x}=(x_{1},\dots,x_{m}),\bm{y}=(y_{1},\dots,y_{m})$.
Moreover, we have obtained in Eq. (97), as a corollary of our generalization
of the MacMahon master theorem:

$\frac{1}{1-\bm{x}^{T}A\bm{y}}=\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}n\in\mathbb{N},\bm{p},\bm{q}\in\mathbb{N}^{m}\\\
|\bm{p}|=|\bm{q}|=n\end{subarray}}n!\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}}).$
(113)

Furthermore, Corollary 1 gives:

$(\bm{x}^{T}A\bm{y})^{n}=\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{p},\bm{q}\in\mathbb{N}^{m}\\\
|\bm{p}|=|\bm{q}|=n\end{subarray}}n!\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}}).$
(114)

In this section we prove Theorem 3, which generalizes all three statements:
for any series $f(z)=\sum_{n=0}^{+\infty}f_{n}z^{n}$,

$f(\bm{x}^{T}A\bm{y})=\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}n\in\mathbb{N},\bm{p},\bm{q}\in\mathbb{N}^{m}\\\
|\bm{p}|=|\bm{q}|=n\end{subarray}}f_{n}n!\frac{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(A_{\bm{p},\bm{q}}).$
(115)

One may prove this statement by linearity, using Eq. (114). In what follows,
we give a direct quantum-inspired proof which is similar to that of Jackson’s
formula from section 5.4.

###### Proof of Theorem 3.

For $\bm{\alpha},\bm{\beta}\in\mathbb{C}^{m}$ and a unitary matrix $U$,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}n\in\mathbb{N},\bm{p},\bm{q}\in\mathbb{N}^{m}\\\
|\bm{p}|=|\bm{q}|=n\end{subarray}}\\!\\!\\!\\!f_{n}n!\frac{\bm{\alpha}^{\bm{p}}\bm{\beta}^{\bm{q}}}{\bm{p}!\bm{q}!}\mathrm{Per}(U_{\bm{p},\bm{q}})$
$\displaystyle=e^{\frac{1}{2}\|\bm{\alpha}\|^{2}+\frac{1}{2}\|\bm{\beta}\|^{2}}\sum_{n=0}^{+\infty}f_{n}n!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{p},\bm{q}\in\mathbb{N}^{m}\\\
|\bm{p}|=|\bm{q}|=n\end{subarray}}\langle\bm{\alpha}^{*}\ket{\bm{p}}\\!\bra{\bm{p}}\hat{U}\ket{\bm{q}}\\!\bra{\bm{q}}\bm{\beta}\rangle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(116)
$\displaystyle=e^{\frac{1}{2}\|\bm{\alpha}\|^{2}+\frac{1}{2}\|\bm{\beta}\|^{2}}\sum_{n=0}^{+\infty}f_{n}n!\bra{\bm{\alpha}^{*}}\Pi_{n}\hat{U}\Pi_{n}\ket{\bm{\beta}}$
$\displaystyle=e^{\frac{1}{2}\|\bm{\alpha}\|^{2}+\frac{1}{2}\|\bm{\beta}\|^{2}}\sum_{n=0}^{+\infty}f_{n}n!\bra{\bm{\alpha}^{*}}\Pi_{n}\hat{U}\ket{\bm{\beta}}$
$\displaystyle=e^{\frac{1}{2}\|\bm{\alpha}\|^{2}+\frac{1}{2}\|\bm{\beta}\|^{2}}\sum_{n=0}^{+\infty}f_{n}n!\bra{\bm{\alpha}^{*}}\Pi_{n}\ket{U\bm{\beta}},$

where we used the Fock basis expansions of coherent states (6) and of
$\hat{U}$ (10) in the first line, the definition of the projector
$\Pi_{n}=\sum_{|\bm{p}|=n}\ket{\bm{p}}\\!\bra{\bm{p}}$ in the second line, the
fact that passive linear operations conserve the total number of photons (12)
in the third line, and the action of $\hat{U}$ on coherent states (13) in the
fourth line. Now,

$\displaystyle
e^{\frac{1}{2}\|\bm{\alpha}\|^{2}+\frac{1}{2}\|\bm{\beta}\|^{2}}\bra{\bm{\alpha}^{*}}\Pi_{n}\ket{U\bm{\beta}}$
$\displaystyle=e^{\frac{1}{2}\|\bm{\alpha}\|^{2}+\frac{1}{2}\|\bm{\beta}\|^{2}}\sum_{|\bm{p}|=n}\langle\bm{\alpha}^{*}\ket{\bm{p}}\\!\bra{\bm{p}}U\bm{\beta}\rangle$
(117)
$\displaystyle=\sum_{|\bm{p}|=n}e^{\frac{1}{2}\|\bm{\beta}\|^{2}-\frac{1}{2}\|U\bm{\beta}\|^{2}}\frac{\bm{\alpha}^{\bm{p}}(U\bm{\beta})^{\bm{p}}}{\bm{p}!}$
$\displaystyle=\sum_{|\bm{p}|=n}\frac{\bm{\alpha}^{\bm{p}}(U\bm{\beta})^{\bm{p}}}{\bm{p}!}$
$\displaystyle=\frac{1}{n!}(\bm{\alpha}^{T}U\bm{\beta})^{n},$

where we used the Fock basis expansion of coherent states (6) in the second
line, the fact that $U$ is unitary in the third line, and the multinomial
theorem in the last line. Combining Eqs. (116) and (117) completes the proof
for unitary matrices. Once again, the relation for a generic nonzero matrix
$A$ is obtained by embedding $\frac{1}{\|A\|}A$ as a submatrix of a unitary
matrix $U$.

As a result, when $f_{n}\neq 0$ we have, for any $m\times m$ matrix $A$, all
$n\in\mathbb{N}$ and all $\bm{p},\bm{q}\in\mathbb{N}^{m}$ such that
$|\bm{p}|=|\bm{q}|=n$,

$\displaystyle\underset{\bm{x},\bm{y}\in\mathbb{T}^{m}}{\mathbb{E}}\left[\frac{\bm{p}!\bm{q}!}{\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}}\frac{f(\bm{x}^{T}A\bm{y})}{\partial^{n}_{z}f(z)\big{|}_{z=0}}\right]$
$\displaystyle=\frac{\bm{p}!\bm{q}!}{\partial^{n}_{z}f(z)\big{|}_{z=0}}\sum_{\begin{subarray}{c}k\in\mathbb{N},\bm{s},\bm{t}\in\mathbb{N}^{m}\\\
|\bm{s}|=|\bm{t}|=k\end{subarray}}f_{k}k!\frac{\mathrm{Per}(A_{\bm{s},\bm{t}})}{\bm{s}!\bm{t}!}$
(118)
$\displaystyle\quad\times\prod_{j=1}^{m}\left(\int_{\theta_{j},\varphi_{j}\in[0,2\pi]}e^{i\theta_{j}(s_{j}-p_{j})}e^{i\varphi_{j}(t_{j}-q_{j})}\frac{d\theta_{j}d\varphi_{j}}{4\pi^{2}}\right)$
$\displaystyle=\frac{1}{\partial^{n}_{z}f(z)\big{|}_{z=0}}f_{n}n!\mathrm{Per}(A_{\bm{p},\bm{q}})$
$\displaystyle=\mathrm{Per}(A_{\bm{p},\bm{q}}),$

where we have set
$\bm{x}=(e^{i\theta_{1}},\dots,e^{i\theta_{m}})\in\mathbb{T}^{m}$ and
$\bm{y}=(e^{i\varphi_{1}},\dots,e^{i\varphi_{m}})\in\mathbb{T}^{m}$ in the
first line, and where we have used
$\int_{0}^{2\pi}e^{i(k-l)\theta}\frac{d\theta}{2\pi}=\delta_{kl}$ in the third
line. ∎

Applications. These generating functions may be used to derive remarkable
identities for the permanent: for instance, for $A$ and $B$ two $m\times m$
matrices, taking $f(z)=e^{z}$ and equating the
$\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}$ coefficients of $e^{\bm{x}^{T}(A+B)\bm{y}}$
and $e^{\bm{x}^{T}A\bm{y}}e^{\bm{x}^{T}B\bm{y}}$ yields the sum formula [1,
2]:

$\mathrm{Per}((A+B)_{\bm{p},\bm{q}})=\sum_{\begin{subarray}{c}\bm{s}+\bm{t}=\bm{p}\\\
\bm{u}+\bm{v}=\bm{q}\end{subarray}}\frac{\bm{p}!\bm{q}!}{\bm{s}!\bm{t}!\bm{u}!\bm{v}!}\mathrm{Per}(A_{\bm{s},\bm{u}})\mathrm{Per}(B_{\bm{t},\bm{v}}),$
(119)

for all $\bm{p},\bm{q}\in\mathbb{N}^{m}$.

Similarly, for $A$ an $m\times m$ matrix, taking $f(z)=z^{p}$ for $p=k,l,k+l$
and equating the $\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}$ coefficients of
$(\bm{x}^{T}A\bm{y})^{k+l}$ and
$(\bm{x}^{T}A\bm{y})^{k}(\bm{x}^{T}A\bm{y})^{l}$ yields the Laplace expansion
formula [1, 2]:

$\mathrm{Per}(A_{\bm{p},\bm{q}})=\frac{k!l!}{(k+l)!}\sum_{\begin{subarray}{c}\bm{s}+\bm{t}=\bm{p}\\\
\bm{u}+\bm{v}=\bm{q}\\\ |\bm{s}|=|\bm{u}|=k\\\
|\bm{t}|=|\bm{v}|=l\end{subarray}}\frac{\bm{p}!\bm{q}!}{\bm{s}!\bm{t}!\bm{u}!\bm{v}!}\mathrm{Per}(A_{\bm{s},\bm{u}})\mathrm{Per}(A_{\bm{t},\bm{v}}),$
(120)

for all $k,l\in\mathbb{N}$ and all $\bm{p},\bm{q}\in\mathbb{N}^{m}$ with
$|\bm{p}|=|\bm{q}|=k+l$.

Finally, for $A$ and $B$ two $m\times m$ matrices, taking $f(z)=-\log(1-z)$
and equating the $\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}$ coefficients of
$-\log(1-\bm{x}^{T}A\bm{y})-\log(1-\bm{x}^{T}B\bm{y})$ and
$-\log[(1-\bm{x}^{T}B\bm{y})(1-\bm{x}^{T}A\bm{y})]$ yields Theorem 4:

$\displaystyle\mathrm{Per}(A_{\bm{p},\bm{q}})+\mathrm{Per}(B_{\bm{p},\bm{q}})$
(121)
$\displaystyle\quad=\sum_{k=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{k}}{\binom{n-1}{k}}\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{a}+\bm{b}+\bm{c}=\bm{p}\\\
\bm{a}^{\prime}+\bm{b}^{\prime}+\bm{c}^{\prime}=\bm{q}\\\
|\bm{a}|=|\bm{b}|=|\bm{a}^{\prime}|=|\bm{b}^{\prime}|=k\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\frac{\bm{p}!\bm{q}!}{\bm{a}!\bm{b}!\bm{c}!\bm{a}^{\prime}!\bm{b}^{\prime}!\bm{c}^{\prime}!}\mathrm{Per}(A_{\bm{a},\bm{a}^{\prime}})\mathrm{Per}(B_{\bm{b},\bm{b}^{\prime}})\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}}),$

after a derivation which we detail below.

###### Proof of Theorem 4.

Let us write

$-\log[(1-\bm{x}^{T}B\bm{y})(1-\bm{x}^{T}A\bm{y})]=-\log[1-\bm{x}^{T}(A+B-B\bm{y}\bm{x}^{T}A)\bm{y})].$
(122)

Applying Theorem 3 for $f(z)=-\log(1-z)$ with
$z=\bm{x}^{T}B\bm{y},\bm{x}^{T}A\bm{y},\bm{x}^{T}(A+B-B\bm{y}\bm{x}^{T}A)\bm{y}$
and considering the $\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}$ coefficient we obtain for
all $n\in\mathbb{N}^{*}$ and all $\bm{p},\bm{q}\in\mathbb{N}$ with
$|\bm{p}|=|\bm{q}|=n$:

$\displaystyle\mathrm{Per}(A_{\bm{p},\bm{q}})+\mathrm{Per}(B_{\bm{p},\bm{q}})$
(123)
$\displaystyle\quad\quad=\frac{\bm{p}!\bm{q}!}{(n-1)!}[\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}]\left(\sum_{\begin{subarray}{c}l\in\mathbb{N}^{*},\bm{s},\bm{t}\in\mathbb{N}^{m}\\\
|\bm{s}|=|\bm{t}|=l\end{subarray}}(l-1)!\frac{\bm{x}^{\bm{s}}\bm{y}^{\bm{t}}}{\bm{s}!\bm{t}!}\mathrm{Per}((A+B-B\bm{y}\bm{x}^{T}A)_{\bm{s},\bm{t}})\right)$
$\displaystyle\quad\quad=\sum_{\begin{subarray}{c}l\in\mathbb{N}^{*},\bm{s},\bm{t}\in\mathbb{N}^{m}\\\
|\bm{s}|=|\bm{t}|=l\end{subarray}}\frac{(l-1)!}{(n-1)!}\frac{\bm{p}!\bm{q}!}{\bm{s}!\bm{t}!}[\bm{x}^{\bm{p}}\bm{y}^{\bm{q}}]\left(\bm{x}^{\bm{s}}\bm{y}^{\bm{t}}\mathrm{Per}((A+B-B\bm{y}\bm{x}^{T}A)_{\bm{s},\bm{t}})\right)$
$\displaystyle\quad\quad=\sum_{\begin{subarray}{c}l\in\mathbb{N}^{*},\bm{s},\bm{t}\in\mathbb{N}^{m}\\\
|\bm{s}|=|\bm{t}|=l\\\
\bm{s}\leq\bm{p},\bm{t}\leq\bm{q}\end{subarray}}\frac{(l-1)!}{(n-1)!}\frac{\bm{p}!\bm{q}!}{\bm{s}!\bm{t}!}[\bm{x}^{\bm{p}-\bm{s}}\bm{y}^{\bm{q}-\bm{t}}]\left(\mathrm{Per}((A+B-B\bm{y}\bm{x}^{T}A)_{\bm{s},\bm{t}})\right)$
$\displaystyle\quad\quad=\sum_{\begin{subarray}{c}l\in\mathbb{N}^{*},\bm{s},\bm{t},\bm{u},\bm{v}\in\mathbb{N}^{m}\\\
|\bm{s}|=|\bm{t}|=l\\\ \bm{s}+\bm{u}=\bm{p}\\\
\bm{t}+\bm{v}=\bm{q}\end{subarray}}\frac{(l-1)!}{(n-1)!}\frac{\bm{p}!\bm{q}!}{\bm{s}!\bm{t}!}[\bm{x}^{\bm{u}}\bm{y}^{\bm{v}}]\left(\mathrm{Per}((A+B-B\bm{y}\bm{x}^{T}A)_{\bm{s},\bm{t}})\right).$

Applying the sum formula from Eq. (119) to the matrices $(A+B)$ and
$(-B\bm{y}\bm{x}^{T}A)$ we obtain, for all $l\in\mathbb{N}^{*}$ and all
$\bm{s},\bm{t}\in\mathbb{N}$ with $|\bm{s}|=|\bm{t}|=l$,

$\displaystyle\mathrm{Per}((A+B-B\bm{y}\bm{x}^{T}A)_{\bm{s},\bm{t}})$
$\displaystyle=\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{c}+\bm{k}=\bm{s}\\\
\bm{c}^{\prime}+\bm{k}^{\prime}=\bm{t}\\\
|\bm{c}|=|\bm{c}^{\prime}|,|\bm{k}|=|\bm{k}^{\prime}|\end{subarray}}\\!\\!\\!\\!\\!\frac{\bm{s}!\bm{t}!}{\bm{c}!\bm{k}!\bm{c}^{\prime}!\bm{k}^{\prime}!}\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}})\mathrm{Per}((-B\bm{y}\bm{x}^{T}A)_{\bm{k},\bm{k}^{\prime}})\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(124)
$\displaystyle=\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{c}+\bm{k}=\bm{s}\\\
\bm{c}^{\prime}+\bm{k}^{\prime}=\bm{t}\\\
|\bm{c}|=|\bm{c}^{\prime}|,|\bm{k}|=|\bm{k}^{\prime}|\end{subarray}}\\!\\!\\!\\!\\!\frac{(-1)^{|\bm{k}|}\bm{s}!\bm{t}!}{\bm{c}!\bm{k}!\bm{c}^{\prime}!\bm{k}^{\prime}!}\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}})\mathrm{Per}((B\bm{y}\bm{x}^{T}A)_{\bm{k},\bm{k}^{\prime}}).\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$

The permanent of the outer product $\bm{v}\bm{w}^{T}$ of two vectors $\bm{v}$
and $\bm{w}$ of the same size $k$ is easily computed as

$\displaystyle\mathrm{Per}(\bm{v}\bm{w}^{T})$
$\displaystyle=\sum_{\sigma\in\mathcal{S}_{k}}\prod_{j=1}^{k}v_{j}w_{\sigma(j)}$
(125) $\displaystyle=\sum_{\sigma\in\mathcal{S}_{k}}v_{1}\dots v_{k}w_{1}\dots
w_{k}$ $\displaystyle=k!\,v_{1}\dots v_{k}w_{1}\dots w_{k}.$

The matrix $(B\bm{y}\bm{x}^{T}A)_{\bm{k},\bm{k}^{\prime}}$ is the outer
product between the vector $(B\bm{y})_{\bm{k}}$ (obtained from the vector
$B\bm{y}$ by repeating $k_{i}$ times its $i^{th}$ entry) and the vector
$(A^{T}\bm{x})_{\bm{k}^{\prime}}$ (obtained from the vector $A^{T}\bm{x}$ by
repeating $k_{i}^{\prime}$ times its $i^{th}$ entry) of size $|\bm{k}|$, so
its permanent is given by

$\mathrm{Per}((B\bm{y}\bm{x}^{T}A)_{\bm{k},\bm{k}^{\prime}})=(|\bm{k}|)!(B\bm{y})^{\bm{k}}(A^{T}\bm{x})^{\bm{k}^{\prime}}.$
(126)

With Eq. (124) we obtain

$\mathrm{Per}((A+B-B\bm{y}\bm{x}^{T}A)_{\bm{s},\bm{t}}=\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{c}+\bm{k}=\bm{s}\\\
\bm{c}^{\prime}+\bm{k}^{\prime}=\bm{t}\\\
|\bm{c}|=|\bm{c}^{\prime}|,|\bm{k}|=|\bm{k}^{\prime}|\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{(-1)^{|\bm{k}|}(|\bm{k}|)!\bm{s}!\bm{t}!}{\bm{c}!\bm{k}!\bm{c}^{\prime}!\bm{k}^{\prime}!}\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}})(B\bm{y})^{\bm{k}}(A^{T}\bm{x})^{\bm{k}^{\prime}}.$
(127)

Plugging this expression into Eq. (123) yields

$\displaystyle\quad\mathrm{Per}(A_{\bm{p},\bm{q}})+\mathrm{Per}(B_{\bm{p},\bm{q}})$
$\displaystyle=\sum_{\begin{subarray}{c}l\in\mathbb{N}^{*},\bm{s},\bm{t},\bm{u},\bm{v}\in\mathbb{N}^{m}\\\
|\bm{s}|=|\bm{t}|=l\\\ \bm{s}+\bm{u}=\bm{p}\\\
\bm{t}+\bm{v}=\bm{q}\end{subarray}}\\!\frac{(l-1)!}{(n-1)!}\frac{\bm{p}!\bm{q}!}{\bm{s}!\bm{t}!}\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{c},\bm{k},\bm{c}^{\prime},\bm{k}^{\prime}\in\mathbb{N}^{m}\\\
\bm{c}+\bm{k}=\bm{s}\\\ \bm{c}^{\prime}+\bm{k}^{\prime}=\bm{t}\\\
|\bm{c}|=|\bm{c}^{\prime}|,|\bm{k}|=|\bm{k}^{\prime}|\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{(-1)^{|\bm{k}|}(|\bm{k}|)!\bm{s}!\bm{t}!}{\bm{c}!\bm{k}!\bm{c}^{\prime}!\bm{k}^{\prime}!}\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}})[\bm{x}^{\bm{u}}\bm{y}^{\bm{v}}](B\bm{y})^{\bm{k}}(A^{T}\bm{x})^{\bm{k}^{\prime}}$
$\displaystyle=\sum_{\begin{subarray}{c}l\in\mathbb{N}^{*},\bm{s},\bm{t},\bm{u},\bm{v}\in\mathbb{N}^{m}\\\
|\bm{s}|=|\bm{t}|=l\\\ \bm{s}+\bm{u}=\bm{p}\\\
\bm{t}+\bm{v}=\bm{q}\end{subarray}}\\!\frac{(l-1)!}{(n-1)!}\bm{p}!\bm{q}!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{c},\bm{k},\bm{c}^{\prime},\bm{k}^{\prime}\in\mathbb{N}^{m}\\\
\bm{c}+\bm{k}=\bm{s}\\\ \bm{c}^{\prime}+\bm{k}^{\prime}=\bm{t}\\\
|\bm{c}|=|\bm{c}^{\prime}|,|\bm{k}|=|\bm{k}^{\prime}|\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\frac{(-1)^{|\bm{k}|}(|\bm{k}|)!}{\bm{c}!\bm{k}!\bm{c}^{\prime}!\bm{k}^{\prime}!}\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}})[\bm{x}^{\bm{u}}](A^{T}\bm{x})^{\bm{k}^{\prime}}[\bm{y}^{\bm{v}}](B\bm{y})^{\bm{k}}$
$\displaystyle=\sum_{\begin{subarray}{c}l\in\mathbb{N}^{*},\bm{s},\bm{t},\bm{u},\bm{v}\in\mathbb{N}^{m}\\\
|\bm{s}|=|\bm{t}|=l\\\ \bm{s}+\bm{u}=\bm{p}\\\
\bm{t}+\bm{v}=\bm{q}\end{subarray}}\\!\frac{(l-1)!}{(n-1)!}\bm{p}!\bm{q}!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{c},\bm{k},\bm{c}^{\prime},\bm{k}^{\prime}\in\mathbb{N}^{m}\\\
\bm{c}+\bm{k}=\bm{s}\\\ \bm{c}^{\prime}+\bm{k}^{\prime}=\bm{t}\\\
|\bm{c}|=|\bm{c}^{\prime}|\\\
|\bm{k}|=|\bm{k}^{\prime}|=|\bm{u}|=|\bm{v}|\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\frac{(-1)^{|\bm{k}|}(|\bm{k}|)!}{\bm{c}!\bm{k}!\bm{c}^{\prime}!\bm{k}^{\prime}!}\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}})\frac{\mathrm{Per}(A^{T}_{\bm{k}^{\prime},\bm{u}})}{\bm{u}!}\frac{\mathrm{Per}(B_{\bm{k},\bm{v}})}{\bm{v}!}$
(128)
$\displaystyle=\sum_{\begin{subarray}{c}l\in\mathbb{N}^{*},\bm{s},\bm{t},\bm{u},\bm{v}\in\mathbb{N}^{m}\\\
|\bm{s}|=|\bm{t}|=l\\\ \bm{s}+\bm{u}=\bm{p}\\\
\bm{t}+\bm{v}=\bm{q}\end{subarray}}\\!\frac{(l-1)!}{(n-1)!}\frac{\bm{p}!\bm{q}!}{\bm{u}!\bm{v}!}\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{c},\bm{k},\bm{c}^{\prime},\bm{k}^{\prime}\in\mathbb{N}^{m}\\\
\bm{c}+\bm{k}=\bm{s}\\\ \bm{c}^{\prime}+\bm{k}^{\prime}=\bm{t}\\\
|\bm{c}|=|\bm{c}^{\prime}|\\\
|\bm{k}|=|\bm{k}^{\prime}|=|\bm{u}|=|\bm{v}|\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\frac{(-1)^{|\bm{k}|}(|\bm{k}|)!}{\bm{c}!\bm{k}!\bm{c}^{\prime}!\bm{k}^{\prime}!}\mathrm{Per}(A_{\bm{u},\bm{k}^{\prime}})\mathrm{Per}(B_{\bm{k},\bm{v}})\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}})$
$\displaystyle=\sum_{k=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{k}k!(n-k-1)!}{(n-1)!}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{a},\bm{b},\bm{c},\bm{a}^{\prime},\bm{b}^{\prime},\bm{c}^{\prime}\in\mathbb{N}^{m}\\\
\bm{a}+\bm{b}+\bm{c}=\bm{p}\\\
\bm{a}^{\prime}+\bm{b}^{\prime}+\bm{c}^{\prime}=\bm{q}\\\
|\bm{a}|=|\bm{a}^{\prime}|=|\bm{b}|=|\bm{b}^{\prime}|=k\end{subarray}}\\!\\!\frac{\bm{p}!\bm{q}!}{\bm{a}!\bm{b}!\bm{c}!\bm{a}^{\prime}!\bm{b}^{\prime}!\bm{c}^{\prime}!}\mathrm{Per}(A_{\bm{a},\bm{a}^{\prime}})\mathrm{Per}(B_{\bm{b},\bm{b}^{\prime}})\mathrm{Per}((A+B)_{\bm{c},\bm{c}^{\prime}}),$

where we used Eq. (61) in the third line,
$\mathrm{Per}(M^{T})=\mathrm{Per}(M)$ in the fourth line, and where we
relabeled $\bm{u}$ as $\bm{a}$, $\bm{k}^{\prime}$ as $\bm{a}^{\prime}$,
$\bm{k}$ as $\bm{b}$, $\bm{v}$ as $\bm{b}^{\prime}$, and
$|\bm{k}|=|\bm{k}^{\prime}|=|\bm{u}|=|\bm{v}|$ as $k=n-l$ in the last line.
The cutoff of the sum at $\lfloor\frac{n}{2}\rfloor$ comes from the fact that
$n=|\bm{p}|=|\bm{a}|+|\bm{b}|+|\bm{c}|=2k+|\bm{c}|\geq 2k$. Writing
$\frac{(-1)^{k}k!(n-k-1)!}{(n-1)!}=\frac{(-1)^{k}}{\binom{n-1}{k}}$ completes
the proof of Theorem 4. ∎

##  7 Boson Sampling with input cat states

In this section, we discuss the classical complexity of sampling from the
output distribution of linear optical computations with input cat states,
which we refer to as Boson Sampling with input cat states.

Since the introduction of Boson Sampling by Aaronson and Arkhipov [6] for the
demonstration of quantum computational advantage using noninteracting bosons,
several variants of this model have been analyzed [48, 49, 50, 51, 52, 53, 54,
55]. These variants were introduced to address two different challenges: on
the one hand, to reduce the experimental burden associated with the
demonstration of quantum advantage; on the other hand, to understand the
resources for this quantum advantage and which computational models are able
to reproduce it. Boson Sampling with input cat states is another such variant
and has been first considered in [21], where its classical hardness was shown
for exact sampling and argued for approximate sampling. We strengthen these
results in the following.

We first give a proof of Lemma 1, which provides a closed form expression for
output amplitudes of Boson Sampling with input cat states.

###### Proof of Lemma 1.

In the quantum-inspired proof of Glynn’s formula in section 5.1 we have
obtained the following identity:

$\langle p_{1}\dots p_{m}|\hat{U}(\ket{\tilde{\mathrm{cat}}_{\alpha}}^{\otimes
n}\otimes\ket{0}^{\otimes
m-n})=\frac{\alpha^{|\bm{p}|-n}e^{-\frac{n}{2}|\alpha|^{2}}}{2^{n}\sqrt{\bm{p}!}}\sum_{\bm{x}\in\\{-1,1\\}^{n}}x_{1}\dots
x_{n}\prod_{i=1}^{m}\left(\sum_{j=1}^{n}u_{ij}x_{j}\right)^{\\!p_{i}}\\!\\!\\!\\!\\!,$
(129)

for all $\alpha\in\mathbb{C}$, all
$\bm{p}=(p_{1},\dots,p_{m})\in\mathbb{N}^{m}$, all $n\leq m$ and any passive
linear operation $\hat{U}$ over $m$ modes with unitary matrix
$U=(u_{ij})_{1\leq i,j\leq m}$, where
$\ket{\tilde{\mathrm{cat}}_{\alpha}}=\frac{1}{2\alpha}(\ket{\alpha}-\ket{-\alpha})$
is an unnormalized cat state. By definition of the cat states (5), these
states are related to normalized cat states $\ket{\mathrm{cat}_{\alpha}}$ by

$\ket{\mathrm{cat}_{\alpha}}=\frac{\alpha
e^{\frac{1}{2}|\alpha|^{2}}}{\sqrt{\sinh(|\alpha|^{2})}}\ket{\tilde{\mathrm{cat}}_{\alpha}}.$
(130)

Setting $|\bm{p}|=n$ and replacing unnormalized cat states by normalized ones
in Eq. (129), we obtain Lemma 1:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\langle p_{1}\dots
p_{m}|\hat{U}(\ket{\mathrm{cat}_{\alpha}}^{\otimes n}\otimes\ket{0}^{\otimes
m-n})$
$\displaystyle=\frac{\alpha^{n}}{\sqrt{\sinh^{n}(|\alpha|^{2})}2^{n}\sqrt{\bm{p}!}}\\!\sum_{\bm{x}\in\\{-1,1\\}^{n}}\\!\\!\\!\\!\\!x_{1}\dots
x_{n}\prod_{i=1}^{m}\\!\left(\sum_{j=1}^{n}u_{ij}x_{j}\right)^{\\!p_{i}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(131)
$\displaystyle=\frac{\alpha^{n}}{\sqrt{\sinh^{n}(|\alpha|^{2})}}\frac{\mathrm{Per}(U_{\bm{p},\bm{1}\oplus\bm{0}})}{\sqrt{\bm{p}!}}$
$\displaystyle=\frac{\alpha^{n}}{\sqrt{\sinh^{n}(|\alpha|^{2})}}\langle
p_{1}\dots p_{m}|\hat{U}(\ket{1}^{\otimes n}\otimes\ket{0}^{\otimes m-n}),$

where we used Eq. (42) in the second line and the Fock basis expansion of
$\hat{U}$ (10) in the last line. ∎

Interestingly, Lemma 1 implies that, up to a global factor, $n$ cat states of
amplitude $\alpha$ reproduce exactly the Boson Sampling statistics of $n$
single photons (however, detection events $\bm{p}$ with $|\bm{p}|>n$ can occur
for input cat states, contrary to single-photons).

This result has two consequences, summarized by Theorem 5. Firstly, it implies
that Boson Sampling with input cat states is hard to sample exactly for $m\geq
2n$ and for all choices of nonzero cat state amplitude unless the polynomial
hierarchy collapses to its third level, since some of its outcome
probabilities are $\\#\textsf{P}$-hard to estimate multiplicatively [6].
Alternative proofs of this statement based either on universality under post-
selection or on rejection sampling can be found in [21]. Secondly, as we show
hereafter, it implies that Boson Sampling with input cat states of small
enough (nonzero) amplitudes is also hard to sample approximately in the same
regime as Boson Sampling.

###### Proof of Theorem 5.

Our proof uses arguments similar to the ones used in [21], extended to the
case of approximate sampling.

At a high level, we show how to convert any classical algorithm for
approximate sampling from the output probability distribution of an $m$-mode
Boson Sampling computation with $n$ input cat states to a classical algorithm
for approximate sampling from the output probability distribution of an
$m$-mode Boson Sampling computation with $n$ input Fock states (Lemma 2). We
do so by performing rejection sampling: we run the first algorithm and only
keep the samples with total photon number equal to $n$. In particular, we show
that if the fraction of samples with correct photon number $n$ is large enough
(namely, at least inverse polynomial in $m$), the rejection sampling
subroutine is efficient and the final classical probability distribution is
close to the ideal distribution of Boson Sampling with Fock state input. Then,
we determine the regime of input cat state amplitude such that the fraction of
samples with correct photon number $n$ is at least inverse polynomial in $m$
(Lemma 3).

Let $P_{\text{cat}}(\bm{p}|n,\alpha):=|\langle p_{1}\dots
p_{m}|\hat{U}(\ket{\mathrm{cat}_{\alpha}}^{\otimes n}\otimes\ket{0}^{\otimes
m-n})|^{2}$ be the probability of detecting $\bm{p}=(p_{1},\dots,p_{m})$
output photons in a Boson Sampling experiment with interferometer $\hat{U}$
and input cat states $\ket{\mathrm{cat}_{\alpha}}^{\otimes
n}\otimes\ket{0}^{\otimes m-n}$, and let $P_{\text{BS}}(\bm{p}|n):=|\langle
p_{1}\dots p_{m}|\hat{U}(\ket{1}^{\otimes n}\otimes\ket{0}^{\otimes
m-n})|^{2}$ be the probability of detecting $\bm{p}=(p_{1},\dots,p_{m})$
output photons in a Boson Sampling experiment with interferometer $\hat{U}$
and input single photons $\ket{1}^{\otimes n}\otimes\ket{0}^{\otimes m-n}$.

We show the following reduction: as long as the fraction of samples with
photon number $n$ is large enough, then any efficient classical algorithm for
approximate sampling from $P_{\text{cat}}$ can be converted to an efficient
classical algorithm for approximate sampling from $P_{\text{BS}}$.

######  Lemma 2.

Suppose there exists an efficient classical algorithm $C$ for approximate
sampling from $P_{\text{cat}}$, i.e., that takes as input the description of
the $m$-mode boson sampler with $n$ input cat states with output probability
distribution $P_{\text{cat}}$ and an error bound $\epsilon$, and that samples
from a distribution $P_{C}$ such that
$\|P_{C}-P_{\text{cat}}\|_{\text{TV}}\leq\epsilon$ in $\text{poly
}(m,\frac{1}{\epsilon})$ time, where $\|\cdot\|_{\text{TV}}$ is the total
variation distance.

Assume further that the fraction of samples from $P_{\text{cat}}$ with photon
number equal to $n$ is at least inverse polynomial in $m$, i.e.
$\sum_{|\bm{p}|=n}P_{\text{cat}}(\bm{p}|n,\alpha)\geq\frac{1}{\text{poly }m}$.

Then, there exists an efficient classical algorithm $\tilde{C}$ for
approximate sampling from $P_{\text{BS}}$, i.e., a classical algorithm that
takes as input the description of the $m$-mode boson sampler with $n$ input
Fock states with output probability distribution $P_{\text{BS}}$ and an error
bound $\epsilon$, and that samples from a distribution $P_{C}^{\text{rej}}$
such that $\|P_{C}^{\text{rej}}-P_{\text{BS}}\|_{\text{TV}}\leq\epsilon$ in
$\text{poly }(m,\frac{1}{\epsilon})$ time.

###### Proof of Lemma 2.

Given a sampling algorithm from a probability distribution $P$ over
$\mathbb{N}^{m}$, we define the following rejection sampling subroutine:
sample $\bm{p}$ from $P$ and compute $|\bm{p}|$; discard the sample if
$|\bm{p}|\neq n$; otherwise, output $\bm{p}$. We denote by $P^{\text{rej}}$
the corresponding output probability distribution. For all
$\bm{p}\in\mathbb{N}^{m}$, we have

$P^{\text{rej}}(\bm{p})=\frac{\delta_{|\bm{p}|,n}P(\bm{p})}{\sum_{|\bm{q}|=n}P(\bm{q})},$
(132)

where $\delta$ is the Kronecker symbol. Applying this subroutine to
$P_{\text{cat}}$ gives a new probability distribution

$P_{\text{cat}}^{\text{rej}}(\bm{p}|n,\alpha)=\frac{\delta_{|\bm{p}|,n}P_{\text{cat}}(\bm{p}|n,\alpha)}{\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)}.$
(133)

By Lemma 1 (see Eq. (131)) we have, for all $\bm{p}$ such that $|\bm{p}|=n$,

$P_{\text{cat}}(\bm{p}|n,\alpha)=\frac{|\alpha|^{2n}}{\sinh^{n}(|\alpha|^{2})}P_{\text{BS}}(\bm{p}|n).$
(134)

Hence

$\displaystyle\sum_{|\bm{p}|=n}P_{\text{cat}}(\bm{p}|n,\alpha)$
$\displaystyle=\frac{|\alpha|^{2n}}{\sinh^{n}(|\alpha|^{2})}\sum_{|\bm{p}|=n}P_{\text{BS}}(\bm{p}|n)$
(135) $\displaystyle=\frac{|\alpha|^{2n}}{\sinh^{n}(|\alpha|^{2})},$

where we used $\sum_{|\bm{p}|=n}P_{\text{BS}}(\bm{p}|n)=1$ in the second line
since $\hat{U}$ conserves the total number of photons (12). Combining Eqs.
(133-135), we thus obtain

$P_{\text{BS}}(\bm{p}|n)=P_{\text{cat}}^{\text{rej}}(\bm{p}|n,\alpha).$ (136)

This reproduces the argument for exact sampling hardness from [21].

To show approximate hardness, we define the classical approximate sampling
algorithm $\tilde{C}$ by running the algorithm $C$ with input error bound
$\frac{\epsilon}{2}\sum_{|\bm{p}|=n}P_{\text{cat}}(\bm{p}|n,\alpha)$ and
applying the rejection sampling subroutine. By assumption, running the
algorithm $C$ can be done in time

$\text{poly}\left(m,\frac{2}{\epsilon\sum_{|\bm{p}|=n}P_{\text{cat}}(\bm{p}|n,\alpha)}\right)=\text{poly}\left(m,\frac{1}{\epsilon}\right),$
(137)

since $\sum_{|\bm{p}|=n}P_{\text{cat}}(\bm{p}|n,\alpha)\geq\frac{1}{\text{poly
}m}$, also by assumption. On the other hand, the rejection sampling subroutine
induces a computational overhead scaling with the inverse of the fraction of
kept samples $\sum_{|\bm{q}|=n}P_{C}(\bm{q})$. We have

$\displaystyle\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)-\sum_{|\bm{q}|=n}P_{C}(\bm{q})$
$\displaystyle\leq\left|\sum_{|\bm{q}|=n}P_{C}(\bm{q})-\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)\right|$
(138)
$\displaystyle\leq\sum_{|\bm{q}|=n}\left|P_{C}(\bm{q})-P_{\text{cat}}(\bm{q}|n,\alpha)\right|$
$\displaystyle\leq\sum_{\bm{q}\in\mathbb{N}^{m}}\left|P_{C}(\bm{q})-P_{\text{cat}}(\bm{q}|n,\alpha)\right|$
$\displaystyle=2\|P_{C}-P_{\text{cat}}\|_{\text{TV}}$
$\displaystyle\leq\epsilon\sum_{|\bm{p}|=n}P_{\text{cat}}(\bm{p}|n,\alpha),$

where we used the triangle inequality in the second line, the definition of
the total variation distance in the fourth line, and the definition of $C$ in
the last line. Hence,

$\displaystyle\sum_{|\bm{q}|=n}P_{C}(\bm{q})$
$\displaystyle\geq(1-\epsilon)\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)$
(139) $\displaystyle\geq\frac{1-\epsilon}{\text{poly }m},$

so that the computational overhead induced by the rejection sampling
subroutine scales as $\text{poly}(m,\frac{1}{\epsilon})$.

Overall, the classical algorithm $\tilde{C}$ outputs samples from the
probability distribution $P_{C}^{\text{rej}}$ in
$\text{poly}(m,\frac{1}{\epsilon})$ time (with probability exponentially close
to $1$). Moreover,

$\displaystyle\|P_{C}^{\text{rej}}-P_{\text{BS}}\|_{\text{TV}}$
$\displaystyle=\frac{1}{2}\sum_{\bm{p}\in\mathbb{N}^{m}}|P_{C}^{\text{rej}}(\bm{p})-P_{\text{BS}}(\bm{p}|n)|$
$\displaystyle=\frac{1}{2}\sum_{|\bm{p}|=n}|P_{C}^{\text{rej}}(\bm{p})-P_{\text{BS}}(\bm{p}|n)|$
$\displaystyle=\frac{1}{2}\sum_{|\bm{p}|=n}\left|P_{C}^{\text{rej}}(\bm{p})-P_{\text{cat}}^{\text{rej}}(\bm{p}|n,\alpha)\right|$
$\displaystyle=\frac{1}{2}\sum_{|\bm{p}|=n}\left|\frac{P_{C}(\bm{p})}{\sum_{|\bm{q}|=n}P_{C}(\bm{q})}-\frac{P_{\text{cat}}(\bm{p}|n,\alpha)}{\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)}\right|$
$\displaystyle\leq\frac{1}{2}\sum_{|\bm{p}|=n}\Bigg{[}\frac{\left|P_{C}(\bm{p})-P_{\text{cat}}(\bm{p}|n,\alpha)\right|}{\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)}$
(140)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad+P_{C}(\bm{p})\left|\frac{1}{\sum_{|\bm{q}|=n}P_{C}(\bm{q})}-\frac{1}{\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)}\right|\Bigg{]}$
$\displaystyle\leq\frac{1}{2\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)}\sum_{|\bm{p}|=n}\Bigg{[}\left|P_{C}(\bm{p})-P_{\text{cat}}(\bm{p}|n,\alpha)\right|$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad+\frac{P_{C}(\bm{p})}{\sum_{|\bm{q}|=n}P_{C}(\bm{q})}\left|\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)-\sum_{|\bm{q}|=n}P_{C}(\bm{q})\right|\Bigg{]}$
$\displaystyle\leq\frac{1}{\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)}\Bigg{[}\frac{1}{2}\sum_{\bm{p}\in\mathbb{N}^{m}}\left|P_{C}(\bm{p})-P_{\text{cat}}(\bm{p}|n,\alpha)\right|$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad+\sum_{|\bm{p}|=n}\frac{P_{C}(\bm{p})}{\sum_{|\bm{q}|=n}P_{C}(\bm{q})}\frac{1}{2}\sum_{\bm{q}\in\mathbb{N}^{m}}\left|P_{\text{cat}}(\bm{q}|n,\alpha)-P_{C}(\bm{q})\right|\Bigg{]}$
$\displaystyle=\frac{2\|P_{C}-P_{\text{cat}}\|_{\text{TV}}}{\sum_{|\bm{q}|=n}P_{\text{cat}}(\bm{q}|n,\alpha)}$
$\displaystyle\leq\epsilon.$

where we used the definition of the total variation distance in the first and
eighth lines, the fact that $P_{C}^{\text{rej}}$ and $P_{\text{BS}}$ are
supported on samples with total photon number $n$ in the second line, Eq.
(136) in the third line, Eq. (132) in the fourth line and the triangle
inequality in the fifth line.

Hence, $\tilde{C}$ is an efficient classical algorithm for approximate
sampling from $P_{\text{BS}}$. ∎

We are left with identifying the regime of input cat state amplitude
$\alpha\in\mathbb{C}$ such that the assumption in Lemma 2 is satisfied, i.e.
$\sum_{|\bm{p}|=n}P_{\text{cat}}(\bm{p}|n,\alpha)\geq\frac{1}{\text{poly }m}$.

######  Lemma 3.

For $0<|\alpha|=\mathcal{O}(n^{-1/4}\log^{1/4}m)$ and $m=\text{poly }n$, the
detection events satisfying $|\bm{p}|=n$ represent at least an inverse-
polynomial fraction of the outcomes for Boson Sampling with input cat states
of amplitude $\alpha$, i.e.
$\sum_{|\bm{p}|=n}P_{\text{cat}}(\bm{p}|n,\alpha)\geq\frac{1}{\text{poly }m}$.

###### Proof of Lemma 3.

Let $\alpha\in\mathbb{C}$ such that
$0<|\alpha|=\mathcal{O}(n^{-1/4}\log^{1/4}m)$. The fraction of outcomes
$\bm{p}\in\mathbb{N}^{m}$ for $m$-mode Boson Sampling with $n$ input cat
states of amplitude $\alpha$ is given by

$\displaystyle\sum_{|\bm{p}|=n}P_{\text{cat}}(\bm{p}|n,\alpha)$
$\displaystyle=\frac{|\alpha|^{2n}}{\sinh^{n}(|\alpha|^{2})}$ (141)
$\displaystyle=e^{n\log(|\alpha|^{2})-n\log(\sinh|\alpha|^{2})}$
$\displaystyle=e^{n\log(|\alpha|^{2})-n\log(|\alpha|^{2}+\frac{1}{6}|\alpha|^{6}+O(|\alpha|^{10}))}$
$\displaystyle=e^{-n\log(1+\frac{1}{6}|\alpha|^{4}+O(|\alpha|^{8}))}$
$\displaystyle=e^{-\frac{1}{6}n|\alpha|^{4}+O(n|\alpha|^{8})}$
$\displaystyle=e^{-O(\log m)+O(n^{-2}\log^{2}m)}$
$\displaystyle\geq\frac{1}{\text{poly }m},$

where we used Eq. (135) in the first line, the scaling of $|\alpha|$ in the
sixth line, and $O(n^{-2}\log^{2}m)=O(1)$ in the last line. ∎

Combining Lemma 2 and Lemma 3, using any efficient classical algorithm for
approximate sampling from the output probability distribution of a Boson
Sampling instance with input cat states with
$0<|\alpha|=\mathcal{O}(n^{-1/4}\log^{1/4}n)$, one may also efficiently sample
approximately from the output probability distribution of the corresponding
Boson Sampling instance with input single photons efficiently, by keeping only
the samples $\bm{p}$ satisfying $|\bm{p}|=n$.

This proves that Boson Sampling with input cat states with
$0<|\alpha|=\mathcal{O}(n^{-1/4}\log^{1/4}n)$ is hard to sample approximately
in the same regime as Boson Sampling with input Fock states, assuming the same
conjectures as in [6]. ∎

## Acknowledgments

U. C. thanks Atul Singh Arora and Pierre-Emmanuel Emeriau for interesting
discussions. U. C, A. D, and S. M. acknowledge funding provided by the
Institute for Quantum Information and Matter, an NSF Physics Frontiers Center
(NSF Grant PHY-1733907). A. D. also acknowledges support from the National
Science Foundation RAISE-TAQS 1839204 and Amazon Web Services, AWS Quantum
Program.

## References

  * [1] H. Minc, “Permanents,”, vol. 6.  Cambridge University Press, 1984.
  * [2] J. K. Percus, “Combinatorial methods,”, vol. 4.  Springer Science & Business Media, 2012.
  * [3] L. G. Valiant, “The complexity of computing the permanent,” Theoretical computer science 8, 189–201 (1979).
  * [4] E. R. Caianiello, “On quantum field theory—I: explicit solution of Dyson’s equation in electrodynamics without use of Feynman graphs,” Il Nuovo Cimento (1943-1954) 10, 1634–1652 (1953).
  * [5] S. Scheel, “Permanents in linear optical networks,” quant-ph/0406127.
  * [6] S. Aaronson and A. Arkhipov, “The computational Complexity of Linear Optics,” Theory of Computing 9, 143 (2013), arXiv:1011.3245.
  * [7] S. Aaronson, “A linear-optical proof that the permanent is# P-hard,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 467, 3393–3405 (2011).
  * [8] S. Rahimi-Keshari, A. P. Lund, and T. C. Ralph, “What can quantum optics say about computational complexity theory?,” Physical review letters 114, 060501 (2015).
  * [9] D. Grier and L. Schaeffer, “New hardness results for the permanent using linear optics,” arXiv:1610.04670.
  * [10] P. P. Rohde, D. W. Berry, K. R. Motes, and J. P. Dowling, “A Quantum Optics Argument for the $\\#$P-hardness of a Class of Multidimensional Integrals,” arXiv:1607.04960.
  * [11] L. Chakhmakhchyan, N. J. Cerf, and R. Garcia-Patron, “Quantum-inspired algorithm for estimating the permanent of positive semidefinite matrices,” Physical Review A 96, 022329 (2017).
  * [12] A. Meiburg, “Inapproximability of Positive Semidefinite Permanents and Quantum State Tomography,” arXiv:2111.03142.
  * [13] P. A. MacMahon, “Combinatory Analysis, Volumes I and II,”, vol. 137.  American Mathematical Soc., 2001.
  * [14] I. Good, “Proofs of some ‘binomial’identities by means of MacMahon’s ‘Master Theorem’,” in Mathematical Proceedings of the Cambridge Philosophical Society, vol. 58, pp. 161–162, Cambridge University Press.  1962\.
  * [15] L. Carlitz, “An application of MacMahon’s master theorem,” SIAM Journal on Applied Mathematics 26, 431–436 (1974).
  * [16] L. Carlitz, “Some expansions and convolution formulas related to MacMahon’s master theorem,” SIAM Journal on Mathematical Analysis 8, 320–336 (1977).
  * [17] H. J. Ryser, “Combinatorial mathematics,”, vol. 14.  American Mathematical Soc., 1963.
  * [18] K. Balasubramanian, Combinatorics and diagonals of matrices.  PhD thesis, Indian Statistical Institute-Kolkata, 1980.
  * [19] E. T. Bax, Finite-difference algorithms for counting problems.  PhD thesis, California Institute of Technology, 1998.
  * [20] D. G. Glynn, “The permanent of a square matrix,” European Journal of Combinatorics 31, 1887–1891 (2010).
  * [21] P. P. Rohde, K. R. Motes, P. A. Knott, J. Fitzsimons, W. J. Munro, and J. P. Dowling, “Evidence for the conjecture that sampling generalized cat states with linear optics is hard,” Physical Review A 91, 012342 (2015).
  * [22] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, “Gaussian quantum information,” Reviews of Modern Physics 84, 621 (2012).
  * [23] A. Leverrier, “$SU(p,q)$ coherent states and a Gaussian de Finetti theorem,” Journal of Mathematical Physics 59, 042202 (2018).
  * [24] T. Jiang and Y. Ma, “Distances between random orthogonal matrices and independent normals,” arXiv:1704.05205.
  * [25] A. C. Dixon, “On the sum of the cubes of the coefficients in a certain expansion by the binomial theorem,” Messenger of mathematics 20, 79–80 (1891).
  * [26] I. Good, “A short proof of MacMahon’s ‘Master Theorem’,” in Mathematical Proceedings of the Cambridge Philosophical Society, vol. 58, pp. 160–160, Cambridge University Press.  1962\.
  * [27] S. Garoufalidis, T. T. Lê, and D. Zeilberger, “The quantum MacMahon master theorem,” Proceedings of the National Academy of Sciences 103, 13928–13931 (2006).
  * [28] M. Konvalinka and I. Pak, “Non-commutative extensions of the MacMahon Master Theorem,” Advances in Mathematics 216, 29–61 (2007).
  * [29] M. P. Tuite, “Some generalizations of the MacMahon Master Theorem,” Journal of Combinatorial Theory, Series A 120, 92–101 (2013).
  * [30] V. V. Kocharovsky, V. V. Kocharovsky, and S. V. Tarasov, “The Hafnian Master Theorem,” Linear Algebra and its Applications 144–161 (2022).
  * [31] W. Y. Chen, H. Galbraith, and J. Louck, “Angular momentum theory, umbral calculus, and combinatorics,” Computers & Mathematics with Applications 41, 1199–1214 (2001).
  * [32] B. M. Terhal and D. P. DiVincenzo, “Classical simulation of noninteracting-fermion quantum circuits,” Physical Review A 65, 032325 (2002).
  * [33] V. Shchesnovich, “Partial indistinguishability theory for multiphoton experiments in multiport devices,” Physical Review A 91, 013844 (2015).
  * [34] D. Spivak, M. Y. Niu, B. C. Sanders, and H. de Guise, “Generalized interference of fermions and bosons,” Physical Review Research 4, 023013 (2022).
  * [35] E.-J. Kuo, Y. Xu, D. Hangleiter, A. Grankin, and M. Hafezi, “Boson Sampling for Generalized Bosons,” arXiv:2204.08389.
  * [36] A. Clément, N. Heurtel, S. Mansfield, S. Perdrix, and B. Valiron, “LO${}_{\text{v}}$-Calculus: A Graphical Language for Linear Optical Quantum Circuits,” arXiv:2204.11787.
  * [37] G. De Felice and B. Coecke, “Quantum Linear Optics via String Diagrams,” arXiv:2204.12985.
  * [38] B. Peropadre, G. G. Guerreschi, J. Huh, and A. Aspuru-Guzik, “Proposal for microwave boson sampling,” Physical review letters 117, 140505 (2016).
  * [39] S. Girvin, “Schrödinger cat states in circuit qed,” arXiv:1710.03179.
  * [40] X. Gu, A. F. Kockum, A. Miranowicz, Y.-x. Liu, and F. Nori, “Microwave photonics with superconducting quantum circuits,” Physics Reports 718, 1–102 (2017).
  * [41] J. Huh, “A fast quantum algorithm for computing matrix permanent,” arXiv:2205.01328.
  * [42] S. Aaronson and T. Hance, “Generalizing and Derandomizing Gurvits’s Approximation Algorithm for the Permanent,” Quantum Info. Comput. 14, 541–559 (2014).
  * [43] S. Chin and J. Huh, “Generalized concurrence in boson sampling,” Scientific reports 8, 1–9 (2018).
  * [44] M.-H. Yung, X. Gao, and J. Huh, “Universal bound on sampling bosons in linear optics and its computational implications,” National science review 6, 719–729 (2019).
  * [45] V. S. Shchesnovich, “On the classical complexity of sampling from quantum interference of indistinguishable bosons,” International Journal of Quantum Information 18, 2050044 (2020).
  * [46] D. M. Jackson, “The unification of certain enumeration problems for sequences,” Journal of Combinatorial Theory, Series A 22, 92–96 (1977).
  * [47] F. R. Cardoso, D. Z. Rossatto, G. P. Fernandes, G. Higgins, and C. J. Villas-Boas, “Superposition of two-mode squeezed states for quantum information processing and quantum sensing,” Physical Review A 103, 062405 (2021).
  * [48] A. P. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L. O’Brien, and T. C. Ralph, “Boson sampling from a Gaussian state,” Physical review letters 113, 100502 (2014).
  * [49] J. P. Olson, K. P. Seshadreesan, K. R. Motes, P. P. Rohde, and J. P. Dowling, “Sampling arbitrary photon-added or photon-subtracted squeezed states is in the same complexity class as boson sampling,” Physical Review A 91, 022317 (2015).
  * [50] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, “Gaussian boson sampling,” Physical review letters 119, 170501 (2017).
  * [51] A. Lund, S. Rahimi-Keshari, and T. Ralph, “Exact boson sampling using Gaussian continuous-variable measurements,” Physical Review A 96, 022301 (2017).
  * [52] L. Chakhmakhchyan and N. J. Cerf, “Boson sampling with Gaussian measurements,” Physical Review A 96, 032326 (2017).
  * [53] U. Chabaud, T. Douce, D. Markham, P. van Loock, E. Kashefi, and G. Ferrini, “Continuous-variable sampling from photon-added or photon-subtracted squeezed states,” Physical Review A 96, 062307 (2017).
  * [54] N. Quesada, J. M. Arrazola, and N. Killoran, “Gaussian boson sampling using threshold detectors,” Physical Review A 98, 062322 (2018).
  * [55] A. Deshpande, A. Mehta, T. Vincent, N. Quesada, M. Hinsche, M. Ioannou, L. Madsen, J. Lavoie, H. Qi, J. Eisert, et al., “Quantum computational advantage via high-dimensional Gaussian boson sampling,” Science advances 8, 7894 (2022).