Gaussian boson sampling with click-counting detectors

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 graphrelated 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 negligi- ble, 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 I 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.

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 † k ] = δ jk .We can then introduce quadrature operators for each mode, defined as q j = (a j + a † j )/ √ 2 and p j = (a j − a † j )/i √ 2 (where we have set ℏ = 1), and arrange them into the following vector r = (q 1 , p 1 , . . ., q M , p M ) ⊺ . ( A Gaussian state ρ is completely characterized by its vector of first moments (displacement) r and its covariance matrix σ, defined respectively as One can show that the s−ordered PQD for a generic M -mode Gaussian state with covariance matrix σ and vector of first moments α is given by [5] Here s = (s 1 , . . ., s M ) ⊺ is the vector of operator orderings, where s j ∈ R, and the matrix s is defined as The conventions used in Eq. ( 4) are such that for a singlemode coherent state |γ⟩ the covariance matrix is the identity matrix σ = I 2 and the vector of first moments reads α = (Re{γ}, Im{γ}).Note that Eq. ( 4) is well-defined iff σ− s ≥ 0, otherwise the 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 s = −I M , s = 0 and s = 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 {Π x } satisfy the conditions Π x ≥ 0 and x Π x = I , where I is the identity operator on the Hilbert space.The probability p(x) of obtaining a given outcome x when measuring a state ρ is given by the Born rule, i.e. p(x) = Tr{ρΠ x }.
The POVM elements of PNR detection are the projectors on Fock states, namely |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 where k = (k 1 , . . ., k M ) ⊺ are the detected photon numbers according to the detection outcome (see Ref. [14] for more details) and The matrix function appearing in Eq. ( 6) is called the Hafnian and is defined as follows where A is a 2n × 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 ∈ {0, 1} and the POVM elements are Π 0 = |0⟩⟨0| and Π 1 = I − Π 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 where is the Torontonian of a 2n × 2n matrix A, P ([n]) is the power set of [n] = {1, . . ., 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 k ′ is compatible with a threshold detection pattern This observation leads to the following identity [15] Haf where O is a 2n × 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.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 clickcounting detector.We remind the reader that the POVM elements of on/off detection are where I is the identity operator, n is the number operator and : • : denotes normal ordering of bosonic operators [27].One can show that the POVM elements corresponding to a clickdetector made up of N on/off detectors are given by 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. Π 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 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 n → ηn + ν, where 0 ≤ η ≤ 1 and ν ≥ 0 are the efficiency and the dark count rate of the threshold detector, respectively, to obtain the following POVM elements One can then prove that the POVM elements of noisy clickcounting detection are given by We note that this operators can be obtained from Eq. ( 13) upon substituting n → ηn + N ν.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 η and dark count rate N ν, whose related POVM elements read In Appendix D we show that − as one might intuitively expect − in the N → ∞ limit and in the absence of noise we retrieve true PNR detection, i.e.
The noisy case requires some additional attention, as the presence of non-zero dark count rate ν 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 ν ≃ 10 −4 , hence in regimes where simultaneously N ≫ 1 and N ν ≪ 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 clickcounting 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 → ∞ 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 k = (k 1 , . . ., k M ), where 0 ≤ k i ≤ N ∀i.Note that the total number of clicks n = 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 where ρ is a generic M -mode Gaussian state with covariance matrix σ and null vector of first moments and is the POVM that characterize the M -mode click detection.The latter simply reads where ki 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 Note that the operator : e − N −k i +ℓ i N n : corresponds to the vacuum element Π (1) 0 of the noisy on/off detection POVM Eq. ( 15), with an effective detection inefficiency given by λ i ≡ (N − k i + ℓ i )/N ∈ [0, 1].This observation is crucial, as it implies that all is needed in order to obtain p(k) are (noisy) vacuum statistics of marginal states of ρ.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 where Q ρ (β) is the Husimi Q-function of ρ and P Π (N ) k (β) is the P-function of the POVM element.Using Eq. ( 4) with ordering parameter s j = −1 ∀j we obtain where we have introduced the matrix Σ = (σ + I)/2 for convenience.On the other hand, the P-function of Π (N ) k is readily obtained once we know how to compute the P-function of : e −λi n : .One can prove that for λ i ̸ = 1 while for λ i = 1 (i.e. for ℓ i = k i ) we have where β 1 and β 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 k reads Here O = I − Σ −1 and is the Kensingtonian of a 2M × 2M matrix A. D Z is a diagonal matrix defined as and the set Z is defined as follows In Eq. ( 27) we have used the subscript notation (I − A) (Z) to denote the matrix obtained from (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 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 → ∞ 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 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 n < 1.Furthermore, tracing out all output modes but one leaves us with a Gaussian marginal state that approximately looks like a thermal state ν th (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 σ = (2n + 1)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.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 n of the probe thermal state ν th (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 ε 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 where ⟨•⟩ U denotes averaging over M × 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 p.We define the set of collision events of an M −mode sampling problem as (32) By definition we have that p(k ∈ 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 clickcounting GBS reads The total variational distance between p and p reads By close inspection of the click-counting detection POVM Eq. ( 13) we notice that Π (N ) 0 =: e −n := |0⟩⟨0| for every value of N .Using this and the fact that the elements of a POVM resolve the identity we can write 1 . ( 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 where k / ∈ C is a collision-less detection event and C k is the set of all possible collision events compatible with k If we now substitute Eq. (36) into Eq.( 34) we obtain Let us now consider a probability distribution π that approximates p arbitrarily well, i.e. it satisfies ||p − π|| 1 = ε ′ .Using the triangle inequality of the L 1 norm we can bound ||p − π|| 1 as follows (39) If we now assume that there exists a polynomial time algorithm that can sample from π, then Eq. (39) tells us that the same algorithm can sample efficiently from an arbitrarily good approximation of 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 = 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 where 0 ≤ k i ≤ N .For N = 1 (threshold detectors) we have 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 (k) at fixed n can be expressed as a sum of a collision-less term and collision term, namely Here C n is the set of collision detection events at a fixed total number of clicks and we use the symbol ⟨•⟩ n to denote averaging at constant n.When there are no collisions we simply have that F (k) = 2 n .We can thus write where ε is the probability of observing a collision at fixed n, namely It can easily be seen that F (k) ≥ (n + 1), hence 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 (k) using the inequality of the arithmetic and geometric means and show that the following chain of inequalities holds 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 noncollisional 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 ≫ 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 clickcounting 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 × 2M matrix.On the other hand, the authors of Ref. [28] model their experimental setup as larger instance of a M N −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 2M N × 2M N 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.

FIG. 1 .
FIG.1.Schematics of a click-counting detector made up of N threshold detectors.The quantum state ρ to be measured enters an N −mode interferometer (represented by a unitary operation 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.
FIG.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 n of the probe thermal state ν th (n).In the energy regime of interest we considered there is a good agreement, with TVD smaller than 0.05.

where c k (n) = N k 1 N=
n ∂ n x [e x − 1] k | x=0 .(E2) Notice how for k > N we can set Π (N ) k = 0 without loss of generality.The total variational distance between p and p reads ||p − p|| 1 A k − B k , where

TABLE I .
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.