General model of photon-pair detection with an image sensor
Hugo Defienne Matthew Reichert Jason W. Fleischer
quant-ph physics.optics


APS/123-QED




defienne@princeton.edu



Department of Electrical Engineering, Princeton University, Princeton, NJ 08544, USA



             




We develop an analytic model that relates intensity correlation measurements performed by an image sensor to the properties of photon pairs illuminating it. Experiments using both an effective single-photon counting (SPC) camera and a linear electron-multiplying charge-coupled device (EMCCD) camera confirm the model. 












Valid PACS appear here
                             

                              
General model of photon-pair detection with an image sensor
    Jason W. Fleischer
    December 30, 2023
===========================================================















Because it may exhibit quantum features at room temperature, light is one of the most promising platforms to investigate quantum mechanics and its applications in quantum computing, communication, and imaging <cit.>. Pairs of photons represent the simplest system showing genuine quantum entanglement in all their degrees of freedom: spatial, spectral, and polarization <cit.>. Demonstrations range from fundamental tests of Bell's inequality with polarization entangled photons <cit.> to the development of new imaging techniques <cit.>. Spatial entanglement between photons is particularly attractive, since its natural high-dimensional structure <cit.> holds promise for powerful information processing algorithms <cit.> and secure cryptographic protocols <cit.>. While generating photon pairs entangled over a large number of spatial positions is now commonly achieved using spontaneous parametric down-conversion (SPDC) <cit.>, full characterization of entangled photon states in high-dimensional Hilbert spaces remains a challenging task. Indeed, the process requires intensity correlation measurements between all pairs of possible positions, and its efficiency strongly depends on the properties of the detection system.

Light intensity correlation is a type of optical measurement used in imaging techniques, such as scintigraphy <cit.> and ghost imaging <cit.>, and in some characterization procedures, such as dynamic light scattering <cit.> and fluorescence correlation spectroscopy <cit.>. In quantum optics, intensity correlation measurements are used to measure coincidences between correlated photons. The detection apparatus generally involves single-photon sensitive devices connected to an electronic coincidence counting circuit. The number of measurements required scales with both the number of correlated photons and the number of optical modes. Typically, correlation measurements of spatially entangled photon pairs are performed with two avalanche photo-diodes (APDs) that are raster-scanned over the different positions. Since pairs generated by a conventional SPDC source may be entangled over a very large number of spatial modes <cit.>, this raster scanning technique is prohibitively time consuming and cannot be used in practice. 

Both electron multiplying (EM) <cit.> and intensified CCD <cit.> cameras have been used to perform high-dimensional measurements. A threshold applied on the measured images allows them to operate effectively as multi-pixel single-photon counter <cit.>. Recent works have revealed some features of entanglement between pairs of photons generated by SPDC <cit.>, but these techniques have not retrieved the full characteristics the photon pairs, i.e., their full joint probability distribution. Moreover, the theoretical analyses associated with these works were carried out under approximations on the form of the correlation <cit.>. In particular, these works assumed a regime of detection in which the two photons never hit the same pixel of the camera, which is counter-intuitive when measuring pairs that are strongly correlated in position. 

In this work, we provide a general theoretical framework for intensity correlation measurements of entangled photon pairs performed with any type of detection system, with no approximation made on the source. We then compare our model to experiments performed with two different detection systems: 1) an APD-like single-photon counter (SPC) camera, implemented using an EMCCD camera with thresholding <cit.>, and b) a linear EMCCD camera with no threshold. Surprisingly, we show that the joint probability distribution of photon pairs can be measured using an EMCCD camera without thresholding, which provides one of the simplest techniques to characterize high-dimensional spatial entanglement of photon pairs.         
 

Figure <ref> shows the general detection scheme considered in our model. It is studied with respect to two assumptions : 


  * Pixels of the image sensor operate independently

  * The input state is a pure two-photon state

In the input, denoted |ϕ⟩, both photons have the same polarization and frequency spectrum. Its associated two-photon wavefunction depends only on the spatial properties of the pairs <cit.> and can be expressed as

    |ϕ⟩ = ∑_i,j ∈ 1,N ϕ_ij|i,j⟩

where |i,j⟩  is a non-symmetric state defining a configuration in which the first photon of a pair is located at pixel i and the second at pixel j, and ϕ_ij is the spatially dependent two-photon wavefunction discretized over the pixels of the sensor. The joint-probability distribution Γ_ij=|ϕ_ij|^2 represents the probability of the first photon of the pair arriving at pixel i and the second at pixel j. 

Each photon falling on the camera has a probability η to be transformed into a photoelectron. In addition, electrons can also be generated from thermal fluctuations (dark noise). As shown in Figure <ref>, input electrons go through a potential amplification and readout process that converts them into detectable signals. The exact operation performed depends only on the internal characteristics of the image sensor, including the specific sensor technology and its noise properties. This process is then fully characterized by a set of conditional probability functions { P(x_i|k_i) }_i ∈ 1,N, in which k_i ∈ℕ is the number of electrons present at pixel i after the screen and x_i the corresponding output value returned by the sensor. Henceforth, we refer to { P(x_i|k_i) }_i ∈ 1,N as the detector response function. 

Two types of images are returned at the output at the l^th acquisition: a direct image, denoted { x_i^(l)}_i ∈ 1,N, composed of output values returned at each pixel, and a correlation image, denoted { x_i^(l) x_j^(l)}_i,j ∈ 1,N, computed by the tensor product of each direct image with itself. When a large number of images M is recorded, averaging over all of them enables estimation of the mean values: 

    ⟨ x_i ⟩   = lim_M →∞1/M∑_l=0^M x_i^(l)
    ⟨ x_i x_j ⟩   = lim_M →∞1/M∑_l=0^M x_i^(l) x_j^(l)

Assuming stationary illumination, ⟨ x_i ⟩ and ⟨ x_i x_j ⟩ can be written in term of their corresponding probability distributions:

    ⟨ x_i ⟩    = ∑_x_i = 0^+ ∞ x_i   P(x_j) 
    ⟨ x_i x_j ⟩    = ∑_x_i = 0^+ ∞∑_x_j = 0^+ ∞ x_i x_j   P(x_i,x_j)

where P(x_i) represents the probability for the sensor to return value x_i at pixel i and P(x_i,x_j) is the joint probability to return values x_i at pixel i and x_j at pixel j, during the acquisition of each frame. Using Bayes' theorem, ⟨ x_i ⟩ can be expressed as:

    ⟨ x_i ⟩ =  ∑_m=0^+ ∞ P(m) ∑_k_i=0^2m I_k_i  P(k_i|m)

where P(m) is the probability for m ∈ℕ pairs to fall on the screen during the exposure time and P(k_i|m) is the conditional probability of generating k_i photoelectrons at pixel i given m pairs. I_k_i is the mean of the detector response function at pixel i, defined as:

    I_k_i = ∑_x_i = 0^+ ∞ x_i   P(x_i|k_i)

P(k_i|m) can then be written in terms of the marginal probability Γ_i = ∑_iΓ_ij and the probability of measuring both photons of a pair at the same pixel Γ_ii. The complete calculation is detailed in Appendix <ref>, and leads to the following general expression for ⟨ x_i ⟩: 


    ⟨ x_i ⟩ = ∑_m=0^+ ∞  P(m) ∑_k_i=0^2 m I_k_i∑_q=0^⌊ k_i/2 ⌋( η^2 Γ_ii)^q ( 2 η Γ_i-2 η^2 Γ_ii)^k_i-2q(1-2 η Γ_i + η^2 Γ_ii)^m-k_i+qk_i-qq mk_i-q


where nk = n!/k!(n-k)!H(n-k), H is the Heaviside (unit step) function. Using a similar approach, the correlation image ⟨ x_i x_j ⟩ for the system to return a value x_i at pixel i and x_j at pixel j ≠ i is written as:

    ⟨ x_i x_j ⟩=∑_m=0^+ ∞ P(m) ∑ _k_i=0^2m∑_k_j=0^2m I_k_iI_k_j P(k_i,k_j |m)

where P(k_i,k_j|m) is the conditional probability of generating k_i and k_j photoelectrons at pixels i and j ≠ i, respectively, given m photon pairs in each frame. Assuming that η is uniform over the screen, the correlation coefficient ⟨ x_i x_j ⟩ can be related to the Γ_ij. The full calculation, given in appendix <ref>, gives:


    ⟨ x_i x_j ⟩   = ∑_m=0^+ ∞ P(m) ∑_k_i=0^2m∑_k_j=0^2m I_k_i I_k_j∑_q=0^⌊ (k_i+k_j)/2 ⌋∑_l=0^q∑_p=0^q-l (1-2 ηΓ_i-2 ηΓ_j +η^2 Γ_ii + η^2 Γ_jj + 2 η^2 Γ_ij)^m-(k_i+k_j-q)
       ( η^2 Γ_jj)^p ( 2 η^2 Γ_ij)^l( η^2 Γ_ii)^q-l-p (2 ηΓ_i - 2 η^2 Γ_ii- 2 η^2 Γ_ij)^k_i+l-2 (q-p)  (2 ηΓ_j - 2 η^2 Γ_jj- 2 η^2 Γ_ij)^k_j-2 p-l
       k_j-l-ppk_i-q+pq-l-pk_i+k_j-q-lk_i-q+pk_i+k_j-qlmk_i+k_j-q



Equations <ref> and <ref> show that knowing the characteristics of the image sensor, namely its quantum efficiency and response function, as well as the number distribution of incident pairs P(m), relates the direct images, ⟨ x_i ⟩ and ⟨ x_j ⟩, and the correlation image, ⟨ x_i x_j ⟩, to the joint probability distribution of the pairs Γ_ij. Note that these results hold only for i ≠ j ; the case i=j is more subtle and is treated separately in appendix <ref>. 

This set of equations provides a general link between measurements performed by any detector and the joint probability distribution of photon pairs illuminating it. We demonstrate the validity of our model by applying it to the case of an SPC camera, mimicked using a thresholded EMCCD camera <cit.>, and to the case of an EMCCD camera operated without threshold.

SPC cameras generally consist of an array of single-photon avalanche diodes (SPADs) or APDs with all electronics incorporated into each pixel. Photon-to-electron conversion is performed at a given quantum efficiency η, and the detector ideally returns a non-null current (value 1) at the output if at least one electron was present at the input of the amplifier and no current if not (value 0). As shown in appendix <ref>, assuming a Poissonian distribution <cit.> for P(m) in this model simplifies Equations <ref> and <ref>, allowing expression of Γ_ij in terms of ⟨ c_i ⟩ and ⟨ c_i c_j ⟩:

    Γ_ij= 1/2 η^2 m̅ln [ 1+⟨ c_i c_j ⟩ - ⟨ c_i ⟩⟨ c_j ⟩/(1-⟨ c_i ⟩)(1-⟨ c_j ⟩) ]

where the general output variable x has been replaced by a binary variable c (counts) that takes only two possible values c∈{0,1}. The mean photon-pair rate m̅, which can be controlled by adjusting the exposure time of the sensor or the power of the pump laser, and the quantum efficiency η act only as scaling factors. 

Experimental results are shown in Figure <ref>. Pairs are generated by type-I SPDC in a β-barium borate (BBO) crystal pumped by a continuous-wave (CW) laser centered at 403 nm, and near-degenerate down-conversion is selected via spectral filters at 806 ± 3 nm. The far field of the output of the BBO crystal is projected onto the screen of an EMCCD camera. As detailed in <cit.>, applying a threshold on each acquired image effectively enables the EMCCD to operate as an SPC camera with quantum efficiency η_eff≈ 0.44 and a noise probability p_10≈ 0.015 (appendix <ref>.c). To facilitate the analysis, the 4-dimensional space of pair positions is reduced to 2 dimensions by fixing the Y-coordinates Y_1=33 and Y_2=45 (arbitrarily) and measuring only ⟨ c_i ⟩ and ⟨ c_i c_j ⟩ along two X-axis pixels of the camera, denoted {X_1} and {X_2} (Figure <ref>.b). Figure <ref>.c shows the measured joint probability distribution Γ_X_1X_2 together with the marginals Γ_X_1 and Γ_X_2 (taken after background subtraction and normalization, appendix <ref>). The intense anti-diagonal reveals an anti-correlated behavior of the pairs, as expected when measuring photons in the far field of the crystal. As shown in Figure <ref>.d, the measured joint probability distribution is well fit by a double-Gaussian model <cit.> Γ^th_X_1X_2 of parameters σ_- = 926.1  μ m and σ_+ =12.1  μ m (appendix <ref>). Selected profiles Γ_X_1|X_2=65 and Γ^th_X_1|X_2=65, highlighted in Figure <ref>.e, show a good match between the experiment and the double-Gaussian fit, confirming the validity of our model.


It is commonly thought that photon counting is necessary to compute the joint probability distribution of pairs of photons. We now demonstrate the surprising result that simple operation of a camera without thresholding also enables measurement of Γ_ij. In this case, the readout process becomes more complex, but an analytic form of P(x|k) can be calculated quantitatively if the sources of noise are known, e.g. those provided in <cit.>. For EMCCD cameras, the mean of the detector response function I_k depends linearly on the number of electrons k at the input, I_k = A   k + x_0, where the amplification parameter A depends on the mean gain and analog-to-digital conversion and x_0 is a constant background (Appendix <ref>). As shown in Appendix <ref>, this response allows expression of Γ_ij as

    Γ_ij = 1/2 A^2m̅η^2[ ⟨ x_i x_j ⟩-⟨ x_i ⟩⟨ x_j ⟩]

where the parameters A,  η and m̅ contribute only a scaling factor, which may be determined by normalization. 

Figure <ref>.a shows the quantity R_X_1X_2≡⟨ x_X_1 x_X_2⟩-⟨ x_X_1⟩⟨ x_X_2⟩ measured by performing an experiment in the same conditions as that for Figure <ref>.a but without thresholding the output images. These results compare favorably with those of Figure <ref>.b. Profiles R_X_1|X_2=65 (blue) and Γ^th_X_1|X_2=65 (red) shown in Figure <ref>.c highlight the very good agreement between the double-Gaussian fit and the measurement without threshold.


The physical interpretation of Equation <ref> can be seen by expanding the expression R_ij over a finite number of images M≫1 using Equations <ref> and <ref>:

    R_ij≈1/M∑_l=0^M x_i^(l) x_j^(l) - 1/M^2∑_l,l',l ≠ l'^M x_i^(l) x_j^(l')

The first term is the average tensor product of each frame with itself. Intensity correlations in this term originate from detections of both real coincidences (two photons from the same entangled pair) and accidental coincidences (two photons from two different entangled pairs). Since there is zero probability for two photons from the same entangled pair to be detected in two different images, intensity correlations in the second term originate only from photons from different entangled pairs (accidental coincidence). A subtraction between these two terms leaves only genuine coincidences, which is proportional to the joint probability distribution Γ_ij. 

These results show that measuring correlation between pairs of photons is not a task exclusive to single-photon sensitive devices such as SPC cameras, SPADs, or APD arrays but can be achieved using any type of image sensor. Using a megapixel image sensor as a highly parallel intensity correlator offers much promise for measuring high-dimensional entangled states, necessary for quantum computing, communication, and imaging. Moreover, the model can be extended readily to states containing more than two entangled photons, in order to study higher degrees or new forms of entanglement.





















§  DERIVATION OF THE GENERAL EXPRESSION FOR ⟨ X_I ⟩ (EQUATIONS <REF>)


This section provides a step-by-step derivation of the general formula of ⟨ x_i ⟩ written in equation <ref>. Starting from the definition of ⟨ x_i ⟩ (equation <ref>) and introducing Bayes' formula gives:


    ⟨ x_i ⟩   = ∑_x_i=0^+ ∞ x_i P(x_i) 
       = ∑_x_i=0^+ ∞ x_i  ∑_m=0^+ ∞ P(m) ∑_k_i=0^2 m P(x_i|k_i) P(k_i|m) 
       = ∑_m=0^+ ∞ P(m) ∑_k_i=0^2 m  I_k_iP(k_i|m)

where P(k_i|m) is the probability of generating k_i photoelectrons at pixel i given that a total of m photon pairs reaches the sensor during the acquisition time. 

Derivation of Equation <ref> relies on expressing P(k_i|m) in function of the Γ_ij. This quantity is first expanded again using Bayes' formula as: 

    P(k_i|m) = ∑_n_i=0^+ ∞ P(k_i|n_i) P(n_i|m)

where P(n_i|m) represents the probability for n_i photons to fall on pixel i given m photon pairs; it is calculated first in section <ref>. P(k_i|m) is then calculated in section <ref>. 



 §.§  Expression of P(n_i|m)


We illustrate our reasoning by first calculating P(n_i|m) for the case m=1,i.e. the case where only a single pair falls on the screen during the exposure time. Three possible events may occur:



  * Both photons reach pixel i: P(2|1)=Γ_ii

  * No photons reach pixel i: P(0|1)=1-2 Γ_i + Γ_ii

  * One photon of the pair reaches pixel i and the other does not: P(1|1)=2   (Γ_i- Γ_ii)


For m>1, we can simply consider each pair following one of these three possibilities. The problem is then reduced to counting all the possible configurations. As shown on Figure <ref>, we can label the number of pairs of each possibility:


  * m_2 is the number of pairs where both photons reach i (green dashed circle)

  * m_1 is the number of pairs where only one photon reaches i (blue dashed circle)

  * m_0 is the number of pairs that are lost (black dashed circle)


As a set, the three numbers {m_0,m_1,m_2} describe all possible arrangements of the pairs at the input. The probability for the pairs to be in a specific configuration is given by:

    P_m_0,m_1,m_2 = P(0|1)^m_0 P(1|1)^m_1 P(2|1)^m_2
       m_1+m_2m_1m_0+m_1+m_2m_1+m_2

To express the probability of detecting n_i photons at pixel i, we need to consider two conservation equations linking  {m_0,m_1,m_2} to {n_i,m}:

    m   = m_0+m_1+m_2
    
    n_i   =2 m_2+m_1

With these constraints, we can write the variables {m_0,m_1} in terms of two fixed parameters {n_i,m} and one free parameter m_2:

    m_1   = n_i-2m_2 
    
    m_0   =m-n_i+m_2

Since numbers are positive, m_2 ∈ 0, ⌊ n_i/2⌋, where ⌊ n_i/2⌋ is the integer part of  n_i/2. Finally, the conditional probability P(n_i|m) of detecting n_i photons at pixel i is obtained by summing the probabilities associated with all the possible arrangements:

    P(   n_i|m)= ∑_m_2=0^⌊ n_i/2⌋Γ_ii^m_2 (2Γ_i-2Γ_ii)^n_i-2 m_2
       (1-2Γ_i+Γ_ii)^m-n_i+m_2n_i-m_2m_2mn_i-m_2



 §.§  Expression of P(k_i|m) 

Equation <ref> links P(k_i|m) to P(n_i|m) and P(k_i|n_i). Considering the finite quantum efficiency η of the photo-sensitive screen, P(k_i|n_i) can be written as:


    P(k_i|n_i)= n_ik_i η^k_i  (1- η)^n_i-k_i

Combining this expression with Equation <ref>, we find the probability for obtaining k_i electrons for m incident pairs:

    P(   k_i|m) = ∑_n_i=0^2m P(k_i|n_i)P(n_i|m) 
       = ∑_n_i=0^2m∑_m_2=0^⌊ n_i/2 ⌋n_ik_i η^k_i  (1- η)^n_i-k_iΓ_ii^m_2  (2 Γ_i-2 Γ_ii)^n_i-2m_2 
        (1-2 Γ_i + Γ_ii)^m-n_i+m_2n_i-m_2m_2 mn_i-m_2
       = ∑_q=0^⌊ k_i/2 ⌋( η^2 Γ_ii)^q ( 2 ηΓ_i-2 η^2 Γ_ii)^k_i-2q 
       ( 1-2 ηΓ_i + η^2 Γ_ii)^m-k_i+qn_i-qq mn_i-q

where the transition from line <ref> to line <ref> follows mathematical induction on m. Introducing this expression in equation <ref> gives the final expression of ⟨ x_i ⟩ in Equation <ref>.









§  DERIVATION OF THE GENERAL EXPRESSION OF ⟨ X_I X_J ⟩ (EQUATIONS <REF>)


This section provides a derivation of ⟨ x_i x_j ⟩ in Equation <ref>. Assuming all pixels have the same properties and are independent, only two different pixels i and j need to be considered. Here, we consider the case i ≠ j ; the case i=j is treated separately in appendix <ref>. Starting from the difinition of ⟨ x_i x_j ⟩ (equation <ref>) and introducing Bayes'formula gives:

    ⟨ x_i x_j ⟩   = ∑_x_i=0^+ ∞∑_x_j=0^+ ∞ x_i x_j P(x_i,x_j) 
       =∑_m=0^+ ∞  P(m) ∑ _k_i=0^2m∑_k_j=0^2m I_k_i I_k_j P(k_i,k_j |m)

where P(k_i,k_j |m) is the probability of generating k_i photoelectrons at pixel i and k_j photoelectrons at pixel j given that a total of m photon reach the sensor during the acquisition time.
Derivation of equation <ref> relies on expressing P(k_i,k_j|m) as a function of Γ_ij. It can be expanded using Bayes' theorem as:

    P(k_i,k_j|m) = ∑_n_i=0^+ ∞∑_n_j=0^+ ∞ P(n_i,n_j|m) P(k_i|n_i) P(k_j|n_j)

where: 


  * P(n_i,n_j|m) is the probability for n_i and n_j photons to fall on pixels i and j, given m pairs. 

  * P(k_i|n_i) P(k_j|n_j) is the joint probability of generating k_i and k_j photoelectrons from n_i and n_j photons. Its factorized form relies on the assumption that pixels of the camera operate independently (no cross-talk).

An analytic form of P(n_i,n_j|m) is first established in section <ref>. Equation <ref> is then simplified to obtain an expression of P(k_i,k_j|n_i) in section <ref>.



 §.§  Expression of P(n_i,n_j|m)


As before, we first calculate P(n_i,n_j|m) for the simple case m=1 and then generalize it to m. For m=1, only one pair falls on the screen during the exposure time. There are six possibilities with the following probabilities:



  * Both photons reach pixel i: P(2,0|1)=Γ_ii

  * Both photons reach pixel j: P(0,2|1)=Γ_jj

  * One photon reaches pixel i and the other reaches pixel j: P(1,1|1)=2Γ_ij

  * Only one photon reaches pixel i, and the other does not reach i or j : P(1,0|1)=2 Γ_i - 2 Γ_ii- 2 Γ_ij

  * Only one photon reaches pixel j and the other does not reach i or j : P(0,1|1)=2 Γ_j - 2 Γ_jj- 2 Γ_ij 


  * No photons reach either pixel i or pixel j: P(0,0|1)=1-2 Γ_i-2 Γ_j +Γ_ii +Γ_jj + 2 Γ_ij. 


For m>1, we can treat each pair as following one of these possibilities. As shown in Figure <ref>, the number of pairs for each possibility can be labelled as follow:


  * m_20 is the number of pairs where both photons are detected at i (upper green circle)

  * m_02 is the number of pairs where both photons are detected at j (lower green circle)

  * m_11 is the number of pairs where one photon is at i and its pair at j (pink circle)

  * m_10 is the number of pairs where one photon is at i and its pair is not at i or j  (upper blue circle)

  * m_01 is the number of pairs where one photon is at j and its pair is not at i or j (lower blue circle) 


  * m_00 is the number of pairs where both photons are lost (black circle)



The six numbers {m_20,m_02,m_11,m_10,m_01,m_00} describe all possible configurations of the pairs at the input of the image sensor. The probability for the pairs to be in a specific configuration is given by:

    P_m_20,m_02,m_11,m_10,m_01,m_00 = P(2,0|1)^m_20 P(0,2|1)^m_02
        P(1,1|1)^m_11 P(1,0|1)^m_10 P(0,1|1)^m_01 P(0,0|1)^m_00
       m_10+m_20m_10m_01+m_02m_01
       m_10+m_01+m_20+m_02m_10+m_01
       m_10+m_01+m_20+m_02+m_11m_11
       m_10+m_01+m_20+m_02+m_11+m_00m_00

Similar to the single-pixel case, to express the probability of n_i photons arriving at pixel i and n_j photons at pixel j given a total number of pairs m, we need to consider three conservation equations linking {m_00,m_10,m_01,m_20,m_02,m_11} to {n_i,n_j,m}, giving the constraints:

    n_i   =2 m_20+m_10+m_11
    
    n_j   =2 m_02+m_01+m_11
    
    m   =m_20+m_02+m_11+m_10+m_01+m_00

Besides the three fixed parameters {n_i,n_j,m}, we also introduce one parameter q defined as:

    q=m_11+m_20+m_02

q represents the number of pairs in which both photons reach at least one of the two pixels. The numbers of pairs can then be expressed as written at:

    m_20   = q-m_02-m_11
    
    m_01   = n_j-2 m_02-m_11
    
    m_10   = n_i+m_11-2 (q-m_02) 
    
    m_00   = m-n_i-n_j+q

Since these numbers are positive, {m_11,m_02,q} take only the following values :

    m_11   ∈ 0, q 
    
    m_02   ∈ 0, q-m_11
    
    q    ∈ 0, ⌊ (n_i+n_j)/2⌋

The conditional probability P(n_i,n_j|m) is obtained by summing the probabilities associated with all the possible configurations:


    P(n_i,n_j|m)    = ∑_q=0^⌊ (n_i+n_j)/2⌋∑_m_11=0^q∑_m_02=0^q-m_11 (1-2 Γ_i-2 Γ_j +Γ_ii +Γ_jj + 2 Γ_ij)^m-(n_i+n_j-q)
        (2 Γ_i - 2 Γ_ii- 2 Γ_ij)^n_i+m_11-2 (q-m_02)  (2 Γ_j - 2 Γ_jj- 2 Γ_ij)^n_j-2 m_02-m_11 Γ_ii^q-m_02-m_11 Γ_jj^m_02 ( 2 Γ_ij)^m_11
       n_i-q+m_02q-m_02-m_11n_j-m_11-m_02m_02n_i+n_j-q-m_11n_i-q+m_02n_i+n_j-qm_11mn_i+n_j-q





 §.§ Expression of P(k_i,k_j|m)



The probability P(k_i,k_j|m) of generating k_i photo-electrons at pixel i and k_j photo-electrons at pixel j given a total of m pairs is expanded starting from equation <ref>:



    P(k_i,k_j|m) = ∑_n_i=0^2m∑_n_j=0^2m P(n_i,n_j|m) P(k_i|n_i) P(k_j|n_j) 
       = ∑_n_i=0^2 m∑_n_j=0^2 m∑_q=0^⌊ (n_i+n_j)/2⌋∑_m_11=0^q∑_m_02=0^q-m_11n_ik_i η^k_i  (1- η)^n_i-k_in_jk_j η^k_j  (1- η)^n_j-k_j  (1-2 Γ_i-2 Γ_j +Γ_ii +Γ_jj + 2 Γ_ij)^m-(n_i+n_j-q)
        (2 Γ_i - 2 Γ_ii- 2 Γ_ij)^n_i+m_11-2 (q-m_02)  (2 Γ_j - 2 Γ_jj- 2 Γ_ij)^n_j-2 m_02-m_11 Γ_ii^q-m_02-m_11 Γ_jj^m_02 ( 2 Γ_ij)^m_11
       n_i-q+m_02q-m_02-m_11n_j-m_11-m_02m_02n_i+n_j-q-m_11n_i-q+m_02n_i+n_j-qm_11mn_i+n_j-q
       = ∑_q=0^⌊ (k_i+k_j)/2 ⌋∑_l=0^q∑_p=0^q-l  (1-2 ηΓ_i-2 ηΓ_j +η^2 Γ_ii + η^2 Γ_jj + 2 η^2 Γ_ij)^m-(k_i+k_j-q)
        (2 ηΓ_i - 2 η^2 Γ_ii- 2 η^2 Γ_ij)^k_i+l-2 (q-p)  (2 ηΓ_j - 2 η^2 Γ_jj- 2 η^2 Γ_ij)^k_j-2 p-l ( η^2 Γ_ii)^q-p-l
       ( η^2 Γ_jj)^p ( 2 η^2 Γ_ij)^lk_i-q+pq-p-lk_j-l-ppk_i+k_j-q-lk_i-q+pk_i+k_j-qlmk_i+k_j-q


where the transition from line <ref> to line <ref> follows from mathematical induction on the variable m. Introducing the expression of P(k_i,k_j|m) in equation <ref> provides the complete expression of ⟨ x_i x_j ⟩ written in Equation <ref>.












§   DERIVATION OF FORMULAS LINKING Γ_IJ TO ⟨ X_I ⟩ AND ⟨ X_I X_J ⟩ IN THE CASE OF AN SPC CAMERA (EQUATIONS <REF>)


The simple readout process performed in SPC camera, together with an assumption on the pair number distribution, allows simplification of Equations <ref>. This operating mode is modeled by substituting the general output variable x of our model by a binary variable c  (counts) that takes only two possible values c∈{0,1}. The corresponding conditional probability functions P(c|k) are shown in Table <ref>.





 §.§   Simplification of sum over k_j and k_i

Using the model of readout process performed by SPC cameras (Table <ref>), I_k takes the following form:

    I_0    = p_10
    
    I_k_i   = 1  k_i > 0

Using this model, Equation <ref> simplifies as:

    ⟨ c_i ⟩ = ∑_m=0^+ ∞  P(m) ∑_k_i=0^2 m I_k_i∑_q=0^⌊ k_i/2 ⌋ ( 2 η Γ_i-2 η^2 Γ_ii)^k_i-2q
       ( η^2 Γ_ii)^q(1-2 η Γ_i + η^2 Γ_ii)^m-k_i+qk_i-qq mk_i-q
       = 1-(1-p_10) ∑_m=0^+ ∞ P(m) (1-2 η Γ_i + η^2 Γ_ii)^m

and Equation <ref> simplifies as:


    ⟨ c_i c_j ⟩   = ∑_m=0^+ ∞ P(m) ∑_k_i=0^2m∑_k_j=0^2m I_k_i I_k_j∑_q=0^⌊ (k_i+k_j)/2 ⌋∑_l=0^q∑_p=0^q-l  (1-2 ηΓ_i-2 ηΓ_j +η^2 Γ_ii + η^2 Γ_jj + 2 η^2 Γ_ij)^m-(k_i+k_j-q)
        (2 ηΓ_i - 2 η^2 Γ_ii- 2 η^2 Γ_ij)^k_i+l-2 (q-p)  (2 ηΓ_j - 2 η^2 Γ_jj- 2 η^2 Γ_ij)^k_j-2 p-l ( η^2 Γ_ii)^q-p-l
       ( η^2 Γ_jj)^p ( 2 η^2 Γ_ij)^lk_i-q+pq-p-lk_j-l-ppk_i+k_j-q-lk_i-q+pk_i+k_j-qlmk_i+k_j-q
       =1+(p_10^2-1-2 p_10) ∑_m=0^+ ∞ P(m)  (1-2 ηΓ_i-2 ηΓ_j +η^2 Γ_ii + η^2 Γ_jj + 2 η^2 Γ_ij)^m 
       +p_10∑_m=0^+ ∞ P(m) [ (1 - η (2 Γ_i- ηΓ_ii) )^m + (1-η (2 Γ_j- ηΓ_jj) )^m ]


where the transition between line <ref> and line <ref> is achieved using the following mathematical results:


  * H(-l-2p) sets all the terms in the summation to zero except those that satisfy l+2p=0 ⇔ l=0  ∧  p=0

  * H(-2(k-p)+l) sets all the terms in the summation to zero except those that satisfy -2(k-p)+l=0 ⇔ l=0  ∧  k=p

  * An extended version of the binomial theorem:

    (a+b+c)^m= 
       ∑_K=0^2m∑_q=0^⌊ K/2 ⌋ a^q   b^K-2q  c^m-K+qK-qKmK-q




 §.§  Simplification of summation over m


In the case of photon pairs generated through an SPDC process in a nonlinear crystal pumped by a weak continuous-wave laser, P(m) can be modeled by a Poisson distribution <cit.>:

    P(m) = m̅^m e^- m̅/ m !


In this case, Equations <ref> and <ref> simplify:

    ⟨ c_i   ⟩ = 1-(1-p_10)e^- m̅η(2 Γ_i- ηΓ_ii)
    ⟨ c_i c_j   ⟩ = 1-(1-p_10) [ e^- m̅η (2 Γ_i- ηΓ_ii)
       + e^- m̅η( 2 Γ_i- ηΓ_jj) ]
       + (1-p_10)^2 e^- m̅η( 2 Γ_i+2 Γ_j - ηΓ_ii - ηΓ_jj - 2 ηΓ_ij)

where the identity e^a=∑_k=0^+ ∞a^k/k! has been used. 


 §.§ Expression of Γ_ij as a function of ⟨ c_i ⟩, ⟨ c_j ⟩ and ⟨ c_i c_j ⟩ (Equation <ref>)

Combining the two previous formulae gives:

    ⟨ c_i c_j ⟩   = -1+⟨ c_i ⟩+⟨ c_j ⟩
       +(⟨ c_i ⟩-1)(⟨ c_j ⟩-1) e^- m̅  2 η^2 Γ_ij

which leads to an expression of Γ_ij in terms of the direct images ⟨ c_i ⟩ and ⟨ c_i ⟩ and the correlation image ⟨ c_i c_j ⟩:

    Γ_ij= 1/2 η^2 m̅ln [ 1+⟨ c_i c_j ⟩ - ⟨ c_i ⟩⟨ c_j ⟩/(1-⟨ c_i ⟩)(1-⟨ c_j ⟩) ]





§  MODEL OF READOUT PROCESS OF AN EMCCD CAMERA


Readout and amplification processes performed by an EMCCD camera  can be modeled using a quantitative model of noise described in <cit.>. The gray value x returned by the camera at a given pixel is modeled by a random variable X decomposed into


    X =  α( X^sig+X^par+X^ser+X^R)


where:




  * α is a scaling operation performed by the analog-to-digital converter.


  * X^sig models the output value returned because of the amplification of k photo-electrons generated by the photo-sensitive screen:


    P_sig(x|k)=x^k-1 e^-x/g/g^k (k-1)!


where g=(1+p_c)^L is the mean gain of the amplifier, where p_c is the probability of duplication of an electron at each cell of the multiplication register and L is the number of cells in the register. 


  * X^R models the readout noise:


    P_R(x)=1/σ_R √(2 π) e^-(x-μ)^2/2 σ_R^2


where μ and σ_R define respectively the mean and standard deviation of the readout noise.


  * X^par models the clock-induced charge (CIC) noise and the dark noise:


    P_cic(x)=p_pare^-x/g/g


where p_par is the probability for a spurious electron to be present at the input of the multiplication register.


  * X^ser models the electronic noise generated at each cell of the serial amplification register:


    P_ser(x)=∑_l=1^L p_sere^-x/(1+p_c)^L-l/(1+p_c)^L-l


where p_ser is the probability for a spurious electron to be generated at any cell of the multiplication register.



This model uses a total of five fitting parameters {α , p_ser ,p_par ,σ_R,μ} that can be estimated from a calibration measurement P(x|k=0) and the internal characteristics of the camera (in our case provided by Andor). Figure <ref>.a shows reconstructed response functions for k ∈ 0, 3 as well as the calibration measurement taken for an EMCCD Andor Ixon Ultra 888 operating at a horizontal shift frequency of 10MHz and vertical shift period of 0.6 μ s, at a controlled temperature of -60^∘C. Values of the fitting parameters are: {L=506, p_c = 1.37 × 10^-2 , α = 1/19 , p_ser = 3.35 × 10^-5  ,p_par = 1.23 × 10^-2 ,σ_R = 12.2 ,μ = 25.54 }. 

Using this model, calculating the mean detector response function I_k corresponds to estimating the mean value of the random variable X defined by equation <ref>:

    I_k =   ∑_x=0^+∞ x P_ccd(x|k)
    
    =   α( g k+ μ +   p_par g +   p_serg-1/p_c) 
     
    =    A  k + x_0
 
where A = 52.6 and x_0 = 569. Linearity of the response of the EMCCD Andor iXon 888 is also confirmed by experimental tests performed by Andor.

Applying a threshold onto the images recorded at the output allows the EMCCD camera to operate as an SPC camera. Details on this particular operating mode are provided in <cit.>. In our experiment, the threshold is set to the value 516 and the camera works effectively as an SPC camera with a probability P_10 = P(x<516|0) = 0.015 and a quantum efficiency η = 44 %. 
























§   DERIVATION OF FORMULAS LINKING Γ_IJ TO ⟨ X_I ⟩ AND ⟨ X_I X_J ⟩ IN THE CASE OF AN EMCCD CAMERA WITHOUT THRESHOLD (EQUATION <REF>) 

As in equation <ref>, we can write:

    I_k = A   k + x_0

where A is an amplification parameter and x_0 is a background (Details for the readout response of an EMCCD camera are given in Appendix <ref>). This result allows simplification of Equations <ref> and <ref>. Derivation of Equation <ref>, assuming a Poisson distribution for the pairs, is finally demonstrated in section <ref>.



 §.§  Simplification of summations over k_i and k_j

Simplification of ⟨ x_i ⟩ starts from Equation <ref>:


    ⟨ x_i ⟩   = ∑_m=0^+ ∞  P(m) ∑_k_i=0^2 m I_k_i∑_q=0^⌊ k_i/2 ⌋( η^2 Γ_ii)^q ( 2 η Γ_i-2 η^2 Γ_ii)^k_i-2q(1-2 η Γ_i + η^2 Γ_ii)^m-k_i+qk_i-qq mk_i-q
       = ∑_m=0^+ ∞  P(m)  [ x_0 ∑_k_i=0^2m∑_k=0^⌊ k_i/2 ⌋( η^2 Γ_ii)^q ( 2 η Γ_i-2 η^2 Γ_ii)^k_i-2q(1-2 η Γ_i + η^2 Γ_ii)^m-k_i+qk_i-qq mk_i-q
       + A  ∑_k_i=0^2m k_i ∑_k=0^⌊ k_i/2 ⌋( η^2 Γ_ii)^q ( 2 η Γ_i-2 η^2 Γ_ii)^k_i-2q(1-2 η Γ_i + η^2 Γ_ii)^m-k_i+qk_i-qq mk_i-q] 
       = x_0 + 2  A  m̅ η Γ_i


where m̅ = ∑_m=0^+ ∞ mP(m). The transition between line <ref> and line <ref> is facilitated by the identities:

    k x^k = x d [x^k]/dx

and

    ∑_k_i=0^2m∑_k_j=0^2m∑_k=0^⌊ (k_i+k_j)/2 ⌋∑_l=0^k a^m-(k_i+k_j-k) b^p  c^l  d^k-p-l e^k_i+l-2 (k-p)
       f^k_j-2 p-lk_i-q+pq-p-lk_j-l-ppk_i+k_j-q-lk_i-q+p
       k_i+k_j-qlmn_i+n_j-q  = (a+b+c+d+e+f)^m


Similarly, simplification of ⟨ x_i x_j ⟩ starts from Equation <ref>:



    ⟨ x_i x_j ⟩   = ∑_m=0^+ ∞ P(m) ∑_k_i=0^2m∑_k_j=0^2m I_k_i I_k_j∑_q=0^⌊ (k_i+k_j)/2 ⌋∑_l=0^q∑_p=0^q-l (1-2 ηΓ_i-2 ηΓ_j +η^2 Γ_ii + η^2 Γ_jj + 2 η^2 Γ_ij)^m-(k_i+k_j-q)
       ( η^2 Γ_jj)^p ( 2 η^2 Γ_ij)^l( η^2 Γ_ii)^q-p-l (2 ηΓ_i - 2 η^2 Γ_ii- 2 η^2 Γ_ij)^k_i+l-2 (q-p)  (2 ηΓ_j - 2 η^2 Γ_jj- 2 η^2 Γ_ij)^k_j-2 p-l
       k_i-q+pq-p-lk_j-l-ppk_i+k_j-q-l k_i-q+p k_i+k_j-qlmk_i+k_j-q
        = ∑_m=0^+ ∞  P(m) [ x_0^2 + 2A   x_0   m  η[ Γ_i + Γ_j ] +  A^2  [ 4m (m-1) η^2 Γ_i Γ_j + 2 m η^2 Γ_ij] ] 
        = x_0^2 + 2 A   x_0  m̅η[ Γ_i + Γ_j ] +  4 A^2 (m̅^2 + σ_m^2 -m̅)η^2Γ_i Γ_j +  2 A^2 m̅η^2Γ_ij


where σ_m^2 = ∑_m=0^+ ∞ m^2 P(m).





 §.§  Simplification of summation over m


Combining Equations <ref> and <ref> allows one to write

    Γ_ij   =1/2 A^2 m̅η^2[ ⟨ x_i x_j ⟩ - ⟨ x_i ⟩⟨ x_j ⟩
         - σ_m^2 - m̅/m̅^2 (⟨ x_i⟩ - x_0)(⟨ x_j⟩ - x_0)  ]

Assuming that P(m) follows a Poissonian distribution <cit.>, σ_m^2 = m̅ and equation <ref> simplifies to:

    Γ_ij= 1/2 A^2m̅η^2  [ ⟨ x_i x_j ⟩ - ⟨ x_i ⟩⟨ x_j ⟩ ]










§  DOUBLE-GAUSSIAN MODEL OF Γ_IJ


As described in <cit.>, the joint probability distribution of photon pairs generated by SPDC (Figure <ref>.a) can be modeled using a double-Gaussian function of the form:

    Γ_ij^th = a   e^- (x_i+x_j)^2/4 σ_+^2 e^-(x_i-x_j)^2/4 σ_-^2

where a is a normalization parameter, and σ_+ =12.06 μ m and σ_- = 926.12 μ m are two correlation lengths associated with the sum (x_i+x_j) and difference (x_i-x_j) coordinates.



§  NORMALIZATION AND BACKGROUND REMOVAL


If Γ_ij is reconstructed directly using Equations <ref> or <ref>, we observe the presence of a non-zero residual background. The presence of this background can be explained by taking into account two new factors in our general model:


  * Pump power fluctuations: intensity fluctuation of the pump implies that the assumption of a Poisson pair number distribution is no longer valid. 

  * Gain fluctuation: fluctuation of the mean gain g with time (due to, e.g., temperature drifts or variations in the high voltage clock amplitude) implies that the conditional distribution P(x|k) also becomes dependent on the image number l.


Fluctuations of g induce fluctuations in both A and x_0 over time. Equations <ref> and <ref> can be generalized to take into account the mean gain fluctuations by introducing the time averaged quantities ⟨ A ⟩, ⟨ x_0 ⟩, ⟨ Ax_0 ⟩, ⟨A^2⟩ and ⟨ x_0^2 ⟩:

    ⟨ x_i ⟩   = 2  ⟨ A ⟩ m̅ η Γ_i + ⟨ x_0 ⟩
    ⟨ x_i x_j ⟩   =2  ⟨A^2⟩m̅η^2  Γ_ij + 4  ⟨A^2⟩  (m̅^2+σ_m^2-m̅) η^2  Γ_i Γ_j  
       + 2 ⟨ A x_0 ⟩m̅η (Γ_i + Γ_j) +  ⟨ x_0^2 ⟩

These lead to

    ⟨ x_i x_j ⟩ - ⟨ x_i ⟩⟨ x_j ⟩   =2⟨ A^2 ⟩m̅η^2  Γ_ij
       +4   [  ⟨ A^2 ⟩(m̅^2+σ_m^2-m̅)- ⟨ A  ⟩^2m̅^2 ] η^2Γ_i Γ_j 
       +2  [⟨ A x_0 ⟩-⟨ A ⟩⟨ x_0 ⟩] m̅η^2  (Γ_i + Γ_j) 
       + ⟨ x_0^2 ⟩ -⟨ x_0 ⟩^2 
       = 2⟨ A^2 ⟩m̅η^2 Γ_ij+ B(Γ_i,Γ_j)

where B(Γ_i,Γ_j) is the residual background. It can be mitigated in two ways:


  * First estimate ⟨ x_i ⟩⟨ x_j ⟩ using only successive frames, rather than from the sum of all of them:

    ⟨ x_i ⟩⟨ x_j ⟩≈lim_M → +∞1/(M-1)^2∑_l=1^M-1 x_i^(l) x_j^(l+1)

The use of successive frames decreases values of the (co)variances ⟨ A^2 ⟩ -⟨ A ⟩^2, ⟨ A x_0 ⟩-⟨ A ⟩⟨ x_0 ⟩ and ⟨ x_0^2 ⟩ -⟨ x_0 ⟩^2, particularly when the fluctuations are relatively slowly varying.

  * Assume that Γ_ij has a higher spatial frequency spectrum than Γ_i and apply a low-pass filter on the reconstructed image to filter out the term B(Γ_i,Γ_j).


A similar procedure is also used to process reconstructed data in the SPC case, as discussed in <cit.>. Finally, the reconstructed Γ_ij is normalized to ensure that ∑_ijΓ_ij = 1.




§  CASE I=J


In the case i=j, Γ_ij=Γ_ii is expressed as a function of ⟨ x_i^2 ⟩ and ⟨ x_i ⟩. Similar to the calculation of ⟨ x_i ⟩,  ⟨ x_i^2 ⟩ can be written as:


    ⟨ x_i^2 ⟩   = ∑_i=0^+ ∞ x_i^2 P(x_i) 
       =∑_m=0^+ ∞ P(m) ∑ _k_i=0^2m J_k_i  P(k_i|m)
       = ∑_m=0^+ ∞  P(m) ∑_k_i=0^2 m J_k_i∑_q=0^⌊ k_i/2 ⌋( η^2 Γ_ii)^q ( 2 η Γ_i-2 η^2 Γ_ii)^k_i-2q(1-2 η Γ_i + η^2 Γ_ii)^m-k_i+qk_i-qq mk_i-k


where:

    J_k_i = ∑_x_i=0^+ ∞ x_i^2 P(x_i|k_i)


In the case of an SPC camera, the form of P(c|k) (Table <ref>) shows that J_k = I_k. Indeed, ⟨ c_i^2  ⟩ = ⟨ c_i  ⟩ and it is not possible to reconstruct Γ_ii. 

In the case of an EMCCD camera, the model of the response function provided in <cit.> shows that J_k is quadratic in k:

    ∑_x=0^+∞ x^2 P(x|k)    = α^2 ( g^2 k+ σ_R^2 +   p_par g^2 +   p_serg^2-1/p_c(p_c+2)) 
       +I_k^2  
       =  [A k + x_0]^2 + A^2 k + σ_0^2
 
where σ_0^2 is the variance of the background. ⟨ x_i^2 ⟩ can be written as

    ⟨ x_i^2 ⟩   = 2 A^2 m̅η^2 Γ_ii + 4 A^2 (m̅ + σ_m^2 -m̅) η^2 Γ_i^2 
       + 4  (A^2+Ax_0) m̅ηΓ_i + σ_0^2+x_0^2

which can be solved for Γ_ii.


