# Gaussian boson sampling with click-counting detectors

Gabriele Bressanini  QOLS, Blackett Laboratory, Imperial College London,
London SW7 2AZ, United Kingdom  Hyukjoon Kwon  Korea Institute for Advanced
Study, Seoul 02455, South Korea  M.S. Kim  QOLS, Blackett Laboratory, Imperial
College London, London SW7 2AZ, United Kingdom  Korea Institute for Advanced
Study, Seoul 02455, South Korea

###### Abstract

Gaussian boson sampling constitutes a prime candidate for an experimental
demonstration of quantum advantage within reach with current technological
capabilities. The original proposal employs photon-number-resolving detectors,
however the latter are not widely available. On the other hand, inexpensive
threshold detectors can be combined into a single click-counting detector to
achieve approximate photon number resolution. We investigate the problem of
sampling from a general multi-mode Gaussian state using click-counting
detectors and show that the probability of obtaining a given outcome is
related to a new matrix function which is dubbed as the Kensingtonian. We show
how the latter relates to the Torontonian and the Hafnian, thus bridging the
gap between known Gaussian boson sampling variants. We then prove that, under
standard complexity-theoretical conjectures, the model can not be simulated
efficiently.

##  I Introduction

Boson sampling, a computational problem introduced by Aaronson and Arkhipov
[1] and conjectured to be hard to simulate on a classical machine [2],
constitutes a prime candidate for an experimental proof of quantum advantage
using photons. In its original formulation, the task consists of sampling from
the output state of a passive linear optical network (LON) fed with single-
photon input states. Several variants of the task that lie in the same
complexity class have been proposed since. These are usually devised by
considering different classes of input states, such as photon-added coherent
states [3], photon-added or photon-subtracted squeezed vacuum states [4] and,
more recently, non-Gaussian input states involving Kerr-type non-linearities
[5]. Most notably, Gaussian boson sampling (GBS) [6] avoids the experimental
hurdles of generating indistinguishable single-photon Fock states by using
squeezed states of light as the non-classical resource needed to show quantum
advantage [7, 8]. On top of being a more experimentally feasible alternative
to prove quantum advantage on photonic platforms, GBS finds application in
simulating molecular vibronic spectra [9], predicting molecular docking
configurations for drug design [10] and in graph-related problems such as
finding dense subgraphs [11] and perfect matchings counting [12].

Boson sampling variants may also be designed by considering different kinds of
detection, such as Gaussian measurements [13] and photo-counting detection
schemes. Focusing on the latter, two classes of GBS experiments have been
investigated thus far: the initial proposal of GBS [14] makes use of photon-
number resolving (PNR) detectors, and only a couple of years later a GBS
implementation utilizing threshold (on/off) detectors was suggested, as a less
experimentally demanding alternative to show quantum advantage [15]. We remind
the reader that on/off detectors can only detect the presence or absence of
quantum light, while PNR detectors are, in principle, able to perfectly
distinguish between any Fock states. Threshold detectors, such as
superconducting nanowires or avalanche photo-diodes that can be operated at
room temperature [16], are widely available and inexpensive. Hence, to lessen
experimental challenges, most GBS experiments up to date employed threshold
detection [17, 18] and only more recently, a 216-mode Gaussian boson sampler
utilizing PNR detectors was used to claim quantum computational advantage
[19].

On the other hand, improved classical algorithms that exploit photon
collisional events to reduce the simulation’s overhead of GBS experiments
employing thresholds detectors [20] as well as the proposal of new classical
spoofing strategies [20, 21, 22] motivate the development of larger scale
experiments of increasing computational complexity. To this end, the
introduction of PNR detectors into the GBS scenario grants access to detection
events with a much larger total photon number (thus increasing the sample
space size exponentially), a regime which remains unachievable using on/off
detectors. As true PNR detectors are not always accessible, multiple on/off
detectors are routinely combined in a multiplexed fashion to achieve
approximate photon number resolution. The underlying idea of click-counting
detection [23] consists in dividing incoming light into weaker signals that
are then measured with threshold detectors, the final measurement output
simply being the number of on/off detectors that registered the presence of
photons. Hence, click-counting detectors can be seen as an intermediate case
between threshold and PNR detectors. We also recall that photon-number
resolution is often needed in GBS applications, e.g. to study higher-energy
molecular vibronic transitions [9, 24].

In this work, we study the problem of sampling from a generic multi-mode
Gaussian state using click-counting detectors, bridging the gap between GBS
experiments utilizing PNR and threshold detectors, which can be seen as
special instances of click-counting GBS. Besides fundamental interest,
Gaussian boson samplers employing click-counting detectors have a clear
experimental appeal, as it provides an easier way to achieve approximate
photon number resolution. In particular, we provide a closed-form expression
for the probability of observing a given click-pattern outcome and show that
it is related to a new matrix function which is dubbed as the Kensingtonian.
The latter plays an analogous role to the Hafnian and Torontonian in GBS
variants employing PNR and on/off detection, respectively. We show that, when
the probability of observing two or more photons in each output mode is
negligible, our model can not be efficiently simulated using a classical
machine under standard complexity-theoretic conjectures, thus making the setup
suitable to prove quantum advantage. In Table 1 we present the matrix
functions needed to compute the output probability distribution obtained when
sampling from a Gaussian state using photo-counting measurements.

This paper is structured as follows. In Sec. II we revise Gaussian states and
their representations through phase-space quasi-probability distributions. In
Sec. III we review the theoretical aspects of GBS experiments employing PNR
and threshold detectors. In Sec. IV the click-counting detection scheme is
introduced. Sections V and VI are dedicated to the main findings of this work:
we investigate the problem of sampling from a Gaussian state using click-
counting detectors, we give the definition of the Kensigtonian and we prove
the computational complexity of the sampling task. Lastly, in Sec. VII we draw
conclusions and give some final remarks.

Detection | zero-mean Gaussian state | displaced Gaussian state
---|---|---
on/off | Torontonian | loop Torontonian
click | Kensingtonian* | loop Kensingtonian*
PNR | Hafnian | loop Hafnian
Table 1: Matrix functions used to compute the output probability distribution
obtained from measuring multi-mode Gaussian states with on/off, click-counting
and PNR detection. The symbol * denotes the functions that are first
introduced in this paper.

##  II Gaussian States

In this section we briefly review the key aspects of Gaussian states and their
representation in terms of phase-space quasi-probability distributions (PQDs)
[25]. Let us consider a continuous variables system made up of $M$-bosonic
modes described by annihilation operators $a_{j}$ that satisfy the standard
commutation relations $[a_{j},a^{\dagger}_{k}]=\delta_{jk}$. We can then
introduce quadrature operators for each mode, defined as
$q_{j}=({a_{j}+a_{j}^{\dagger}})/{\sqrt{2}}$ and
$p_{j}=({a_{j}-a_{j}^{\dagger}})/{i\sqrt{2}}$ (where we have set $\hbar=1$),
and arrange them into the following vector

$\bm{r}=(q_{1},p_{1},\dots,q_{M},p_{M})^{\intercal}\,.$ (1)

A Gaussian state $\rho$ is completely characterized by its vector of first
moments (displacement) $\bm{\overline{r}}$ and its covariance matrix $\sigma$,
defined respectively as

$\overline{\bm{r}}=\Tr{\rho\bm{r}}\,,$ (2)
$\sigma=\Tr{\rho\\{\bm{r}-\bm{\overline{r}},(\bm{r}-\bm{\overline{r}})^{\intercal}\\}}\,.$
(3)

One can show that the $\bm{s}-$ordered PQD for a generic $M$-mode Gaussian
state with covariance matrix $\sigma$ and vector of first moments
$\bm{\alpha}$ is given by [5]

$W_{\rho}^{(\bm{s})}(\bm{\beta})=\frac{2^{M}}{\pi^{M}\sqrt{\det{\bm{\sigma}-\bm{\tilde{s}}}}}e^{-2(\bm{\beta}-\bm{\alpha})^{\intercal}(\bm{\sigma}-\bm{\tilde{s}})^{-1}(\bm{\beta}-\bm{\alpha})}\,.$
(4)

Here $\bm{s}=(s_{1},\dots,s_{M})^{\intercal}$ is the vector of operator
orderings, where $s_{j}\in\mathbb{R}$, and the matrix $\tilde{\bm{s}}$ is
defined as

$\bm{\tilde{s}}=\bigoplus_{j=1}^{M}s_{j}\mathbb{I}_{2}\,.$ (5)

The conventions used in Eq. (4) are such that for a single-mode coherent state
$\ket{\gamma}$ the covariance matrix is the identity matrix
$\sigma=\mathbb{I}_{2}$ and the vector of first moments reads
$\bm{\alpha}=(\Re{\gamma},\Im{\gamma})$. Note that Eq. (4) is well-defined iff
$\bm{\sigma}-\bm{\tilde{s}}\geq 0$, otherwise the $\bm{s}-$PQD becomes more
singular than a delta function. The well-known Husimi Q-function, Wigner
function and Glauber-Sudarshan P-function are retrieved respectively for
$\bm{s}=-\mathbb{I}_{M}$, $\bm{s}=0$ and $\bm{s}=\mathbb{I}_{M}$.

##  III GBS with PNR and threshold detectors

In this section we briefly review the two most widely studied and
experimentally relevant instances of GBS, i.e. the original proposal, which
consists of sampling from a Gaussian state using PNR detectors, and its
variant that employs more readily available threshold detectors. Both
implementations are thought to have the potential to show quantum advantage
(at least in some regimes), as sampling from their output probability
distribution using classical algorithms can be shown to be a computationally
hard task.

We remind the reader that a quantum measurement can be described in terms of a
positive operator-valued measure (POVM), whose elements $\\{\Pi_{x}\\}$
satisfy the conditions $\Pi_{x}\geq 0$ and $\sum_{x}\Pi_{x}=\mathcal{I}$ ,
where $\mathcal{I}$ is the identity operator on the Hilbert space. The
probability $p(x)$ of obtaining a given outcome $x$ when measuring a state
$\rho$ is given by the Born rule, i.e. $p(x)=\Tr{\rho\Pi_{x}}$.

The POVM elements of PNR detection are the projectors on Fock states, namely
$\outerproduct{k}{k}$ with $k$ non-negative integer. This kind of detector
provides perfect photon number resolution and can thus distinguish between any
Fock states with no uncertainty. It was shown [6] that the photon-number
statistics obtained from an $M-$mode Gaussian state with a null vector of
first moments reads

$p(\bm{k})=\frac{\text{Haf}[XO_{(S)}]}{\sqrt{\det{\Sigma}}\,k_{1}!\dots
k_{M}!}\,,$ (6)

where $\bm{k}=(k_{1},\dots,k_{M})^{\intercal}$ are the detected photon numbers
$k_{i}\geq 0$, $\Sigma$ is the covariance matrix of the state’s Q-function,
$O_{(S)}$ is a matrix constructed from $O=1-\Sigma^{-1}$ according to the
detection outcome (see Ref. [14] for more details) and

$X=\begin{bmatrix}0&\mathbb{I}\\\ \mathbb{I}&0\end{bmatrix}\,.$ (7)

The matrix function appearing in Eq. (6) is called the Hafnian and is defined
as follows

$\text{Haf}[A]=\sum_{\mu\in\text{PMP}}\prod_{j=1}^{n}A_{\mu(2j-1),\mu(2j)}\,\,,$
(8)

where $A$ is a $2n\times 2n$ complex matrix, and PMP is the set of perfect
matching permutations. Eq. (6) can be generalized to the case of an $M-$mode
Gaussian state with non-zero displacement by introducing the loop Hafnian
matrix function [14].

On the other hand, threshold detectors can only distinguish between the
absence or presence of light, without being able to resolve the photon content
of the state. In this case $k_{i}\in\\{0,1\\}$ and the POVM elements are
$\Pi_{0}=\outerproduct{0}{0}$ and $\Pi_{1}=\mathcal{I}-\Pi_{0}$. It was shown
(see Ref. [15] for details) that the probability distribution obtained from
measuring a $M-$mode Gaussian state with threshold detectors reads

$p(\bm{k})=\frac{\text{Tor}[O_{(S)}]}{\sqrt{\det{\Sigma}}}\,,$ (9)

where

$\text{Tor}[A]=\sum_{Z\in
P([n])}(-1)^{|Z|}\frac{1}{\sqrt{\det{\mathbb{I}-A_{(Z)}}}}$ (10)

is the Torontonian of a $2n\times 2n$ matrix $A$, $P([n])$ is the power set of
$[n]=\\{1,\dots,n\\}$ and $A_{(Z)}$ denotes a matrix constructed from $A$ by
eliminating rows and columns according to the set $Z$. Once again, it is
possible to lift the zero-displacement constraint and generalize Eq. (9) by
introducing the loop Torontonian matrix function [26].

A relation between the Torontonian and the Hafnian can be established by
noticing that the probability of obtaining a specific output pattern in a GBS
experiment employing threshold detectors can also be computed by summing the
probabilities of all detection events of a GBS experiment using PNR detectors
that are compatible with that specific click pattern. In particular, we say a
PNR detection pattern $\bm{k}^{\prime}$ is compatible with a threshold
detection pattern $\bm{k}$ if $k_{i}=0\implies k^{\prime}_{i}=0$ and
$k_{i}=1\implies k^{\prime}_{i}\geq 1$. This observation leads to the
following identity [15]

$\text{Haf}[XO]=\frac{1}{n!}\frac{d^{n}}{d\eta^{n}}\text{Tor}[\eta
O]\Big{|}_{\eta=0}\,,$ (11)

where $O$ is a $2n\times 2n$ matrix.

##  IV Click counting detection

In this section we describe the click-detection scheme introduced in Ref.
[23], where a multiplexing setup and on/off detectors are employed to achieve
approximate photon-number resolution. The main idea consists in splitting the
incoming state into weaker signals (distributing the intensity uniformly among
the output ports of the interferometer) that will then be measured using
threshold detectors.

Figure 1: Schematics of a click-counting detector made up of $N$ threshold
detectors. The quantum state $\rho$ to be measured enters an $N-$mode
interferometer (represented by a unitary operation $\mathcal{U}(N)$) that
splits the intensity equally among its output ports. $N$ on/off detectors then
measure the output state and the detection results are summed into the final
signal.

Note how we are only interested in the total number of recorded clicks $k$,
and not in the specific output pattern obtained in a given measurement. In
Fig. (1) we display a schematic representation of a click-counting detector.
We remind the reader that the POVM elements of on/off detection are

$\Pi_{0}^{(1)}=\outerproduct{0}{0}=:e^{-\hat{n}}:\quad\quad\Pi_{1}^{(1)}=\mathcal{I}-\Pi_{0}^{(1)}\,,$
(12)

where $\mathcal{I}$ is the identity operator, $\hat{n}$ is the number operator
and $:\bullet:$ denotes normal ordering of bosonic operators [27]. One can
show that the POVM elements corresponding to a click-detector made up of $N$
on/off detectors are given by

$\Pi_{k}^{(N)}=:\binom{N}{k}e^{-\frac{N-k}{N}\hat{n}}(1-e^{-\frac{\hat{n}}{N}})^{k}:\,,$
(13)

where $k$, i.e. the number of recorded clicks, may vary between $0$ and $N$.
Note how we always refer to clicks and not to photons. As expected, for $N=1$
we retrieve the POVM elements of threshold detection, i.e.

$\Pi_{k}^{(1)}=:(e^{\hat{n}}-\mathcal{I})^{k}e^{-\hat{n}}:$ (14)

with $k=0,1$. This formalism also allows us to easily introduce imperfections,
such as sub-unit efficiency and finite dark count rate. To do so, we simply
have to consider a click detector made up of noisy on/off detectors, whose
POVM is obtained from Eq. (12) by substitution of the number operator
$\hat{n}$ with a suitable response function whose functional form depends on
the specific noise model considered. In this paper we consider the simple and
widely used substitution $\hat{n}\mapsto\eta\hat{n}+\nu$, where $0\leq\eta\leq
1$ and $\nu\geq 0$ are the efficiency and the dark count rate of the threshold
detector, respectively, to obtain the following POVM elements

$\Pi_{0}^{(1)}=:e^{-(\eta\hat{n}+\nu)}:\quad\quad\Pi_{1}^{(1)}=\mathcal{I}-\Pi_{0}^{(1)}\,.$
(15)

One can then prove that the POVM elements of noisy click-counting detection
are given by

$\Pi_{k}^{(N)}=\,:\binom{N}{k}\left[e^{-\left(\eta\frac{\hat{n}}{N}+\nu\right)}\right]^{N-k}\left[1-e^{-\left(\eta\frac{\hat{n}}{N}+\nu\right)}\right]^{k}:\,.$
(16)

We note that this operators can be obtained from Eq. (13) upon substituting
$\hat{n}\mapsto\eta\hat{n}+N\nu$. This is intuitively clear: we expect the
click-counting detector to “inherit” the inefficiency of the threshold
detectors, however the dark count rates from each on/off detector will add up
to the “total” dark count rate of the click-counting detector. The reason for
this is that dark counts coming from different threshold detectors can be
thought as independent Poisson variables and it is then well known that the
sum of Poisson random variables is still a Poisson variable, whose mean value
is the sum of the addends’ mean values. This also means that the performance
of a noisy click-counting detector made up of $N$ threshold detectors
characterized by Eq. (15) should be compared to that of a PNR detector with
sub-unit efficiency $\eta$ and dark count rate $N\nu$, whose related POVM
elements read

$\tilde{\Pi}_{k}=:\frac{(\eta\hat{n}+N\nu)^{k}}{k!}e^{-(\eta\hat{n}+N\nu)}:$
(17)

In Appendix D we show that $-$ as one might intuitively expect $-$ in the
$N\rightarrow\infty$ limit and in the absence of noise we retrieve true PNR
detection, i.e.

$\lim_{N\rightarrow\infty}\Pi_{k}^{(N)}=\,:\frac{\hat{n}^{k}}{k!}e^{-\hat{n}}:\,\equiv\outerproduct{k}{k}\,.$
(18)

The noisy case requires some additional attention, as the presence of non-zero
dark count rate $\nu$ causes the expression to diverge if we are to take the
formal limit. In practice, typical values of a threshold detector’s dark count
rate are of the order $\nu\simeq 10^{-4}$, hence in regimes where
simultaneously $N\gg 1$ and $N\nu\ll 1$ we still expect good convergence of
Eq. (16) to Eq. (17).

##  V The Kensingtonian

In this section we derive the outcome probability distribution of an $M-$mode
GBS experiment employing click-counting detectors, each made up of $N$ ideal
threshold detectors. Note that for $N=1$ we are describing GBS with threshold
detectors, and the probability distribution of outcomes is given by Eq. (9).
On the other hand, in the $N\rightarrow\infty$ limit we retrieve photon-
number-resolving GBS with PNR detection, and the probability distribution
converges to Eq. (6). In this section we interpolate between these two special
cases by deriving a closed formula valid for general $N$.

In click-counting GBS, a single detection event is denoted by
$\bm{k}=(k_{1},\dots,k_{M})$, where $0\leq k_{i}\leq N$ $\forall i$. Note that
the total number of clicks $n=\sum_{i}k_{i}$ is not fixed, as the number of
photons entering the interferometer is not determined in the first place. We
want to compute

$p(\bm{k})=\Tr{\rho\Pi_{\bm{k}}^{(N)}}\,,$ (19)

where $\rho$ is a generic $M$-mode Gaussian state with covariance matrix
$\sigma$ and null vector of first moments and $\Pi_{\bm{k}}^{(N)}$ is the POVM
that characterize the $M$-mode click detection. The latter simply reads

$\Pi_{\bm{k}}^{(N)}=\bigotimes_{i=1}^{M}\Pi_{{k_{i}}}^{(N)}\,,$ (20)

where $\Pi_{{k_{i}}}^{(N)}$ is the POVM of a single click-counting detector,
given by Eq. (13). The latter can also be expressed, by virtue of the binomial
theorem, as

$\Pi_{k_{i}}^{(N)}=\binom{N}{k_{i}}\sum_{\ell_{i}=0}^{k_{i}}\binom{k_{i}}{\ell_{i}}(-1)^{\ell_{i}}:e^{-\frac{N-k_{i}+\ell_{i}}{N}\hat{n}}:\,.$
(21)

Note that the operator $:e^{-\frac{N-k_{i}+\ell_{i}}{N}\hat{n}}:$ corresponds
to the vacuum element $\Pi_{0}^{(1)}$ of the noisy on/off detection POVM Eq.
(15), with an effective detection inefficiency given by
$\lambda_{i}\equiv{(N-k_{i}+\ell_{i})}/{N}\in[0,1]$. This observation is
crucial, as it implies that all is needed in order to obtain $p(\bm{k})$ are
(noisy) vacuum statistics of marginal states of $\rho$. Each of these
contributions can be computed efficiently using the Gaussian formalism, hence
$-$ as we will see $-$ the complexity of the sampling task arises from the
exponential number of terms appearing in the expression of the probability. We
can express the latter as follows

$p(\bm{k})=\pi^{M}\int
d^{2M}\bm{\beta}\,Q_{\rho}(\bm{\beta})\,P_{\Pi_{\bm{k}}^{(N)}}(\bm{\beta})\,,$
(22)

where $Q_{\rho}(\bm{\beta})$ is the Husimi Q-function of $\rho$ and
$P_{\Pi_{\bm{k}}^{(N)}}(\bm{\beta})$ is the P-function of the POVM element.
Using Eq. (4) with ordering parameter $s_{j}=-1\,\forall j$ we obtain

$Q_{\rho}(\bm{\beta})=\frac{1}{\pi^{M}\sqrt{\det{{\Sigma}}}}e^{-\bm{\beta}^{\intercal}{\Sigma}^{-1}\bm{\beta}}\,,$
(23)

where we have introduced the matrix ${\Sigma}=({{\sigma}+\mathbb{I}})/{2}$ for
convenience. On the other hand, the P-function of $\Pi_{\bm{k}}^{(N)}$ is
readily obtained once we know how to compute the P-function of
$:e^{-\lambda_{i}\hat{n}}:$ . One can prove that for $\lambda_{i}\neq 1$

$P_{:e^{-\lambda_{i}\hat{n}}:}(\beta_{1},\beta_{2})=\frac{1}{\pi(1-\lambda_{i})}e^{-\frac{\lambda_{i}}{1-\lambda_{i}}(\beta_{1}^{2}+\beta_{2}^{2})}\,.$
(24)

while for $\lambda_{i}=1$ (i.e. for $\ell_{i}=k_{i}$) we have

$P_{:e^{-\hat{n}}:}(\beta_{1},\beta_{2})=\delta(\beta_{1})\delta(\beta_{2})\,,$
(25)

where $\beta_{1}$ and $\beta_{2}$ are the two (real) Cartesian coordinates of
the complex plane. After cumbersome calculations (see Appendix A for a
detailed derivation) it is possible to show that the probability of observing
a given click pattern $\bm{k}$ reads

$p(\bm{k})=\frac{\text{Ken}[O]}{\sqrt{\det{\Sigma}}}\,.$ (26)

Here $O=\mathbb{I}-\Sigma^{-1}$ and

$\begin{split}\text{Ken}[A]&=\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}(-1)^{k_{i}-d_{i}}\right]\\\
&\prod_{j\notin
Z}\left[\frac{N}{d_{j}}\right]\frac{1}{\sqrt{\det{(\mathbb{I}-A)_{(Z)}+D_{Z}}}}\,\end{split}$
(27)

is the Kensingtonian of a $2M\times 2M$ matrix $A$. $D_{Z}$ is a diagonal
matrix defined as

$D_{Z}=\bigoplus_{i\notin
Z}\left(\frac{N-d_{i}}{d_{i}}\right)\mathbb{I}_{2}\,,$ (28)

and the set $Z$ is defined as follows

$Z=\\{i|1\leq i\leq M,\,d_{i}=0\\}\,.$ (29)

In Eq. (27) we have used the subscript notation $(\mathbb{I}-A)_{(Z)}$ to
denote the matrix obtained from $(\mathbb{I}-A)$ by eliminating the
corresponding rows/columns according to the set $Z$. In particular, if
$Z=\\{a,b\\}$ we remove the rows(columns) numbered $2a-1,2a,2b-1$ and $2b$. In
Appendix C we generalize Eq. (26) to displaced Gaussian input states by
introducing the loop Kensingtonian matrix function.

For $N=1$, Eq. (26) must coincide with the probability distribution of a GBS
experiment employing threshold detectors Eq. (9). In Appendix B we explicitly
prove the following identity that links the Kensingtonian and the Torontonian

$\text{Ken}[A]\stackrel{{\scriptstyle\small{N=1}}}{{=}}\text{Tor}[A_{(\mathcal{K})}]\,.$
(30)

On the other hand, the outcome probability distribution of click-counting GBS
Eq. (26) converges to that of PNR GBS Eq. (6) in the $N\rightarrow\infty$
limit, as a consequence of the fact that in that same limit the click-counting
detection POVM converges to that of true PNR detection. This in turn implies
that the Hafnian can be retrieved as a limiting case of the Kensingtonian. We
leave it as an open question whether this can be exploited to develop faster
algorithms to approximate the Hafnian of a matrix.

One might also be interested in knowing the regimes where the photon counting
statistics coming from a GBS experiment employing click-counting detectors
well approximates Eq. (6). Of course, the convergence of the click-counting
detector’s POVM to that of the PNR detector implies that we can always
increase $N$ to make the total variational distance (TVD) between the two
distributions arbitrarily small. However, in what follows we want to give a
more quantitative and practical indication of the number of threshold
detectors needed in order to have a good agreement between Eq. (26) and Eq.
(6). This is particularly relevant for GBS experiments aimed at applications
(rather than those aimed at proving quantum advantage), where the Hafnian’s
specific functional form and properties are crucial to map the sampling task
onto problems in graph theory and chemistry.

We remind the reader that, on average, a Haar-random LON equally distributes
the intensity among its output modes, and denote with $\overline{n}$ the
average photon density per mode. Most GBS experiments up to date, especially
those targeting applications, have modest values of photon densities per
output mode, usually $\overline{n}<1$. Furthermore, tracing out all output
modes but one leaves us with a Gaussian marginal state that approximately
looks like a thermal state $\nu_{th}(\overline{n})$. We thus restrict
ourselves to this single-mode case and numerically compute the TVD to give a
rough estimate of how many threshold detectors making up a click-detector are
needed so that the latter well approximates an ideal PNR detector in the
energy regime of interest. The two probability distributions needed to compute
this TVD are given by Eq. (6) and Eq. (26), with $M=1$ and
$\sigma=(2\overline{n}+1)\mathbb{I}_{2}$. In Fig. (2) we show how using just
$N=8$ threshold detectors (corresponding to 3 multiplexing steps, if we are
using the simple scheme where intensity is halved at each layer made up of
balanced beamsplitters) already provides a good agreement between the two
distributions.

Figure 2: The total variational distance (TVD) between the probability
distributions coming from a PNR and click-counting detector made up of $N=8$
threshold detectors, as a function of the mean photon number $\overline{n}$ of
the probe thermal state $\nu_{th}(\overline{n})$. In the energy regime of
interest we considered there is a good agreement, with TVD smaller than 0.05.

##  VI Complexity of click-GBS

In this section we prove that sampling from the output probability
distribution of a Gaussian boson sampler utilizing click-counting detectors
constitutes a hard problem to solve on a classical computer. Proofs of
hardness of the Boson sampling task and its variants rely on the assumption
that the probability of observing a collision (i.e. two or more photons) in
any of the interferometer’s output modes is negligible. This is also known as
the non-collisional regime. In Ref. [1] the authors proved that the
probability $\varepsilon$ of observing at least a collision at the output of
an $M-$mode Haar-random LON fed with $n$ single photons can be bounded as
follows

$\expectationvalue{\varepsilon}_{\mathcal{U}}\leq\frac{2n^{2}}{M}\,,$ (31)

where $\expectationvalue{\cdot}_{\mathcal{U}}$ denotes averaging over $M\times
M$ Haar-random unitary matrices. In GBS the number of photons is not fixed,
hence the squared number of photons $n^{2}$ on right-hand side of Eq. (31)
needs to be further averaged according to the photon number probability
distribution of the Gaussian state. This means that it is possible to set the
collision probability to be a small constant by choosing the scaling of the
number of modes appropriately.

Intuitively, in this regime, the statistics coming from PNR, threshold and
click-counting detectors should all behave similarly. In the following we
rigorously formalize this idea, by adjusting the arguments used in Ref. [15]
to prove that sampling from a Gaussian boson sampler with on/off detectors is
a computationally-hard problem. Roughly speaking, the authors proved that $-$
when collision probability is small $-$ an approximation of threshold GBS
would also constitute a good approximation of GBS with PNR detectors, an event
considered to be unlikely as it would imply the collapse of the polynomial
hierarchy to the third level. Hence, it is concluded that $-$ under standard
complexity-theoretic arguments $-$ threshold GBS is in the same complexity
class as the original GBS proposal.

It thus suffices to show that in the non-collisional regime the output
probability distribution $p$ of click-counting GBS is arbitrarily close to
that of threshold GBS, which we denote with the symbol $\tilde{p}$. We define
the set of collision events of an $M-$mode sampling problem as

$\mathcal{C}=\\{\bm{k}=(k_{1},\dots,k_{M})\,|\,k_{i}>1\text{ for some
}i\in\\{1,\dots,M\\}\\}\,.$ (32)

By definition we have that $\tilde{p}(\bm{k}\in\mathcal{C})=0$, as threshold
detectors can either click once or not click at all. On the other hand, the
probability of observing a collision event for click-counting GBS reads

$\varepsilon=\sum_{\bm{k}\in\mathcal{C}}p(\bm{k})\,.$ (33)

The total variational distance between $p$ and $\tilde{p}$ reads

$\begin{split}||p-\tilde{p}||_{1}&=\frac{1}{2}\sum_{\bm{k}}|p(\bm{k})-\tilde{p}(\bm{k})|\\\
&=\frac{1}{2}\sum_{\bm{k}\in\mathcal{C}}|p(\bm{k})-\tilde{p}(\bm{k})|+\frac{1}{2}\sum_{\bm{k}\notin\mathcal{C}}|p(\bm{k})-\tilde{p}(\bm{k})|\\\
&=\frac{1}{2}\sum_{\bm{k}\in\mathcal{C}}p(\bm{k})+\frac{1}{2}\sum_{\bm{k}\notin\mathcal{C}}|p(\bm{k})-\tilde{p}(\bm{k})|\\\
&=\frac{\varepsilon}{2}+\frac{1}{2}\sum_{\bm{k}\notin\mathcal{C}}|p(\bm{k})-\tilde{p}(\bm{k})|\,.\end{split}$
(34)

By close inspection of the click-counting detection POVM Eq. (13) we notice
that $\Pi^{(N)}_{0}=:e^{-\hat{n}}:=\outerproduct{0}{0}\quad$ for every value
of $N$. Using this and the fact that the elements of a POVM resolve the
identity we can write

$\sum_{k=1}^{N}\Pi_{k}^{(N)}=\mathcal{I}-\outerproduct{0}{0}=\Pi_{1}^{(1)}\,.$
(35)

This implies that the probability of a threshold detector clicking is equal to
the probability of a click detector detecting any number of clicks between $1$
and $N$. Generalizing this to $M$ modes we obtain

$\tilde{p}(\bm{k})=p(\bm{k})+\sum_{\bm{k}^{\prime}\in\mathcal{C}_{\bm{k}}}p(\bm{k}^{\prime})\,,$
(36)

where $\bm{k}\notin\mathcal{C}$ is a collision-less detection event and
$\mathcal{C}_{\bm{k}}$ is the set of all possible collision events compatible
with $\bm{k}$

$\mathcal{C}_{\bm{k}}=\\{\bm{k}^{\prime}\in\mathcal{C}\,|\,\text{if
}k_{i}=0\implies k_{i}^{\prime}=0\\}\,.$ (37)

If we now substitute Eq. (36) into Eq. (34) we obtain

$\begin{split}||p-\tilde{p}||_{1}&=\frac{\varepsilon}{2}+\frac{1}{2}\sum_{\bm{k}\notin\mathcal{C}}|p(\bm{k})-\tilde{p}(\bm{k})|\\\
&=\frac{\varepsilon}{2}+\frac{1}{2}\sum_{\bm{k}\notin\mathcal{C}}\sum_{\bm{k}^{\prime}\in\mathcal{C}_{\bm{k}}}p(\bm{k}^{\prime})\\\
&=\frac{\varepsilon}{2}+\frac{1}{2}\sum_{\bm{k}^{\prime}\in\mathcal{C}}p(\bm{k}^{\prime})\\\
&=\varepsilon\,.\end{split}$ (38)

Let us now consider a probability distribution $\pi$ that approximates $p$
arbitrarily well, i.e. it satisfies $||p-\pi||_{1}=\varepsilon^{\prime}$.
Using the triangle inequality of the $L^{1}$ norm we can bound
$||\tilde{p}-\pi||_{1}$ as follows

$||\tilde{p}-\pi||_{1}=||\tilde{p}-p+p-\pi||_{1}\leq||\tilde{p}-p||_{1}+||p-\pi||_{1}=\varepsilon+\varepsilon^{\prime}\,.$
(39)

If we now assume that there exists a polynomial time algorithm that can sample
from $\pi$, then Eq. (39) tells us that the same algorithm can sample
efficiently from an arbitrarily good approximation of $\tilde{p}$, thus
causing the collapse of the polynomial hierarchy to the third level and
concluding our proof of hardness for click-counting GBS. In Appendix E we
outline an alternative proof of hardness that does not make use of the output
probability distribution of a Gaussian boson sampler employing threshold
detectors.

We can now discuss the time complexity of computing the Kensingtonian
according to its definition Eq. (27). Let us first call $n=\sum_{i}k_{i}$ the
total number of clicks. It is clear that, similarly to the Torontonian, the
complexity of computing the Kensingtonian arises from the number of
determinant contributions one needs to evaluate, which is given by

$\sum_{0\leq\bm{d}\leq\bm{k}}=\prod_{i=1}^{M}(k_{i}+1)\equiv F(\bm{k})\,,$
(40)

where $0\leq k_{i}\leq N$. For $N=1$ (threshold detectors) we have
$k_{i}=0,1$, hence $\prod_{i=1}^{M}(k_{i}+1)=2^{n}$. In this special case the
result depends solely on the total number of clicks $n$, however one can
easily see that this is not the case for $N>1$, as we would need and
increasing number of functions of all the $k_{i}$ to completely characterize
the expression, making it impractical. An exact, simple, closed formula is
thus out of reach, however in the following we show that, on average, we still
obtain an exponential scaling in $n$.

The average value of $F(\bm{k})$ at fixed $n$ can be expressed as a sum of a
collision-less term and collision term, namely

$\langle
F\rangle_{n}=\sum_{\bm{k}\notin\mathcal{C}_{n}}p(\bm{k})F({\bm{k}})+\sum_{\bm{k}\in\mathcal{C}_{n}}p(\bm{k})F({\bm{k}})\,.$
(41)

Here $\mathcal{C}_{n}$ is the set of collision detection events at a fixed
total number of clicks and we use the symbol $\expectationvalue{\cdot}_{n}$ to
denote averaging at constant $n$. When there are no collisions we simply have
that $F(\bm{k})=2^{n}$. We can thus write

$\langle
F\rangle_{n}=(1-\tilde{\varepsilon})2^{n}+\sum_{\bm{k}\in\mathcal{C}_{n}}p(\bm{k})F({\bm{k}})\,,$
(42)

where $\tilde{\varepsilon}$ is the probability of observing a collision at
fixed $n$, namely

$\tilde{\varepsilon}=\sum_{\bm{k}\in\mathcal{C}_{n}}p(\bm{k})\,.$ (43)

It can easily be seen that $F(\bm{k})\geq(n+1)$, hence

$\langle
F\rangle_{n}\geq(1-\tilde{\varepsilon})2^{n}+\tilde{\varepsilon}(n+1)\,,$ (44)

which confirms the exponential scaling of the number of terms to be evaluated
in Eq. (27) with the total number of clicks. Analogously, we can upper bound
$F(\bm{k})$ using the inequality of the arithmetic and geometric means and
show that the following chain of inequalities holds

$(1-\tilde{\varepsilon})2^{n}+\tilde{\varepsilon}(n+1)\,\leq\langle
F\rangle_{n}\leq(1-\tilde{\varepsilon})2^{n}+\tilde{\varepsilon}e^{n}\,.$ (45)

We recall that the determinants present in the definition of the Kensingtonian
can be computed efficiently using the standard algorithm based on the Cholesky
decomposition, whose time complexity scales with the cube of the matrix
dimension. Hence it follows that a direct evaluation of Eq. (27), in the non-
collisional regime, leads to a complexity upper bounded by $O(n^{3}2^{n})$. We
emphasize that we do not claim optimality, and we leave it as an open question
the possibility of exploiting the structure of the Kensingtonian to find
faster algorithms for its evaluation.

##  VII Conclusions

In this paper we have investigated the problem of sampling from Gaussian
states with click-counting detectors and found a closed-form expression for
the probability distribution. The latter is related to the Kensingtonian, a
new matrix function that plays an analogous role to the Hafnian and the
Torontonian in GBS experiments employing PNR and threshold detectors,
respectively. We then proved that, in the non-collisional regime, the problem
at study still gives rise to a computationally hard problem, intractable using
classical sampling algorithms, and showed how the Kensingtonian is related to
known matrix functions in limiting cases of interest.

Our work leaves some open questions. We recall that Eq. (26) converges to Eq.
(6) in the $N\gg 1$ limit, hence it would be interesting to investigate
whether the Kensingtonian’s structure could be exploited to design new
algorithms to approximate the Hafnian. Future efforts will also focus on
studying the classical simulability of a GBS experiment employing noisy click-
counting detection, where we envision the existence of a trade-off relation
between $N$, i.e. the number of threshold detectors making up a single click-
counting detector, and the noise parameters that characterize each on/off
detector, for the system to enter a regime where achieving quantum advantage
is not ruled out.

Note added $-$ As we are finalizing the manuscript we became aware of a recent
experiment [28] reporting an implementation of a Gaussian boson sampling
device employing unbalanced click-counting detectors, i.e. the intensity of
the incoming light is split unevenly among the threshold detectors that make
up a single click-counting detector. Our theoretical modeling and that
presented in Ref. [28] are consistent with each other, as they originate from
the same POVM describing a click-counting detector. In this paper, we have
exploited the structure of the POVM to obtain a closed analytical formula for
the outcome probability distribution, that can be readily applied to an
$M-$mode GBS task employing balanced click-counting detectors. In particular,
we showed that the probability of a particular detection outcome may be
obtained by computing the Kensingtonian of a $2M\times 2M$ matrix. On the
other hand, the authors of Ref. [28] model their experimental setup as larger
instance of a $MN-$mode GBS experiment employing on/off detectors. As a
result, they sum over the probabilities of (combinatorially-many) threshold-
detection outcomes that give rise to the same click-detection pattern. Each of
these probabilities requires the evaluation of the Torontonian of matrices
that are up to $2MN\times 2MN$ in dimension, and may therefore lead to a less
efficient computation of a click-detection outcome probability with respect to
the approach presented in this paper.

##  VIII Acknowledgments

G.B. is part of the AppQInfo MSCA ITN which received funding from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Sklodowska-Curie grant agreement No 956071. G.B. thanks Jan Sperling for
introducing him to click-counting detection and Francesco Viganò for helpful
discussions. The authors thank Changhun Oh for providing useful comments. H.
K. is supported by the KIAS Individual Grant No. CG085301 at Korea Institute
for Advanced Study. MSK acknowledges the KIST Open Research Programme, Samsung
GRC programme and the KIAS visiting professorship.

## References

  * 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_, STOC ’11 (Association for Computing Machinery, New York, NY, USA, 2011) p. 333–342.
  * Scheel  _et al._ [2003] S. Scheel, K. Nemoto, W. J. Munro, and P. L. Knight, Measurement-induced nonlinearity in linear optics, Phys. Rev. A 68, 032310 (2003).
  * Seshadreesan  _et al._ [2015] K. P. Seshadreesan, J. P. Olson, K. R. Motes, P. P. Rohde, and J. P. Dowling, Boson sampling with displaced single-photon fock states versus single-photon-added coherent states: The quantum-classical divide and computational-complexity transitions in linear optics, Phys. Rev. A 91, 022334 (2015).
  * Olson  _et al._ [2015] 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, Phys. Rev. A 91, 022317 (2015).
  * Bressanini  _et al._ [2022] G. Bressanini, H. Kwon, and M. S. Kim, Noise thresholds for classical simulability of nonlinear boson sampling, Phys. Rev. A 106, 042413 (2022).
  * Hamilton  _et al._ [2017] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian boson sampling, Phys. Rev. Lett. 119, 170501 (2017).
  * Deshpande  _et al._ [2022] A. Deshpande, A. Mehta, T. Vincent, N. Quesada, M. Hinsche, M. Ioannou, L. Madsen, J. Lavoie, H. Qi, J. Eisert, D. Hangleiter, B. Fefferman, and I. Dhand, Quantum computational advantage via high-dimensional gaussian boson sampling, Science Advances 8, eabi7894 (2022), https://www.science.org/doi/pdf/10.1126/sciadv.abi7894 .
  * Grier  _et al._ [2022] D. Grier, D. J. Brod, J. M. Arrazola, M. B. d. A. Alonso, and N. Quesada, The Complexity of Bipartite Gaussian Boson Sampling, Quantum 6, 863 (2022).
  * 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).
  * Banchi  _et al._ [2020] L. Banchi, M. Fingerhuth, T. Babej, C. Ing, and J. M. Arrazola, Molecular docking with gaussian boson sampling, Science Advances 6, eaax1950 (2020), https://www.science.org/doi/pdf/10.1126/sciadv.aax1950 .
  * Arrazola and Bromley [2018] J. M. Arrazola and T. R. Bromley, Using gaussian boson sampling to find dense subgraphs, Phys. Rev. Lett. 121, 030503 (2018).
  * Brádler  _et al._ [2018] K. Brádler, P.-L. Dallaire-Demers, P. Rebentrost, D. Su, and C. Weedbrook, Gaussian boson sampling for perfect matchings of arbitrary graphs, Phys. Rev. A 98, 032310 (2018).
  * Chakhmakhchyan and Cerf [2017] L. Chakhmakhchyan and N. J. Cerf, Boson sampling with gaussian measurements, Phys. Rev. A 96, 032326 (2017).
  * 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).
  * Quesada  _et al._ [2018] N. Quesada, J. M. Arrazola, and N. Killoran, Gaussian boson sampling using threshold detectors, Phys. Rev. A 98, 062322 (2018).
  * Hadfield [2009] R. H. Hadfield, Single-photon detectors for optical quantum information applications, Nature photonics 3, 696 (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).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, _et al._ , Quantum computational advantage using photons, Science 370, 1460 (2020).
  * Madsen  _et al._ [2022] L. S. Madsen, F. Laudenbach, M. F. Askarani, F. Rortais, T. Vincent, J. F. Bulmer, F. M. Miatto, L. Neuhaus, L. G. Helt, M. J. Collins, _et al._ , Quantum computational advantage with a programmable photonic processor, Nature 606, 75 (2022).
  * Bulmer  _et al._ [2022a] J. F. Bulmer, B. A. Bell, R. S. Chadwick, A. E. Jones, D. Moise, A. Rigazzi, J. Thorbecke, U.-U. Haus, T. Van Vaerenbergh, R. B. Patel, _et al._ , The boundary for quantum advantage in gaussian boson sampling, Science advances 8, eabl9236 (2022a).
  * Villalonga  _et al._ [2021] B. Villalonga, M. Y. Niu, L. Li, H. Neven, J. C. Platt, V. N. Smelyanskiy, and S. Boixo, Efficient approximation of experimental gaussian boson sampling, arXiv preprint arXiv:2109.11525 (2021).
  * Oh  _et al._ [2022] C. Oh, L. Jiang, and B. Fefferman, Spoofing cross entropy measure in boson sampling, arXiv preprint arXiv:2210.15021 (2022).
  * Sperling  _et al._ [2012] J. Sperling, W. Vogel, and G. S. Agarwal, True photocounting statistics of multiple on-off detectors, Phys. Rev. A 85, 023820 (2012).
  * Wang  _et al._ [2020] C. S. Wang, J. C. Curtis, B. J. Lester, Y. Zhang, Y. Y. Gao, J. Freeze, V. S. Batista, P. H. Vaccaro, I. L. Chuang, L. Frunzio, L. Jiang, S. M. Girvin, and R. J. Schoelkopf, Efficient multiphoton sampling of molecular vibronic spectra on a superconducting bosonic processor, Phys. Rev. X 10, 021060 (2020).
  * Serafini [2017] A. Serafini,  _Quantum continuous variables: a primer of theoretical methods_ (CRC press, 2017).
  * Bulmer  _et al._ [2022b] J. F. F. Bulmer, S. Paesani, R. S. Chadwick, and N. Quesada, Threshold detection statistics of bosonic states, Phys. Rev. A 106, 043712 (2022b).
  * Cahill and Glauber [1969] K. E. Cahill and R. J. Glauber, Ordered expansions in boson amplitude operators, Phys. Rev. 177, 1857 (1969).
  * Deng  _et al._ [2023] Y.-H. Deng, Y.-C. Gu, H.-L. Liu, S.-Q. Gong, H. Su, Z.-J. Zhang, H.-Y. Tang, M.-H. Jia, J.-M. Xu, M.-C. Chen, _et al._ , Gaussian boson sampling with pseudo-photon-number resolving detectors and quantum computational advantage, arXiv preprint arXiv:2304.12240 (2023).

##  Appendix A Click-counting probability

In this section we derive the probability of observing a given click-counting
pattern $\bm{k}=(k_{1},\dots,k_{M})$ when a generic $M-$mode Gaussian state
$\rho$ with covariance matrix $\sigma$ and null vector of first moments is
sampled using $M$ click-counting detectors, each composed of $N$ threshold
detectors. We remind the reader that each click-detector can measure at most
$N$ clicks. The probability reads

$p(\bm{k})=\Tr{\rho\Pi_{\bm{k}}^{(N)}}\,,$ (46)

where $\Pi_{\bm{k}}^{(N)}$ is the POVM element that characterize the $M$-mode
click detection

$\Pi_{\bm{k}}^{(N)}=\bigotimes_{i=1}^{M}\Pi_{{k_{i}}}^{(N)}\,.$ (47)

Here $\Pi_{{k_{i}}}^{(N)}$ is the POVM of a single click-counting detector and
is given by

$\Pi_{k_{i}}^{(N)}=:\binom{N}{k_{i}}e^{-\frac{N-k_{i}}{N}\hat{n}}(1-e^{-\frac{\hat{n}}{N}})^{k_{i}}:=\binom{N}{k_{i}}\sum_{\ell_{i}=0}^{k_{i}}\binom{k_{i}}{\ell_{i}}(-1)^{\ell_{i}}:e^{-\frac{N-k_{i}+\ell_{i}}{N}\hat{n}}:\,.$
(48)

Note that the operator $:e^{-\frac{N-k_{i}+\ell_{i}}{N}\hat{n}}:$ corresponds
to the vacuum element $\Pi_{0}^{(1)}$ of the POVM associated with noisy
threshold detection with effective detection inefficiency given by
$\lambda_{i}\equiv({N-k_{i}+\ell_{i}})/{N}\in[0,1]$. In particular,
$:e^{-\lambda_{i}\hat{n}}:$ represents an unnormalized thermal state for
$\lambda_{i}\in(0,1)$, the vacuum state for $\lambda_{i}=1$ and the identity
operator for $\lambda_{i}=0$, respectively. The probability can be expressed
as

$p(\bm{k})=\pi^{M}\int
d^{2M}\bm{\beta}\,Q_{\rho}(\bm{\beta})\,P_{\Pi_{\bm{k}}^{(N)}}(\bm{\beta})\,,$
(49)

where $Q_{\rho}(\bm{\beta})$ is the Husimi Q-function of $\rho$ and
$P_{\Pi_{\bm{k}}^{(N)}}(\bm{\beta})$ is the P-function of
$\Pi_{\bm{k}}^{(N)}$. Using Eq. (4) with ordering parameters
$s_{j}=-1\,\forall j$ we obtain

$Q_{\rho}(\bm{\beta})=\frac{1}{\pi^{M}\sqrt{\det{{\Sigma}}}}e^{-\bm{\beta}^{\intercal}{\Sigma}^{-1}\bm{\beta}}\,,$
(50)

where we have introduced the matrix ${\Sigma}=({{\sigma}+\mathbb{I}})/{2}$.
In order to compute the P-function of $\Pi_{\bm{k}}^{(N)}$ we just need the
P-function of the unnormalized thermal state $:e^{-\lambda_{i}\hat{n}}:$ .
Using Eq. (4) one can prove that for $\lambda_{i}\neq 1$

$P_{:e^{-\lambda_{i}\hat{n}}:}(\beta_{1},\beta_{2})=\frac{1}{\pi(1-\lambda_{i})}e^{-\frac{\lambda_{i}}{1-\lambda_{i}}(\beta_{1}^{2}+\beta_{2}^{2})}=\frac{N}{\pi(k_{i}-\ell_{i})}e^{-\frac{N-k_{i}+\ell_{i}}{k_{i}-\ell_{i}}(\beta_{1}^{2}+\beta_{2}^{2})}\,.$
(51)

For $\lambda_{i}=1$ (i.e. for $\ell_{i}=k_{i}$) we have
$P_{:e^{-\hat{n}}:}(\beta_{1},\beta_{2})=\delta(\beta_{1})\delta(\beta_{2})$.
In particular, in the following we will always understand the term

$\frac{N}{\pi(k_{i}-\ell_{i})}e^{-\frac{N-k_{i}+\ell_{i}}{k_{i}-\ell_{i}}(\beta_{1}^{2}+\beta_{2}^{2})}$
(52)

to be the delta function $\delta(\beta_{1})\delta(\beta_{2})$ whenever
$\ell_{i}=k_{i}$. In fact, recall that

$\lim_{\varepsilon\rightarrow
0}\frac{1}{\pi\varepsilon^{2}}e^{-(x^{2}+y^{2})/{\varepsilon^{2}}}=\delta(x)\delta(y)\,.$
(53)

Note that throughout this work we will use Cartesian coordinates of the
complex plane, meaning that $\beta_{1}$ and $\beta_{2}$ are real variables.
Now recall that the P-function of a tensor product is simply the product of
the P-functions, i.e.

$P_{\Pi_{\bm{k}}}^{(N)}(\bm{\beta})=\prod_{i=1}^{M}P_{\Pi_{k_{i}}}^{(N)}(\beta_{2i-1},\beta_{2i})\,.$
(54)

After some calculations, one obtains

$P_{\Pi_{k_{i}}^{(N)}}(\beta_{2i-1},\beta_{2i})=\binom{N}{k_{i}}\sum_{\ell_{i}=0}^{k_{i}}\binom{k_{i}}{\ell_{i}}(-1)^{\ell_{i}}\frac{N}{\pi(k_{i}-\ell_{i})}e^{-\frac{N-k_{i}+\ell_{i}}{k_{i}-\ell_{i}}(\beta_{2i-1}^{2}+\beta_{2i}^{2})}\,,$
(55)

i.e. a linear combination of Gaussians. Note that, in the previous expression,
for $\ell_{i}=k_{i}$ we obtain a delta function, which we will need to take
care of separately. Consequently, the probability reads

$p(\bm{k})=\int
d^{2M}\bm{\beta}\,\frac{e^{-\bm{\beta}^{\intercal}{\Sigma}^{-1}\bm{\beta}}}{\sqrt{\det{{\Sigma}}}}\prod_{i=1}^{M}\binom{N}{k_{i}}\sum_{\ell_{i}=0}^{k_{i}}\binom{k_{i}}{\ell_{i}}\frac{N(-1)^{\ell_{i}}}{\pi(k_{i}-\ell_{i})}e^{-\frac{N-k_{i}+\ell_{i}}{k_{i}-\ell_{i}}(\beta_{2i-1}^{2}+\beta_{2i}^{2})}\,.$
(56)

First, we want to invert the product and the sum that appear in the previous
expression. This is easily done by using the following property

$\prod_{i=1}^{M}\sum_{\ell_{i}=0}^{k_{i}}f(\ell_{i},k_{i})=\sum_{0\leq\ell_{1}\leq
k_{1}}\cdots\sum_{0\leq\ell_{M}\leq
k_{M}}\prod_{i=1}^{M}f(\ell_{i},k_{i})\equiv\sum_{0\leq\bm{\ell}\leq\bm{k}}\prod_{i=1}^{M}f(\ell_{i},k_{i})\,,$
(57)

where $f(\ell_{i},k_{i})$ is a generic function of $\ell_{i}$ and $k_{i}$. The
probability then reads

$\begin{split}p(\bm{k})&=\sum_{0\leq\bm{\ell}\leq\bm{k}}\int
d^{2M}\bm{\beta}\,\frac{e^{-\bm{\beta}^{\intercal}{\Sigma}^{-1}\bm{\beta}}}{\sqrt{\det{{\Sigma}}}}\prod_{i=1}^{M}\left[\binom{N}{k_{i}}\binom{k_{i}}{\ell_{i}}\frac{N(-1)^{\ell_{i}}}{\pi(k_{i}-\ell_{i})}e^{-\frac{N-k_{i}+\ell_{i}}{k_{i}-\ell_{i}}(\beta_{2i-1}^{2}+\beta_{2i}^{2})}\right]\\\
&=\sum_{0\leq\bm{\ell}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{k_{i}}\binom{k_{i}}{\ell_{i}}\frac{N(-1)^{\ell_{i}}}{\pi(k_{i}-\ell_{i})}\right]\int
d^{2M}\bm{\beta}\,\frac{e^{-\bm{\beta}^{\intercal}{\Sigma}^{-1}\bm{\beta}}}{\sqrt{\det{{\Sigma}}}}e^{{\sum_{i=1}^{M}-\frac{N-k_{i}+\ell_{i}}{k_{i}-\ell_{i}}(\beta_{2i-1}^{2}+\beta_{2i}^{2})}}\,.\end{split}$
(58)

To further ease the notation we introduce the multinomial coefficient

$\binom{N}{\ell_{i}}\binom{\ell_{i}}{k_{i}}=\frac{N!}{(N-\ell_{i})!(\ell_{i}-k_{i})!k_{i}!}\equiv\binom{N}{N-\ell_{i},\ell_{i}-k_{i},k_{i}}$
(59)

and make the variable change $d_{i}=k_{i}-\ell_{i}$.

$\begin{split}p(\bm{k})=\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}\frac{N(-1)^{k_{i}-d_{i}}}{\pi
d_{i}}\right]\int
d^{2M}\beta\,\frac{e^{-\bm{\beta}^{\intercal}{\Sigma}^{-1}\bm{\beta}}}{\sqrt{\det{{\Sigma}}}}e^{{\sum_{i=1}^{M}-\frac{N-d_{i}}{d_{i}}(\beta_{2i-1}^{2}+\beta_{2i}^{2})}}\end{split}$
(60)

Let us also define the set

$Z=\\{i|1\leq i\leq M,\,d_{i}=0\\}\,,$ (61)

which identifies the presence of delta functions in the integrand. Integrating
over $\delta(\beta_{2i-1})\delta(\beta_{2i})$ has the effect of setting
$\beta_{2i-1}=\beta_{2i}=0$ which, in turn, is equivalent to deleting the
corresponding rows and columns from the matrix $\Sigma^{-1}$. We will denote
this new matrix with $(\Sigma^{-1})_{(i)}$. In general, we will use the
notation $(\Sigma^{-1})_{(Z)}$ to denote the matrix obtained from
$\Sigma^{-1}$ by eliminating the corresponding rows/columns according to the
set $Z$. In particular, if $Z=\\{a,b\\}$ we will eliminate the rows(columns)
numbered $2a-1,2a,2b-1$ and $2b$. Note that, if we delete all the $2M$ rows
and columns from $\Sigma^{-1}$ we are left with 1 by definition. After
integrating over all the delta functions we are left with

$p(\bm{k})=\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}(-1)^{k_{i}-d_{i}}\right]\prod_{i\notin
Z}\left[\frac{N}{\pi d_{i}}\right]\int
d^{2(M-|Z|)}\bm{\beta}\frac{1}{\sqrt{\det{\Sigma}}}e^{-\bm{\beta}^{\intercal}[(\Sigma^{-1})_{(Z)}+D_{Z}]\bm{\beta}}\,,$
(62)

where $|Z|$ is the cardinality of $Z$ and $D_{Z}$ is a diagonal matrix defined
as

$D_{Z}=\bigoplus_{i\notin
Z}\left(\frac{N-d_{i}}{d_{i}}\right)\mathbb{I}_{2}\,.$ (63)

The remaining integrals in Eq. (62) are multi-dimensional Gaussian integrals
that we can evaluate straightforwardly

$\int
d^{2(M-|Z|)}\bm{\beta}\,e^{-\bm{\beta}^{\intercal}[(\Sigma^{-1})_{(Z)}+D_{(Z)}]\bm{\beta}}=\frac{\pi^{M-|Z|}}{\sqrt{\det{(\Sigma^{-1})_{(Z)}+D_{(Z)}}}}\,.$
(64)

Putting everything together we obtain

$p(\bm{k})=\frac{1}{\sqrt{\det{\Sigma}}}\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}(-1)^{k_{i}-d_{i}}\right]\prod_{i\notin
Z}\left[\frac{N}{d_{i}}\right]\frac{1}{\sqrt{\det{(\Sigma^{-1})_{(Z)}+D_{Z}}}}\,.$
(65)

To conclude the derivation we define the Kensingtonian, a new matrix function
whose action is specified once $N$ and $\bm{k}$ have been fixed. This allows
us to express the probability as

$p(\bm{k})=\frac{\text{Ken}[O]}{\sqrt{\det{\Sigma}}}\,,$ (66)

where $O=\mathbb{I}-\Sigma^{-1}$ and

$\text{Ken}[A]=\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}(-1)^{k_{i}-d_{i}}\right]\prod_{j\notin
Z}\left[\frac{N}{d_{j}}\right]\frac{1}{\sqrt{\det{(\mathbb{I}-A)_{(Z)}+D_{Z}}}}\,.$
(67)

##  Appendix B Relation between the Kensingtonian and the Torontonian

In this section we show that, starting from the expression for the click-
counting probability distribution for a GBS experiment Eq. (26), we can
retrieve the results of Ref. [15] by setting $N=1$. In this scenario, there
are only two possible measurement outcomes for each detector, namely
$k_{i}=0,1$. This in turn implies that the summation variable $d_{i}$ can
either be equal to $0$ or $1$ and consequently that the matrix $D_{Z}$
vanishes, as can readily be seen from its definition Eq. (28). Eq. (27) then
simplifies to

$\text{Ken}[A]=\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[(-1)^{k_{i}-d_{i}}\right]\frac{1}{\sqrt{\det{(\mathbb{I}-A)_{(Z)}}}}=\sum_{0\leq\bm{d}\leq\bm{k}}(-1)^{n-\sum_{i}d_{i}}\frac{1}{\sqrt{\det{(\mathbb{I}-A)_{(Z)}}}}$
(68)

where $n=\sum_{i}k_{i}$ is the total number of clicks. Let us also define two
new sets for future convenience

$\mathcal{K}=\\{i|1\leq i\leq M,\,k_{i}=0\\}\,,$ (69) $\mathcal{X}=\\{i|1\leq
i\leq M,\,d_{i}=0\,\text{and}\,k_{i}=1\\}\,.$ (70)

By noticing that $n-\sum_{i}d_{i}$ corresponds to the cardinality of
$\mathcal{X}$ and that $Z=\mathcal{X}\cup\mathcal{K}$ we can rewrite Eq. (68)
as

$\begin{split}\text{Ken}[A]=\sum_{0\leq\bm{d}\leq\bm{k}}(-1)^{|\mathcal{X}|}\frac{1}{\sqrt{\det{(\mathbb{I}-A)_{(\mathcal{K}\cup\mathcal{X})}}}}=\sum_{Y\in
P([n])}(-1)^{|\mathcal{Y}|}\frac{1}{\sqrt{\det{(\mathbb{I}-A_{(\mathcal{K})})_{(\mathcal{Y})}}}}\equiv\text{Tor}[A_{(\mathcal{K})}]\,,\end{split}$
(71)

where $P([n])$ is the power set of $[n]=\\{1,\dots,n\\}$ and Tor is the
Torontonian of a matrix. Note how, strictly speaking, the Torontonian is
defined slightly differently in Ref. [15]. The reason lies in the fact that
the authors used a different quadrature operators ordering and this choice is
reflected in different rules for deleting rows and column from the matrix $A$.
Consequently, we should more accurately say that we have found an equivalent
expression of the Torontonian valid for the operator ordering set by Eq. (1).
We have thus shown that for $N=1$ we correctly retrieve the probability
distribution obtained in Ref. [15]

$p(\bm{k})=\frac{\text{Ken}[O]}{\sqrt{\det{\Sigma}}}=\frac{\text{Tor}[O_{(\mathcal{K})}]}{\sqrt{\det{\Sigma}}}\,.$
(72)

Lastly, note that the $\Sigma$ matrix used in Ref. [15] differs from the one
used in this paper by a unitary transformation. This transformation, however,
does not affect the determinants contained in the probability formula, leaving
the latter unchanged.

##  Appendix C The loop-Kensingtonian

In this section we generalize the formula for the click-counting detection
probability distribution of an $M-$mode Gaussian state Eq. (26), when we lift
the zero-displacement constraint. This is achieved by introducing a new matrix
function that we name the loop Kensingtonian. The latter plays an analogous
role to the loop Torontonian Ref. [26] and the loop Hafnian Ref. [14]
functions in GBS setups employing threshold detection and PNR detection,
respectively. Apart from fundamental interest, this generalization is also
motivated by the fact that some of the applications of GBS, like computing
molecular vibronic spectra, require displacement in order to encode the
problem into a Gaussian boson sampler.
The derivation proceeds similarly to what we have seen in the zero-
displacement case, the only difference being that the $Q$ function of the
Gaussian state $\rho$ to be substituted in Eq. (49) now reads

$Q_{\rho}(\bm{\beta})=\frac{1}{\pi^{M}\sqrt{\det{{\Sigma}}}}e^{-(\bm{\beta}-\bm{\alpha})^{\intercal}{\Sigma}^{-1}(\bm{\beta}-\bm{\alpha})}\,,$
(73)

where $\bm{\alpha}$ is the displacement of $\rho$. Recall that our conventions
are such that a single mode coherent state $\ket{\gamma}$ has displacement
vector equal to $\bm{\alpha}=(\Re{\gamma},\Im{\gamma})^{\intercal}$. Following
the same steps of Appendix A we get to

$\begin{split}p(\bm{k})&=\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}\frac{N(-1)^{k_{i}-d_{i}}}{\pi
d_{i}}\right]\int
d^{2M}\bm{\beta}\,\frac{e^{-(\bm{\beta}-\bm{\alpha})^{\intercal}{\Sigma}^{-1}(\bm{\beta}-\bm{\alpha})}}{\sqrt{\det{{\Sigma}}}}e^{{\sum_{i=1}^{M}-\frac{N-d_{i}}{d_{i}}(\beta_{2i-1}^{2}+\beta_{2i}^{2})}}\\\
&=\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}\frac{N(-1)^{k_{i}-d_{i}}}{\pi
d_{i}}\right]\frac{e^{-\bm{\alpha}^{\intercal}{\Sigma}^{-1}\bm{\alpha}}}{\sqrt{\det{\Sigma}}}\int
d^{2M}\bm{\beta}\,e^{-\bm{\beta}^{\intercal}{\Sigma}^{-1}\bm{\beta}+2\bm{\alpha}^{\intercal}\Sigma^{-1}\bm{\beta}}e^{{\sum_{i=1}^{M}-\frac{N-d_{i}}{d_{i}}(\beta_{2i-1}^{2}+\beta_{2i}^{2})}}\,,\end{split}$
(74)

where we have used the fact that $\Sigma^{-1}$ is a symmetric matrix.
Integrating over the delta function contributions (i.e. $d_{i}=0$) we obtain

$p(\bm{k})=\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}(-1)^{k_{i}-d_{i}}\right]\prod_{i\notin
Z}\left[\frac{N}{\pi
d_{i}}\right]\frac{{e^{-\bm{\alpha}^{\intercal}{\Sigma}^{-1}\bm{\alpha}}}}{\sqrt{\det{\Sigma}}}\int
d^{2(M-|Z|)}\bm{\beta}\,e^{-\bm{\beta}^{\intercal}[(\Sigma^{-1})_{(Z)}+D_{Z}]\bm{\beta}+\bm{\beta}^{\intercal}(2\Sigma^{-1}\bm{\alpha})_{(Z)}}\,,$
(75)

where $Z$ and $D_{Z}$ are defined by Eq. (61) and Eq. (28), respectively. Note
that the subscript notation $(2\Sigma^{-1}\bm{\alpha})_{(Z)}$ applied to a
vector denotes the elimination of its elements according to the set $Z$. We
can now easily evaluate the remaining Gaussian integrals in the previous
expression

$\int
d^{2(M-|Z|)}\bm{\beta}\,e^{-\bm{\beta}^{\intercal}[(\Sigma^{-1})_{(Z)}+D_{Z}]\bm{\beta}+\bm{\beta}^{\intercal}(2\Sigma^{-1}\bm{\alpha})_{(Z)}}=\frac{\pi^{M-|Z|}}{\sqrt{\det{(\Sigma^{-1})_{(Z)}+D_{Z}}}}e^{(\Sigma^{-1}\bm{\alpha})^{\intercal}_{(Z)}[(\Sigma^{-1})_{(Z)}+D_{Z}]^{-1}(\Sigma^{-1}\bm{\alpha})_{(Z)}}\,.$
(76)

Putting everything together we obtain

$p(\bm{k})=\frac{e^{-\bm{\alpha}^{\intercal}\Sigma^{-1}\bm{\alpha}}}{\sqrt{\det{\Sigma}}}\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}(-1)^{k_{i}-d_{i}}\right]\prod_{j\notin
Z}\left[\frac{N}{d_{j}}\right]\frac{e^{(\Sigma^{-1}\bm{\alpha})^{\intercal}_{(Z)}[(\Sigma^{-1})_{(Z)}+D_{Z}]^{-1}(\Sigma^{-1}\bm{\alpha})_{(Z)}}}{\sqrt{\det{(\Sigma^{-1})_{(Z)}+D_{Z}}}}\,.$
(77)

Finally, we can rewrite the previous expression as

$p(\bm{k})=p(\bm{0})\text{lken}[O,\bm{\alpha}]\,,$ (78)

where
$p(\bm{0})=e^{-\bm{\alpha}^{\intercal}\Sigma^{-1}\bm{\alpha}}/\sqrt{\det{\Sigma}}$,
$O=\mathbb{I}-\Sigma^{-1}$ and

$\text{lken}[A,\bm{\alpha}]=\sum_{0\leq\bm{d}\leq\bm{k}}\prod_{i=1}^{M}\left[\binom{N}{N-k_{i},k_{i}-d_{i},d_{i}}(-1)^{k_{i}-d_{i}}\right]\prod_{j\notin
Z}\left[\frac{N}{d_{j}}\right]\frac{{e^{((\mathbb{I}-A)\bm{\alpha})^{\intercal}_{(Z)}[(\mathbb{I}-A)_{(Z)}+D_{Z}]^{-1}((\mathbb{I}-A)\bm{\alpha})_{(Z)}}}}{\sqrt{\det{(\mathbb{I}-A)_{(Z)}+D_{Z}}}}\,.$
(79)

is the loop Kensingtonian of a matrix $2M\times 2M$ matrix $A$.

##  Appendix D Click-counting detection POVM in the $N\rightarrow\infty$ limit

In this section we prove that in the $N\rightarrow\infty$ limit, ideal click-
counting detection converges PNR detection. To do so, we simply Taylor-expand
the POVM element Eq. (13)

$\begin{split}\Pi_{k}^{(N)}&=\,:\binom{N}{k}e^{-\hat{n}}(e^{\frac{\hat{n}}{N}}-1)^{k}\,:=\,:\frac{N!}{k!(N-k)!}\left(\frac{\hat{n}}{N}\right)^{k}e^{-\hat{n}}:+\,O\left(\frac{1}{N}\right)\\\
&=\frac{N(N-1)\cdots(N-k+1)}{N^{k}}:\frac{\hat{n}^{k}}{k!}e^{-\hat{n}}:+\,O\left(\frac{1}{N}\right)\stackrel{{\scriptstyle\small{N\rightarrow\infty}}}{{\rightarrow}}\,:\frac{\hat{n}^{k}}{k!}e^{-\hat{n}}:\,\equiv\outerproduct{k}{k}\,.\end{split}$
(80)

##  Appendix E Complexity of click-GBS (alternative proof)

In what follows we present an alternative proof of hardness for click-counting
GBS that does not rely on using the output probability distribution of a GBS
task employing threshold detectors. For simplicity of exposition, we consider
a single mode scenario and let $\rho$ be the output state of the LON. We show
that, in the non-collisional regime, the output probability distribution $p$
of click-counting GBS is arbitrarily close (in total variational distance) to
that of a Gaussian boson sampler employing PNR detectors, which we denote with
$\tilde{p}$. The latter reads $\tilde{p}(k)=\Tr\\{\rho\tilde{\Pi}_{k}\\}$ for
$k\in\\{0,1,\dots\\}$, and the POVM element is a projector on a Fock state,
i.e. $\tilde{\Pi}_{k}=\outerproduct{k}{k}$. For click-counting detection we
have $p(k)=\Tr\\{\rho\Pi^{(N)}_{k}\\}$ with $k\in\\{0,\dots,N\\}$, and the
POVM element $\Pi_{k}^{(N)}$ is given by Eq. (13). The latter may be
alternatively expressed as [23]

$\Pi^{(N)}_{k}=\binom{N}{k}\sum_{n=k}^{\infty}\frac{1}{N^{n}}\partial_{x}^{n}[e^{x}-1]^{k}|_{x=0}\outerproduct{n}{n}\equiv\sum_{n=k}^{\infty}c_{k}(n)\outerproduct{n}{n}=c_{k}(k)\outerproduct{k}{k}+\sum_{n=k+1}^{\infty}c_{k}(n)\outerproduct{n}{n}\,,$
(81)

where

$c_{k}(n)=\binom{N}{k}\frac{1}{N^{n}}\partial_{x}^{n}[e^{x}-1]^{k}|_{x=0}\,.$
(82)

Notice how for $k>N$ we can set $\Pi_{k}^{(N)}=0$ without loss of generality.
The total variational distance between $p$ and $\tilde{p}$ reads

$||\tilde{p}-p||_{1}=\frac{1}{2}\sum_{k=0}^{\infty}|\tilde{p}(k)-p(k)|=\frac{1}{2}\sum_{k=0}^{\infty}|\Tr\\{\rho(\tilde{\Pi}_{k}-\Pi_{k}^{(N)})\\}|\,.$
(83)

We can then write $\tilde{\Pi}_{k}-\Pi_{k}^{(N)}=A_{k}-B_{k}$, where

$A_{k}=\begin{cases}(1-c_{k}(k))\outerproduct{k}{k}&0\leq k\leq N\\\
\outerproduct{k}{k}&k>N\end{cases}\,,\quad\quad
B_{k}=\begin{cases}\sum_{n=k+1}^{\infty}c_{k}(n)\outerproduct{n}{n}&0\leq
k\leq N\\\ 0&k>N\end{cases}\,,$ (84)

and

$c_{k}(k)=\frac{N!}{(N-k)!N^{k}}\,.$ (85)

Since the POVM elements sum to the identity
$\sum_{k=0}^{\infty}\tilde{\Pi}_{k}=\sum_{k=0}^{\infty}{\Pi}^{(N)}_{k}=\mathcal{I}$,
we have that

$\sum_{k=0}^{\infty}(\tilde{\Pi}_{k}-\Pi_{k}^{(N)})=\sum_{k=0}^{\infty}(A_{k}-B_{k})=0\,,$
(86)

which in turn implies that
$\sum_{k=0}^{\infty}A_{k}=\sum_{k=0}^{\infty}B_{k}$. Using this and the fact
that both $A_{k}$ and $B_{k}$ are positive semi-definite operators, we can
write

$\begin{split}||\tilde{p}-p||_{1}&=\frac{1}{2}\sum_{k=0}^{\infty}|\Tr{\rho
A_{k}}-\Tr{\rho B_{k}}|\leq\frac{1}{2}\sum_{k=0}^{\infty}\left(\Tr{\rho
A_{k}}+\Tr{\rho B_{k}}\right)=\sum_{k=0}^{\infty}\Tr{\rho A_{k}}\\\
&=\sum_{k=0}^{N}(1-c_{k}(k))\bra{k}\rho\ket{k}+\sum_{k=N+1}^{\infty}\bra{k}\rho\ket{k}=\sum_{k=0}^{N}(1-c_{k}(k))\tilde{p}(k)+\sum_{k=N+1}^{\infty}\tilde{p}(k)\,.\end{split}$
(87)

Notice how the $k=0$ and $k=1$ terms are actually null, since Eq. (85) implies
that $c_{0}(0)=c_{1}(1)=1$. In the non-collisional regime, the probability
$\varepsilon$ of observing two or more photons

$\varepsilon=\sum_{k=2}^{\infty}\tilde{p}(k)$ (88)

becomes negligible. Hence, we can finally write

$||\tilde{p}-p||_{1}\leq\varepsilon\,,$ (89)

ans use the same arguments presented in the main text to claim the
computational complexity of click-counting GBS. The generalization to the
$m-$mode, although cumbersome, may be easily obtained following the steps
outlined above.