Current affiliation: ]PsiQuantum, Palo Alto, CA, United States

# Threshold detection statistics of bosonic states

J. F. F. Bulmer<EMAIL_ADDRESS>[  Quantum Engineering
Technology Laboratories, University of Bristol, Bristol, BS8 1FD, UK  S.
Paesani  Center for Hybrid Quantum Networks (Hy-Q), Niels Bohr Institute,
University of Copenhagen, Blegdamsvej 17, DK-2100 Copenhagen, Denmark  R. S.
Chadwick  Quantum Engineering Technology Laboratories, University of Bristol,
Bristol, BS8 1FD, UK  Quantum Engineering Centre for Doctoral Training,
University of Bristol, Bristol BS8 1FD, UK  N. Quesada  Department of
Engineering Physics, École Polytechnique de Montréal, Montréal, QC, H3T 1JK,
Canada

###### Abstract

In quantum photonics, threshold detectors, distinguishing between vacuum and
one or more photons, such as superconducting nanowires and avalanche
photodiodes, are routinely used to measure Fock and Gaussian states of light.
Despite being the standard measurement scheme, there is no general closed form
expression for measurement probabilities with threshold detectors, unless
accepting coarse approximations or combinatorially scaling summations. Here,
we present new matrix functions to fill this gap. We develop the Bristolian
and the loop Torontonian functions for threshold detection of Fock and
displaced Gaussian states, respectively, and connect them to each other and to
existing matrix functions. By providing a unified picture of bosonic
statistics for most quantum states of light, we provide novel tools for the
design and analysis of photonic quantum technologies.

##  I Introduction

Quantum photonic experiments can generally be described as preparing quantum
states of light, evolving them through linear optical interferometers, and
detecting the output photons. While the most common types of photonic states,
Fock states and Gaussian states, can be routinely prepared via spontaneous
processes in optical non-linearities or (artificial) atomic systems and
processed with high-fidelity linear optical components [1], photon number
detection is typically approximated via the use of threshold photon detectors.
Threshold detection, i.e. a measurement distinguishing only between vacuum and
the presence of one or more photons, is widely available, e.g. via high-
efficiency superconducting nanowires [2] or room-temperature avalanche
photodiodes [3], making it the standard measurement apparatus in quantum
photonics. Its use in experiments, represented in Fig. 1, encompasses many
areas of quantum research, including demonstrations of quantum advantages,
e.g. in computation [4], measurement sensitivity [5], and loophole-free tests
of non-locality [6]. However, threshold detection only provides a meaningful
approximation of the desired Fock basis measurement projectors in the regime
of low mean photon numbers per mode. On the other hand, as technology
progresses, mean photon numbers increase and higher fidelities are demanded
[4], making this approximation less appropriate.

Figure 1: Types of typical quantum photonic experiments which can be modelled
using the results of this work. Threshold detection statistics of Gaussian
states (highlighted in red) formed by squeezing, displacement and linear
optics are captured by the loop Torontonian matrix function. Displacement can
be generated using coherent states from lasers, squeezing can come from
nonlinear processes such as Spontaneous Parametric Down Conversion (SPDC) or
Spontaneous Four-Wave Mixing (SFWM). Linear optics can be implemented in a
variety of platforms, including bulk and integrated optics. States created by
the linear optical interference of Fock states (highlighted in blue), as
generated by (artificial) atoms, lead to threshold detection probabilities
given by the Bristolian matrix function.

To circumvent this issue, experiments can be described directly using the
output statistics of threshold detection instead of its photon number
resolving approximation. However, despite the wide adoption of such systems,
there are in general no closed form expressions in the literature for
computing measurement probabilities of threshold detectors. In fact, while an
expression for the threshold detection of zero-displaced Gaussian states is
known, given by the Torontonian matrix function [7], no analogous expressions
exist for other commonly used states, e.g. Fock or displaced Gaussian states.
For example, for Fock states with fixed photon number, threshold probabilities
could be exactly calculated by summing over all possible output states which
lead to the given threshold detector outcome. However, this method requires
calculating a number of probabilities scaling combinatorially with the number
of clicked detectors, rendering it impractical already for smaller-scale
experiments [8, 9, 10]. New methods are required to describe quantum photonic
technologies which use threshold detection.

Here, we provide such methods by developing a unified picture to compute
threshold statistics for most quantum photonic states of experimental
interest. As described in Table 1, this is achieved by introducing two new
matrix functions, the Bristolian and the loop Torontonian, for threshold
statistics with Fock and displaced Gaussian states, respectively, and
demonstrating close connections between them and to other existing matrix
functions. The developed tools provide exact simulation, design, and analysis
methods for current [11, 8, 9, 12, 4, 10] and future quantum photonic systems
that use threshold detection. We wish to highlight the different challenges
between computing probabilities, known as strong simulation, which we focus on
in this work, and drawing samples from a probability distribution, known as
weak simulation [13]. These tasks often have very different complexity. For
example, using methods from Ref. [14], we can sample threshold detector
outcomes without ever calculating a threshold detection probability.

State | Detector
---|---
number resolving | threshold
Fock | permanent | Bristolian*
zero-mean Gaussian | Hafnian | Torontonian
displaced Gaussian | loop Hafnian | loop Torontonian*
Table 1: Matrix functions for the calculation of measurement probabilities in
quantum photonics. (*)-symbol is used to indicate functions which are
introduced in this work.

##  II Threshold detection statistics from vacuum statistics

Threshold detectors are described by the measurement operators

$\displaystyle\hat{\Pi}_{j}^{(0)}$ $\displaystyle=\ket{0_{j}}\\!\bra{0_{j}},$
(1a) $\displaystyle\hat{\Pi}_{j}^{(1)}$
$\displaystyle=\sum_{k=1}^{\infty}\ket{k_{j}}\\!\bra{k_{j}}=\mathbb{I}-\ket{0_{j}}\\!\bra{0_{j}},$
(1b)

for vacuum (0) and click (1) outcomes on a mode described by label $j$. We use
$\ket{0}$ to denote the vacuum state of an optical mode,
$\ket{k}=(\hat{a}^{\dagger})^{k}\ket{0}/\sqrt{k!}$ for Fock states of the
optical mode, and $\mathbb{I}$ is the identity operator (we will always assume
its dimension to be the same as the other operators appearing in the
equation).

We write the outcome of $M$ threshold detectors, labelled with
$j\in[M]=\\{1,2,\dots,M\\}$, using a length-$M$ bit-string $\vec{d}$, where
the $j$th element gives the measurement outcome of the $j$th mode. Defining a
set of modes which clicked, $C=\\{j\in[M]\ |\ d_{j}=1\\}$, and a set for modes
with the vacuum outcome, $V=\\{j\in[M]\ |\ d_{j}=0\\}$, we can write the
multimode measurement operator as

$\displaystyle\hat{\Pi}^{(\vec{d})}$
$\displaystyle=\bigotimes_{j=1}^{M}\hat{\Pi}^{(d_{j})}_{j}=\bigotimes_{j\in
C}\left(\mathbb{I}-\ket{0_{j}}\\!\bra{0_{j}}\right)\bigotimes_{k\in
V}\ket{0_{k}}\\!\bra{0_{k}},$ (2)

which can be rearranged to give

$\displaystyle\hat{\Pi}^{(\vec{d})}=\ket{\vec{0}_{V}}\\!\bra{\vec{0}_{V}}\sum_{Z\in
P(C)}(-1)^{|Z|}\ket{\vec{0}_{Z}}\\!\bra{\vec{0}_{Z}}.$ (3)

Here, we use $P(C)$ to denote the powerset of $C$ and $|Z|$ for the number of
elements in a set $Z$. $\ket{\vec{0}_{V}}\\!\bra{\vec{0}_{V}}$ describes the
vacuum projector in all the vacuum outcome modes and
$\ket{\vec{0}_{Z}}\\!\bra{\vec{0}_{Z}}$ describes the vacuum projector in all
modes in a subset $Z\subseteq C$.

Eq. (3) indicates that to calculate the threshold detection probabilities for
any state it is sufficient to calculate marginal vacuum probabilities, which
are used in an inclusion/exclusion sum as described by Eq. (3). Using this
measurement operator and the Born rule
$p(\vec{d})=\operatorname{tr}(\hat{\Pi}^{(\vec{d})}\rho)$ on some state
$\rho$, we find:

$\displaystyle p(\vec{d})=\sum_{Z\in
P(C)}(-1)^{|Z|}p(\vec{d}_{V}=\vec{0},\vec{d}_{Z}=\vec{0}),$ (4)

where
$p(\vec{d}_{V}=\vec{0},\vec{d}_{Z}=\vec{0})=\operatorname{tr}(\ket{\vec{0}_{V}}\\!\bra{\vec{0}_{V}}\otimes\ket{\vec{0}_{Z}}\\!\bra{\vec{0}_{Z}}\rho)$.
This formula provides our starting point for deriving general expressions for
threshold detection statistics.

##  III Marginal vacuum probabilities from the photon number probability
generating function

For an $M$-mode linear optical interferometer, described by an $M\times M$
matrix $U$ and the operator $\hat{\mathcal{U}}$, the creation operators are
transformed as

$\displaystyle\hat{\mathcal{U}}\hat{a}^{\dagger}_{j}\hat{\mathcal{U}}^{\dagger}=\sum_{k=1}^{M}U_{kj}\hat{a}^{\dagger}_{k}.$
(5)

Considering an input state $\ket{\Phi_{0}}$, the output photon number
probability distribution is then

$\displaystyle
p(\vec{m})=\left|\bra{\vec{m}}\hat{\mathcal{U}}\ket{\Phi_{0}}\right|^{2}$ (6)

where $\vec{m}$ is a length-$M$ list describing the photon number in each mode
at the output of the interferometer, and
$\ket{\vec{m}}=\bigotimes_{j=1}^{M}\left((\hat{a}^{\dagger}_{j})^{m_{j}}/\sqrt{m_{j}!}\right)\ket{0}$.

Following Ref. [15], considering the Fourier transform of the probability
distribution of photon number basis measurements we define the characteristic
function

$\displaystyle\chi(\vec{\phi})=\sum_{\vec{m}}\exp\left(\mathrm{i}\sum_{j=1}^{M}\phi_{j}m_{j}\right)p(\vec{m})$
(7)

which, with some manipulation (see Appendix B), can be expressed as

$\displaystyle\chi(\vec{\phi})=\bra{\Phi_{0}}\hat{\mathcal{U}}^{\dagger}\
\hat{\mathcal{U}}_{\vec{\phi}}\ \hat{\mathcal{U}}\ket{\Phi_{0}}$ (8)

where $\hat{\mathcal{U}}_{\vec{\phi}}$ is the operator given by the evolution
due to the linear optical transformation

$U_{\vec{\phi}}=\bigoplus_{j=1}^{M}\exp(\mathrm{i}\phi_{j}).$ (9)

We can transform this into a probability generating function, $G$, using the
substitution $x_{j}=\exp(\mathrm{i}\phi_{j})$:

$\displaystyle
G(\vec{x})=\sum_{\vec{m}}\left(\prod_{j=1}^{M}x_{j}^{m_{j}}\right)p(\vec{m}).$
(10)

The function $G(\vec{x})$ has the following useful properties. To marginalise
the $j$th mode, we simply set $x_{j}=1$. If we set $x_{j}=0$, this gives us
the probability for $n_{j}=0$. Therefore, if we want to calculate the
probability that some subset of the modes, $V$, measure vacuum and we
marginalise over all other modes, $B$, we can evaluate

$\displaystyle
p(\vec{m}_{V}=\vec{0})=G(\vec{x}_{V}=\vec{0},\vec{x}_{B}=\vec{1}),$ (11)

where $\vec{0}$ ($\vec{1}$) is a vector with 0 (1) in all entries.
$G(\vec{x})$ is the probability distribution generating function, so by taking
derivatives of $G(\vec{x})$, we can find information about the photon number
basis probability distribution [15].

By using the expression for the characteristic function in Eq. (8), we can see
that this amounts to calculating the scattering amplitude of $\ket{\Phi_{0}}$
to itself, through a linear optical interferometer described by the
transformation $U^{\dagger}U_{\vec{\phi}}U$. From again using the substitution
$x_{j}=\exp(\mathrm{i}\phi_{j})$ in Eq. (9), we see that $U_{\vec{\phi}}$
physically corresponds to either zero transmission for modes in $V$, or unit
transmission for modes in $B$. As we show in the next sections, we can use
this, in conjunction with Eq. (4) to calculate threshold detection
probabilities for all the experimental scenarios outlined in Fig. 1.

##  IV Fock state inputs

Recall that the scattering amplitudes of Fock states evolved through a
lossless interferometer are given by the permanent matrix function [16]

$\displaystyle\bra{\vec{m}}\hat{\mathcal{U}}\ket{\vec{n}}=\frac{\operatorname{per}(U_{\vec{m},\vec{n}})}{\sqrt{\prod_{j=1}^{M}n_{j}!m_{j}!}}$
(12)

where $U_{\vec{m},\vec{n}}$ is constructed from $U$ by repeating its $j$th row
$m_{j}$ times and its $j$th column $n_{j}$ times for all $j\in[M]$.

Therefore if we have an $N$-photon input Fock state,
$\ket{\Phi_{0}}=\ket{\vec{n}}$, we can use Eq. (8), Eq. (10) and Eq. (12) to
write

$\displaystyle
G(\vec{x})=\frac{\operatorname{per}\left([U^{\dagger}U_{\vec{x}}U]_{\vec{n},\vec{n}}\right)}{\prod_{j=1}^{M}n_{j}!}$
(13)

where $U_{\vec{x}}$ is formed like $U_{\vec{\phi}}$, but with diagonal matrix
elements: $[U_{\vec{x}}]_{jj}=x_{j}$. For this equation to be valid, we must
have a lossless unitary transformation. However, we are free to marginalise
over modes by allowing elements of $\vec{x}$ to be set to 1 for any mode we
wish to marginalise over, including any loss modes, as shown in Eq. (11).

Because Eq. (13) provides us with a closed form expression for marginal vacuum
probabilities, and because Eq. (4) shows us that marginal vacuum probabilities
are sufficient to calculate threshold detection probabilities, we can use this
to derive a matrix function for calculating threshold detection probabilities
of Fock states. For more generality, we first consider a linear optical
transformation with losses, described by an $M_{\text{in}}\times
M_{\text{out}}$ matrix $T$, with singular values upper bounded by $1$ [17]. In
Appendix C we show that if the input state is an $M_{\text{in}}$-mode Fock
state, $\vec{n}$, then by combining Eq. (4) with Eq. (13), we can calculate
the threshold detection probability of the outcome described by an
$M_{\text{out}}$-length bit-string $\vec{d}$ as:

$\displaystyle
p(\vec{d})=\frac{\operatorname{brs}\left(T_{\vec{d},\vec{n}},E(T)_{\vec{n},\vec{n}}\right)}{\prod_{j=1}^{M_{\text{in
}}}n_{j}!}.$ (14)

Here we have introduced a matrix function, the Bristolian, defined as

$\operatorname{brs}\left(A,E\right)=\\\ \sum_{Y\in
P([m])}(-1)^{m-|Y|}\operatorname{per}\left([A_{Y}]^{\dagger}A_{Y}+E\right),$
(15)

where $A$ is an $m\times n$ matrix and $E$ is an $n\times n$ matrix. $A_{Y}$
denotes selecting the rows of $A$ according to the elements of $Y$, and
$[m]=\\{1,2,\dots,m\\}$. We have also defined a matrix which accounts for the
mixing with vacuum in the environment modes

$\displaystyle E(T)$ $\displaystyle=\mathbb{I}-T^{\dagger}T.$ (16)

Our naming of the Bristolian is inspired by the convention established by the
Hafnian and Torontonian matrix functions, which are named after the cities of
their discovery. By noticing that $\mathbb{I}-T^{\dagger}T$ gives a zero
matrix when $T$ is unitary, the $\operatorname{brs}$ function can be
simplified when $T$ is a unitary matrix $U$, only requiring the rows of $U$
which correspond to modes with a detector click, providing

$\displaystyle
p(\vec{d})=\frac{\operatorname{ubrs}\left(U_{\vec{d},\vec{n}}\right)}{\prod_{j=1}^{M}n_{j}!}.$
(17)

Here we defined the Unitary Bristolian acting on an $m\times m$ matrix, $A$,
as

$\displaystyle\operatorname{ubrs}(A)=\sum_{Y\in
P([m])}(-1)^{m-|Y|}\operatorname{per}\left([A_{Y}]^{\dagger}A_{Y}\right).$
(18)

##  V Displaced Gaussian state inputs

Gaussian states are the set of states that have a Gaussian characteristic
function. A Gaussian state $\rho$ is uniquely characterized by its vector of
means with entries

$\displaystyle\vec{\alpha}_{i}=\text{tr}\left[\rho\hat{\vec{\zeta}}_{i}\right],$
(19)

and its Husimi covariance matrix with entries

$\displaystyle\Sigma_{i,j}=\tfrac{1}{2}\text{tr}\left(\left[\hat{\vec{\zeta}}_{i}\hat{\vec{\zeta}}_{j}^{\dagger}+\hat{\vec{\zeta}}_{j}^{\dagger}\hat{\vec{\zeta}}_{i}\right]\rho\right)-\vec{\alpha}_{i}\vec{\alpha}_{j}^{*}+\tfrac{1}{2}\delta_{i,j},$
(20)

where we have used a vector of creation and annihilation operators

$\displaystyle\hat{\vec{\zeta}}=\left(\hat{a}_{1},\ldots,\hat{a}_{M},\hat{a}_{1}^{\dagger},\ldots,\hat{a}_{M}^{\dagger}\right).$
(21)

The Husimi function $Q(\vec{r})=\bra{\vec{r}}\rho\ket{\vec{r}}$ maps
displacement vectors, $\vec{r}$, to probabilities, so to calculate vacuum
probabilities we can evaluate the Husimi function at the origin. Noting that
we can marginalise over modes by deleting all the corresponding elements of
$\Sigma$ and $\vec{\alpha}$, we obtain [18]

$\displaystyle p(\vec{m}_{V}=\vec{0})$ $\displaystyle=Q(\vec{r}_{V}=\vec{0})$
(22)
$\displaystyle=\frac{\exp\left(-\frac{1}{2}\vec{\alpha}_{V}^{\dagger}[\Sigma_{VV}]^{-1}\vec{\alpha}_{V}\right)}{\sqrt{\det(\Sigma_{VV})}}.$
(23)

The notation $\Sigma_{VV}$ and $\vec{\alpha}_{V}$ differs slightly here from
the previous section, as now there are two basis vectors for each mode of our
system, corresponding to each mode’s $\hat{a}$ and $\hat{a}^{\dagger}$
operator. We form $\Sigma_{VV}$ by selecting both rows/columns of $\Sigma$
which correspond to each element of $V$ and we form $\vec{\alpha}_{V}$ by
selecting both elements of $\vec{\alpha}$ corresponding to each element of
$V$.

We can use this to immediately arrive at a threshold detection probability for
displaced Gaussian states using Eq. (4). However, here we must invert and
compute determinants for square matrices of size $2(|V|+|Z|)$. It would be
preferable if we could reduce these to matrices of size $2|Z|$. It would also
be helpful conceptually to have a formula which can be connected to other
relevant matrix functions, the Torontonian [7] and the loop Hafnian [19].
Therefore, it is of interest to write this probability in terms of

$\displaystyle O=\mathbb{I}-\Sigma^{-1}\text{ and
}\vec{\gamma}=(\Sigma^{-1}\vec{\alpha})^{*}.$ (24)

In Appendix F, we show how Eq. (4) can be rearranged into the following:

$\displaystyle
p(\vec{d})=p(\vec{0})\operatorname{ltor}\left(O_{CC},\vec{\gamma}_{C}\right),$
(25)

where $C$ is given by the index of the elements of $\vec{d}$ where $d_{j}=1$,
so $O_{CC}$ and $\vec{\gamma}_{C}$ are the matrix and vector formed by
selecting the rows/columns of $O$ and elements of $\vec{\gamma}$ which
correspond to modes which see a detector click. $p(\vec{0})$ is the
probability of detecting vacuum in all modes, and can be calculated using Eq.
(23). We introduce the loop Torontonian, which is defined as

$\operatorname{ltor}\left(O,\vec{\gamma}\right)=\\\ \sum_{Y\in
P([m])}(-1)^{m-|Y|}\frac{\exp\left[\tfrac{1}{2}\vec{\gamma}_{Y}^{t}[\mathbb{I}-O_{YY}]^{-1}\vec{\gamma}_{Y}^{*}\right]}{\sqrt{\det(\mathbb{I}-O_{YY})}},$
(26)

where $O$ is a $2m\times 2m$ matrix and $\vec{\gamma}$ is a $2m$-length
vector.

##  VI Connections between matrix functions

In the limit of no displacement $\vec{\alpha}=\vec{\gamma}=\vec{0}$, the
exponential terms in the numerator of Eq. (26) becomes 1 and thus
$\operatorname{ltor}\left(O,\vec{0}\right)=\operatorname{tor}\left(O\right)$,
where $\operatorname{tor}$ is the Torontonian function from Ref. [7]. One can
show, using the scattershot construction [20], that the Torontonian and
Bristolian are related via the following limit

$\operatorname{brs}\left(T_{\vec{d},\vec{n}},E(T)_{\vec{n},\vec{n}}\right)=\\\
\lim_{\varepsilon\to
0}(\varepsilon^{-2}-1)^{N}\operatorname{tor}\left(O(\varepsilon)_{CC}\right),$
(27) $O(\varepsilon)=\varepsilon\begin{pmatrix}0&0&0&T\\\ 0&\varepsilon
E(T)^{*}&T^{t}&0\\\ 0&T^{*}&0&0\\\ T^{\dagger}&0&0&\varepsilon
E(T)\end{pmatrix}.$ (28)

where $\vec{n}$ is a bitstring (implying that this identity is only valid for
single-photon or vacuum inputs), $N=\sum_{i}\vec{n}_{i}$ and $C$ is the union
of the labels of the modes in which single photons were input into the
interferometer and the labels of the modes in which clicks are registered.
This relation is proven in Appendix H.

As we show in Appendix G, the loop Torontonian can also be used as a
generating function for the loop Hafnian,

$\operatorname{lhaf}(XO_{CC},\vec{\gamma}_{C})=\left.\frac{1}{\ell!}\frac{d^{\ell}}{d\eta^{\ell}}\operatorname{ltor}\left(\eta
O_{CC},\sqrt{\eta}\vec{\gamma}_{C}\right)\right|_{\eta=0}$ (29)

where $X=\left[\begin{smallmatrix}0&\mathbb{I}\\\
\mathbb{I}&0\end{smallmatrix}\right]$ and $\ell=|C|$. We use this to derive
the trace formula for the loop Hafnian, the fastest known method for computing
photon number resolved measurement probabilities on displaced Gaussian states
[21, 19, 22]. Because the loop Hafnian of a bipartite graph is given by the
matrix permanent [21], all the matrix functions in Table 1 can be derived from
the loop Torontonian.

We also see a connection between the Bristolian and the permanent when an $N$
photon Fock state results in $N$ threshold detector clicks. In this case, each
threshold detector must have seen exactly 1 photon, so we can describe the
measurement operator of each threshold detector click as a single photon
projector, which leads to describing the event with permanents, as given by
Eq. (12). In Appendix D, we show this link directly by first describing the
Unitary Bristolian for $N$ photon, $N$ click events as the permanent of an
$N\times N\times N$ 3-tensor [23].

##  VII Time complexities

In Appendix I, we discuss the time complexities for the Bristolian and the
loop Torontonian. We find that, using the formulae presented in this work, the
Bristolian, $\operatorname{brs}(A,E)$, has a time complexity of
$\mathcal{O}(n2^{n+m})$ for an $m\times n$ matrix $A$ and $n\times n$ matrix
$E$ and the loop Torontonian, $\operatorname{ltor}(O,\vec{\gamma})$, has time
complexity of $\mathcal{O}(m^{3}2^{m})$ for a $2m\times 2m$ matrix $O$ and
$2m$-length vector $\vec{\gamma}$. For the loop Torontonian, this complexity
can be reduced using a recursive strategy which exploits Cholesky
decomposition [24]. We also believe that the Bristolian’s time complexity can
likely be reduced, and we leave this as an open problem.

##  VIII Improved accuracy of a threshold detection model

To assess the improvements offered by using the correct description of
threshold detection over the common approximation of single photon projective
measurement, we present two representative examples. By simulating the
probability distribution for 100 different Haar random unitaries in lossy 4
input photon Fock state Boson sampling experiments on mode numbers from 4 to
12, we evaluate the total variation distance (TVD) between probability
distributions from the exact model, which uses the Bristolian, and an
approximate model, which uses a sum over matrix permanents, as discussed in
Appendix E. Although the TVD is reduced for higher numbers of modes, the
approximation is always 5% - 12% removed from the correct distribution. To
test the loop Torontonian, we use experimental data from Ref. [10]. We see
that for the 2 photon distribution for different levels of displacement, the
loop Torontonian offers a better match to the experiment of up to 16%. See
Appendix J for more detail.

##  IX Conclusion

The new methods we have derived, in particular the Bristolian and the loop
Torontonian functions, are useful tools to model and analyse a wide variety of
quantum photonic experiments and applications. For example, the Bristolian is
relevant to applications including linear-optical quantum computing [25, 26,
27], Boson Sampling [28, 9] and quantum communications [29], commonly based on
threshold detection. The loop Torontonian can be applied to applications
including Gaussian state reconstruction [10], measuring graph similarity [30],
calculations of vibronic spectra of molecules [31], and quantum metrology
[32], and has already been applied for evaluating proposed quantum
communication protocols 111The initial inspiration for us to derive the loop
Torontonian came from the need to calculate threshold detection statistics for
the quantum communication protocols proposed in Ref. [47]. To facilitate their
use, we provide example calculations of common experimental scenarios in
Appendix A using the Bristolian and loop Torontonian, and have made available
implementations in the open-source Python package The Walrus [34]. Details for
the software implementation are provided in Appendix K. The connections that
we have shown between the Bristolian and the permanent (Appendix D), the loop
Torontonian and the loop Hafnian (Appendix G), and the Bristolian and the
Torontonian (Appendix H) indicate that these functions can provide a useful
mathematical and conceptual tool for a deeper understanding of bosonic
statistics in photonic experiments.

## Acknowledgements

JFFB and RSC acknowledge support from EPSRC (EP/N509711/1, EP/LO15730/1). NQ
acknowledges support from the Ministère de l’Économie et de l’Innovation du
Québec and the Natural Sciences and Engineering Research Council of Canada. SP
acknowledges funding from the Cisco University Research Program Fund nr.
2021-234494. We thank G. S. Thekkadath for useful discussions and sharing
experimental data from Ref. [10]. NQ thanks S. Duque Mesa, B. Lanthier, D.
Leclerc, B. Turcotte, and J. Zhao for valuable discussions. We thank G. Morse
for implementing the generalisation of the recursive Torontonian formula [24]
to the loop Torontonian, see pull request (332) to The Walrus [34].

##  Appendix A Examples

We present some examples of how to apply the Unitary Bristolian, Bristolian
and loop Torontonian to some representative situations.

###  A.1 Lossless Hong-Ou-Mandel

Our first example is Hong-Ou-Mandel interference [35] of single photons on a
lossless 50/50 beam splitter. Here, the input state is $\vec{n}=(1,1)$ and

$\displaystyle U=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\ 1&-1\end{pmatrix}.$
(30)

The probability of detecting a coincidence is given by

$\displaystyle p(\vec{d}=(1,1))=$ $\displaystyle\operatorname{ubrs}(U).$ (31)

So, $C=\\{1,2\\}$. We expand the $\operatorname{ubrs}$ function as

$\displaystyle\operatorname{ubrs}(U)=$
$\displaystyle\operatorname{per}\left(U^{\dagger}U\right)$ (32)
$\displaystyle-\operatorname{per}\left([U_{1}]^{\dagger}U_{1}\right)-\operatorname{per}\left([U_{2}]^{\dagger}U_{2}\right),$
(33)

with $U_{1}=(1,1)/\sqrt{2}$ and $U_{2}=(1,-1)/\sqrt{2}$ being the 1st and 2nd
rows of $U$ respectively. We also use that
$\operatorname{per}\left(U^{\dagger}U\right)=\operatorname{per}(\mathbb{I})=1$,
and calculate

$\displaystyle\operatorname{per}\left([U_{1}]^{\dagger}U_{1}\right)$
$\displaystyle=\operatorname{per}\left[\frac{1}{2}\begin{pmatrix}1&1\\\
1&1\end{pmatrix}\right]$ $\displaystyle=\frac{1}{2},$ (34)
$\displaystyle\operatorname{per}\left([U_{2}]^{\dagger}U_{2}\right)$
$\displaystyle=\operatorname{per}\left[\frac{1}{2}\begin{pmatrix}1&-1\\\
-1&1\end{pmatrix}\right]$ $\displaystyle=\frac{1}{2}$ , (35)

giving $p(\vec{d}=(1,1))=0$, as expected. Here, we have not included the term
where $Y$ is the empty set. In this case, we are considering the permanent of
the all zeros matrix, which is zero and so does not contribute.

We can also see that
$p(\vec{d}=(1,0))=\operatorname{per}\left([U_{1}]^{\dagger}U_{1}\right)=1/2$,
and similarly $p(\vec{d}=(0,1))=1/2$. These are due to the $(2,0)$ and $(0,2)$
photon number output terms.

###  A.2 3-mode zero transmission law

The zero transmission law (ZTL) tells us that many output states of a Fourier
transform interferometer are suppressed due to multi-photon interference.

Using $\omega=\exp(-2i\pi/3)$, the three mode Fourier transform interferometer
is given by

$\displaystyle U=\frac{1}{\sqrt{3}}\begin{pmatrix}1&1&1\\\
1&\omega&\omega^{2}\\\ 1&\omega^{2}&\omega\end{pmatrix}.$ (36)

We consider $\vec{n}=(1,1,1)$. According to the ZTL, all permutations of the
output $\vec{m}=(2,1,0)$ should be suppressed, and so $\vec{d}=(1,1,0)$, due
to $\vec{m}=(2,1,0)$ and $\vec{m}=(1,2,0)$, is also expected to be suppressed.
The probability of this threshold detector outcome is given by

$\displaystyle p(\vec{d}=(1,1,0))=$
$\displaystyle\operatorname{ubrs}\left(U_{\vec{d},\vec{n}}\right),$ (37)

where

$\displaystyle\operatorname{ubrs}\left(U_{\vec{d},\vec{n}}\right)=$
$\displaystyle\operatorname{per}\left([U_{(1,2)}]^{\dagger}U_{(1,2)}\right)$
(38)
$\displaystyle-\operatorname{per}\left([U_{1}]^{\dagger}U_{1}\right)-\operatorname{per}\left([U_{2}]^{\dagger}U_{2}\right).$

By evaluating these permanents, we find

$\displaystyle\operatorname{ubrs}\left(U_{\vec{d},\vec{n}}\right)=\frac{4}{9}-\frac{2}{9}-\frac{2}{9}=0,$
(39)

showing a suppression as expected. We can also calculate

$\displaystyle p(\vec{d}=(1,1,1))=$
$\displaystyle\operatorname{ubrs}(U)=\frac{1}{3},$ (40)

which agrees with the prediction that
$\operatorname{ubrs}(U)=\left|\operatorname{per}(U)\right|^{2}$ when $U$ is
square, since $\operatorname{per}(U)=-1/\sqrt{3}$.

###  A.3 Lossy Hong-Ou-Mandel

The Hong-Ou-Mandel effect is preserved under balanced loss. In contrast to the
previous examples, in this example example, the loss means that we need to use
the Bristolian, as the Unitary Bristolian is no longer valid. We consider a
transmission matrix defined like $U$ in Eq. (30), but with transmission
$\eta$, giving $T=\sqrt{\eta}U$. In this case

$\displaystyle p(\vec{d}=(1,1))=$ $\displaystyle\operatorname{brs}(T,E)$ (41)

where $E=(1-\eta)\mathbb{I}$. We evaluate the Bristolian to find

$\displaystyle\operatorname{brs}(T,E)=$
$\displaystyle\operatorname{per}\left(T^{\dagger}T+E\right)$ (42)
$\displaystyle-\operatorname{per}\left([T_{1}]^{\dagger}T_{1}+E\right)-\operatorname{per}\left([T_{2}]^{\dagger}T_{2}+E\right)$
$\displaystyle+\operatorname{per}(E)$ $\displaystyle=$ $\displaystyle\
1-2\left(1-\eta+\eta^{2}/2\right)+(1-2\eta+\eta^{2})$ $\displaystyle=$
$\displaystyle\ 0$

which confirms that the coincidence event is still suppressed under balanced
loss.

###  A.4 Lossy zero transmission law

When there is loss, we start to witness threshold detector outcomes which were
suppressed in the lossless case. We repeat the example in section A.2 but
adding a transmission $\eta\leq 1$. So, $T=\sqrt{\eta}U$ for $U$ defined in
Eq. (36).

Using the Bristolian, we find

$\displaystyle
p(\vec{d}=(1,1,0))=\operatorname{brs}(T_{\vec{d},\vec{n}},E(T)_{\vec{n},\vec{n}}),$
(43)

where $E(T)=(1-\eta)\mathbb{I}$, and this gives

$\displaystyle\operatorname{brs}(T_{\vec{d},\vec{n}},E(T)_{\vec{n},\vec{n}})=$
$\displaystyle\operatorname{per}\left([T_{(1,2)}]^{\dagger}T_{(1,2)}+E(T)\right)$
(44)
$\displaystyle-\operatorname{per}\left([T_{1}]^{\dagger}T_{1}+E(T)_{\vec{n},\vec{n}}\right)$
$\displaystyle-\operatorname{per}\left([T_{2}]^{\dagger}T_{2}+E(T)_{\vec{n},\vec{n}}\right)$
$\displaystyle+\operatorname{per}(E(T)_{\vec{n},\vec{n}})$ $\displaystyle=$
$\displaystyle\ \eta^{2}(1-\eta)/3.$ (45)

We can check this problem intuitively by considering all the losses to be
applied just before the measurement. In this picture, the outcome
$\vec{d}=(1,1,0)$ can only occur when we have the state $\vec{m}=(1,1,1)$
before the losses, then the first two photons are transmitted and the last
photon is lost. Therefore we expect the probability found in Eq. (45). The
factor $1/3$ comes from section A.2, where $p(\vec{m}=(1,1,1))=1/3$ before any
loss is applied.

For our final Bristolian example calculation, we consider a case where we have
to repeat columns more than once. Consider the same $T$ as above, but with the
input state $\vec{n}=(1,2,0)$. We calculate the probability of
$\vec{d}=(0,1,1)$ (so $C=\\{2,3\\}$), denoting the elements of $T$ as
$t_{jk}$.

$\displaystyle
p(\vec{d}=(0,1,1))=\frac{1}{2}\operatorname{brs}\left(T_{\vec{d},\vec{n}},E(T)_{\vec{n},\vec{n}}\right),$
(46)

where,

$\displaystyle T_{\vec{d},\vec{n}}=$
$\displaystyle\begin{pmatrix}t_{22}&t_{23}&t_{23}\\\
t_{32}&t_{33}&t_{33}\end{pmatrix}$ (47) $\displaystyle=$ $\displaystyle\
\sqrt{\frac{\eta}{3}}\begin{pmatrix}\omega&\omega^{2}&\omega^{2}\\\
\omega^{2}&\omega&\omega\end{pmatrix}$ $\displaystyle E(T)=$ $\displaystyle\
(1-\eta)\begin{pmatrix}1&0&0\\\ 0&1&0\\\ 0&0&1\end{pmatrix}$ (48)
$\displaystyle E(T)_{\vec{n},\vec{n}}=$ $\displaystyle\
(1-\eta)\begin{pmatrix}1&0&0\\\ 0&1&1\\\ 0&1&1\end{pmatrix},$ (49)

and finally,

$\displaystyle\operatorname{brs}\left(T_{\vec{d},\vec{n}},E(T)_{\vec{n},\vec{n}}\right)=$
(50)
$\displaystyle\operatorname{per}([T_{\vec{d},\vec{n}}]^{\dagger}T_{\vec{d},\vec{n}}+E(T)_{\vec{n},\vec{n}})$
$\displaystyle-\operatorname{per}\left([T_{2,\vec{n}}]^{\dagger}T_{2,\vec{n}}+E(T)_{\vec{n},\vec{n}}\right)$
$\displaystyle-\operatorname{per}\left([T_{3,\vec{n}}]^{\dagger}T_{3,\vec{n}}+E(T)_{\vec{n},\vec{n}}\right)$
$\displaystyle+\operatorname{per}(E(T)_{\vec{n},\vec{n}}).$

Here: $T_{2,\vec{n}}=(t_{22},t_{23},t_{23})$ and
$T_{3,\vec{n}}=(t_{32},t_{33},t_{33})$. We will not symbolically evaluate this
expression, but in the interest of providing simple test cases for future
software implementations, we see numerically that $p(\vec{d})=0.222\dots$ for
$\eta=1$, $p(\vec{d})=0.189$ for $\eta=0.9$ and $p(\vec{d})=0.069444\dots$ for
$\eta=0.5$.

###  A.5 Using the loop Torontonian

When performing calculations using the loop Torontonian, we need to know the
matrix $O$ and the vector $\vec{\gamma}$ for our state. Tools such as
Strawberry Fields [36] allow for conveniently computing the real means vector,
$\vec{\mu}$, and covariance matrix, $\sigma$, of a Gaussian state’s Wigner
function. To use $\sigma$ and $\vec{\mu}$ to find $O$ and $\vec{\gamma}$, we
can use functionality from The Walrus [34]. We can convert from $\sigma$ to
$\Sigma$ by using the Qmat function (The Walrus uses a different ordering for
$\vec{\zeta}$ and so we also apply a complex conjugate to match our definition
of $\Sigma$), and similarly convert from $\vec{\mu}$ to $\vec{\alpha}$ using
complex_to_real_displacements. These can be used to find
$O=\mathbb{I}-\Sigma^{-1}$ and $\vec{\gamma}=(\Sigma^{-1}\vec{\alpha})^{*}$.

To show how to select the appropriate rows/columns for calculating threshold
detector outcome probabilities, we will consider a 5-mode experiment, with the
outcome $\vec{d}=(1,0,1,0,0)$, so $C=\\{1,3\\}$. This state can be represented
by a $10\times 10$ matrix $O$, with elements $o_{jk}$, and a 10-element vector
$\vec{\gamma}$ with elements $\gamma_{j}$. To evaluate the loop Torontonian,
we form $O_{CC}$ and $\vec{\gamma}_{C}$:

$\displaystyle O_{CC}=\begin{pmatrix}o_{11}&o_{13}&o_{16}&o_{18}\\\
o_{31}&o_{33}&o_{36}&o_{38}\\\ o_{61}&o_{63}&o_{66}&o_{68}\\\
o_{81}&o_{83}&o_{86}&o_{88}\\\
\end{pmatrix}\quad\vec{\gamma}_{C}=\begin{pmatrix}\gamma_{1}\\\ \gamma_{3}\\\
\gamma_{6}\\\ \gamma_{8}\end{pmatrix}.$ (51)

This assumed the basis vector ordering convention specified by Eq. 21, where
mode $j$ corresponds to basis vectors $j$ and $j+M$, as is used in The Walrus
[34].

##  Appendix B Derivation of the characteristic function

We start from the photon number probability distribution in Eq. 6 and the
definitions of the characteristic function in Eq. 7 and rearrange terms

$\displaystyle\chi(\vec{\phi})$
$\displaystyle=\sum_{\vec{n}}\exp\left(\mathrm{i}\sum_{i=1}^{M}\phi_{i}n_{i}\right)\left|\bra{\vec{n}}\hat{\mathcal{U}}\ket{\Phi_{0}}\right|^{2}$
(52)
$\displaystyle=\sum_{\vec{n}}\exp\left(\mathrm{i}\sum_{i=1}^{M}\phi_{i}n_{i}\right)\bra{\Phi_{0}}\hat{\mathcal{U}}^{\dagger}\ket{\vec{n}}\bra{\vec{n}}\hat{\mathcal{U}}\ket{\Phi_{0}}$
(53)
$\displaystyle=\bra{\Phi_{0}}\hat{\mathcal{U}}^{\dagger}\sum_{\vec{n}}\exp\left(\mathrm{i}\sum_{i=1}^{M}\phi_{i}n_{i}\right)\ket{\vec{n}}\bra{\vec{n}}\hat{\mathcal{U}}\ket{\Phi_{0}}.$
(54)

Notice that we can define an operator, $\hat{\mathcal{U}}_{\vec{\phi}}$, which
acts like

$\displaystyle\hat{\mathcal{U}}_{\vec{\phi}}\ket{\vec{n}}=\exp\left(\mathrm{i}\sum_{i=1}^{M}\phi_{i}n_{i}\right)\ket{\vec{n}},$
(55)

as the operator given by the linear optical transformation $U_{\phi}$

$\displaystyle U_{\phi}=\bigoplus_{i=1}^{M}\exp(\mathrm{i}\phi_{i}).$ (56)

If we include this in our expression above, we find

$\displaystyle\chi(\vec{\phi})$
$\displaystyle=\bra{\Phi_{0}}\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{U}}_{\vec{\phi}}\sum_{\vec{n}}\ket{\vec{n}}\bra{\vec{n}}\hat{\mathcal{U}}\ket{\Phi_{0}}$
(57)
$\displaystyle=\bra{\Phi_{0}}\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{U}}_{\vec{\phi}}\hat{\mathcal{U}}\ket{\Phi_{0}},$
(58)

where we have used the resolution of the identity to arrive at the answer.

##  Appendix C Derivation of the Bristolian

The Unitary Bristolian has a slightly simpler derivation than the Bristolian,
so we will begin by considering this case. We start by combining Eq. 13 with
Eq. 4:

$\displaystyle p(\vec{d})$ $\displaystyle=$ (59)
$\displaystyle\left(\prod_{j=1}^{M}n_{j}!\right)^{-1}\sum_{Z\in
P(C)}(-1)^{|Z|}\operatorname{per}([U^{\dagger}U_{\vec{x}(V,Z)}U]_{\vec{n},\vec{n}}).$

$U_{\vec{x}(V,Z)}$ is formed by defining $U_{\vec{x}}$ with $x_{j}=0$ if $j\in
V$ or $j\in Z$, and $x_{j}=1$ otherwise. Since $U_{\vec{x}(V,Z)}$ contains
only zeros on the rows and columns given by the elements of $V$ and $Z$, this
is equivalent to deleting the rows and columns of $U$ and $U^{\dagger}$
respectively according to the elements of $V$ and $Z$. Therefore, we can write
the sum above as

$\displaystyle\sum_{Y\in
P(C)}(-1)^{|C|-|Y|}\operatorname{per}\left([(U_{Y})^{\dagger}U_{Y}]_{\vec{n},\vec{n}}\right).$
(60)

Here, we sum over the modes which are marginalised, instead of summing over
the modes being projected into the vacuum state. This corresponds to using the
substitution $Z$ for $Y=C\setminus Z$, and noticing that the sum over
$Y=C\setminus Z:\ Z\in P(C)$ is the same as the sum over $Z\in P(C)$. We use
$U_{Y}$ to denote selecting only the rows of $U$ according to the elements in
the set $Y$. Note that since we delete rows when all their elements are set to
all zeros, the permanent when $Y$ is the empty set should be zero, as it
corresponds to the permanent of an all zeros matrix, rather an empty matrix.

Because Eq. (60) contains an inclusion/exclusion formula, like that of Ryser’s
permanent formula [37], it could be viewed as a 3-dimensional permanent,
similar to those which appear elsewhere in quantum photonics [23, 38].
However, we note that to compute the probability for input state
$\ket{\vec{n}}$, interferometer transformation $U$ and threshold detector
pattern $\vec{d}$, it is sufficient to know only the rows of $U$ which
correspond to nonzero elements of $\vec{d}$ and columns of $U$ given by
$\vec{n}$. Therefore the input to this problem is a matrix, so we chose to
write this probability in terms of a new matrix function, the Unitary
Bristolian

$\displaystyle
p(\vec{d})=\frac{\operatorname{ubrs}(U_{\vec{d},\vec{n}})}{\prod_{j=1}^{M}n_{j}!}.$
(61)

We construct $U_{\vec{d},\vec{n}}$ from $U$ as described under Eq. 12. The
Unitary Bristolian, $\operatorname{ubrs}$, is a matrix function which acts on
some $m\times n$ matrix, $A$:

$\displaystyle\operatorname{ubrs}(A)=\sum_{Y\in
P([m])}(-1)^{m-|Y|}\operatorname{per}([A_{Y}]^{\dagger}A_{Y}),$ (62)

which is the form this is reported in Eq. 18.

Now, we are ready to derive the more general formula for the Bristolian.
Unlike the permanent, the Bristolian can be generalised for calculating
marginal detection probabilities. Consider that we wish to calculate the
probability of observing detector clicks for modes in $C$, vacuum in modes
$V$, and marginalise over modes in $B$, so the union of $C$, $V$ and $B$ is
$[M]$. We can write this probability as

$\displaystyle p(\vec{d}_{C}=\vec{1},\vec{d}_{V}=\vec{0})=$ (63)
$\displaystyle\left(\prod_{i=1}^{M}n_{i}!\right)^{-1}\sum_{Y\in
P(C)}(-1)^{|C|-|Y|}\operatorname{per}\left([(U_{Y\cup B})^{\dagger}U_{Y\cup
B}]_{\vec{n},\vec{n}}\right).$

Here $U_{Y\cup B}$ is constructed by selecting rows of $U$ according to the
elements of $Y$ and $B$. We could also write the summation in Eq. (63) as

$\displaystyle\sum_{Y\in
P(C)}(-1)^{|C|-|Y|}\operatorname{per}\left[\begin{pmatrix}U_{Y,\vec{n}}\\\
U_{B,\vec{n}}\end{pmatrix}^{\dagger}\begin{pmatrix}U_{Y,\vec{n}}\\\
U_{B,\vec{n}}\end{pmatrix}\right],$ (64)

where $\begin{pmatrix}U_{Y,\vec{n}}\\\ U_{B,\vec{n}}\end{pmatrix}$ is the
augmented $(|Y|+|B|)\times N$ matrix formed by stacking the matrices
$U_{Y,\vec{n}}$ and $U_{B,\vec{n}}$, where $N=\sum_{j}n_{j}$. The notation
$U_{C,\vec{n}}$ is used to show that we take rows of $U$ according to the set
$C$ and repeat the columns of $U$ according to $\vec{n}$, and equivalently for
$U_{B,\vec{n}}$.

A counter-intuitive feature of this formula is that it depends on matrix
elements of the linear transformation that are ignored by our measurements. In
an experiment, changing the elements of $U_{B}$ should have no impact on
$p(\vec{d}_{C}=\vec{1},\vec{d}_{V}=\vec{0})$. Following this argument, we
propose that we only need to know $U_{C}$, and we can construct $U_{B}$ by
performing a unitary dilation. This is particularly helpful for lossy
experiments, where we typically do not have an understanding of the full
unitary transformation acting on the both the experiment’s (lossy) modes and
the loss modes of its environment.

Any open quantum dynamics can be expressed as unitary evolution on of a larger
system via unitary dilation. For a non-unitary transformation given by a
rectangular matrix, $T$, with singular values all $\leq 1$, we can write the
unitary dilation

$\displaystyle
U(T)=\begin{pmatrix}T&\left(\mathbb{I}-TT^{\dagger}\right)^{1/2}\\\
\left(\mathbb{I}-T^{\dagger}T\right)^{1/2}&-T^{\dagger}\end{pmatrix}.$ (65)

With this, we are now ready to write down the probability for the general case
of a Fock state input, linear optical experiment. We consider that we have
$M_{\text{in}}$ input modes, initialised in the state $\vec{n}$. These
propagate through a nonunitary linear transformation, $T$, before being
detected by $M_{\text{out}}$ threshold detectors, which give an outcome
$\vec{d}$, which is a length-$M_{\text{out}}$ bit-string. So $T$ is given by
an $M_{\text{out}}\times M_{\text{in}}$ matrix.

First, we dilate $T$ according to Eq. (65), giving an
$(M_{\text{in}}+M_{\text{out}})$-dimension unitary matrix. In this
construction, $U_{C}=T$ and $U_{B}$ is given by
$\left(\mathbb{I}-T^{\dagger}T\right)^{1/2}$. We also notice that we can
explicitly write out the multiplication of the augmented matrices:

$\displaystyle\begin{pmatrix}A\\\ R\end{pmatrix}^{\dagger}\begin{pmatrix}A\\\
R\end{pmatrix}=A^{\dagger}A+R^{\dagger}R.$ (66)

Combining these observations with Eq. (63) and Eq. (64), we can write down the
probability of measuring a click pattern $\vec{d}$ on a Fock state $\vec{n}$
evolving through a nonunitary transformation $T$ as

$\displaystyle
p(\vec{d})=\frac{\operatorname{brs}\left(T_{\vec{d},\vec{n}},E(T)_{\vec{n},\vec{n}}\right)}{\prod_{j=1}^{M}n_{j}!}.$
(67)

Where we have introduced a new matrix function, the Bristolian,
$\operatorname{brs}$

$\operatorname{brs}\left(A,E\right)=\\\ \sum_{Y\in
P([m])}(-1)^{m-|Y|}\operatorname{per}\left([A_{Y}]^{\dagger}A_{Y}+E\right),$
(68)

where $m=\sum_{j}d_{j}=|C|$, is the total number of clicks. We also define the
matrix

$\displaystyle E(T)=\mathbb{I}-T^{\dagger}T.$ (69)

To arrive at this equation, we have used

$\displaystyle R(T)_{:,\vec{n}}$
$\displaystyle=\left[\left(\mathbb{I}-T^{\dagger}T\right)^{1/2}\right]_{:,\vec{n}}$
(70) $\displaystyle[R(T)_{:,\vec{n}}]^{\dagger}R(T)_{:,\vec{n}}$
$\displaystyle=E(T)_{\vec{n},\vec{n}},$ (71)

where the subscript notation $:,\vec{n}$ means that we select the columns
according to $\vec{n}$.

##  Appendix D Unitary Bristolian reduction to the permanent for clicks equal
to photons

When the number of input photons is equal to the number of threshold detector
clicks, the Unitary Bristolian has a square matrix as an input. To see $N$
clicks for an $N$ photon input state, we know that each threshold detector
must have detected exactly 1 photon, and so this event could also be modelled
by using photon number projectors instead of the click measurement operator,
and we can use Eq. 12 to calculate its amplitude. Therefore, we expect that
the Bristolian of a square matrix should reduce to the absolute square of the
permanent of the same matrix.

For a square $N\times N$ matrix, $A$, with elements $a_{jk}$, we can write the
Unitary Bristolian as

$\displaystyle\operatorname{ubrs}(A)=$ $\displaystyle\sum_{Y\in
P([N])}(-1)^{N-|Y|}\operatorname{per}([A_{Y}]^{\dagger}A_{Y}).$ (72)

Then we can expand the permanent using Ryser’s formula

$\operatorname{per}([A_{Y}]^{\dagger}A_{Y})=\\\ \sum_{Z\in
P([N])}(-1)^{N-|Z|}\prod_{i=1}^{N}\sum_{j\in
Z}\left[(A_{Y})^{\dagger}A_{Y}\right]_{ij}$ (73) $\displaystyle=\sum_{Z\in
P([N])}(-1)^{N-|Z|}\prod_{i=1}^{N}\sum_{j\in Z}\sum_{k\in Y}a^{*}_{ki}a_{kj}.$
(74)

This gives

$\displaystyle\operatorname{ubrs}(A)=$ $\displaystyle\sum_{Y,Z\in
P([N])}(-1)^{|Y|+|Z|}\prod_{i=1}^{N}\sum_{\begin{subarray}{c}j\in Z\\\ k\in
Y\end{subarray}}a^{*}_{ki}a_{kj},$ (75)

which is the Ryser-style formula for the 3-tensor permanent, as defined in
Ref. [23]

$\displaystyle\operatorname{ubrs}(A)=\operatorname{per}(B)=\sum_{\sigma,\rho\in
S_{N}}\prod_{i=1}^{N}b_{i\sigma(i)\rho(i)},$ (76)

with $\sigma$ and $\rho$ being elements of the permutation group $S_{N}$.
However, this 3-tensor, $B$, with elements $b_{ijk}=a^{*}_{ki}a_{kj}$ is very
structured. This kind of structure is discussed in Ref. [23], where it is
shown that it allows us to factor this expression into

$\displaystyle\operatorname{ubrs}(A)=$
$\displaystyle\operatorname{per}(A^{*})\operatorname{per}(A)=\operatorname{per}(A)^{*}\operatorname{per}(A)$
(77) $\displaystyle=$ $\displaystyle\ |\operatorname{per}(A)|^{2}.$ (78)

##  Appendix E Computing the Bristolian as a sum over Fock state probabilities

In the absence of our expression for the Bristolian in Eq. 15, the only known
way in the literature to compute probabilities for Fock states measured with
threshold detectors is to consider all possible events which could lead to the
witnessed outcome, and sum all their corresponding probabilities.

To provide an example for how the complexity of this method compares to the
Bristolian, we consider an $M$ mode experiment with uniform transmission,
$\eta$, $N$ input photons and $n$ clicks (with $n\leq N$). For interferometers
with imbalanced losses, we must consider losses at both the input and output
of the interferometer. However, here we are considering balanced loss which
allows us to consider that any photon loss occurs before the photons reach the
interferometer, as balanced loss commutes with linear optics. Because any
number of photons between $n$ and $N$ can be transmitted through this loss
channel and lead to an $n$ click event, we must consider all
$\sum_{j=n}^{N}\binom{N}{j}$ possible configurations for how these photons
could have been transmitted. For each input configuration, we then need to
consider all the ways that these photons can bunch within the detectors. We
know that at least 1 photon must arrive in each detector, which leaves $j-n$
photons left, which can be configured in any arrangement. There are
$\binom{X+Y-1}{Y}$ ways of arranging $Y$ photons into $X$ modes, and here we
have $j-n$ photons which can be arrive in $n$ modes, meaning we need to
calculate $\binom{j-n+n-1}{j-n}=\binom{j-1}{j-n}$ permanents for each input
configuration. This results in a total of
$\sum_{j=n}^{N}\binom{N}{j}\binom{j-1}{j-n}$ permanents.

Each permanent for a $j$ photon configuration has a time complexity of
$j2^{j}$, giving an overall time complexity which is lower bounded by
$\binom{N}{n}n2^{n}$. By comparison, the Bristolian complexity of $n2^{2n}$
provides a superexponential speedup when $N$ and $N-n$ are large.

##  Appendix F Derivation of the loop Torontonian

We start with some definitions. For a Gaussian state with complex Husimi
covariance matrix $\Sigma$ and complex vector of means $\alpha$ we define

$\displaystyle O$ $\displaystyle=\mathbb{I}-\Sigma^{-1},$ (79)
$\displaystyle\vec{\gamma}$ $\displaystyle=(\Sigma^{-1}\vec{\alpha})^{*},$
(80)

which uniquely specify the _photon number_ statistics of the Gaussian state
[22] via loop Hafnians [21]. Given a photon number outcome
$\vec{n}=(n_{1},\ldots,n_{M})$ its probability is given by

$\displaystyle
p(\vec{n})=p(\vec{0})\operatorname{lhaf}(XO_{\vec{n},\vec{n}},\vec{\gamma}_{\vec{n}}),$
(81)

where $p(\vec{0})$ is the vacuum probability,
$X=\left[\begin{smallmatrix}0&\mathbb{I}\\\
\mathbb{I}&0\end{smallmatrix}\right]$ and $O_{\vec{n},\vec{n}}$ and
$\vec{\gamma}_{\vec{n}}$ are submatrices of $O$ and $\vec{\gamma}$, found by
using $n_{j}$ repetitions of the rows/columns corresponding to mode $j$
(recalling that each mode $j$ corresponds to 2 rows/columns of $O$).

We will now show that the threshold probabilities can also be written in terms
of the quantities defined in Eq. (79) and Eq. (80).

By ordering our basis vectors such that modes which see a click ($C$) are
arranged to be before modes which see vacuum ($V$), the matrix $O$ can be
written in block form

$\displaystyle O=\begin{pmatrix}O_{CC}&O_{CV}\\\ O_{VC}&O_{VV}\end{pmatrix},$
(82)

then, using Schur complements, we can see that

$\displaystyle
O_{CC}=\mathbb{I}-\left[\Sigma_{CC}-\Sigma_{CV}[\Sigma_{VV}]^{-1}\Sigma_{VC}\right]^{-1}.$
(83)

We are now ready to investigate threshold probabilities. We start with Eq. 4
which we write as

$\displaystyle p(\vec{d})=p(\vec{0})\sum_{Z\in
P(C)}(-1)^{|Z|}\frac{p(\vec{m}_{Z}=0,\vec{m}_{V}=0)}{p(\vec{0})}.$ (84)

For a given $Z$, the term inside the sum can be written as in Eq. 23. Note
that the argument inside the exponential in said equation can be rewritten as

$\displaystyle\vec{\alpha}_{W}^{\dagger}[\Sigma_{WW}]^{-1}\vec{\alpha}_{W}=$
$\displaystyle\vec{\gamma}_{W}^{t}\Sigma_{WW}\vec{\gamma}_{W}^{*}$ (85)
$\displaystyle+\vec{\gamma}_{W}^{t}\Sigma_{WY}\vec{\gamma}_{Y}^{*}+\vec{\gamma}_{Y}^{t}\Sigma_{YW}\vec{\gamma}_{W}^{*}$
$\displaystyle+\vec{\gamma}_{Y}^{t}\Sigma_{YW}[\Sigma_{WW}]^{-1}\Sigma_{WY}\vec{\gamma}_{Y}^{*},$

where $W$ denotes the union of sets $Z$ and $V$, and $Y$ denotes the modes not
included in $W$. Here we have used that
$\vec{\alpha}_{W}=\Sigma_{WW}\vec{\gamma}_{W}^{*}+\Sigma_{WY}\vec{\gamma}_{Y}^{*}$,
$\Sigma_{WW}^{\dagger}=\Sigma_{WW}$ and $\Sigma_{WY}^{\dagger}=\Sigma_{YW}$.
We can similarly write the argument of the exponential in $p(\vec{0})$ as

$\displaystyle\vec{\alpha}^{\dagger}[\Sigma]^{-1}\vec{\alpha}=\vec{\gamma}^{t}\Sigma\vec{\gamma}^{*}=$
$\displaystyle\vec{\gamma}_{W}^{t}\Sigma_{WW}\vec{\gamma}_{W}^{*}+\vec{\gamma}_{W}^{t}\Sigma_{WY}\vec{\gamma}_{Y}^{*}$
(86)
$\displaystyle+\vec{\gamma}_{Y}^{t}\Sigma_{YW}\vec{\gamma}_{W}^{*}+\vec{\gamma}_{Y}^{t}\Sigma_{YY}\vec{\gamma}_{Y}^{*}.$

With these two expressions we can then write

$\frac{p(\vec{m}_{W}=0)}{p(\vec{0})}=\sqrt{\frac{\det(\Sigma)}{\det(\Sigma_{WW})}}\\\
\times\exp\left[\tfrac{1}{2}\vec{\gamma}_{Y}^{t}[\Sigma_{YY}-\Sigma_{YW}[\Sigma_{WW}]^{-1}\Sigma_{WY}]\vec{\gamma}_{Y}^{*}\right]$
(87)
$\displaystyle=\frac{\exp\left[\tfrac{1}{2}\vec{\gamma}_{Y}^{t}[\mathbb{I}-O_{YY}]^{-1}\vec{\gamma}_{Y}^{*}\right]}{\sqrt{\det(\mathbb{I}-O_{YY})}}.$
(88)

In the last equation we used the result in Eq. (83), together with factorising
the determinant

$\det(\Sigma)=\\\
\det(\Sigma_{WW}){\det(\Sigma_{YY}-\Sigma_{YW}[\Sigma_{WW}]^{-1}\Sigma_{WY})},$
(89)

to show

$\displaystyle\frac{\det(\Sigma)}{\det(\Sigma_{WW})}=$
$\displaystyle\det(\Sigma_{YY}-\Sigma_{YW}[\Sigma_{WW}]^{-1}\Sigma_{WY})$ (90)
$\displaystyle=$ $\displaystyle\det\left([\mathbb{I}-O_{YY}]^{-1}\right)$ (91)
$\displaystyle=$ $\displaystyle\frac{1}{\det(\mathbb{I}-O_{YY})}.$ (92)

This corresponds to the identity
$\det(\Sigma_{WW})=\det(\Sigma)\det([\Sigma^{-1}]_{YY})$ from Ref. [39]. With
these observations we can write the sought after probability as

$\displaystyle
p(\vec{d})=p(\vec{0})\operatorname{ltor}\left(O_{CC},\vec{\gamma}_{C}\right),$
(93)

where $C$ is the set of modes with threshold detector clicks, and
$\operatorname{ltor}$ stands for the loop Torontonian defined as

$\operatorname{ltor}\left(O,\vec{\gamma}\right)=\\\ \sum_{Y\in
P([m])}(-1)^{m-|Y|}\frac{\exp\left[\tfrac{1}{2}\vec{\gamma}_{Y}^{t}[\mathbb{I}-O_{YY}]^{-1}\vec{\gamma}_{Y}^{*}\right]}{\sqrt{\det(\mathbb{I}-O_{YY})}}$
(94)

for a $2m\times 2m$ matrix $O$, and a $2m$-length vector $\vec{\gamma}$.

In arriving at this formula, we have swapped the summation over $Z$ for a
summation over $Y$, as we did in Eq. (60). This is the form that the equation
appears in the main text.

##  Appendix G Generating the loop Hafnian from the loop Torontonian

Just like the Torontonian [7], the loop Torontonian has an interesting
interpretation of being a generating function for photon number probabilities
of Gaussian states (now with nonzero displacement).

To this end, recall that a threshold probability can be obtained as a sum
(coarse graining) of many photon number events. For any $\vec{d}$, we can
write

$\displaystyle\frac{p(\vec{d})}{p(\vec{0})}=$
$\displaystyle\frac{p(\vec{n}=\vec{d})}{p(\vec{0})}+\sum_{\vec{k}\in\mathcal{C}_{\vec{d}}}\frac{p(\vec{n}=\vec{d}+\vec{k})}{p(\vec{0})}$
(95) $\displaystyle\operatorname{ltor}\left(O_{CC},\vec{\gamma}_{C}\right)=$
$\displaystyle\operatorname{lhaf}(XO_{\vec{d},\vec{d}},\vec{\gamma}_{\vec{d}})$
(96)
$\displaystyle+\sum_{\vec{k}\in\mathcal{C}_{\vec{d}}}\operatorname{lhaf}(XO_{(\vec{d}+\vec{k}),(\vec{d}+\vec{k})},\vec{\gamma}_{\vec{d}+\vec{k}}),$

where $\mathcal{C}_{\vec{d}}$ is the set of all vectors of integers that have
zero in the positions where $\vec{d}$ has zero and strictly positive integers
in all the other positions. Note that any element $\vec{k}$ in the set
satisfies $K=\sum_{i=1}^{M}k_{i}>0$. We now recall that loop hafnians satisfy
the following scaling property [21]

$\displaystyle\operatorname{lhaf}((\eta{A})_{\vec{n},\vec{n}},(\sqrt{\eta}{\gamma})_{\vec{n}})=\eta^{|N|}\operatorname{lhaf}({A}_{\vec{n},\vec{n}},{\gamma}_{\vec{n}}),$
(97)

with $N=\sum_{j}n_{j}$. This allows us to write

$\operatorname{ltor}\left(\eta
O_{CC},\sqrt{\eta}\vec{\gamma}_{C}\right)=\eta^{|C|}\operatorname{lhaf}(XO_{\vec{d},\vec{d}},\vec{\gamma}_{\vec{d}})\\\
+\sum_{\vec{k}\in\mathcal{C}_{\vec{d}}}\eta^{|C|+K}\operatorname{lhaf}(XO_{(\vec{d}+\vec{k}),(\vec{d}+\vec{k})},\vec{\gamma}_{(\vec{d}+\vec{k})}).$
(98)

If we set $\ell=|C|$ we can derive a formula for the loop hafnian

$\displaystyle\operatorname{lhaf}(XO_{\vec{d},\vec{d}},\vec{\gamma}_{\vec{d}})=$
$\displaystyle\left.\frac{1}{\ell!}\frac{d^{\ell}}{d\eta^{\ell}}\operatorname{ltor}\left(\eta
O_{\vec{d},\vec{d}},\sqrt{\eta}\vec{\gamma}_{\vec{d}}\right)\right|_{\eta=0}$
(99) $\displaystyle=$ $\displaystyle\sum_{Y\in
P([m])}(-1)^{m-|Y|}f\left(O_{YY},\vec{\gamma}_{Y}\right)$ (100)

where we define

$f(O,\vec{\gamma})=\left.\frac{1}{\ell!}\frac{d^{\ell}}{d\eta^{\ell}}q(O,\vec{\gamma})\right|_{\eta=0}$
(101)
$q(O,\vec{\gamma})=\frac{\exp\left[\tfrac{1}{2}\sqrt{\eta}\vec{\gamma}^{t}[\mathbb{I}-\eta
O]^{-1}\sqrt{\eta}\vec{\gamma}^{*}\right]}{\sqrt{\det(\mathbb{I}-\eta O)}}.$
(102)

One way of interpreting $f(O,\vec{\gamma})$ is that it selects the $\ell$th
coefficient of the polynomial expansion of $q$. Therefore, we do not require
knowledge of $q(O,\vec{\gamma})$ beyond order-$\ell$. Using the Mercator
series, we expand the denominator in Eq. (101) as

$\displaystyle\det(\mathbb{I}-\eta
O)^{-1/2}=\exp\left[\sum_{k=1}^{\infty}\frac{\operatorname{tr}\left(\eta^{k}O^{k}\right)}{2k}\right],$
(103)

so we can combine the denominator into the exponent. We also notice that we
can use a binomial expansion, $(1-x)^{-1}=\sum_{j=0}^{\infty}x^{j}$, to
substitute $[\mathbb{I}-\eta O]^{-1}=\sum_{k=1}^{\infty}\left(\eta
O\right)^{k-1}$. This gives us an expression

$q(O,\vec{\gamma})=\\\
\exp\left(\sum_{k=1}^{\infty}\left[\frac{\operatorname{tr}\left(O^{k}\right)}{2k}+\frac{\vec{\gamma}^{t}\left(O^{k-1}\right)\vec{\gamma}^{*}}{2}\right]\eta^{k}\right),$
(104)

which we Taylor-expand up to order $\ell$ to obtain

$q(O,\vec{\gamma})=\\\
\sum_{j=0}^{\ell}\frac{1}{j!}\left(\sum_{k=1}^{\ell}\left[\frac{\operatorname{tr}\left(O^{k}\right)}{2k}+\frac{\vec{\gamma}^{t}\left(O^{k-1}X\right)\vec{\gamma}}{2}\right]\eta^{k}\right)^{j}.$
(105)

We have used that $\vec{\gamma}^{*}=X\vec{\gamma}$. This gives us the exact
form of the trace formula algorithm for the loop Hafnian as presented in Ref.
[21] (bar some typos in the referenced paper). It is interesting that this
derivation uses various arguments which rely specifically on $O$ and
$\vec{\gamma}$ being formed by a Gaussian state covariance matrix (in
particular Eq. (95) and Eq. (96)), however the trace formula algorithm that we
are able to derive here is applicable to arbitrary symmetric matrices.

Figure 2:  Scattershot construction. Gates between modes, symbolised by
vertical lines with dots at the ends, show weak two-mode squeezing. The two-
mode squeezing operator between modes $i$ and $j$ is given
$\hat{S}_{i,j}(t)=\exp\left(r\left[\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}-\hat{a}_{i}\hat{a}_{j}\right]\right)$.

##  Appendix H Connecting the Bristolian and the Torontonian

In this Appendix we provide a formal link between the Bristolian and the
Torontonian. Concretely, we show that the Bristolian associated with the
threshold detection of a multimode Fock state with single photon or vacuum
inputs can be evaluated as a certain limit of a Torontonian.

The starting point of our derivation is the scattershot Boson sampling
construction shown in Fig. 2 and introduced in Ref. [20]. The diagram
represents $M$ two-mode squeezed vacuum states where the first half of the
modes (the heralded modes) are sent into an interferometer with transmission
matrix $T$, and the second half of the modes (the heralding modes) are sent
into threshold detectors.

For each individual two-mode squeezed vacuum with squeezing parameter $t$, it
is straightforward to show that conditioned on a click in the heralding mode,
the state of the heralded mode collapses to [40]

$\displaystyle\rho_{\checkmark}=\frac{1}{1+\bar{n}}|1\rangle\langle
1|+\sum_{m=2}^{\infty}\frac{1}{\bar{n}}\left(\frac{\bar{n}}{1+\bar{n}}\right)^{m}|m\rangle\langle
m|,$ (106)

where $\bar{n}=\sinh^{2}t$ is the mean photon number of either mode of the
two-mode squeezed vacuum. The probability of heralding the state
$\rho_{\checkmark}$ is given by

$\displaystyle p_{\checkmark}=\frac{\bar{n}}{1+\bar{n}}=\varepsilon^{2}\text{
with }\varepsilon=\tanh t.$ (107)

Note that as $\varepsilon\to 0$ the fidelity between $\rho_{\checkmark}$ and a
single photon Fock state approaches one, but at the same time the probability
of heralding the state approaches zero.

We can now study the covariance matrix $\Sigma$ of the Gaussian circuit in
Fig. 2. After some algebra, it can be shown that the matrix
$O=\mathbb{I}-\Sigma^{-1}$ dictating the threshold probabilities of the
Gaussian state is given by

$\displaystyle O(\varepsilon)=\varepsilon\begin{pmatrix}0&0&0&T\\\
0&\varepsilon E(T)^{*}&T^{t}&0\\\ 0&T^{*}&0&0\\\ T^{\dagger}&0&0&\varepsilon
E(T)\end{pmatrix}.$ (108)

We now want to calculate the probability that a subset of the first $M$ modes
clicks conditioned on all the modes on the second half clicking. We specify
the modes that clicked in the first half by $\vec{d}$ and then can write the
click pattern for the $2M$ modes to be $\vec{e}=\vec{d}\oplus\vec{1}$ where,
recall, $\vec{1}$ is the all ones vector of length $M$. We can now write the
conditional probability as

$\displaystyle
p(\vec{d}|\vec{1}_{M})=\frac{p(\vec{d}\oplus\vec{1})}{p_{\checkmark}^{M}}=\frac{p(\vec{e})}{p_{\checkmark}^{M}}.$
(109)

For the probability in the numerator of the right-hand side in the last
equation we can use Eq. 25. The vacuum probability of the $2M$-mode state is
given by

$\displaystyle
p(\vec{0})=\sqrt{\det(\mathbb{I}-O(\varepsilon))}=(1-\varepsilon^{2})^{M}.$
(110)

We can now write

$\displaystyle
p(\vec{d}|\vec{1}_{M})=\frac{(1-\varepsilon^{2})^{M}}{\varepsilon^{2M}}\operatorname{tor}(O(\varepsilon)_{CC})$
(111)

where $C$ is the union of the $\\{M+1,\ldots,2M\\}$ heralding modes and the
modes that have a one in the vector $\vec{d}$, i.e., the modes that click in
the first half.

As explained at the beginning of this section, in the limit where
$\varepsilon\to 0$ we know that the conditional state of the input heralded
modes becomes a product of single photons in each mode, and thus in this same
limit we can write

$\operatorname{brs}\left(T_{\vec{d},\vec{1}},E(T)_{\vec{1},\vec{1}}\right)=\\\
\lim_{\varepsilon\to
0}(\varepsilon^{-2}-1)^{M}\operatorname{tor}\left(O(\varepsilon)_{CC}\right).$
(112)

Note that we have so far only considered the case where single photons are
input in all the modes of the interferometer. The more general case where
vacuum is fed into some of the modes can be dealt with by applying a loss
channel with zero-transmission to the relevant modes. This is equivalent to
setting to zero the columns of $T$ where vacuum is fed. Using this argument we
find that for cases where $\vec{n}$ is a bistring, i.e., we only allow single
photons inputs, and setting $N=\sum_{i}\vec{n}_{i}$ we can write

$\operatorname{brs}\left(T_{\vec{d},\vec{n}},E(T)_{\vec{n},\vec{n}}\right)=\\\
\lim_{\varepsilon\to
0}(\varepsilon^{-2}-1)^{N}\operatorname{tor}\left(O(\varepsilon)_{CC}\right),$
(113)

where now $C$ is the union of the labels of the modes in which single photons
where input into the interferometers and the labels of the modes in which
clicks are registered. It is interesting to consider that one could
potentially derive the form of the Bristolian in terms of sums of permanents
by using the connection between permanents and determinants provided by the
MacMahon Master theorem [41].

Finally, note that one can also write circuits to herald multi-photon Fock
states using only threshold detectors as shown in Appendix D of Ref. [40].

##  Appendix I Time complexities

To calculate a marginal vacuum probability for a Fock state evolved through a
linear interferometer, we can compute a permanent, as given in Eq. 13 If we
are detecting vacuum in modes given by $V$ and marginalising over all other
modes, given by $B$, then the matrix
$[U^{\dagger}U_{\vec{x}}U]_{\vec{n},\vec{n}}$ has rank $\leq|B|$, so its
permanent can be computed in time $\mathcal{O}(N^{\mathcal{O}(|B|)})$ using
the algorithm introduced in Sec. III of Ref. [42]. If we are marginalising
over many modes, and therefore $|V|$ is small, we can instead consider using
the algorithm introduced in Sec. IV. A of Ref. [15] to compute the permanent
in time $\mathcal{O}(N^{2|V|+1})$.

For a general permanent of an $n\times n$ matrix, Ryser’s algorithm [37] has
the best known complexity of $\mathcal{O}(n2^{n})$. For computing the
Bristolian, to calculate the permanents inside the summation in Eq. 15, we may
sometimes be able to use the faster algorithms above, however the dominant
complexity for this formula would still come from computing the intermediate
cases, when neither of the faster algorithms are applicable, where Ryser’s
algorithm may be the fastest option. This is upper bounded by
$\mathcal{O}(n2^{n})$ and there are $\mathcal{O}(2^{m})$ terms in the sum.
Therefore, a Bristolian of an $m\times n$ matrix $A$, and $n\times n$ matrix
$E$ has a time complexity of $\mathcal{O}(n2^{m+n})$.

We do not claim that these complexities are optimal. The structure of this
matrix function may be exploited to reduce the complexity, for example by
using methods similar to those for low rank permanents [42], exploiting
recursion [24] and using Laplace expansions [43]. However, we leave it as an
open problem to find faster algorithms for the Bristolian.

For the loop Torontonian, we find a comparable complexity to the original
algorithm for the Torontonian, with complexity $\mathcal{O}(m^{3}2^{m})$ for a
$2m\times 2m$ matrix $O$ and $2m$-length vector $\vec{\gamma}$. In each step,
we must compute a matrix inverse and a matrix determinant, both having
$m^{3}$-time algorithms. However, these steps can make use of the Cholesky
decomposition of $O$, so we can improve the polynomial prefactor, following
the methods described in Ref. [24].

Both of these methods see a quadratic penalty as compared to the fastest
methods for calculating photon number probabilities of pure states. This can
be understood as being caused by the threshold detection operators in Eq. 1b
having high rank, whereas the photon number operators are rank-1 projectors.
For Gaussian state calculations, we also see a quadratic cost for calculating
photon number probabilities when the state is mixed [44, 22], so we can also
understand this quadratic penalty as being a result of introducing mixture
into the projected state. This differs from the case of sampling, where it was
shown that sampling threshold detector clicks can be simulated with the same
complexity as single photon measurements [14]. If we accept approximate
expressions, accurate to additive error, we can efficiently compute the
probabilities presented in this work using Monte-Carlo phase space methods
[45]. However, for events with small probability, these methods can quickly
become impractical due to a large relative error.

##  Appendix J Accuracy improvements of a threshold detection model over a
single photon projection approximation

To provide a quantitative demonstration of the accuracy improvements of an
exact model of threshold detection versus the typical approximation of using
single photon projection, we calculate the full probability distribution for a
4 photon Fock state boson sampling experiment [28] using the Bristolian with
the number of modes ranging from 4 to 12, choosing a transitivity of
$\eta=0.6$. We also find the probability distribution given by 0 or 1 photon
Fock state projective measurement. We calculate the total variation distance
(TVD) between these distributions for 100 different Haar random linear optical
interferometers for each number of modes. The results, shown in Fig. 3, show
that even as we approach the $M=N^{2}$ regime, we do not see a convergence
between these two distributions. In particular, we see that the TVD for all
the experiments lies typically within the range of 5%-12%. This highlights the
importance of using the correct mathematical description of the experiment in
order to best understand the results.

We also performed an analysis of data from a recent experiment of displaced
Gaussian boson sampling [10], which used a two-mode squeezed vacuum and a
single coherent state as input states. Here, we look at the 2 photon
probability distribution as estimated using the experiment for different
levels of displacement, as labelled by the mean photon number of the input
coherent state, $\langle n_{\alpha}\rangle$, in Fig. 4. For each probability
distribution, we compare the TVD to a model which uses threshold detection,
using the loop Torontonian, and a model which assumes Fock state projections,
using the loop Hafnian. The experiment of Ref. [10] uses threshold detection,
and we see that the model which uses loop Torontonians provides a more
accurate model of the experiment. For the largest displacements, when looking
at the ratio of the TVD, the loop Torontonian gives a probability distribution
which is 16% closer to the experimental data.

Figure 3: Total variation distance (TVD) between a simulation of a threshold
detection based, 4 photon, Fock state boson sampling experiment and an
approximate model which uses Fock state projection. For each number of modes,
100 Haar random unitary matrices are sampled and the full probability
distribution is calculated. Threshold detection calculations are performed
using the Bristolian, Fock state projection calculation uses a formula based
on matrix permanents, as discussed in Appendix E. Figure 4: Total variation
distance (TVD) of 2 photon displaced Gaussian boson sampling distributions
from the experiment reported in Ref. [10] against a threshold detector model,
using the loop Torontonian ($\operatorname{ltor}$, dotted line), and a Fock
state projection model, using the loop Hafnian ($\operatorname{lhaf}$, dashed
line). This is calculated for different mean photon numbers of the input
coherent state $\langle n_{\alpha}\rangle$.

##  Appendix K Software implementation

For an efficient and parallelisable implementation of the Bristolian and loop
Torontonian, we use just-in-time compilation provided by Numba [46]. Our code
is available in the open-source Python package The Walrus [34]
(https://github.com/XanaduAI/thewalrus) in releases from 0.19.0 onward, and
were contributed in pull requests (316) and (317).

## References

  * Flamini  _et al._ [2018] F. Flamini, N. Spagnolo, and F. Sciarrino, Photonic quantum information processing: a review, Reports on Progress in Physics 82, 016001 (2018).
  * Reddy  _et al._ [2020] D. V. Reddy, R. R. Nerem, S. W. Nam, R. P. Mirin, and V. B. Verma, Superconducting nanowire single-photon detectors with 98% system detection efficiency at 1550 nm, Optica 7, 1649 (2020).
  * Warburton  _et al._ [2009] R. E. Warburton, M. Itzler, and G. S. Buller, Free-running, room temperature operation of an ingaas/inp single-photon avalanche diode, Applied Physics Lett. 94, 071116 (2009).
  * Zhong  _et al._ [2021] H.-S. Zhong, Y.-H. Deng, J. Qin, H. Wang, M.-C. Chen, L.-C. Peng, Y.-H. Luo, D. Wu, S.-Q. Gong, H. Su, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, J. J. Renema, C.-Y. Lu, and J.-W. Pan, Phase-programmable Gaussian boson sampling using stimulated squeezed light, Phys. Rev. Lett. 127, 180502 (2021).
  * Slussarenko  _et al._ [2017] S. Slussarenko, M. M. Weston, H. M. Chrzanowski, L. K. Shalm, V. B. Verma, S. W. Nam, and G. J. Pryde, Unconditional violation of the shot-noise limit in photonic quantum metrology, Nature Photonics 11, 700 (2017).
  * Shalm  _et al._ [2015] L. K. Shalm, E. Meyer-Scott, B. G. Christensen, P. Bierhorst, M. A. Wayne, M. J. Stevens, T. Gerrits, S. Glancy, D. R. Hamel, M. S. Allman, K. J. Coakley, S. D. Dyer, C. Hodge, A. E. Lita, V. B. Verma, C. Lambrocco, E. Tortorici, A. L. Migdall, Y. Zhang, D. R. Kumor, W. H. Farr, F. Marsili, M. D. Shaw, J. A. Stern, C. Abellán, W. Amaya, V. Pruneri, T. Jennewein, M. W. Mitchell, P. G. Kwiat, J. C. Bienfang, R. P. Mirin, E. Knill, and S. W. Nam, Strong loophole-free test of local realism, Phys. Rev. Lett. 115, 250402 (2015).
  * Quesada  _et al._ [2018] N. Quesada, J. M. Arrazola, and N. Killoran, Gaussian boson sampling using threshold detectors, Phys. Rev. A 98, 062322 (2018).
  * Wang  _et al._ [2018] H. Wang, W. Li, X. Jiang, Y.-M. He, Y.-H. Li, X. Ding, M.-C. Chen, J. Qin, C.-Z. Peng, C. Schneider, M. Kamp, W.-J. Zhang, H. Li, L.-X. You, Z. Wang, J. P. Dowling, S. Höfling, C.-Y. Lu, and J.-W. Pan, Toward scalable boson sampling with photon loss, Physical review letters 120, 230502 (2018).
  * Wang  _et al._ [2019] H. Wang, J. Qin, X. Ding, M.-C. Chen, S. Chen, X. You, Y.-M. He, X. Jiang, L. You, Z. Wang, C. Schneider, J. J. Renema, S. Höfling, C.-Y. Lu, and J.-W. Pan, Boson sampling with 20 input photons and a 60-mode interferometer in a $1{0}^{14}$-dimensional Hilbert space, Phys. Rev. Lett. 123, 250503 (2019).
  * Thekkadath  _et al._ [2022] G. Thekkadath, S. Sempere-Llagostera, B. Bell, R. Patel, M. Kim, and I. Walmsley, Experimental demonstration of gaussian boson sampling with displacement, PRX Quantum 3, 020336 (2022).
  * Bentivegna  _et al._ [2015] M. Bentivegna, N. Spagnolo, C. Vitelli, F. Flamini, N. Viggianiello, L. Latmiral, P. Mataloni, D. J. Brod, E. F. Galvão, A. Crespi, _et al._ , Experimental scattershot boson sampling, Science advances 1, e1400255 (2015).
  * Paesani  _et al._ [2019] S. Paesani, Y. Ding, R. Santagati, L. Chakhmakhchyan, C. Vigliar, K. Rottwitt, L. K. Oxenløwe, J. Wang, M. G. Thompson, and A. Laing, Generation and sampling of quantum states of light in a silicon chip, Nature Physics 15, 925 (2019).
  * Van Den Nes [2010] M. Van Den Nes, Classical simulation of quantum computation, the Gottesman-Knill theorem, and slightly beyond, Quantum Information & Computation 10, 258 (2010).
  * Bulmer  _et al._ [2022] J. F. F. Bulmer, B. A. Bell, R. S. Chadwick, A. E. Jones, D. Moise, A. Rigazzi, J. Thorbecke, U.-U. Haus, T. V. Vaerenbergh, R. B. Patel, I. A. Walmsley, and A. Laing, The boundary for quantum advantage in gaussian boson sampling, Science Advances 8, eabl9236 (2022).
  * Ivanov and Gurvits [2020] D. A. Ivanov and L. Gurvits, Complexity of full counting statistics of free quantum particles in product states, Phys. Rev. A 101, 012303 (2020).
  * Scheel [2004] S. Scheel, Permanents in linear optical networks, arXiv preprint quant-ph/0406127 (2004).
  * García-Patrón  _et al._ [2019] R. García-Patrón, J. J. Renema, and V. Shchesnovich, Simulating boson sampling in lossy architectures, Quantum 3, 169 (2019).
  * Serafini [2017] A. Serafini,  _Quantum Continuous Variables, A Primer of Theoretical Methods_, 1st ed. (CRC Press, 2017).
  * Quesada [2019] N. Quesada, Franck-Condon factors by counting perfect matchings of graphs with loops, J. Chem. Phys. 150, 164113 (2019).
  * Lund  _et al._ [2014] A. P. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L. O’Brien, and T. C. Ralph, Boson sampling from a Gaussian state, Phys. Rev. Lett. 113, 100502 (2014).
  * Björklund  _et al._ [2019] A. Björklund, B. Gupt, and N. Quesada, A faster hafnian formula for complex matrices and its benchmarking on a supercomputer, Journal of Experimental Algorithmics (JEA) 24, 1 (2019).
  * Quesada  _et al._ [2019] N. Quesada, L. G. Helt, J. Izaac, J. M. Arrazola, R. Shahrokhshahi, C. R. Myers, and K. K. Sabapathy, Simulating realistic non-Gaussian state preparation, Phys. Rev. A 100, 022341 (2019).
  * Tichy [2015] M. C. Tichy, Sampling of partially distinguishable bosons and the relation to the multidimensional permanent, Phys. Rev. A 91, 022316 (2015).
  * Kaposi  _et al._ [2021] Á. Kaposi, Z. Kolarovszki, T. Kozsik, Z. Zimborás, and P. Rakyta, Polynomial speedup in Torontonian calculation by a scalable recursive algorithm, arXiv preprint arXiv:2109.04528 (2021).
  * Knill  _et al._ [2001] E. Knill, R. Laflamme, and G. J. Milburn, A scheme for efficient quantum computation with linear optics, Nature 409, 46 (2001).
  * Kok  _et al._ [2007] P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, Linear optical quantum computing with photonic qubits, Rev. Mod. Phys. 79, 135 (2007).
  * Rudolph [2017] T. Rudolph, Why I am optimistic about the silicon-photonic route to quantum computing, APL Photonics 2, 030901 (2017).
  * Aaronson and Arkhipov [2011] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, in  _Proceedings of the forty-third annual ACM symposium on Theory of computing_ (2011) pp. 333–342.
  * You  _et al._ [2021] X. You, M.-Y. Zheng, S. Chen, R.-Z. Liu, J. Qin, M.-C. Xu, Z.-X. Ge, T.-H. Chung, Y.-K. Qiao, Y.-F. Jiang, _et al._ , Quantum interference between independent solid-state single-photon sources separated by 300 km fiber, arXiv preprint arXiv:2106.15545 (2021).
  * Schuld  _et al._ [2020] M. Schuld, K. Brádler, R. Israel, D. Su, and B. Gupt, Measuring the similarity of graphs with a Gaussian boson sampler, Phys. Rev. A 101, 032314 (2020).
  * Huh  _et al._ [2015] J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean, and A. Aspuru-Guzik, Boson sampling for molecular vibronic spectra, Nature Photonics 9, 615 (2015).
  * Afek  _et al._ [2010] I. Afek, O. Ambar, and Y. Silberberg, High-NOON states by mixing quantum and classical light, Science 328, 879 (2010).
  * Note [1] The initial inspiration for us to derive the loop Torontonian came from the need to calculate threshold detection statistics for the quantum communication protocols proposed in Ref. [47].
  * Gupt  _et al._ [2019] B. Gupt, J. Izaac, and N. Quesada, The Walrus: a library for the calculation of hafnians, hermite polynomials and Gaussian boson sampling, Journal of Open Source Software 4, 1705 (2019).
  * Hong  _et al._ [1987] C.-K. Hong, Z.-Y. Ou, and L. Mandel, Measurement of subpicosecond time intervals between two photons by interference, Phys. Rev. Lett. 59, 2044 (1987).
  * Killoran  _et al._ [2019] N. Killoran, J. Izaac, N. Quesada, V. Bergholm, M. Amy, and C. Weedbrook, Strawberry fields: A software platform for photonic quantum computing, Quantum 3, 129 (2019).
  * Ryser [1963] H. J. Ryser,  _Combinatorial Mathematics_, Vol. 14 (Mathematical Association of America, 1963).
  * Rudolph [2021] T. Rudolph, Perhaps they are everywhere? Undetectable distributed quantum computation and communication for alien civilizations can be established using thermal light from stars, arXiv preprint arXiv:2107.13023 (2021).
  * Shi and Byrnes [2022] J. Shi and T. Byrnes, Effect of partial distinguishability on quantum supremacy in gaussian boson sampling, npj Quantum Information 8, 1 (2022).
  * Bourassa  _et al._ [2021] J. E. Bourassa, N. Quesada, I. Tzitrin, A. Száva, T. Isacsson, J. Izaac, K. K. Sabapathy, G. Dauphinais, and I. Dhand, Fast simulation of bosonic qubits via Gaussian functions in phase space, PRX Quantum 2, 040315 (2021).
  * MacMahon [2001] P. A. MacMahon,  _Combinatory Analysis, Volumes I and II_ , Vol. 137 (American Mathematical Soc., 2001).
  * Barvinok [1996] A. I. Barvinok, Two algorithmic results for the traveling salesman problem, Mathematics of Operations Research 21, 65 (1996).
  * Clifford and Clifford [2018] P. Clifford and R. Clifford, The classical complexity of boson sampling, in  _Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms_ (SIAM, 2018) pp. 146–155.
  * Kruse  _et al._ [2019] R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Detailed study of Gaussian boson sampling, Phys. Rev. A 100, 032326 (2019).
  * Drummond  _et al._ [2022] P. D. Drummond, B. Opanchuk, A. Dellios, and M. D. Reid, Simulating complex networks in phase space: Gaussian boson sampling, Physical Review A 105, 012427 (2022).
  * Lam  _et al._ [2015] S. K. Lam, A. Pitrou, and S. Seibert, Numba: A LLVM-based python JIT compiler, in  _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_ (2015) pp. 1–6.
  * Bacco  _et al._ [2021] D. Bacco, J. F. F. Bulmer, M. Erhard, M. Huber, and S. Paesani, Proposal for practical multidimensional quantum networks, Phys. Rev. A 104, 052618 (2021).