Characterization of space-momentum entangled photon with a time resolving CMOS SPAD array

Single photon avalanche diode arrays can provide both the spatial and temporal information of each detected photon. We present here the characterization of entangled light with a sensor specifically designed for quantum imaging applications. The sensors is time-tagging each detection events at the pixel level with sub-nanosecond accuracy, within frames of 50 ns. The spatial correlations between any number of detections in a defined temporal window can thus be directly extracted from the data. We show the ability of the sensor to characterize space-momentum entangled photon pairs emitted by spontaneous parametric downconversion. Their entanglement is demonstrated by violating an EPR-type inequality.

Single photon avalanche diode arrays can provide both the spatial and temporal information of each detected photon. We present here the characterization of entangled light with a sensor specifically designed for quantum imaging applications. The sensors is time-tagging each detection events at the pixel level with sub-nanosecond accuracy, within frames of 50 ns. The spatial correlations between any number of detections in a defined temporal window can thus be directly extracted from the data. We show the ability of the sensor to characterize space-momentum entangled photon pairs emitted by spontaneous parametric downconversion. Their entanglement is demonstrated by violating an EPR-type inequality.

I. INTRODUCTION
Quantum states of light are a fundamental tool to implement quantum information processing protocols and quantum metrology methods in the optical domain, as quantum imaging [1]. They can be described either by discrete variables, the photons, or by continuous fields. Correspondingly measurements of light can be mainly split into two classes: measurements of fields, by homodyne measurements for example, and intensity measurements, that correspond fundamentally to photon counting. These are the tool of choice to characterize discrete quantum optical states; for example to measure correlations in quantum states with few photons. An ideal detector should be able to detect the position and time of arrival of all impinging photons. Such universal detector do not yet exist, but there are various type of single photon sensitive detectors, that are optimized with respect to specific characteristics. Most of the experiments based on correlation between photons are realized with single photon detectors with high temporal resolution but no spatial resolution [2], including photomultipliers, single photon avalanche diodes (SPAD) or superconductor detectors. On the other side, imaging experiments are usually implemented, even at the single photon level, with low noise cameras as EMCCD or scientific CMOS cameras [3,4], that are able to localize single photons among many pixels, but have low temporal resolution. Spatial entanglement in photon pairs have been demonstrated with EMCCD [5] A way to combine high temporal and spatial resolution is fast gating of imaging sensors. Detection of spatial entangled photon pairs have been realized with ICCD cameras as in [6] and commercial gated single photon avalanche array sensors [7]. However, while gating gives access to short times, the full spatio-temporal correlations of light across a wide range of detection times cannot be measured. In addition, gating is well suited for pulsed operations, but a continuous source of correlated photons would ideally be measured with a non-synchronous detector. This is easily implemented with single pixel detectors, that are always ready to trigger, up to their deadtime. However when operating many detectors in parallel, as in detector arrays, the readout mechanism usually requires a frame-based operation. In that case one of the criterion to be optimized is the sensor duty cycle, given by the frame rate of the detector times the duration of each frame. For gated detectors, the requirement to have short gates to improve the temporal discrimination of the correlated photons contradicts the need of long frames to improve the photons collection efficiency. The solution is to be able time-stamping with high temporal resolution each detection event within a long frame. Combining single photon detectors and external time-stamping electronics can be realized only for a limited number of detectors [8]. This cannot however be extended to imaging with thousands of pixels without integrating light sensitive devices with digital electronics. CMOS offers the possibility of full integration for silicon sensors. Here we demonstrate the capability of a recently developed SPAD array sensor [9,10] to characterize spatial entanglement from spontaneous parametric downconversion (SPDC) by measuring both near-and far-field correlations [5,11]. We demonstrate entanglement of the state by testing an Einstein Podolsky Rosen (EPR)-type inequality. We estimate the violation of the inequality both directly from the measured join probabilities, without a priori assumption on the quantum state, and from fitted data assuming a gaussian model for the SPDC emission. In section II we recall the quantum state of the photon pairs emitted by SPDC and of the expected join probability distribution that have to be measured in order to verify their spatial entanglement, by violating an EPR-type inequality. Section III describes the experimental setup and the processing of the data acquired with the SPAD array with corrections of the accidental events and cross-talk. Finally, the results are presented in section IV, where the different estimations of the violation of the inequality are presented.

II. THEORY OF SPATIAL ENTANGLEMENT IN SPDC
Spontaneous Parametric Down-Conversion (SPDC) is nowadays a common source of entangled photon pair, where a pump photon of frequency ω p is annihilated inside a non-linear crystal (NLC) with a non-vanishing second order susceptibility χ (2) and two photons with frequencies ω 1 and ω 2 and wavevectors k 1 and k 2 (historically called signal and idler ) are simultaneously created [12]. Here we consider the case of type-0 quasi-phase-matching in a periodically poled crystal of length L with poling period G and refractive index n(ω), where all fields have the same polarization.
Under the approximation of a monochromatic pump with a frequency ω cp , the envelop of the pump field is given by The quantum state of the generated signal and idler photons is then given by with the one photon Fock states |q j , ω j =â † (q j , ω j ) |0 j and where q i = (k xi , k yi ) is the transverse component of the wave vector k i . The two-photon wavefunction, or joint spectral amplitude (JSA) is with the phase matching relation is and the normalization function We note that the JSA Λ(q 1 , q 2 , ω 2 ) given by (3) cannot in general be factorized into functions Λ i (q 1 , ω 2 ) and Λ s (q 2 , ω 2 ) and thereby the state is entangled. Assuming narrow-band frequency filtering of the down-converted photons and thus only considering the spatial dependency, the JSA factorizes into momentum and frequency parts [13] Equivalently the JSA can be expressed in term of transverse positions (ρ 1 , ρ 2 ) and times (t 1 , t 2 ) through Fourier transform FΛ In the monochromatic approximation, the temporal part only depends on the time difference τ = t 1 −t 2 . Coincidence detections measure events around τ = 0. Therefore the temporal part of the JSA only account for a proportionality factor and thus will be further omitted.
The position ρ i = (x i , y i ) and momentum q i = ( k xi , k yi ) operators do not commute and therefore a criterion for EPR-type correlations between two systems 1 and 2 (that correspond in our case to the signal and idler photons) can be established [14]. In the SPDC emission, the x and y components are decoupled. We thus consider the one dimensional case with x i and p i = k xi . The results of correlation measurements are therefore described by the join probabilities P(ρ 1 , ρ 2 ) and P(q 1 , q 2 ), or P(x 1 , x 2 ) and P(p 1 , p 2 ) in the 1D case. The minimal inferred variance of (x 1 , x 2 ) is defined by where ∆ 2 (x 1 |x 2 ) is the variance of the conditional probability P(x 1 |x 2 ) and P(x 2 ) is the marginal probability of system 2. The same definitions apply for the variables (p 1 , p 2 ). According to [14], the fulfillment of the following inequality on the product of minimum inferred variances indicates the EPR-type correlations or expressed as a dimensionless quantity Experimentally the join probabilities are derived from correlation measurements, that correspond to second order coherence functions of the electric field. They are given at transverse positions ρ 1 and ρ 2 and times t 1 and t 2 by Correspondingly, in the monochromatic approximation the join probability is directly proportional to the second order coherence function and is explicitly related at the crystal position to the spatial JSA of the state (2) by In general, the imaging setup from the crystal to the sensor plan determines the field operators that have to be introduced in (11). In the following we performed the two types of measurements relevant for inequality (10): nearand far-field imaging.

A. Second Order Near-Field Correlations
When imaging the crystal with a magnification M , the correlation function at the image plan using Fourier optics can be expressed as This is essentially the Fourier transform of the propagated two-photon wave function.

B. Second Order Far-Field Correlations
In a far-field imaging setup, using a lens of focal length f , the relation between the transverse position at the imaging plan ρ and the transverse momentum q at the object is where k is the light wave-vector. Hence, the second-order correlation function in the far-field reads The spatial dependency of the correlation function thus directly reflects the JSA expressed in the momentum space.

III. EXPERIMENTAL EVALUATION OF CORRELATIONS
The join probability distribution that leads to the evaluation of Eq. (10) has to be estimated from the raw correlations acquired with the sensor. Ideally it should be estimated directly from the non-processed data, in order to achieve a non-conditional violation of the inequality. However due to the sensor imperfections: reduced efficiency, dark counts and cross-talk between pixels, some assumptions have to be introduced in order to process and correct the raw data. In a second step, the corrected data can be either directly used to numerically evaluate Eq.(10) or, by introducing additional assumptions, fitted with a model of the expected correlations as for instance in [5]. The experimental setup is shown in Fig.1. A continuous wave (CW) laser (Toptica DL PRO HP 405) with a maximal power output of 30 mW at a wavelength of 405 nm and a spectral bandwidth of 80 MHz serves as the pump. It is slightly focuses onto the crystal's center plane Σ 0 by a two lens system (f 1 = f 2 = 200 mm). The beam waist of the pump at this plane is w 0x = 250 µm in the x direction and and w 0y = 300 µm in the y direction. λ/4-and a λ/2-plates are used to achieve the desired horizontal polarization in the crystal. A 1 × 2 × 12 mm 3 periodically poled KTiOPO 4 (PPKTP) non-linear crystal with poling period G 0 = 3.510 43 µm is embedded into a temperature controlled oven and provides the source for the down-converted photons. The oven is maintained at a temperature of 26.0 • C for an almost collinear phase-matching, that maximize the photon flux onto the sensor in the far field. The down-converted photons are separated from the pump by a bandpass filter (BP) centered at 810 nm and with FWHM 10 nm. Given the source parameters, waists of the pump beam and the length of the SPDC crystal, we can estimate the minimal inferred standard deviations to be and the violation of the Heisenberg-inferred inequality to be Two lenses (f 3 = f 4 = 50 mm) form a 4f imaging system such that the electric field at the plane Σ N F is an exact replica of the field in the crystal at Σ 0 . Either the near-or far-field of this plan is then imaged onto the sensor, depending on the selected lens. For near-field, a lens f N F = 25.4 mm images the plane Σ N F onto the sensor with a magnification factor M = 9. For far-field, a lense f F F = 150 mm images the far-field such that the area covered by the sensor corresponds in k-space to ±36.2 mm −1 , as q = x k f F F . The coincidence detection is performed with a single photon avalanche detector (SPAD) array recently developed [9]. It is a fully digital 32×32 pixels sensor array based on CMOS technology. The total sensitive area is 1.4×1.4 mm 2 and the pixel pitch ∆L is 44.67 µm with an overall fill-factor of 19.48 %. A Time-to-digital converter (TDC) integrated in each pixel allows for pixel-wise time stamping the first detection event within a frame. The time resolution (timebin length) of the TDC is 205 ps and the frame length is 255 time-bins that corresponds to about 50 ns. A maximal frame observation rate of ∼ 1 MHz can be reached, which corresponds to a duty cycle of 4.5 % . The total photon detection efficiency is 5 % at 400 nm and 0.8 % at 810 nm. The dark count rate is below 1 kHz at room temperature on average. The data from the sensor's FPGA are transfered to a computer through a USB 3.0 interface. A Labview interface controls the acquisition, allowing for either on-line simple processing of the data, or acquisition of all raw sensor data. For the following measurements, we acquired the raw data and extracted the correlations off-line, with Matlab and C programs.

B. Data processing
The coordinate and time of all detection events within one frame, together with the frame identification number, constitute the acquired raw data. In the following we associate the spatial coordinate ρ i = (x i , y i ) with the pixel coordinate p i = (p xi , p yi ) using ρ i = ∆Lp i . An estimation of the correlation G (2) (p 1 , p 2 ) between a pair of pixels is obtained by counting all coincidence events between those pixels that occur within a defined temporal window. The histogram of all time differences between pairs of detection events across the sensor is shown on Fig. 2. While the temporal jitter of each pixel is of the order of 200 ps, in practice temporal shifts across the sensor widen the coincidence peak. In principle, a pixel-wise temporal calibration would allow to narrow the window. Here we define the coincidence window to be 10 TDC steps (∼ 2 ns) in order to catch all coincidence events. The obtained histogram is further normalized to a number of coincidence per frames or million of frames (MFrames) [13]. The second order correlation in a transverse plane is a function of four variables G (2)  In order to be able to plot a 2D representation of the full correlations, we introduce the pixel indexes numerating each pixel from 1 . . . 1024p i = p xi + 32 · (p yi − 1).
(21) Figure 3 shows the full correlation matrix between every pair of pixels after accidentals subtraction (see bellow) in the case of far-field measurement. The photon pairs exhibit anti-correlation in their detection position, that reflects into the ≈ 8 anti-diagonal lines. The exact diagonal cancels, as the sensor cannot measure self-correlations on one pixel. The nearest neighbors correlations (separated by pixel indexes ±1 and ±32) are affected from cross-talk, forming the four diagonal lines. In order to apply the inequality (10) for only one dimension, the correlation functions are projected onto either x or y dimnesions by Alternatively, one can use the centroid and difference coordinates to project the 4D correlation function on 2D plans defined by fixed value of ρ 1 + ρ 2 or ρ 1 − ρ 2 , respectively. The projected correlation functions read .

Removing Accidentals
Accidental coincidences are coincidences events that occur in the defined coincidence window but are neither temporally nor spatially correlated. They steam from detections triggered by background light, dark counts, and mostly from SPDC photons that do not belong to the same photon pair. The total measured raw correlation signal is thus given by These accidentals events contributes significantly to the raw signal as we can observe on Fig. 4 (a). They can be however estimated and corrected for, either by measuring the coincidence signal in a shifted time window or by using the fact that uncorrelated light obeys The proportionally factor depends on the coincidence window length, and is estimated by comparing event rates in G (2) (ρ 1 , ρ 2 ) and G (1) (ρ 1 ) ⊗ G (1) (ρ 2 ) for pixels for which no correlations are expected. G (1) (ρ) is the first order correlation function, or corresponding to the intensity. Figure 4 shows the raw (a) and corrected data after removing accidental coincidences (b). The anti-correlations from the photon pairs is clearly visible, but also correlations due to cross talk, that are further removed in Figure 4 (c) (see bellow).

Cross-talk
A single detection event, triggered by a dark count or a photon can trigger nearby pixels. This leads to undesirable detection events called cross-talk and are an artifact intrinsic to the detector. The physical process behind optical crosstalk are photons, created by the charge avalanche of the first triggering SPAD. These photons may reach and trigger neighboring pixels, in a very short time scale. Cross-talk is especially undesired in the near-field correlations, as it overlays the real signal. The cross-talk probability, averaged over all pixels, can be extracted from the correlationpeak in the far-field, as in that case only anti-correlations are expected. Hence, the mean probability that a pixel at distance (∆x, ∆y) is triggered from cross-talk is given by where the sum over G (1) F F normalizes with the total number of counts. To take into account that we cannot measure all cross-talk events of the pixels close to the sensor's edges, only an inner window of the total sensor area is used to estimate the cross-talk behavior. G (2) F F in (28) is reduced to cover only an inner window of 29×29 pixels. Consequently, the normalization is also adjusted to the total counts within this window. The second order correlation measurements can then be corrected [15] according to where G (2) m is the measured correlation function. Fig. 5 the anti-correlation peak of the far-field signal (a) is shown with corresponding estimated cross-talk probability map (b), which is point-symmetric with respect to (∆x, ∆y) = (0, 0). Cross-talk probability extracted from the correlation peak and normalized with the total number of counts. The pixel at (∆x, ∆y) = 0 is the emitter of the cross-talk.

IV. MEASUREMENT RESULTS
The full correlation map results of far-field measurements are shown on Fig. 4. The same data can then be represented in the other coordinates previously introduced. Fig. 6 shows the anti-correlation (a) and correlation (b) peaks in the far-field after cross-talk removal. As far-field measurements map the momentum space onto the sensor, the plot are labeled with the coordinates q + and q − , but still in units of pixel. We clearly observe the presence of an anti-correlation peak and the almost disappearance of the correlation peak.
The same processing can be applied for the near-field measurements. The full correlation matrix are shown in Fig. 7 without (a) and with (b) cross-talk correction. We observe that the cross-talk along the diagonal overlays the actual signal, even after correction. As the cross-talk contribution is of the same order as the signal, a correction is very sensitive to the estimation of the cross-talk probabilities. A full calibration at the pixel level will be performed in the future. This is why in the present work, the next neighbor correlations will not be taken into account in the quantitative estimation of the correlations. The partial suppression of the cross-talk can also be see on the correlation peaks of Fig. 7 (c) and (d), where the bright central peak is due to cross-talk, while the broad peak indicates correlation between photons. Note that the magnification factor has been chosen in order to obtain a correlation peak broader than the cross-talk. While the identification of EPR-type correlations should be ideally done on the raw data, usually some additional assumption are used to correct for the imperfections of the sensor. Here, accidentals and cross-talk corrections are needed, as the detector sensitivity is limited at the photons wavelength. In a second step, the data from the second order correlations can be either directly numerically evaluated, or fitted with a model of the SPDC emission. The results of the various evaluations are summarized on table I.

Numerical Evaluation of the EPR inequality
First, we present in detail the numerical evaluation for the x coordinates, i.e. G (2) (x 1 , x 2 ) and G (2) (q x1 , q x2 ). The treatment for the y coordinate is identical. First the conditional and unconditional probability density functions are estimated from the corrected experimental data Examples of experimentally obtained probability distributions are shown on Fig. 9. Next, expectation values are given by and the conditional variances are Finally, the minimum inferred variances are computed We introduce the quantity and EPR-type correlations are identified by violating the inequality The same procedure is applied for the y coordinate and the result is denoted as The direct numerical evaluation of V min from the corrected measured data is sensitive to noise and to the accidental and cross-talk corrections. The measured conditional probability distributions take non-zero values even far away from their expected peak value, as can be seen in Fig. 9. This effect shifts the expectation value µ x1,qx1 away from the true one and also increases the conditional variances. The second challenge arises from the cross-talk. Setting to 0 the pixels mostly affected in the near-field correlations increases the variance of the conditioned probability distributions. As a consequence the numerically determined values of V min given in Table I  An alternative approach to direct numerical evaluation is to constrain the data into a set of "reasonable" data. The shape of the JSA from SPDC emission can be, in some conditions, approximated by Gaussian functions. Therefore fitting Gaussian functions to the data is a way to to obtain P (q x2 ), P (x 2 ), P (q x1 |q x2 ) and P (x 1 |x 2 ). Figures 9 and 10 show example of fitted data. One can then extract the variances directly from these fitted probability density functions and calculate the EPR-criterion according to the definition of Eqs. (8) and (10). The results are shown in the second line of Table I and are in agreement with the expected values from the experimental parameters. . 9: Example of the measured conditional probabilities for the x variable in the far-field P(q x1 |q x2 ) (a) and near-field P(x 1 |x 2 ) (b). The variance is either extracted numerically or using a Gaussian fit. The pixels affected from cross-talk are skipped (indicated by the red asterisks) FIG. 10: Example of the measured unconditional probabilities for the x variable in the far-field (a) and near-field (b).

2D Gaussian Fitting
Instead of fitting multiple Gaussians by slicing the data sets G (2) , one can also directly fit 2D Gaussians over G (2) (x 1 , x 2 ) and G (2) (q x1 , q x2 ) (and similarly for y), as shown on Fig. 11. Note that we only assume the correlations to be well-approximated by 2D Gaussians with variances ∆ 2 x+,x− and ∆ 2 qy+,qx− , but we don't impose the quantum state to follow the double gaussian model [16]. The conditional variances are then related to the fitting parameters by [17]

(Anti-) Correlation Peaks
The projections onto x and y of the (anti-) correlation peaks G (2) (q + ) and G (2) (ρ − ) are shown on Fig. 12. The width of the fitted Gaussian are respectively σ x+ , σ y+ and σ qx+ , σ qy+ . They can be related to the conditional variances if we assume the correlations to appear along the + and − coordinates. Next we require the correlations in the nearfield to be expanded very far along the + coordinate (i.e. σ x+ , σ y+ → ∞) and similarly in the far-field for the − coordinate (σ qx+ , σ qy+ → ∞). We can then estimate the unconditioned variances from G (2) (x − ) and G (2) (q x+ ) (see Fig. 12). The conditional variances (or its square roots) are then given by and similarly for the y coordinate.

V. DISCUSSION
We make use of the capacity of a recently developed 32x32 pixels CMOS SPAD array sensor to time-tag each individual photons with high temporal resolution in order to measure the emission of spatially entangled photon pairs from SPDC. Anti-correlations in the far-field and correlations in the near-field are observed. From the experimentally measured probability distributions, the violation of an EPR-type inequality is demonstrated. Interestingly the inequality is violated even without fitting the data with a specific model. When assuming Gaussian shape for the probability distributions, the obtained value for the variances are in good agreement with the expected ones derived from the parameters of the photon source. Those results demonstrate the ability of CMOS SPAD arrays to effectively measure simultaneously spatial and temporal correlations between photons. They are a promising tool for practical implementations of quantum imaging schemes, that will require further developments towards higher number of pixels.