Simulating the photon statistics of multimode Gaussian states by automatic differentiation of generating functions

Recent 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.¹ 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,3 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–7 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.⁸ Transition edge sensors combine high detection efficiencies and PNR capabilities but can only be operated at comparatively low count rates.⁹ 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 photons¹⁰ or by subdividing the light-sensitive area into multiple pixels.^11–14 SNSPD detectors achieving photon-number resolution with eight pixels are commercially available.¹⁵ Not only SNSPDs but also single-photon avalanche detectors have been multiplexed, in time^16,17 or space,¹⁸ 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,19

The task of finding the detection probability p(n) = p(n₁, …, n_S) for n₁ photons in mode 1, n₂ 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.²⁰ Experimentally, GBS has been realized, e.g., using networks of beam splitters on photonic chips.^21,22 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–25 Reference 26 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–34 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,34} 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.³⁴ 

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. II. 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. II B, 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. II C, 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. II D and provide an implementation example in Ref. 36.

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. III.

As a more complex application, we present in Sec. IV a simulation of an entanglement-based QKD system we recently presented in Ref. 37. 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.

Our 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. II A, 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.

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 by^38–41 

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–41 This method has recently been used by Nauth to express the PND of biphoton states in terms of the covariance matrix.⁴² 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,

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 p_noise(n) = e^−ννⁿ/n! and noise parameter ν. The probability generating function of the noise is given by h_noise(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.⁴³ 

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. 44 and 45, the operator $ĝ$ has been extended to

and has been related to the density matrix elements of $ϱ̂$ by introducing u and v.

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. 46, which we use in Sec. IV. 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 by

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 y₁…y_D for the detectors d = 1…D, detecting M_d modes with additional Poissonian noise¹¹⁰ ν_d. The total number of modes is S = ∑_dM_d, and the index s runs over all mode indices, i.e., enumerates m_d = 1_d…M_d for all D detectors in the order 1₁, 2₁, …, M₁, 1₂, …M₂, …, M_D. Abbreviating G(w) = G(0, 0, w), we obtain by employing a multi-index notation,¹¹¹ 

with

with

with

with

with

Equations (6)–(11) connect the covariance formalism and the photon statistics by repeated derivatives and are therefore crucial for our simulation method. Alternative expressions to calculate p(n), restricted to the noiseless case, can be found in the literature and are discussed and compared to Eq. (7) in Sec. V. Our generating-function approach additionally yields Eqs. (8)–(11) for the cumulative probabilities, raw moments (for μ = 0), central moments (for $μ=⟨N̂⟩$⁠), and rising and falling factorial moments directly incorporating the detection efficiency and noise. For general multi-mode GSs, to the best of our knowledge, expressions for these quantities have so far not been presented in the literature.

The cumulative probabilities are useful when ranges of photon numbers are to be analyzed jointly, especially when detection events in multiple detectors are considered. For example, the calculation of the coincidence probability to detect $≤n1$ photons in one detector and $≤n2$ in a second detector from the PND requires the evaluation of the probabilities for (n₁ + 1)(n₂ + 1) combinations of photon numbers, but only one evaluation of Eq. (8). The mean and variance are two of the most important moments of the PND and can be directly calculated from Eq. (9). The factorial moments are useful for the calculation of photon statistics of photon-added and photon-subtracted states (cf. Sec. II C).

By using a formulation with detection operators $Π̂$⁠, more complex detection events involving multiple detectors can be calculated. For example, the probability for the detection of n_A photons in detector A and n_B photons in detector B is given by $⟨Π̂NA=nB,NB=nB⟩=⟨Π̂NA=nAΠ̂NB=nB⟩$⁠, and the probability for n₁ or n₂ photons in the same detector with n₁ ≠ n₂ is given by $⟨Π̂N=n1 or N=n2⟩=⟨Π̂N=n1+Π̂N=n2⟩$⁠. Operators for complementary events can be defined as well, enabling the calculation of the probability to detect any photon number except for n or for more than n photons with $Π̂N≠n=1̂−Π̂N=n$ and $Π̂N>n=1̂−Π̂N≤n$⁠. The detection operators for different modes are then joined before the generating function is evaluated. This procedure was used, for example, in Ref. 46 to relate the detection probabilities for non-PNR detectors to the detection of vacuum by considering operators $Π̂N>0=1̂−Π̂N=0$⁠. We extend this method to PNR detection. As an example, consider the probability to not detect n_A photons in detector A and to detect more than n_B photons in detector B,

The last term, for example, is given by

where the modes for A and B are kept in Λ and d.

Finally, we derive multivariate expressions for matrix elements of $ϱ̂$⁠. In Appendix B 1, we present a simple derivation of the matrix elements in terms of $g(u,v,w)=⟨ĝ(u,v,w)⟩$⁠. We generalize the single-mode expressions from Eqs. (B3) and (B5) to the multimode case by using Eq. (6),

In Eq. (15), we define l by l_s = min(n_s, m_s) as well as Δn = n − l and Δm = m − l. An equation similar to Eq. (15) with l = 0, Δn = n,and Δm = m can be derived differently⁴⁷ and can be related to an expression involving the loop Hafnian function. Compared to this expression, our Eq. (15) has the advantage that for each mode the number of derivatives to be evaluated is reduced from n_s + m_s to max(n_s, m_s). When the derivatives are evaluated numerically, for example, by automatic differentiation (AD), our expression is especially advantageous for close-by values of n_s and m_s.

The range of possible applications for our simulation method of the photon statistics can be extended by adapting it to apply to certain non-Gaussian states.

A special class of states derived from GSs $ϱ̂$ are photon-added GSs $ϱ̂+k$ and photon-subtracted GSs $ϱ̂−k$⁠,⁴⁸ 

with k = k₁, k₂, …, etc. photons added or subtracted from the modes 1, 2, …, etc. Here, the normally and anti-normally ordered moments⁴⁴ $m(k)=tr(âkϱ̂â†k)=⟨â†kâk⟩$ and $m(k)=tr(â†kϱ̂âk)=⟨âkâ†k⟩$ ensure normalization. They can be obtained from the falling and rising factorial moments of the GS by setting η = 1 and ν = 0,

Similarly, we directly obtain m^(k) = n^(k)(ν = 0, η = 1). Here, our Eqs. (10) and (11) for factorial moments are helpful because they allow to calculate n_(k) and n^(k) directly.

Photon addition and photon subtraction can have non-trivial effects on the PND of quantum states. An example is photon subtraction from a thermal state, which increases the mean photon number of the state.^48,49 Photon-subtracted states can be approximately realized by inserting a beam splitter with high transmission into the beam and by conditioning the detection of the transmitted quantum state on the detection of a reflected photon.^48–50 Photon addition has been realized by seeding an SPDC process with low gain so that the SPDC signal modes are emitted into the mode of the incident state. The generation of the photon-added state is then conditioned on the detection of an SPDC idler photon.^51,52 Photon addition and subtraction can be used, for example, to enhance the signal-to-noise ratio in quantum ghost imaging.⁵³ 

The matrix elements of photon-added and photon-subtracted GSs can be directly calculated by using $âk|l〉=l!/(l−k)!|l−k〉$⁠, given that l ≥ k, and $â†k|l〉=(k+l)!/l|k+l〉$⁠,

The expressions for $⟨n+k|ϱ̂|l+k⟩$ and $⟨n−k|ϱ̂|l−k⟩$ can be calculated from the GS’s matrix elements by using Eq. (15).

To derive explicit expressions for the photon statistics, we first consider the single mode case of $ϱ̂−k$⁠. The operator $â†kĝâk$ is in normal order and we obtain from the optical equivalence theorem in Eq. (B1),

By comparison to the phase-space representation of $ĝ$ from Eq. (B2), the expectation value of $ĝ$ with respect to $ϱ̂−k$ can be expressed in terms of the generating function of the underlying GS $ϱ̂$⁠,

For the photon-added GS, we consider for simplicity only G(w), as u and v are essentially only required for the calculation of the matrix elements discussed earlier. In Appendix B 3, we derive the expression for the single-mode case , which generalizes to

with $w̃s=1−[(1−ws)−1−rs]−1$⁠.

This means that all the quantities of the photon statistics for which we have derived expressions above can also be calculated for multi-mode photon-added and photon-subtracted GSs. In this context, $G(w)=tr(ϱ̂ĝ(w))$ in Eqs. (7)–(11) is replaced by the expressions for $tr(ϱ̂−kĝ(w))$ and $tr(ϱ̂+kĝ(w))$ from Eqs. (21) and (22).

To our knowledge, expressions for the PND of photon-added and photon-subtracted GSs have so far only been reported for single-mode states.^54–58 Our expressions enable the calculation of the photon statistics for the general multi-mode case and additionally enable the calculation of moments, factorial moments, and matrix elements. Moreover, the expressions are evaluated by repeated differentiation and can, therefore, be calculated with the same tools as the expressions for GSs.

Certain non-Gaussian states $ϱ̂′$⁠, such as NOON states and cat states, can be generated from larger GSs $ϱ̂$ by PNR detection.^47,59–61 The effects of the optical setup on such states can be modeled using the covariance formalism as described in Ref. 60. Consider a GS $ϱ̂$ with (M′ + M) modes, of which M modes are detected with PNR detectors so that a state $ϱ̂′$ with M′ modes remains. By conditioning the further analysis of $ϱ̂′$ on the detected PND, the non-Gaussian states are produced probabilistically.⁶⁰ Unitary transformations U′ acting on $ϱ̂′$ can simply be modeled by applying the unitary transformation $U=U′⊗1̂M×M$ to $ϱ̂$⁠.⁶⁰ For the simulation, the required PNR detection pattern on the M modes is fixed and the photon detection pattern of interest is applied to the remaining M′ modes. Unitary transformations that are linear or quadratic in $â$ and $â†$ can be applied to $ϱ̂$ in the covariance matrix formalism before the detection.⁶⁰ By using this procedure from Ref. 60, our method can be directly applied to simulate states that can be obtained from GSs via PNR detection.

To summarize, we note that our method to simulate the photon statistics of GSs can also be applied to photon-added and photon-subtracted GSs as well as to non-Gaussian states derived from larger GSs by PNR detection. The class of states that can be simulated and the range of possible applications of the simulation method are thereby greatly extended.

Multiple options exist to evaluate the derivatives of the generating functions. One option used in Ref. 62 is to use finite-difference approximations, but this method can accumulate numerical inaccuracies. Another option to retrieve probabilities from the generating function is to approximate Cauchy’s integral formula on a circle around the origin in the complex plane. In Ref. 63, this method has been applied to probability generating functions, and error bounds for the approximations have been given. It is extensively discussed in Ref. 64, and expressions for moments and multi-dimensional distributions are given in Refs. 65–67.

A versatile method to evaluate derivatives of functions numerically is automatic differentiation (AD).^68,69 AD breaks down a function to be differentiated into elementary operations such as addition, multiplication, or sin(x). The differentiation rules for these operations are then applied and combined according to the chain rule to the derivative of the overall function, which is finally evaluated. For example, instead of approximating the derivative of sin(x) numerically, AD uses the fact that the derivative is given by cos(x) and evaluates the cosine function. Therefore, in contrast to finite-difference approximations, the derivatives obtained are accurate to their working precision.

Applying AD to the generating functions allows us to readily calculate all the quantities we have derived. Equations (7)–(11) as well as Eqs. (21) and (22) only involve derivatives and basic functions as well as linear algebra that are easily implemented. With modern software tools, multivariate higher-order derivatives can be conveniently obtained with a few lines of code. From a practical point of view, this is convenient because the user is not confronted with the complexity of the calculation, which is hidden in the repeated derivatives. AD thereby greatly improves the practical usability of the generating functions we derived. However, as mentioned in the introduction, the computational complexity to calculate the PND of GBS problems with state-of-the-art algorithms scales exponentially with the number of photons, so that the numerical evaluation of the derivatives of the generating function via AD will become infeasible for large photon numbers.

A variety of AD software exists, and an overview of available tools is provided in Ref. 70. AD is widely used in machine learning, for example, for training artificial neural networks.⁷¹ Therefore, AD functionalities are part of popular machine learning libraries such as PyTorch^72,73 and TensorFlow.⁷⁴ Machine learning has gained considerable attention over the past few years, and the range of its applications is consistently growing. Consequently, software for machine learning will be further developed, and it can, therefore, be expected that AD software will gain further capabilities during the next few years, which of course is beneficial for our application. State-of-the-art machine learning libraries provide numerous possibilities to optimize their performance, for example, by using parallel computing on graphics processing units (GPUs).

Throughout this paper, we use the autograd function from the machine learning framework PyTorch 1.11.0 for AD. We use PyTorch in a very basic configuration, i.e., we do not use the option for acceleration by GPU computations, and the only setting we change is the default precision, which we increase from float32 to float64. All simulations are run on a regular desktop computer to show that for small photon numbers, the generating functions can be evaluated without much effort, demonstrating the practicability of our approach. In Ref. 36, we present a straightforward code example to demonstrate the use of PyTorch for computing the PND.

In this section, we apply our simulation method to two common GSs with well-known photon statistics: the single-mode displaced squeezed thermal state and the two-mode squeezed vacuum with multiple squeezers. These practical examples demonstrate that our simulation method can easily take into account effects such as non-unity detection efficiencies, noise, and multiple modes entering the same detector.

As a first example, we consider the single-mode displaced squeezed thermal state,⁶ 

with displacement operator $D(α)=exp(αâ†−α*â)$ and squeezing operator $S(χ)=exp[(χâ†2−χ*â2)/2]$⁠. Here, χ = re^iθ is the squeezing parameter and

is a thermal state with mean photon number μ_th.⁶ 

By using Eq. (6), we derive an analytic expression for the generating function,¹¹² 

Here, we use the abbreviations

and

Meanwhile, an explicit expression for the PND can be derived from this generating function; we do not present it here as it can be found in Ref. 56, where it was derived using a different approach. The photon statistics of the single-mode squeezed thermal state, the displaced squeezed thermal state, and the squeezed displaced thermal state have been studied in Refs. 75 and 76.

As a simplification, we now consider a single-mode displaced squeezed state |ψ⟩ = D(α)S(χ)|0⟩, i.e., with μ_th = 0, for which the PND is given by⁵⁶ 

where H_n is the n-th Hermite polynomial.

In Fig. 1, we verify that deriving the PND of the state from the generating function Eqs. (7) and (25) for η_1/2 = 1 and ν_1/2 = 0 by AD yields the same results as the analytic expression from Eq. (28) within a tolerance of a few units of least precision (ULP) of the analytical value. The maximum absolute value of the deviation from Eqs. (7) and (25) to the analytical result from Eq. (28) is 118 ULP, i.e., the relative error of the numerical results is less than 2 × 10⁻¹⁴ for both methods. This is remarkable because for n = 11, computing the result via Eq. (7) already took a few minutes and required multiple GB of RAM, indicating a considerable number of numerical operations. The fact that the result is, nevertheless, numerically accurate underlines that the AD capabilities of PyTorch well suit our simulation method.

The analytic formula has the advantage that it can be quickly evaluated even for high photon numbers. However, an advantage of our formulation with a generating function in terms of Λ and d is that effects such as noise and detection efficiencies can easily be taken into account (cf. Fig. 1).

A strength of our approach is that the detection statistics of states with multiple modes entering the same detector can be easily evaluated. For example, the spectral degree of freedom can be incorporated into a simulation by discretizing the frequency spectrum and treating each frequency as a separate mode. For the simulation in Sec. IV, we want to take into account the photon statistics, but we do not resolve the frequency spectrum of the SPDC process. Because the spectrum and the photon statistics are related to each other, we recap some basic properties of the multi-mode SPDC photon statistics.

SPDC with one signal mode s and one idler mode i generates the two-mode squeezed vacuum (TMSV) state $|ψ〉TMSV=exp(reiϑâs†âi†−re−iϑâsâi)|0〉$ with mean photon number μ = sinh²r , which can be written as⁴ 

The photon statistics of either the signal or idler mode alone are, therefore, given by p(n) = sech²(r) tanh²ⁿr, i.e., resemble those of a thermal state [cf. Eq. (24)], with the probability generating function h(y) = sech²(r)/(1 − y tanh²r).⁷⁸ In practice, the bandwidths of the parametric process and the pump light result in two-mode squeezing between a number M of pairs of orthogonal Schmidt modes determined by the Schmidt decomposition of the joint spectral amplitude of signal and idler photons.^62,78 The resulting state is produced by applying M independent two-mode squeezers to vacuum, each with thermal statistics of photon pairs. Therefore, the generating function for the number of photon pairs is given by⁷⁸ 

with $rieiθi=Cλi/(iℏ)$⁠. The constant C is determined by the properties of the nonlinear medium and the pump light and $λi$ is the i-th Schmidt coefficient. For an increasing number of equally strong squeezers, the PND changes from thermal to Poissonian statistics,⁷⁸ as shown in Fig. 2.

As a further example of our simulation method, we consider the bivariate photon statistics of SPDC while taking into account noise and detection efficiencies. First, we calculate the two-dimensional joint distributions p(n_A, n_B) for different parameters of the detection efficiencies and noise, as well as the cumulative distribution. The results obtained via our approach using the covariance matrix are shown in Fig. 3. It can be seen that, as expected, decreasing the detection efficiency and adding noise blurs the line of delta functions that are observed for the ideal state along the diagonal n_A = n_B.

Second, we use Eq. (9) to analyze mixed moments of the PND of a TMSV state in Fig. 4. For comparison, analytic values for the mean values and the first mixed moment are shown. However, if the photon statistics are more complicated, e.g., due to loss, simultaneous detection of multiple modes, or noise, analytic expressions for the moments may not be easy to find but can be easily calculated by applying AD to the moment-generating function Eq. (9). In Fig. 4, we show the influence of loss and noise on the moments.

To investigate the scaling of our method, we have measured the run time for the PND computation from Fig. 3(c) as a function of n₂ with n₁ = 0 and for different numbers of modes. As expected for a GBS computation method, we observed an approximately exponential increase in the run time with the photon number. On average, each additional photon prolonged the run time by a factor of 2.8 for 256 two-mode squeezers and 3.5 for four two-mode squeezers. In comparison, Hafnian-based algorithms achieve a scaling of $O(N32N/2)$ for N photons³³ and even faster methods have recently been demonstrated.³⁴ It is important to note that the run time of our simulation method also grows with the number of modes. However, computing p(n₁ = 0, n₂ = 6) for a 1024 × 1024 covariance matrix representing 256 two-mode squeezers still took less than 13 s.

For states and setups more complex than the examples considered above, analytical solutions can become infeasible and numerical simulations are required. In this section, we consider an application demonstrating that by employing our method, the photon statistics can be simulated in a relatively simple way, even for complex optical setups. We model a multi-user QKD system recently published based on the BBM92 protocol.³⁷ Although the system uses non-PNR detectors, the photon statistics are relevant as multi-photon pair emission is an important effect limiting the achievable quantum key rates. Detection efficiencies and noise are taken into account. Imperfect interference is modeled by using the mode mismatch model from Ref. 46. All these effects are readily included by using our simulation method.

QKD enables the distribution of secure keys between users based on information-theoretic principles from quantum physics.^2,79,80 The entanglement-based QKD system for four users named Alice (A), Bob (B), Charlie (C), and Diana (D)³⁷ implements the BBM92 time-bin protocol.^81–86 Figure 5(a) shows the setup for the users Alice and Bob. The setup for Charlie and Diana is identical, but with different values for parameters such as insertion losses and coupling ratios as well as different detector properties.

The central component of the QKD system is a source of time-bin entangled photon pairs. Laser pulses are sent through an imbalanced interferometer, resulting in double pulses with a time and phase difference given by the interferometer’s path length imbalance. The double pulses produce photon pairs by type-0 SPDC in a nonlinear crystal. The mean photon pair number per double pulse is μ ≈ 0.034. The spectral width of the pump is much narrower than the photon frequency spectrum, so that the photons are frequency-correlated; for a pump frequency 2ω₀, energy conservation requires that the frequencies of the signal and idler photon sum up to ω_s + ω_i = 2ω₀. The photons are separated by wavelength demultiplexing (WDM). The phase matching function is almost constant over several dozen nanometers, so that photons are produced in multiple pairs of 100 GHz-wide WDM channels. A pair of users can only obtain entangled photons if their wavelength channels are symmetrically arranged around the center frequency ω₀. Multiple user pairs can thereby exchange quantum keys simultaneously and independently, and the user combinations can be reconfigured by changing the WDM channel assignments.

In the receiver station of each user, the photons pass through an interferometer with a path length difference matched to that of the source’s interferometer. They are, therefore, detected in one of three different time bins, depending on the path combination of the laser pulse and the photons in the interferometers. The quantum bit values in the time basis are obtained from the arrival time of the photons. When both photons are detected in the central time bin, two-photon Franson interference⁸⁷ leads to emission from correlated interferometer outputs, and the bit values in the phase basis are obtained from the detector numbers 0 and 1.^37,82,83

However, various imperfections in the real setup lead to bit errors. For example, the interferometer delays may not be perfectly matched, the detectors have a non-zero dark count probability, and the beam splitters are not ideal 50/50 splitters. The ratio of the number of bit errors to the number of measured bits is called the quantum bit error rate (QBER). The higher the QBER, the lower the secure key rate that can be extracted. A simulation of the system can help to quantify the contribution of the different effects to the QBER and to optimize the setup. A higher value of μ, for example, yields a higher emission probability for photon pairs per laser pulse, resulting in a higher bit rate. However, due to the statistical nature of SPDC, increasing μ increases the probability of multi-pair emission and thereby leads to a higher QBER. Therefore, a simulation inevitably needs to take into account the effect of multi-pair emission.

For the simulation, we treat the three time bins as separate modes. The setup from Fig. 5(a) can thus be unfolded as depicted in Fig. 5(b). It is divided into four parts: state generation, propagation in fibers and interferometers, imperfect interference with mode-mismatch at the beam splitters A₂ and B₂, and detection in three time bins. The state after the source interferometer is modeled by initializing two TMSV states with covariance matrices Γ_TMSV(χ_S) and Γ_TMSV(χ_L), where the beam splitter asymmetries and losses in the arms of the source interferometer are absorbed into the squeezing parameters,

and

The coefficient r₀ is chosen such that both states together have a certain mean photon pair number μ. The propagation in the fibers and interferometers is represented by the state transformations for beam splitters, phase shifts, and losses, and the interference described by the mode-mismatch model from Ref. 46 is also modeled by beam splitters. Although the setup is complex, the construction of the final covariance matrix is straightforward: it is a combination of only four basic building blocks, namely, the covariance matrix of a TMSV state and the Gaussian state transformations for beam splitters, phase shifts, and losses. Combining the building blocks is simply achieved by multiplying the matrices representing them. The representations of the most important Gaussian states and their transformations can be found, e.g., in Refs. 4–7 and 46, so that constructing a setup’s model becomes a relatively simple task.

In Fig. 2, we have shown the effect of the number of two-mode squeezers on the photon statistics. To obtain results representing the correct photon statistics, we would have to replicate the setup from Fig. 5(b) in the simulation for the correct number of two-mode squeezers, each with its own individual mean photon pair number. The Schmidt number K can be used to estimate the number of two-mode squeezers based on the shape of the joint spectral amplitude of the SPDC process.^78,88,89 For the simulation, we roughly estimate that K is in the order of magnitude of 100, and we assume equally strong two-mode squeezers resulting in almost Poissonian photon pair statistics. As the beam splitters and losses are applied to all Schmidt modes, the covariance matrix is block-diagonal with K identical blocks. Instead of directly calculating the determinant of this large matrix, we apply the rule for block determinants and raise the determinant of one block to the power of K = 100.

In the experiment, we used single-photon avalanche detectors, as thoroughly analyzed in Ref. 90. Relevant parameters are the dead time τ_dead ≈ 10 µs, the efficiency η ≈ 20%, the after-pulse probability of a few percent, and the dark count rates in the range from 300 cps to more than 3000 cps depending on the detector. The flexibility of our simulation method allows us to include all these effects in the model as described in  Appendix C.

For the simulation, we have measured the values of all relevant experimental parameters, i.e., for the beam splitter transmissions, propagation losses in the fiber links, insertion losses of the components, and interferometer mode mismatches. Then, we used these values to determine the parameters of the simulation model in Fig. 5(b). For the detection efficiencies, dark count rates, and afterpulse probabilities of our detectors, we have included values obtained from tomographic measurements⁹⁰ in the detection operators as described in  Appendix C.

We have simulated the sifted key rates and quantum bit error rates by multiplying the source repetition rate with the detection probabilities for key bits and quantum bit errors in the time basis and phase basis. We compare the results with the measurements we have reported in Ref. 37 to check the validity of our simulation. For these experiments, the mean photon pair number had been set to μ = 0.034. For comparison, we consider the time basis QBER instead of the overall QBER because in the experiment the interferometer phases fluctuate and are readjusted by tuning the interferometer temperatures so as to minimize the phase basis QBER. Therefore, the phase basis QBER is not necessarily always at the minimum. Although it is possible to adjust the phases in the simulation to mimic phase fluctuations, we have chosen the phases for optimal interference and compare only the time basis QBER, which is unaffected by the phase fluctuations.

A comparison between the measured and simulated sifted key rates and QBERs in the time basis is shown in Fig. 6. The simulation represents both the dependence of the key rates and of the time basis QBER on the user combination and distance very well. The key rate decreases with increasing distance between the participants due to the increasing propagation losses, as expected. Since, for the distances involved and the chosen pulse durations, dispersion is not significant, the time QBER is almost independent of the transmission distance and transmission losses. However, the QBER variations can be attributed to the individual dark count rates and afterpulse probabilities of the detectors. Multiple effects lead to quantum bit errors in the time basis: uncorrelated clicks from two noise counts can lead to coincidence. Or, one photon from a pair may be lost, and the other photon can be detected in coincidence with a noise count or a photon from a different pair. Investigating the different contributions helps to determine the effects of limiting the secure key rate and optimizing the setup. Therefore, PNR simulations enable a detailed analysis of the different effects contributing to the QBER, as we demonstrate in Sec. IV C.

The process to derive the initial quantum state from the measurement outcome by applying Bayes’ theorem is known as quantum retrodiction^91–94 and has been applied, for example, to imperfect detectors with non-unity quantum efficiency and dark counts⁹⁵ as well as to beam splitters and amplifiers.⁹⁶ Bayes’ theorem states that the conditional probability p(E₁|E₂) for event E₁, given that event E₂ occurred, can be calculated from

The ability to simulate the photon number statistics enables such retrodictive analysis of the quantum state, i.e., it enables us to answer the question: Given that a certain measurement result was observed, what is the probability that it was produced by an initial quantum state with a specific photon number?

PNR simulations with GSs as presented earlier can be used to obtain p(E₁|E₂) when the probabilities p(E₁) and p(E₂) can be obtained via the simulation and p(E₂|E₁) is also accessible. Such retrodictive analysis enables the investigation of effects that may not be easily accessible by measurements. For example, the influence of multi-photon pair emission on the QBER can be compared to the contribution from noise counts. Instead of simulating the probability for multi-photon pair production, it can be calculated from the photon statistics of the pair source or it can be estimated by measuring the heralded g⁽²⁾ auto-correlation function of the source.⁹⁷ Not all events in which multiple photon pairs are produced lead to quantum bit errors, and conversely, not all quantum bit errors originate from multi-photon pair emission. Using PNR simulation, these effects can be included systematically. However, for example, the quantum bit errors originating from dark counts could also be estimated without setting up such a simulation.

When each of the two pump pulses produces a photon pair, Alice can detect a photon of the first pair in the early time bin and Bob can detect a photon of the second pair in the late time bin, causing a quantum bit error in the time basis. Here, we only want to briefly illustrate the method, and for simplicity, we neglect the dependency of the time bins in the same repetition cycle due to the dead time among each other. Furthermore, we only consider the contribution from up to one photon pair produced by each pump pulse as well as raw coincidences between the detectors A₀ and B₀ instead of all exclusive coincidences that would lead to error bits in the time basis.

Two photon pairs can only lead to an error in the time basis when one pair was produced by the first pump pulse happening with probability p(f) and the other one is produced by the second pulse occurring with probability p(s). Both values are obtained by simulating the probability of obtaining one photon pair directly after the SPDC process. From Bayes’ theorem in Eq. (33), the probability can be calculated that when a coincidence between A_0,E and B_0,L is measured, actually f and s pairs were produced by the first and second pulse, respectively,

The emission events are independent of each other, so that p(f, s) = p(f)p(s). The detections are independent as well, so that in Eq. (34) we can factorize $pA0,E,B0,L|f,s=p(A0,E|f)p(B0,L|s)$⁠. Here, p(A_0,E|f) is the probability that detector A₀ is triggered in the early time bin, given that f pairs have been produced in the first pump pulse. For f = 1, it can be simply obtained by tracing the path of the photon through the setup to calculate the total photon transmission probability $TA0,E$ to the detector, including the beam splitter transmissions, propagation losses, and the detector efficiency. The value for p(B_0,L|s) is similarly obtained,

For f = 0 and s = 0, the detection probability is given by the probability for a noise count, which is equivalent to setting T = 0 in Eqs. (35) and (36). Finally, the raw time error coincidence probability p(A_0,E, B_0,L) is also obtained from the simulation. Thereby, all probabilities can be calculated that are required to obtain p(f, s|A_0,E, B_0,L) from Eq. (34). In Fig. 7, we show p(f, s|A_0,E, B_0,L) as a function of the produced total mean number of photon pairs per repetition cycle μ and for different combinations of f and s.

For small values of μ below 10⁻³, the largest contribution to the time basis error coincidence probability p(A_0,E, B_0,L) comes from coincidences between dark counts. Coincidences involving one noise count and one photon pair are relevant for a wide range of μ-values from μ < 10⁻⁴ to μ > 1. The difference between the cases f, s = 1, 0 and f, s = 0, 1 is mainly due to different transmission losses from the fiber lengths of 26.9 km to Alice and 50.0 km to Bob and due to the much higher dark count probabilities for A₀ of 3045 cps compared to B₀ with 606 cps. The complementary probability in Fig. 7 shows that in the range around μ ≈ 0.1 and above, effects from more than one photon pair become relevant. In the range of μ = 0.034 where the QKD setup is operated, a majority of 53% of the A_0,E, B_0,L-coincidences occur when each pump pulse produces one photon pair. This underpins the fact that a simulation of the QKD system should take into account effects from multiple photon pairs to produce accurate results for the QBERs.

In the previous sections, we have shown the versatile applications of our expressions for the calculation of the photon statistics of GSs. To the best of our knowledge, we are the first to derive expressions for cumulative probabilities, moments, and factorial moments for general multi-mode GSs. For the PND p(n), alternative calculation methods can be found in the literature. One approach requires 2n derivatives to evaluate the photon number n. For that, Γ and d are expressed in the complex $(â,â†)$-basis instead of the real $(x̂,p̂)$-basis via $σQ=Ω†ΓΩ/2+1/2$ and τ = Ω^†d. In our notation, the expression reads^20,31

with $γT=(αT,α†)$ and $A=01N1N0(12N−σQ−1)$⁠.

The exponential to be differentiated has been recognized as the generating function for multivariate Hermite polynomials and, therefore, the PND can be expressed in terms of these polynomials.^98–100 However, these polynomials are given by complex recursion formulas, and the expressions are, therefore, considered difficult to evaluate for higher numbers of photons.^101–103 

For GSs with d = 0, Eq. (37) can be rewritten in terms of the Hafnian function as

where A_S is derived from A by repeating rows and columns, depending on the number of photons to be detected in a particular mode.^20,30,31,33 For states with d ≠ 0, a similar expression involves the loop Hafnian function.^32,33,35,47 For non-PNR detection, the probabilities involve the Torontonian and loop Torontonian function.^30,34 The software library The Walrus¹⁰⁴ provides software related to GBS and has implemented algorithms to evaluate these functions. In contrast to such highly specialized algorithms, we use the general-purpose tool of AD for our computations.

Generating operators for the photon statistics are treated in textbooks such as Refs. 40 and 105. In Refs. 78 and 88, Mauerer et al. have applied AD to the generating function from Eq. (30) to obtain the PND resulting from multiple two-mode squeezers. In Ref. 42, a probability-generating function for TMSV states in terms of the covariance matrix has been presented. However, generating functions in combination with AD, have rarely been applied for practical calculations in the context of photon statistics. A possible reason is that the evaluation of higher-order derivatives is in general a resource-intensive computational task so that analytical expressions are preferred. However, for the calculation of the photon statistics, this is not an a priori disadvantage, as the required computational resources for GBS computations scale exponentially even with optimized state-of-the-art algorithms.^33,34 It can, nevertheless, be expected that specially tailored algorithms yield performance benefits becoming relevant for higher photon numbers.

Often only the lowest few moments of probability distributions, such as order 3 or 4, are considered, which are easily evaluated by our method. For the PND, we have been able to calculate photon numbers up to 10–12 without any optimizations. The simulations took a few minutes and required multiple GBs of RAM. By tuning the parameters of PyTorch and using different AD software or more powerful hardware, even higher numbers of photons can be simulated. In the field of GBS complexity research, PNDs for considerably higher numbers of photons have been computed on much more powerful computers. For example, GBS probabilities have been computed on a workstation with 96 CPUs for 50 photons in Ref. 32 or for up to 92 photons on a 100 000-core supercomputer in Ref. 34.

However, quantum optics applications often consider states with low mean photon numbers, for example, when coincidences are considered. An example is entanglement-based QKD, where μ ≪ 1 is required to avoid quantum bit errors from coincidences between photons from different pairs. The QKD system we have modeled in Sec. IV is operated with a mean photon pair number per pulse of around μ = 0.034 (cf. Ref. 37). For such low mean photon numbers, only the first few photon numbers are of practical relevance as the probabilities of finding higher photon numbers decrease rapidly. Therefore, for many applications, a computational limitation to low photon numbers is not a significant restriction.

Another method to calculate the PND of GSs has recently been presented in Ref. 62, where one of the major results written in our notation is the expression

with Y = diag(w) ⊕ diag(w) and $ws=yj$⁠. The authors of Ref. 62 derive it by applying a fanning-out transformation to infinitely many virtual non-PNR detectors, motivated by the intuition that in this limiting case, the probability that any of these virtual detectors receives more than one photon vanishes. Therefore, this approach formalizes the experimental approach to realize PNR detection by multiplexing multiple non-PNR detectors mentioned in Sec. I. Equation (39) can be rewritten as a special case of Eq. (7), which generalizes the result from Ref. 62 in multiple ways: it already incorporates detection efficiencies and noise counts and it can be applied to states with non-zero displacement vector d. Besides Eq. (7), our derivation yields generating functions that enable the calculation of the moments, cumulative probabilities, density matrix elements, and factorial moments of the PND via Eqs. (8)–(11), (14), and (15). Furthermore, our expressions require, in contrast to Eq. (37), only one derivative per photon number facilitating the numerical evaluation. Meanwhile, analytic expressions for the PND of single-mode photon-added and photon-subtracted GSs have been reported in Refs. 54–56, we have presented the more general Eqs. (21) and (22) that cover the multi-mode case and which provide, besides the PND, also the cumulative probabilities, moments, and factorial moments.

Conceptually, our generating-function approach provides a new point of view on photon statistics that is different from the Hafnian approach closely related to graph theory, and also different from the approach using the fanning-out transformation. Therefore, it can be a starting point for further theoretical investigation. For example, the analytic generating functions could also be analyzed for the complexity of their computation. Another example is the analytical generating function we have presented for the single-mode displaced squeezed thermal state.

We have shown that the flexibility of our approach enables us to easily incorporate imperfections such as noise of arbitrary statistics, provided its generating functions are known. In Ref. 62, it is noted that in comparison to other formulas for PNR detection, such as Eq. (37) or Eq. (38), an advantage of Eq. (39) is that modes entering the same detector are treated by the same derivative, whereas the alternative formulas treat all modes separately, so that terms for all the possible distribution patterns of the photons over the modes in a detector have to be calculated separately. In Ref. 62, this fact has been used to incorporate the spectral modes into calculations of the Hong-Ou-Mandel visibility of SPDC sources. The same argument holds for our Eqs. (7)–(11), which involve one differentiation variable per detector, independent of the number of modes entering it. This is because our formulas obey the usual rules for generating functions for probabilities and moments and, therefore, the convolution of probability distributions is simply represented by the multiplication of the generating functions.

We have presented various generating functions for the photon statistics of multi-mode GSs in terms of the covariance matrix and displacement vector, from which the probabilities, moments, and factorial moments are retrieved by repeated AD. Common effects in experiments, such as noise, non-unity detection efficiency, and the joint detection of multiple modes, can be easily taken into account. Furthermore, expressions for the generating functions of multi-mode photon-added and photon-subtracted GSs are derived, relating them to the generating functions of the underlying GSs, thus greatly extending the range of their possible applications.

We have implemented AD by using the software framework PyTorch to effortlessly retrieve probabilities for photon numbers up to values of 12 on a desktop computer without any optimization for performance. The code for a simple example demonstrating the calculation of a PND using PyTorch can be found in our technical report (Ref. 36).

As an application, we have modeled an entanglement-based system for QKD that we have recently developed, and we have observed very good agreement between the simulated and measured sifted key rates and quantum bit error rates. Finally, we have used our PNR simulation of the setup to analyze the fraction of quantum bit errors due to the production of up to two photon pairs in the same repetition cycle.

Our expressions for the generating functions in Eqs. (7)–(11) are formulated in terms of simple operations on the covariance matrix and displacement vector of GSs. Together with AD for the differentiation of the generating functions, this enables practical, easy-to-implement, and versatile simulations of quantum-optical setups.

This research has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), under Grant No. SFB 1119–236615297. We thank Philipp Kleinpaß for fruitful discussions about the simulation of the QKD system. We also thank Julian K. Nauth, Vladimir M. Stojanović, and Gernot Alber for their valuable comments.

The authors have no conflicts to disclose.

Erik Fitzke: Conceptualization (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Validation (lead); Visualization (supporting); Writing – original draft (lead); Writing – review & editing (equal). Florian Niederschuh: Conceptualization (supporting); Methodology (supporting); Software (lead); Validation (supporting); Visualization (lead). Thomas Walther: Funding acquisition (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon request.

In this section, we briefly summarize the basic properties of general multi-mode GSs and of generating functions for probability distributions.

A photonic state with S modes can be described by creation and annihilation operators $âs†$ and $âs$ with s = 1…S. The quadrature operators $x̂s=(âs+âs†)/2$ and $p̂s=(âs−âs†)/i2$ can be expressed in matrix notation as $q̂=Ωâ$ with

with the basis-changing matrix $Ω=1211−i1i1$ and with $1$ denoting the S × S identity matrix.

The characteristic function of an operator $Ô$⁠,

has 2S arguments that can be collected into the vector $ξT=(ξxT,ξpT)=ξx1…ξxSξp1…ξpS$⁠. The expectation value $⟨Ô⟩$ w.r.t. a quantum state $ϱ̂$ can be written as^4,46

GSs are the states described by a Gaussian characteristic function^4–6,46

with a real, symmetric 2S × 2S covariance matrix Γ and a displacement vector d defined by¹¹³ 

Hamiltonians linear or quadratic in $â$ and $â†$ map GSs to GSs⁶ and can be expressed in matrix notation as

and