I. INTRODUCTION

Section:

ChooseTop of pageABSTRACTI. INTRODUCTION <<II. GENERATING FUNCTIONS ...III. APPLICATION TO COMMO...IV. DETECTION PROBABILITI...V. DISCUSSIONVI. CONCLUSIONSREFERENCESRecent progress in the generation, manipulation, and detection of photonic quantum states has led to new applications in the field of photonic quantum information processing and sensing. An example is photonic quantum computing, which can be realized by using only single-photon sources, beam splitters, phase shifters, and photon detectors.11. E. Knill, R. Laflamme, and G. J. Milburn, “A scheme for efficient quantum computation with linear optics,” Nature 409, 46–52 (2001). https://doi.org/10.1038/35051009 Other applications such as quantum key distribution (QKD) or quantum imaging have gained considerable attention over the last few years and become commercially relevant, requiring the development of sophisticated optical setups working at the single-photon level.2,32. F. Xu, X. Ma, Q. Zhang, H.-K. Lo, and J.-W. Pan, “Secure quantum key distribution with realistic devices,” Rev. Mod. Phys. 92, 025002 (2020). https://doi.org/10.1103/revmodphys.92.0250023. P.-A. Moreau, E. Toninelli, T. Gregory, and M. J. Padgett, “Imaging with quantum states of light,” Nat. Rev. Phys. 1, 367–380 (2019). https://doi.org/10.1038/s42254-019-0056-0 Photon-number resolved (PNR) detection of quantum states opens new pathways to experiments and applications requiring the simulation of such experiments.The strength of nonlinear optical interactions between light and matter typically decreases rapidly with the order of the nonlinear effect. Therefore, many common photonic quantum states are described by Hamiltonians that are at most quadratic in the creation and annihilation operators. States with such Hamiltonians are called Gaussian states (GSs) and include, e.g., vacuum, coherent states, squeezed states, thermal states, or states generated by spontaneous parametric down-conversion (SPDC). Transformations by optical setups described by such Hamiltonians, introduced, e.g., by phase shifters or beam splitters, map GSs to other GSs. Hence, the capability to simulate the photon statistics of GSs is of great practical relevance. The covariance formalism describes GSs using a covariance matrix and a displacement vector and allows the modeling of many common effects on GSs in experiments, such as losses, phase shifts, and interference at beam splitters, using relatively simple matrix transformations of the covariance matrix and displacement vector.4–74. X.-B. Wang, T. Hiroshima, A. Tomita, and M. Hayashi, “Quantum information with Gaussian states,” Phys. Rep. 448, 1–111 (2007). https://doi.org/10.1016/j.physrep.2007.04.0055. C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, “Gaussian quantum information,” Rev. Mod. Phys. 84, 621–669 (2012). https://doi.org/10.1103/revmodphys.84.6216. S. Olivares, “Quantum optics in the phase space,” Eur. Phys. J.: Spec. Top. 203, 3–24 (2012). https://doi.org/10.1140/epjst/e2012-01532-47. G. Adesso, S. Ragy, and A. R. Lee, “Continuous variable quantum information: Gaussian states and beyond,” Open Syst. Inf. Dyn. 21, 1440001 (2014). https://doi.org/10.1142/s1230161214400010 Therefore, it lends itself to the implementation of simulations of quantum optical experiments. Importantly, the covariance formalism, described briefly in Appendix A 1, carries the full information about the photon statistics of the state.While non-PNR detectors such as single-photon avalanche photodiodes are common, a variety of detector types capable of PNR detection have been developed as well.88. A. Divochiy, F. Marsili, D. Bitauld, A. Gaggero, R. Leoni, F. Mattioli, A. Korneev, V. Seleznev, N. Kaurova, O. Minaeva, G. Gol’tsman, K. G. Lagoudakis, M. Benkhaoul, F. Lévy, and A. Fiore, “Superconducting nanowire photon-number-resolving detector at telecommunication wavelengths,” Nat. Photonics 2, 302–306 (2008). https://doi.org/10.1038/nphoton.2008.51 Transition edge sensors combine high detection efficiencies and PNR capabilities but can only be operated at comparatively low count rates.99. M. Schmidt, M. von Helversen, M. López, F. Gericke, E. Schlottmann, T. Heindel, S. Kück, S. Reitzenstein, and J. Beyer, “Photon-number-resolving transition-edge sensors for the metrology of quantum light sources,” J. Low Temp. Phys. 193, 1243–1250 (2018). https://doi.org/10.1007/s10909-018-1932-1 Superconducting nanowire single-photon detectors (SNSPDs) can be operated at higher count rates and allow to realize PNR detection by evaluating the electronic output pulse shape depending on the number of incident photons1010. C. Cahall, K. L. Nicolich, N. T. Islam, G. P. Lafyatis, A. J. Miller, D. J. Gauthier, and J. Kim, “Multi-photon detection using a conventional superconducting nanowire single-photon detector,” Optica 4, 1534–1535 (2017). https://doi.org/10.1364/optica.4.001534 or by subdividing the light-sensitive area into multiple pixels.11–1411. X. Tao, S. Chen, Y. Chen, L. Wang, X. Li, X. Tu, X. Jia, Q. Zhao, L. Zhang, L. Kang, and P. Wu, “A high speed and high efficiency superconducting photon number resolving detector,” Supercond. Sci. Technol. 32, 064002 (2019). https://doi.org/10.1088/1361-6668/ab079912. S. Steinhauer, S. Gyger, and V. Zwiller, “Progress on large-scale superconducting nanowire single-photon detectors,” Appl. Phys. Lett. 118, 100501 (2021). https://doi.org/10.1063/5.004405713. G. He, H. Li, R. Yin, L. Zhang, D. Dong, J. Lv, Y. Fei, X. Wang, Q. Chen, F. Li, H. Li, H. Wang, X. Tu, Q. Zhao, X. Jia, J. Chen, L. Kang, and P. Wu, “Simultaneous resolution of photon numbers and positions with series-connected superconducting nanowires,” Appl. Phys. Lett. 120, 124001 (2022). https://doi.org/10.1063/5.008474414. R. Cheng, Y. Zhou, S. Wang, M. Shen, T. Taher, and H. X. Tang, “A 100-pixel photon-number-resolving detector unveiling photon statistics,” Nat. Photonics 17, 112–119 (2022). https://doi.org/10.1038/s41566-022-01119-3 SNSPD detectors achieving photon-number resolution with eight pixels are commercially available.1515. ID Quantique SA, id281 superconducting nanowire series, 2022, product Brochure. Not only SNSPDs but also single-photon avalanche detectors have been multiplexed, in time16,1716. D. Achilles, C. Silberhorn, C. Sliwa, K. Banaszek, I. A. Walmsley, M. J. Fitch, B. C. Jacobs, T. B. Pittman, and J. D. Franson, “Photon-number-resolving detection using time-multiplexing,” J. Mod. Opt. 51, 1499–1515 (2004). https://doi.org/10.1080/0950034040823528817. R. Kruse, J. Tiedau, T. J. Bartley, S. Barkhofen, and C. Silberhorn, “Limits of the time-multiplexed photon-counting method,” Phys. Rev. A 95, 023815 (2017). https://doi.org/10.1103/physreva.95.023815 or space,1818. R. Heilmann, J. Sperling, A. Perez-Leija, M. Gräfe, M. Heinrich, S. Nolte, W. Vogel, and A. Szameit, “Harnessing click detectors for the genuine characterization of light states,” Sci. Rep. 6, 19489 (2016). https://doi.org/10.1038/srep19489 to realize PNR detectors. In such multiplexing setups, multiple photons can hit the same sub-detector so that the photons are not resolved. However, with an increasing number of detectors, the counting statistics converge on the actual PND.18,1918. R. Heilmann, J. Sperling, A. Perez-Leija, M. Gräfe, M. Heinrich, S. Nolte, W. Vogel, and A. Szameit, “Harnessing click detectors for the genuine characterization of light states,” Sci. Rep. 6, 19489 (2016). https://doi.org/10.1038/srep1948919. J. Sperling, W. Vogel, and G. S. Agarwal, “True photocounting statistics of multiple on-off detectors,” Phys. Rev. A 85, 023820 (2012). https://doi.org/10.1103/physreva.85.023820The task of finding the detection probability p(n) = p(n1, …, nS) for n1 photons in mode 1, n2 photons in mode 2, etc. of a GS with S modes, i.e., its photon number distribution (PND), is known as the Gaussian boson sampling (GBS) problem.2020. C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, “Gaussian boson sampling,” Phys. Rev. Lett. 119, 170501 (2017). https://doi.org/10.1103/physrevlett.119.170501 Experimentally, GBS has been realized, e.g., using networks of beam splitters on photonic chips.21,2221. M. Tillmann, B. Dakić, R. Heilmann, S. Nolte, A. Szameit, and P. Walther, “Experimental boson sampling,” Nat. Photonics 7, 540–544 (2013). https://doi.org/10.1038/nphoton.2013.10222. M. A. Broome, A. Fedrizzi, S. Rahimi-Keshari, J. Dove, S. Aaronson, T. C. Ralph, and A. G. White, “Photonic boson sampling in a tunable circuit,” Science 339, 794–798 (2013). https://doi.org/10.1126/science.1231440 Large-scale GBS experiments have been realized, for example, in the context of quantum computing and to pursue the demonstration of the computational advantage of quantum computers over classical computers.23–2523. D. J. Brod, E. F. Galvão, A. Crespi, R. Osellame, N. Spagnolo, and F. Sciarrino, “Photonic implementation of boson sampling: A review,” Adv. Photonics 1, 034001 (2019). https://doi.org/10.1117/1.AP.1.3.03400124. 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, 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, C.-Y. Lu, and J.-W. Pan, “Quantum computational advantage using photons,” Science 370, 1460–1463 (2020). https://doi.org/10.1126/science.abe877025. 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). https://doi.org/10.1103/physrevlett.127.180502 Reference 2626. T. R. Bromley, J. M. Arrazola, S. Jahangiri, J. Izaac, N. Quesada, A. D. Gran, M. Schuld, J. Swinarton, Z. Zabaneh, and N. Killoran, “Applications of near-term photonic quantum computers: Software and algorithms,” Quantum Sci. Technol. 5, 034010 (2020). https://doi.org/10.1088/2058-9565/ab8504 discusses possible applications for GBS-based quantum computing. The computational complexity of simulating GBS has been investigated as a benchmark for optical quantum computing comparing it to GBS simulations on classical computers.20,27–3420. C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, “Gaussian boson sampling,” Phys. Rev. Lett. 119, 170501 (2017). https://doi.org/10.1103/physrevlett.119.17050127. 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, 2011), pp. 333–342.28. A. P. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L. O’Brien, and T. C. Ralph, “Boson sampling from a Gaussian state,” Phys. Rev. Lett. 113, 100502 (2014). https://doi.org/10.1103/physrevlett.113.10050229. S. Rahimi-Keshari, A. P. Lund, and T. C. Ralph, “What can quantum optics say about computational complexity theory?,” Phys. Rev. Lett. 114, 060501 (2015). https://doi.org/10.1103/PhysRevLett.114.06050130. N. Quesada, J. M. Arrazola, and N. Killoran, “Gaussian boson sampling using threshold detectors,” Phys. Rev. A 98, 062322 (2018). https://doi.org/10.1103/physreva.98.06232231. 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). https://doi.org/10.1103/physreva.100.03232632. N. Quesada, R. S. Chadwick, B. A. Bell, J. M. Arrazola, T. Vincent, H. Qi, and R. García−Patrón, “Quadratic speed-up for simulating Gaussian boson sampling,” PRX Quantum 3, 010306 (2022). https://doi.org/10.1103/prxquantum.3.01030633. A. Björklund, B. Gupt, and N. Quesada, “A faster hafnian formula for complex matrices and its benchmarking on a supercomputer,” ACM J. Exp. Algorithmics 24, 1–17 (2019). https://doi.org/10.1145/332511134. J. F. 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, I. A. Walmsley, and A. Laing, “The boundary for quantum advantage in Gaussian boson sampling,” Sci. Adv. 8, eabl9236 (2022). https://doi.org/10.1126/sciadv.abl9236 The fact that GBS is investigated in the context of computational quantum supremacy indicates that calculating the solution to GBS problems can require substantial computational resources. The PND can be calculated by evaluating expressions involving matrix functions called the Hafnian and loop Hafnian for PNR detection as well as functions called the Torontonian and loop Torontonian for non-PNR detection.20,30,33,35,3420. C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, “Gaussian boson sampling,” Phys. Rev. Lett. 119, 170501 (2017). https://doi.org/10.1103/physrevlett.119.17050130. N. Quesada, J. M. Arrazola, and N. Killoran, “Gaussian boson sampling using threshold detectors,” Phys. Rev. A 98, 062322 (2018). https://doi.org/10.1103/physreva.98.06232233. A. Björklund, B. Gupt, and N. Quesada, “A faster hafnian formula for complex matrices and its benchmarking on a supercomputer,” ACM J. Exp. Algorithmics 24, 1–17 (2019). https://doi.org/10.1145/332511135. N. Quesada, “Franck-Condon factors by counting perfect matchings of graphs with loops,” J. Chem. Phys. 150, 164113 (2019). https://doi.org/10.1063/1.508638734. J. F. 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, I. A. Walmsley, and A. Laing, “The boundary for quantum advantage in Gaussian boson sampling,” Sci. Adv. 8, eabl9236 (2022). https://doi.org/10.1126/sciadv.abl9236 The operation count for the evaluation of these functions scales exponentially with the number of detected photons, and optimization of the algorithms is an active field of research.3434. J. F. 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, I. A. Walmsley, and A. Laing, “The boundary for quantum advantage in Gaussian boson sampling,” Sci. Adv. 8, eabl9236 (2022). https://doi.org/10.1126/sciadv.abl9236

The previously mentioned methods for GBS computations have been designed to efficiently calculate PNDs with many photons but do not offer much flexibility. Typically, quantum optical setups not designed for the particular task of GBS quantum computation are operated at low photon numbers. Examples are applications in quantum ghost imaging and QKD, where coincidences involving only a relatively small number of detectors are relevant. If such a setup is to be simulated, minimizing the required computational resources is not necessarily the main concern. However, often a simulation method is required that is flexible enough to include the imperfections that are inevitably present in real setups. In this paper, we present such a simulation method. The purpose of our method is not to compete with specialized GBS algorithms in terms of performance, but rather to provide a flexible and simple way to compute the photon statistics for imperfect setups.

Two relevant effects in experimental setups are the simultaneous detection of multiple modes by the same detector and noise in the detection process. These effects require the convolution of probability distributions, which can conveniently be expressed by the multiplication of the corresponding probability-generating functions. Therefore, we derive generating functions of the photon statistics in terms of the covariance matrix and displacement vector in Sec. . Our generating-function approach is an alternative to the established expressions for the PND involving Hafnian-type functions. This different perspective on the GBS problem allows for deriving several related expressions for the generating functions of cumulative probabilities, raw or central moments, and rising or falling factorial moments of the detection statistics in Sec. , for which so far no systematic method existed. Furthermore, our method allows us to calculate the same quantities for certain non-Gaussian states called photon-added and photon-subtracted GSs, as shown in Sec. , considerably extending the range of possible applications.Generating functions need to be differentiated repeatedly to retrieve probabilities or moments. We show that these derivatives can be evaluated numerically by automatic differentiation (AD) without much effort. An advantage of AD is that it provides accurate numerical results while hiding the whole complexity of the calculation from the user. We discuss the usage of AD for our application in Sec. and provide an implementation example in Ref. 3636. E. Fitzke, F. Niederschuh, and T. Walther, “Simulating the photon statistics of Gaussian states employing automatic differentiation from PyTorch,” Technical Report, 2022..Our simulation method consists of two steps, connected by the generating function. First, the quantum state and the optical setup are modeled in the covariance formalism. Second, the photon statistics are obtained by applying the AD method to the generating function expressed in terms of the covariance matrix and displacement vector. Both steps can be easily implemented with widely available software, making the method very practical. However, it is expected that optimized algorithms for the evaluation of Hafnian-type functions will outperform general-purpose AD algorithms for the particular task of calculating the PND. To show that our method can nevertheless be used to simulate various aspects of Gaussian photon statistics efficiently for small photon numbers, we apply it to multiple examples with common GSs in Sec. .As a more complex application, we present in Sec. a simulation of an entanglement-based QKD system we recently presented in Ref. 3737. E. Fitzke, L. Bialowons, T. Dolejsky, M. Tippmann, O. Nikiforov, T. Walther, F. Wissel, and M. Gunkel, “Scalable network for simultaneous pairwise quantum key distribution via entanglement-based time-bin coding,” PRX Quantum 3, 020341 (2022). https://doi.org/10.1103/prxquantum.3.020341. It demonstrates the strengths of our method to readily incorporate relevant imperfections such as noise, detection efficiencies, and simultaneous detection of multiple modes. Although the experimental setup is complex, the simulation is straightforward because it can be broken down into elementary operations described by basic Gaussian state transformations. The simulation results are in very good agreement with the experimental values. Furthermore, we show that the contribution of multi-photon pair emission to the quantum bit error rate of the QKD system can be estimated by applying Bayes’ theorem to PNR simulation results.

II. GENERATING FUNCTIONS FOR THE DETECTION STATISTICS OF GAUSSIAN STATES

Section:

ChooseTop of pageABSTRACTI. INTRODUCTIONII. GENERATING FUNCTIONS ... <<III. APPLICATION TO COMMO...IV. DETECTION PROBABILITI...V. DISCUSSIONVI. CONCLUSIONSREFERENCESOur method to simulate the photon statistics of GSs uses generating functions for the photon statistics in terms of the covariance matrix Γ and the displacement vector d. They connect the covariance formalism to the photon statistics and are therefore essential for our simulation method. In Appendix A 2, we briefly summarize the relevant properties of generating functions for probability distributions. In the following Sec. , we briefly recap the formulation of the PND in terms of a generating operator and extend it to generating operators for moments and factorial moments.

A. Generating operators for photon detection probabilities and moments

Assuming a photon detection process where photons in a single mode are detected independently of each other with an efficiency η, the probability to detect n photons from a quantum state ϱ is given by38–4138. P. L. Kelley and W. H. Kleiner, “Theory of electromagnetic field measurement and photoelectron counting,” Phys. Rev. 136, A316–A334 (1964). https://doi.org/10.1103/physrev.136.a31639. M. Lax and M. Zwanziger, “Exact photocount statistics: Lasers near threshold,” Phys. Rev. A 7, 750–771 (1973). https://doi.org/10.1103/physreva.7.75040. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, Cambridge, 1995).41. D. Walls and G. J. Milburn, Quantum Optics, 2nd ed. (Springer, 2008).p(n)=:(ηN̂)nn!e−ηN̂:.(1)Here, normal order is indicated by :: and the photon number operator is N̂=â†â. By inserting Eq. (1) into the defining equation of a probability generating function h(y) [cf. Eq. (A14)], the PND can be obtained from the expectation value of a generating operator ĥ(y).38–4138. P. L. Kelley and W. H. Kleiner, “Theory of electromagnetic field measurement and photoelectron counting,” Phys. Rev. 136, A316–A334 (1964). https://doi.org/10.1103/physrev.136.a31639. M. Lax and M. Zwanziger, “Exact photocount statistics: Lasers near threshold,” Phys. Rev. A 7, 750–771 (1973). https://doi.org/10.1103/physreva.7.75040. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, Cambridge, 1995).41. D. Walls and G. J. Milburn, Quantum Optics, 2nd ed. (Springer, 2008). This method has recently been used by Nauth to express the PND of biphoton states in terms of the covariance matrix.4242. J. K. Nauth, “Full time-dependent counting statistics of highly entangled biphoton states,” Phys. Rev. A 106, 053716 (2022). https://doi.org/10.1103/physreva.106.053716 We extend this procedure to obtain similar generating operators for the moment generating function M(μ, y) and the rising factorial moment generating function R(y) as described in Appendix A 2,ĥ(y)=:eη(y−1)N̂:,(2)M̂(μ,y)=e−yμ:expη(ey−1)N̂:,(3)R̂(y)=11−y:eηyN̂/(1−y):.(4)These operators can be modified to take into account noise in the detection process. For that, according to the multiplication rule for generating functions (cf. Appendix A 2), the generating operator is simply multiplied by the generating function of the noise process. As an example, we consider noise with Poissonian statistics pnoise(n) = e−ννn/n! and noise parameter ν. The probability generating function of the noise is given by hnoise(y) = eν(y−1), so that the generating operators from Eqs. (2)–(4) including noise read eν(y−1)ĥ(y), expν(ey−1)M̂(μ,y) and eνy/(1−y)R̂(y).Similarly, noise with different statistics could be taken into account by multiplying the generating operators from Eqs. (2)–(4) with the generating function of the respective noise process. Another option is to include noise in the covariance formalism. Thermal noise can be modeled by coupling in a thermal state with a beam splitter or a matrix can be added to the covariance matrix to represent, e.g., classical Gaussian noise or noise from amplification.4343. J. Eisert and M. M. Wolf, Gaussian Quantum Channels (Imperial College Press, distributed by Wold Scientific Publishing Co., 2007), pp. 23–42.Equations (2)–(4) have in common that they involve special cases of the generating operator ĝ(w(y))=:exp(−w(y)â†â): for different functions w(y).In Refs. 4444. G. Adam, “Density matrix elements and moments for generalized Gaussian state fields,” J. Mod. Opt. 42, 1311–1328 (1995). https://doi.org/10.1080/09500349514551141 and 4545. J. Perina, Quantum Statistics of Linear and Nonlinear Optical Phenomena, 2nd ed. (Springer Science & Business Media, 1991)., the operator ĝ has been extended toĝ(u,v,w)=:exp(uâ+vâ†−wâ†â):(5)and has been related to the density matrix elements of ϱ̂ by introducing u and v.

B. Generating functions in terms of the covariance matrix and displacement vector

In experiments, often multiple modes, e.g., different polarization directions or frequency modes, enter the same detector. For example, calculating the simultaneous detection of multiple modes is required for modeling imperfect interference due to a mode mismatch with the model from Ref. 4646. M. Takeoka, R.-B. Jin, and M. Sasaki, “Full analysis of multi-photon pair effects in spontaneous parametric down conversion based photonic quantum information processing,” New J. Phys. 17, 043030 (2015). https://doi.org/10.1088/1367-2630/17/4/043030, which we use in Sec. . Furthermore, often joint detection probabilities between multiple detectors such as coincidence probabilities are of interest. Therefore, we generalize the calculation to multi-mode states.In Appendix B 2, we show that for multiple modes with vectors u, v, and w, G(u,v,w)=⟨ĝ(u,v,w)⟩ can be expressed in terms of the covariance matrix Γ and the displacement vector d byG(u,v,w)=exp(−zTΛ−1Wz/2+Z)detΛ,(6)with the diagonal matrix W = diag(w) ⊕ diag(w), Λ=1+W(Γ−1)/2, Z=∑susvsws−1 and z = d + ζ with ζ as defined in Eq. (B9). Equation (6) is the generating function for the multivariate photon statistics, from which the probabilities and moments are retrieved by repeated differentiation w.r.t. D parameters y1…yD for the detectors d = 1…D, detecting Md modes with additional Poissonian noise110110. Similarly, the generating function of noise processes with a different, non-Poissonian statistics could be taken into account by multiplying with the corresponding generating function. νd. The total number of modes is S = ∑dMd, and the index s runs over all mode indices, i.e., enumerates md = 1d…Md for all D detectors in the order 11, 21, …, M1, 12, …M2, …, MD. Abbreviating G(w) = G(0, 0, w), we obtain by employing a multi-index notation,111111. For tuples x and k we write xk=∏ixiki, so that e.g. w−1=∏sws−1. We also use n! = ∏dnd! and ∂yn=∏d∂nd/∂ydnd.p(n,ν,η)=⟨Π̂N=n⟩=1n!∂ynexp∑d(yd−1)νdG(w)y=0,withp(N≤n,ν,η)=⟨Π̂N≤n⟩=1n!∂yn(1−y)−1⁡×exp∑d(yd−1)νdG(w)y=0,withM(μ,k,ν,η)=⟨(N̂−μ)k⟩=∂ykexp∑d(eyd−1)νd−μdyd×G(w)y=0,withn(k)(ν,η)=∏dN̂d(kd)=∂ykexp∑dydνdG(w)y=0,withn(k)