{'fraction_non_alphanumeric': 0.03898324667822068, 'fraction_numerical': 0.011808203350664357, 'mean_word_length': 4.3348126232741615, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Threshold detection statistics of bosonic states

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


I. INTRODUCTION
Quantum photonic experiments can generally be described as preparing quantum states of light, evolving them through linear optical interferometers, and detecting the output photons. While the most common types of photonic states, Fock states and Gaussian states, can be routinely prepared via spontaneous processes in optical non-linearities or (artificial) atomic systems and processed with high-fidelity linear optical components [1], photon number detection is typically approximated via the use of threshold photon detectors. Threshold detection, i.e. a measurement distinguishing only between vacuum and the presence of one or more photons, is widely available, e.g. via high-efficiency superconducting nanowires [2] or room-temperature avalanche photodiodes [3], making it the standard measurement apparatus in quantum photonics. Its use in experiments, represented in Fig. 1, encompasses many areas of quantum research, including demonstrations of quantum advantages, e.g. in computation [4], measurement sensitivity [5], and loophole-free tests of non-locality [6]. However, threshold detection only provides a meaningful approximation of the desired Fock basis measurement projectors in the regime of low mean photon numbers per mode. On the other hand, as technology progresses, mean photon numbers increase and higher fidelities are demanded [4], making this approximation less appropriate.
To circumvent this issue, experiments can be described directly using the output statistics of threshold detection instead of its photon number resolving approximation. However, despite the wide adoption of such systems, there are in general no closed form expressions in the literature for computing measurement probabilities of threshold detectors. In fact, while an expression for the threshold detection of zero-displaced Gaussian states is known, given by the Torontonian matrix function [7], no analogous expressions exist for other commonly used states, e.g. Fock or displaced Gaussian states. For example, for Fock states with fixed photon number, thresh- formed by squeezing, displacement and linear optics are captured by the loop Torontonian matrix function. Displacement can be generated using coherent states from lasers, squeezing can come from nonlinear processes such as Spontaneous Parametric Down Conversion (SPDC) or Spontaneous Four-Wave Mixing (SFWM). Linear optics can be implemented in a variety of platforms, including bulk and integrated optics. States created by the linear optical interference of Fock states (highlighted in blue), as generated by (artificial) atoms, lead to threshold detection probabilities given by the Bristolian matrix function. old probabilities could be exactly calculated by summing over all possible output states which lead to the given threshold detector outcome. However, this method requires calculating a number of probabilities scaling combinatorially with the number of clicked detectors, rendering it impractical already for smaller-scale experiments [8][9][10]. New methods are required to describe quantum photonic technologies which use threshold detection.
Here, we provide such methods by developing a unified picture to compute threshold statistics for most quan- tum photonic states of experimental interest. As described in Table I, this is achieved by introducing two new matrix functions, the Bristolian and the loop Torontonian, for threshold statistics with Fock and displaced Gaussian states, respectively, and demonstrating close connections between them and to other existing matrix functions. The developed tools provide exact simulation, design, and analysis methods for current [4,[8][9][10][11][12] and future quantum photonic systems that use threshold detection. We wish to highlight the different challenges between computing probabilities, known as strong simulation, which we focus on in this work, and drawing samples from a probability distribution, known as weak simulation [13]. These tasks often have very different complexity. For example, using methods from Ref. [14], we can sample threshold detector outcomes without ever calculating a threshold detection probability.

II. THRESHOLD DETECTION STATISTICS FROM VACUUM STATISTICS
Threshold detectors are described by the measurement operatorsΠ for vacuum (0) and click (1) outcomes on a mode described by label j. We use |0 to denote the vacuum state of an optical mode, |k = (â † ) k |0 / √ k! for Fock states of the optical mode, and I is the identity operator (we will always assume its dimension to be the same as the other operators appearing in the equation).
We write the outcome of M threshold detectors, labelled with j ∈ [M ] = {1, 2, . . . , M }, using a length-M bit-string d, where the jth element gives the measurement outcome of the jth mode. Defining a set of modes which clicked, C = {j ∈ [M ] | d j = 1}, and a set for modes with the vacuum outcome, V = {j ∈ [M ] | d j = 0}, we can write the multimode measurement operator asΠ which can be rearranged to givê Here, we use P (C) to denote the powerset of C and |Z| for the number of elements in a set Z. | 0 V 0 V | describes the vacuum projector in all the vacuum outcome modes and | 0 Z 0 Z | describes the vacuum projector in all modes in a subset Z ⊆ C. Eq. (3) indicates that to calculate the threshold detection probabilities for any state it is sufficient to calculate marginal vacuum probabilities, which are used in an inclusion/exclusion sum as described by Eq. (3). Using this measurement operator and the Born rule p( d) = tr(Π ( d) ρ) on some state ρ, we find: where This formula provides our starting point for deriving general expressions for threshold detection statistics.

III. MARGINAL VACUUM PROBABILITIES FROM THE PHOTON NUMBER PROBABILITY GENERATING FUNCTION
For an M -mode linear optical interferometer, described by an M × M matrix U and the operatorÛ, the creation operators are transformed aŝ Considering an input state |Φ 0 , the output photon number probability distribution is then where m is a length-M list describing the photon number in each mode at the output of the interferometer, and [15], considering the Fourier transform of the probability distribution of photon number basis measurements we define the characteristic function which, with some manipulation (see Appendix B), can be expressed as whereÛ φ is the operator given by the evolution due to the linear optical transformation We can transform this into a probability generating function, G, using the substitution x j = exp(iφ j ): The function G( x) has the following useful properties. To marginalise the jth mode, we simply set x j = 1. If we set x j = 0, this gives us the probability for n j = 0. Therefore, if we want to calculate the probability that some subset of the modes, V , measure vacuum and we marginalise over all other modes, B, we can evaluate where 0 ( 1) is a vector with 0 (1) in all entries. G( x) is the probability distribution generating function, so by taking derivatives of G( x), we can find information about the photon number basis probability distribution [15]. By using the expression for the characteristic function in Eq. (8), we can see that this amounts to calculating the scattering amplitude of |Φ 0 to itself, through a linear optical interferometer described by the transformation U † U φ U . From again using the substitution x j = exp(iφ j ) in Eq. (9), we see that U φ physically corresponds to either zero transmission for modes in V , or unit transmission for modes in B. As we show in the next sections, we can use this, in conjunction with Eq. (4) to calculate threshold detection probabilities for all the experimental scenarios outlined in Fig. 1.

IV. FOCK STATE INPUTS
Recall that the scattering amplitudes of Fock states evolved through a lossless interferometer are given by the permanent matrix function [16] where U m, n is constructed from U by repeating its jth row m j times and its jth column n j times for all j ∈ [M ]. Therefore if we have an N -photon input Fock state, |Φ 0 = | n , we can use Eq. (8), Eq. (10) and Eq. (12) to write where U x is formed like U φ , but with diagonal matrix elements: [U x ] jj = x j . For this equation to be valid, we must have a lossless unitary transformation. However, we are free to marginalise over modes by allowing elements of x to be set to 1 for any mode we wish to marginalise over, including any loss modes, as shown in Eq. (11). Because Eq. (13) provides us with a closed form expression for marginal vacuum probabilities, and because Eq. (4) shows us that marginal vacuum probabilities are sufficient to calculate threshold detection probabilities, we can use this to derive a matrix function for calculating threshold detection probabilities of Fock states. For more generality, we first consider a linear optical transformation with losses, described by an M in ×M out matrix T , with singular values upper bounded by 1 [17]. In Appendix C we show that if the input state is an M in -mode Fock state, n, then by combining Eq. (4) with Eq. (13), we can calculate the threshold detection probability of the outcome described by an M out -length bit-string d as: Here we have introduced a matrix function, the Bristolian, defined as where A is an m×n matrix and E is an n×n matrix. A Y denotes selecting the rows of A according to the elements of Y , and [m] = {1, 2, . . . , m}. We have also defined a matrix which accounts for the mixing with vacuum in the environment modes Our naming of the Bristolian is inspired by the convention established by the Hafnian and Torontonian matrix functions, which are named after the cities of their discovery. By noticing that I−T † T gives a zero matrix when T is unitary, the brs function can be simplified when T is a unitary matrix U , only requiring the rows of U which correspond to modes with a detector click, providing Here we defined the Unitary Bristolian acting on an m × m matrix, A, as Gaussian states are the set of states that have a Gaussian characteristic function. A Gaussian state ρ is uniquely characterized by its vector of means with entries and its Husimi covariance matrix with entries where we have used a vector of creation and annihilation operatorsˆ The Husimi function Q( r) = r| ρ | r maps displacement vectors, r, to probabilities, so to calculate vacuum probabilities we can evaluate the Husimi function at the origin. Noting that we can marginalise over modes by deleting all the corresponding elements of Σ and α, we obtain [18] .
The notation Σ V V and α V differs slightly here from the previous section, as now there are two basis vectors for each mode of our system, corresponding to each mode'ŝ a andâ † operator. We form Σ V V by selecting both rows/columns of Σ which correspond to each element of V and we form α V by selecting both elements of α corresponding to each element of V . We can use this to immediately arrive at a threshold detection probability for displaced Gaussian states using Eq. (4). However, here we must invert and compute determinants for square matrices of size 2(|V | + |Z|). It would be preferable if we could reduce these to matrices of size 2|Z|. It would also be helpful conceptually to have a formula which can be connected to other relevant matrix functions, the Torontonian [7] and the loop Hafnian [19]. Therefore, it is of interest to write this probability in terms of In Appendix F, we show how Eq. (4) can be rearranged into the following: where C is given by the index of the elements of d where d j = 1, so O CC and γ C are the matrix and vector formed by selecting the rows/columns of O and elements of γ which correspond to modes which see a detector click.
p( 0) is the probability of detecting vacuum in all modes, and can be calculated using Eq. (23). We introduce the loop Torontonian, which is defined as where O is a 2m×2m matrix and γ is a 2m-length vector.

VI. CONNECTIONS BETWEEN MATRIX FUNCTIONS
In the limit of no displacement α = γ = 0, the exponential terms in the numerator of Eq. (26) becomes 1 and thus ltor O, 0 = tor (O), where tor is the Torontonian function from Ref. [7]. One can show, using the scattershot construction [20], that the Torontonian and Bristolian are related via the following limit where n is a bitstring (implying that this identity is only valid for single-photon or vacuum inputs), N = i n i and C is the union of the labels of the modes in which single photons were input into the interferometer and the labels of the modes in which clicks are registered. This relation is proven in Appendix H.
As we show in Appendix G, the loop Torontonian can also be used as a generating function for the loop Hafnian, where X = 0 I I 0 and = |C|. We use this to derive the trace formula for the loop Hafnian, the fastest known method for computing photon number resolved measurement probabilities on displaced Gaussian states [19,21,22]. Because the loop Hafnian of a bipartite graph is given by the matrix permanent [21], all the matrix functions in Table I can be derived from the loop Torontonian. We also see a connection between the Bristolian and the permanent when an N photon Fock state results in N threshold detector clicks. In this case, each threshold detector must have seen exactly 1 photon, so we can describe the measurement operator of each threshold detector click as a single photon projector, which leads to describing the event with permanents, as given by Eq. (12).
In Appendix D, we show this link directly by first describing the Unitary Bristolian for N photon, N click events as the permanent of an N × N × N 3-tensor [23].

VII. TIME COMPLEXITIES
In Appendix I, we discuss the time complexities for the Bristolian and the loop Torontonian. We find that, using the formulae presented in this work, the Bristolian, brs(A, E), has a time complexity of O(n2 n+m ) for an m × n matrix A and n × n matrix E and the loop Torontonian, ltor(O, γ), has time complexity of O(m 3 2 m ) for a 2m × 2m matrix O and 2m-length vector γ. For the loop Torontonian, this complexity can be reduced using a recursive strategy which exploits Cholesky decomposition [24]. We also believe that the Bristolian's time complexity can likely be reduced, and we leave this as an open problem.

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

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

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

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

Lossless Hong-Ou-Mandel
Our first example is Hong-Ou-Mandel interference [35] of single photons on a lossless 50/50 beam splitter. Here, the input state is n = (1, 1) and The probability of detecting a coincidence is given by So, C = {1, 2}. We expand the ubrs function as with U 1 = (1, 1)/ √ 2 and U 2 = (1, −1)/ √ 2 being the 1st and 2nd rows of U respectively. We also use that per U † U = per(I) = 1, and calculate giving p( d = (1, 1)) = 0, as expected. Here, we have not included the term where Y is the empty set. In this case, we are considering the permanent of the all zeros matrix, which is zero and so does not contribute. We can also see that p( d = (1, 0)) = per [U 1 ] † U 1 = 1/2, and similarly p( d = (0, 1)) = 1/2. These are due to the (2, 0) and (0, 2) photon number output terms.

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

Lossy Hong-Ou-Mandel
The Hong-Ou-Mandel effect is preserved under balanced loss. In contrast to the previous examples, in this example example, the loss means that we need to use the Bristolian, as the Unitary Bristolian is no longer valid. We consider a transmission matrix defined like U in Eq. (A1), but with transmission η, giving T = √ ηU .
In this case where E = (1 − η)I. We evaluate the Bristolian to find which confirms that the coincidence event is still suppressed under balanced loss.

Lossy zero transmission law
When there is loss, we start to witness threshold detector outcomes which were suppressed in the lossless case. We repeat the example in section A 2 but adding a transmission η ≤ 1. So, T = √ ηU for U defined in Eq. (A7).
Using the Bristolian, we find p( d = (1, 1, 0)) = brs(T d, n , E(T ) n, n ), where E(T ) = (1 − η)I, and this gives brs(T d, n , E(T ) n, n ) = per [T (1,2) We can check this problem intuitively by considering all the losses to be applied just before the measurement. In this picture, the outcome d = (1, 1, 0) can only occur when we have the state m = (1, 1, 1) before the losses, then the first two photons are transmitted and the last photon is lost. Therefore we expect the probability found in Eq. (A16). The factor 1/3 comes from section A 2, where p( m = (1, 1, 1)) = 1/3 before any loss is applied. For our final Bristolian example calculation, we consider a case where we have to repeat columns more than once. Consider the same T as above, but with the input state n = (1, 2, 0). We calculate the probability of d = (0, 1, 1) (so C = {2, 3}), denoting the elements of T as t jk . where, and finally, Here: T 2, n = (t 22 , t 23 , t 23 ) and T 3, n = (t 32 , t 33 , t 33 ). We will not symbolically evaluate this expression, but in the interest of providing simple test cases for future software implementations, we see numerically that p( d) = 0.222 . . . for η = 1, p( d) = 0.189 for η = 0.9 and p( d) = 0.069444 . . . for η = 0.5.

Using the loop Torontonian
When performing calculations using the loop Torontonian, we need to know the matrix O and the vector γ for our state. Tools such as Strawberry Fields [36] allow for conveniently computing the real means vector, µ, and covariance matrix, σ, of a Gaussian state's Wigner function. To use σ and µ to find O and γ, we can use functionality from The Walrus [34]. We can convert from σ to Σ by using the Qmat function (The Walrus uses a different ordering for ζ and so we also apply a complex conjugate to match our definition of Σ), and similarly convert from µ to α using complex to real displacements. These can be used to find O = I − Σ −1 and γ = (Σ −1 α) * .
To show how to select the appropriate rows/columns for calculating threshold detector outcome probabilities, we will consider a 5-mode experiment, with the outcome d = (1, 0, 1, 0, 0 This assumed the basis vector ordering convention specified by Eq. 21, where mode j corresponds to basis vectors j and j + M , as is used in The Walrus [34].

Appendix B: Derivation of the characteristic function
We start from the photon number probability distribution in Eq. 6 and the definitions of the characteristic function in Eq. 7 and rearrange terms Notice that we can define an operator,Û φ , which acts likeÛ as the operator given by the linear optical transformation If we include this in our expression above, we find where we have used the resolution of the identity to arrive at the answer.

Appendix C: Derivation of the Bristolian
The Unitary Bristolian has a slightly simpler derivation than the Bristolian, so we will begin by considering this case. We start by combining Eq. 13 with Eq. 4: , and x j = 1 otherwise. Since U x(V,Z) contains only zeros on the rows and columns given by the elements of V and Z, this is equivalent to deleting the rows and columns of U and U † respectively according to the elements of V and Z. Therefore, we can write the sum above as Here, we sum over the modes which are marginalised, instead of summing over the modes being projected into the vacuum state. This corresponds to using the substitution Z for Y = C \ Z, and noticing that the sum over Y = C \ Z : Z ∈ P (C) is the same as the sum over Z ∈ P (C). We use U Y to denote selecting only the rows of U according to the elements in the set Y . Note that since we delete rows when all their elements are set to all zeros, the permanent when Y is the empty set should be zero, as it corresponds to the permanent of an all zeros matrix, rather an empty matrix. Because Eq. (C2) contains an inclusion/exclusion formula, like that of Ryser's permanent formula [37], it could be viewed as a 3-dimensional permanent, similar to those which appear elsewhere in quantum photonics [23,38]. However, we note that to compute the probability for input state | n , interferometer transformation U and threshold detector pattern d, it is sufficient to know only the rows of U which correspond to nonzero elements of d and columns of U given by n. Therefore the input to this problem is a matrix, so we chose to write this probability in terms of a new matrix function, the Unitary Bristolian We construct U d, n from U as described under Eq. 12. The Unitary Bristolian, ubrs, is a matrix function which acts on some m × n matrix, A: which is the form this is reported in Eq. 18. Now, we are ready to derive the more general formula for the Bristolian. Unlike the permanent, the Bristolian can be generalised for calculating marginal detection probabilities. Consider that we wish to calculate the probability of observing detector clicks for modes in C, vacuum in modes V , and marginalise over modes in B, so the union of C, V and B is [M ]. We can write this probability as Here U Y ∪B is constructed by selecting rows of U according to the elements of Y and B. We could also write the summation in Eq. (C5) as where U Y, n U B, n is the augmented (|Y | + |B|) × N matrix formed by stacking the matrices U Y, n and U B, n , where N = j n j . The notation U C, n is used to show that we take rows of U according to the set C and repeat the columns of U according to n, and equivalently for U B, n . A counter-intuitive feature of this formula is that it depends on matrix elements of the linear transformation that are ignored by our measurements. In an experiment, changing the elements of U B should have no impact on p( d C = 1, d V = 0). Following this argument, we propose that we only need to know U C , and we can construct U B by performing a unitary dilation. This is particularly helpful for lossy experiments, where we typically do not have an understanding of the full unitary transformation acting on the both the experiment's (lossy) modes and the loss modes of its environment.
Any open quantum dynamics can be expressed as unitary evolution on of a larger system via unitary dilation. For a non-unitary transformation given by a rectangular matrix, T , with singular values all ≤ 1, we can write the unitary dilation (C7) With this, we are now ready to write down the probability for the general case of a Fock state input, linear optical experiment. We consider that we have M in input modes, initialised in the state n. These propagate through a nonunitary linear transformation, T , before being detected by M out threshold detectors, which give an outcome d, which is a length-M out bit-string. So T is given by an M out × M in matrix.
First, we dilate T according to Eq. (C7), giving an (M in +M out )-dimension unitary matrix. In this construction, U C = T and U B is given by I − T † T 1/2 . We also notice that we can explicitly write out the multiplication of the augmented matrices: where m = j d j = |C|, is the total number of clicks. We also define the matrix To arrive at this equation, we have used where the subscript notation :, n means that we select the columns according to n.

Appendix D: Unitary Bristolian reduction to the permanent for clicks equal to photons
When the number of input photons is equal to the number of threshold detector clicks, the Unitary Bristolian has a square matrix as an input. To see N clicks for an N photon input state, we know that each threshold detector must have detected exactly 1 photon, and so this event could also be modelled by using photon number projectors instead of the click measurement operator, and we can use Eq. 12 to calculate its amplitude. Therefore, we expect that the Bristolian of a square matrix should reduce to the absolute square of the permanent of the same matrix.
For a square N × N matrix, A, with elements a jk , we can write the Unitary Bristolian as Then we can expand the permanent using Ryser's formula This gives which is the Ryser-style formula for the 3-tensor permanent, as defined in Ref. [23] ubrs with σ and ρ being elements of the permutation group S N . However, this 3-tensor, B, with elements b ijk = a * ki a kj is very structured. This kind of structure is discussed in Ref. [23], where it is shown that it allows us to factor this expression into In the absence of our expression for the Bristolian in Eq. 15, the only known way in the literature to compute probabilities for Fock states measured with threshold detectors is to consider all possible events which could lead to the witnessed outcome, and sum all their corresponding probabilities.
To provide an example for how the complexity of this method compares to the Bristolian, we consider an M mode experiment with uniform transmission, η, N input photons and n clicks (with n ≤ N ). For interferometers with imbalanced losses, we must consider losses at both the input and output of the interferometer. However, here we are considering balanced loss which allows us to consider that any photon loss occurs before the photons reach the interferometer, as balanced loss commutes with linear optics. Because any number of photons between n and N can be transmitted through this loss channel and lead to an n click event, we must consider all N j=n N j possible configurations for how these photons could have been transmitted. For each input configuration, we then need to consider all the ways that these photons can bunch within the detectors. We know that at least 1 photon must arrive in each detector, which leaves j − n photons left, which can be configured in any arrangement. There are X+Y −1 Y ways of arranging Y photons into X modes, and here we have j − n photons which can be arrive in n modes, meaning we need to calculate j−n+n−1 j−n = j−1 j−n permanents for each input configuration. This results in a total of N j=n N j j−1 j−n permanents. Each permanent for a j photon configuration has a time complexity of j2 j , giving an overall time complexity which is lower bounded by N n n2 n . By comparison, the Bristolian complexity of n2 2n provides a superexponential speedup when N and N − n are large.

Appendix F: Derivation of the loop Torontonian
We start with some definitions. For a Gaussian state with complex Husimi covariance matrix Σ and complex vector of means α we define which uniquely specify the photon number statistics of the Gaussian state [22] via loop Hafnians [21]. Given a photon number outcome n = (n 1 , . . . , n M ) its probability is given by where p( 0) is the vacuum probability, X = 0 I I 0 and O n, n and γ n are submatrices of O and γ, found by using n j repetitions of the rows/columns corresponding to mode j (recalling that each mode j corresponds to 2 rows/columns of O).
We will now show that the threshold probabilities can also be written in terms of the quantities defined in Eq. (F1) and Eq. (F2).
By ordering our basis vectors such that modes which see a click (C) are arranged to be before modes which see vacuum (V ), the matrix O can be written in block form then, using Schur complements, we can see that We are now ready to investigate threshold probabilities. We start with Eq. 4 which we write as For a given Z, the term inside the sum can be written as in Eq. 23. Note that the argument inside the exponential in said equation can be rewritten as where W denotes the union of sets Z and V , and Y denotes the modes not included in W . Here we have used that We can similarly write the argument of the exponential in p( 0) as With these two expressions we can then write In the last equation we used the result in Eq. (F5), together with factorising the determinant This corresponds to the identity det(Σ W W ) = det(Σ) det([Σ −1 ] Y Y ) from Ref. [39]. With these observations we can write the sought after probability as where C is the set of modes with threshold detector clicks, and ltor stands for the loop Torontonian defined as for a 2m × 2m matrix O, and a 2m-length vector γ.
In arriving at this formula, we have swapped the summation over Z for a summation over Y , as we did in Eq. (C2). This is the form that the equation appears in the main text. Just like the Torontonian [7], the loop Torontonian has an interesting interpretation of being a generating function for photon number probabilities of Gaussian states (now with nonzero displacement).
To this end, recall that a threshold probability can be obtained as a sum (coarse graining) of many photon number events. For any d, we can write where C d is the set of all vectors of integers that have zero in the positions where d has zero and strictly positive integers in all the other positions. Note that any element k in the set satisfies K = M i=1 k i > 0. We now recall that loop hafnians satisfy the following scaling property [21] lhaf((ηA) n, n , ( with N = j n j . This allows us to write If we set = |C| we can derive a formula for the loop hafnian where we define One way of interpreting f (O, γ) is that it selects the th coefficient of the polynomial expansion of q. Therefore, we do not require knowledge of q(O, γ) beyond order-.
Using the Mercator series, we expand the denominator in Eq. (G7) as so we can combine the denominator into the exponent.
We also notice that we can use a binomial expansion, This gives us an expression which we Taylor-expand up to order to obtain We have used that γ * = X γ. This gives us the exact form of the trace formula algorithm for the loop Hafnian as presented in Ref. [21] (bar some typos in the referenced paper). It is interesting that this derivation uses various arguments which rely specifically on O and γ being formed by a Gaussian state covariance matrix (in particular Eq. (G1) and Eq. (G2)), however the trace formula algorithm that we are able to derive here is applicable to arbitrary symmetric matrices.

FIG. 2.
Scattershot construction. Gates between modes, symbolised by vertical lines with dots at the ends, show weak two-mode squeezing. The two-mode squeezing operator between modes i and j is givenŜi,

Appendix H: Connecting the Bristolian and the Torontonian
In this Appendix we provide a formal link between the Bristolian and the Torontonian. Concretely, we show that the Bristolian associated with the threshold detection of a multimode Fock state with single photon or vacuum inputs can be evaluated as a certain limit of a Torontonian.
The starting point of our derivation is the scattershot Boson sampling construction shown in Fig. 2 and introduced in Ref. [20]. The diagram represents M two-mode squeezed vacuum states where the first half of the modes (the heralded modes) are sent into an interferometer with transmission matrix T , and the second half of the modes (the heralding modes) are sent into threshold detectors.
For each individual two-mode squeezed vacuum with squeezing parameter t, it is straightforward to show that conditioned on a click in the heralding mode, the state of the heralded mode collapses to [40] ρ = 1 1 +n |1 1| + ∞ m=2 1 n n 1 +n m |m m|, (H1) wheren = sinh 2 t is the mean photon number of either mode of the two-mode squeezed vacuum. The probability of heralding the state ρ is given by Note that as ε → 0 the fidelity between ρ and a single photon Fock state approaches one, but at the same time the probability of heralding the state approaches zero.
We can now study the covariance matrix Σ of the Gaussian circuit in Fig. 2. After some algebra, it can be shown that the matrix O = I−Σ −1 dictating the threshold probabilities of the Gaussian state is given by We now want to calculate the probability that a subset of the first M modes clicks conditioned on all the modes on the second half clicking. We specify the modes that clicked in the first half by d and then can write the click pattern for the 2M modes to be e = d ⊕ 1 where, recall, 1 is the all ones vector of length M . We can now write the conditional probability as For the probability in the numerator of the right-hand side in the last equation we can use Eq. 25. The vacuum probability of the 2M -mode state is given by We can now write where C is the union of the {M + 1, . . . , 2M } heralding modes and the modes that have a one in the vector d, i.e., the modes that click in the first half.
As explained at the beginning of this section, in the limit where ε → 0 we know that the conditional state of the input heralded modes becomes a product of single photons in each mode, and thus in this same limit we can write Note that we have so far only considered the case where single photons are input in all the modes of the interferometer. The more general case where vacuum is fed into some of the modes can be dealt with by applying a loss channel with zero-transmission to the relevant modes. This is equivalent to setting to zero the columns of T where vacuum is fed. Using this argument we find that for cases where n is a bistring, i.e., we only allow single photons inputs, and setting N = i n i we can write where now C is the union of the labels of the modes in which single photons where input into the interferometers and the labels of the modes in which clicks are registered. It is interesting to consider that one could potentially derive the form of the Bristolian in terms of sums of permanents by using the connection between permanents and determinants provided by the MacMahon Master theorem [41].
Finally, note that one can also write circuits to herald multi-photon Fock states using only threshold detectors as shown in Appendix D of Ref. [40].

Appendix I: Time complexities
To calculate a marginal vacuum probability for a Fock state evolved through a linear interferometer, we can compute a permanent, as given in Eq. 13 If we are detecting vacuum in modes given by V and marginalising over all other modes, given by B, then the matrix [U † U x U ] n, n has rank ≤ |B|, so its permanent can be computed in time O(N O(|B|) ) using the algorithm introduced in Sec. III of Ref. [42]. If we are marginalising over many modes, and therefore |V | is small, we can instead consider using the algorithm introduced in Sec. IV. A of Ref. [15] to compute the permanent in time O(N 2|V |+1 ).
For a general permanent of an n × n matrix, Ryser's algorithm [37] has the best known complexity of O(n2 n ). For computing the Bristolian, to calculate the permanents inside the summation in Eq. 15, we may sometimes be able to use the faster algorithms above, however the dominant complexity for this formula would still come from computing the intermediate cases, when neither of the faster algorithms are applicable, where Ryser's algorithm may be the fastest option. This is upper bounded by O(n2 n ) and there are O(2 m ) terms in the sum. Therefore, a Bristolian of an m × n matrix A, and n × n matrix E has a time complexity of O(n2 m+n ).
We do not claim that these complexities are optimal. The structure of this matrix function may be exploited to reduce the complexity, for example by using methods similar to those for low rank permanents [42], exploiting recursion [24] and using Laplace expansions [43]. However, we leave it as an open problem to find faster algorithms for the Bristolian.
For the loop Torontonian, we find a comparable complexity to the original algorithm for the Torontonian, with complexity O(m 3 2 m ) for a 2m × 2m matrix O and 2m-length vector γ. In each step, we must compute a matrix inverse and a matrix determinant, both having m 3time algorithms. However, these steps can make use of the Cholesky decomposition of O, so we can improve the polynomial prefactor, following the methods described in Ref. [24].
Both of these methods see a quadratic penalty as compared to the fastest methods for calculating photon number probabilities of pure states. This can be understood as being caused by the threshold detection operators in Eq. 1b having high rank, whereas the photon number operators are rank-1 projectors. For Gaussian state calculations, we also see a quadratic cost for calculating photon number probabilities when the state is mixed [22,44], so we can also understand this quadratic penalty as being a result of introducing mixture into the projected state. This differs from the case of sampling, where it was shown that sampling threshold detector clicks can be simulated with the same complexity as single photon measurements [14]. If we accept approximate expressions, accurate to additive error, we can efficiently compute the probabilities presented in this work using Monte-Carlo phase space methods [45]. However, for events with small probability, these methods can quickly become impractical due to a large relative error.
Appendix J: Accuracy improvements of a threshold detection model over a single photon projection approximation To provide a quantitative demonstration of the accuracy improvements of an exact model of threshold detection versus the typical approximation of using single photon projection, we calculate the full probability distribution for a 4 photon Fock state boson sampling experiment [28] using the Bristolian with the number of modes ranging from 4 to 12, choosing a transitivity of η = 0.6. We also find the probability distribution given by 0 or 1 photon Fock state projective measurement. We calculate the total variation distance (TVD) between these distributions for 100 different Haar random linear optical interferometers for each number of modes. The results, shown in Fig. 3, show that even as we approach the M = N 2 regime, we do not see a convergence between these two distributions. In particular, we see that the TVD for all the experiments lies typically within the range of 5%-12%. This highlights the importance of using the correct mathematical description of the experiment in order to best understand the results. We also performed an analysis of data from a recent experiment of displaced Gaussian boson sampling [10], which used a two-mode squeezed vacuum and a single coherent state as input states. Here, we look at the 2 photon probability distribution as estimated using the experiment for different levels of displacement, as labelled by the mean photon number of the input coherent state, n α , in Fig. 4. For each probability distribution, we compare the TVD to a model which uses threshold detection, using the loop Torontonian, and a model which assumes Fock state projections, using the loop Hafnian. The experiment of Ref. [10] uses threshold detection, and we see that the model which uses loop Torontonians provides a more accurate model of the experiment. For the largest displacements, when looking at the ratio of the TVD, the loop Torontonian gives a probability distribution which is 16% closer to the experimental data.

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