Imaging and certifying high-dimensional entanglement with a single-photon avalanche diode camera

Spatial correlations between two photons are the key resource in realising many quantum imaging schemes. Measurement of the bi-photon correlation map is typically performed using single-point scanning detectors or single-photon cameras based on charged coupled device (CCD) technology. However, both approaches are limited in speed due to the slow scanning and the low frame rate of CCD-based cameras, resulting in data acquisition times on the order of many hours. Here, we employ a high frame rate, single-photon avalanche diode (SPAD) camera, to measure the spatial joint probability distribution of a bi-photon state produced by spontaneous parametric down-conversion, with statistics taken over 10 7 frames. Through violation of an Einstein – Podolsky – Rosen criterion by 227 sigmas, we con ﬁ rm the presence of spatial entanglement between our photon pairs. Furthermore, we certify, in just 140 s, an entanglement dimensionality of 48. Our work demonstrates the potential of SPAD cameras in the rapid characterisation of photonic entanglement, leading the way towards real-time quantum imaging and quantum information processing.


INTRODUCTION
Individual single-photon avalanche diodes (SPADs) have long been the workhorse of many quantum optics experiments 1,2 . This is due largely to their single-photon level sensitivity, and also to the Geiger mode operation, which allows for straightforward methods of single-photon discrimination and counting, provided the detector operates in the photon-starved regime. Furthermore, precise timing electronics results in an impulse response function that can be as short as 20 ps 3 , which is ideal for measuring temporal correlations between multiple photons while reducing the influence of background radiation and dark counts. These properties make SPADs one of the leading technologies for measuring photon-photon correlations and entanglement.
Arrays of SPADs, or SPAD cameras, fabricated with standard CMOS (complementary metal oxide semiconductor) technology, have been produced in recent years and are now commercially available (e.g. from Photon Force, Micro Photon Devices). The maturity of the technology has enabled the production of compact arrays 4,5 , as well as the reduction of the cost per device through bulk manufacturing processes 6 . Thus far, imaging devices based on SPADs have demonstrated their capabilities in fluorescence lifetime imaging [6][7][8] , LiDAR (light detection and ranging) [9][10][11][12] , non-line-of-sight imaging 13,14 , imaging through strongly scattering media 15 and time-resolved correlation measurements 16,17 . However, SPAD cameras have yet to make their mark due to their overall efficiency and resolution; the fill factor of the earliest available SPAD cameras was on the order of a few percent 6,18 , which, despite the quantum efficiency of the single SPAD pixels being on par with commercial single element SPADs, equates to a large overall loss. This high loss is particularly detrimental to the detection of quantum states formed of multiple photons as it scales with the power of the photon number.
Charged coupled device (CCD)-based single-photon cameras have, therefore, been typically preferred for quantum imaging applications, where one of the routine tasks involves measuring spatial correlations between photons to build an image. This, in turn, is a result of the photon source of choice that is usually spontaneous parametric down-conversion (SPDC) in non-linear crystals, where a pump photon is converted with a given probability, into a pair photons. The governing law for this process is momentum conservation that ensures correlations between the photons in the pair 19,20 . These correlations can be exploited for ghost imaging [21][22][23] , imaging Bell-type non-local behaviour 24 , imaging at enhanced spatial resolution 25,26 , quantum-enhanced target detection 27 and to distil an image encoded in quantum states in the presence of classical background radiation 28,29 . We note that recently, Ianzano et al. 30 have demonstrated the measurement of polarisation entanglement using a camera, owing to its high-temporal resolution (1.5 ns). However, this demonstration did not take full advantage of the spatial resolution of the camera to measure the spatial correlations in the photon pairs; these were spatially filtered using single-mode fibres before being measured with the camera. While this ensured mode matching and spatial indistinguishability that are necessary for polarisation entanglement, the approach comes at the cost of reducing the number of modes to only one, insufficient for imaging, and limits the dimensionality of entanglement to two.
A quantity of particular interest in quantum imaging is the spatial bi-photon joint probability distribution (JPD) describing the correlations between photon pairs. The reconstruction of the JPD can be achieved through statistical averaging over a large number of intensity images 31 -typically on the order of 10 6 to 10 7 images -of identically prepared photon pairs. Given that CCD-based detectors provide frame rates on the order of 100 frames s −1 , the total acquisition time can vary from a few hours to over a day. Such long acquisition times constitute a hindrance to the widespread adoption of quantum imaging schemes for practical applications.
In the following, we show the reconstruction of the JPD in both position and momentum using a commercially available SPAD 1 School of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, UK. 2  camera (SPC3, Micro Photon Devices). Correlation measurements in position and momentum space are averaged over a total of 10 7 images, amounting to an acquisition time of 140 s. These measurements are of a quality that allows to demonstrate spatial entanglement through violation of an Einstein-Podolsky-Rosen (EPR) criterion 32,33 with a confidence of 227 sigmas. We show an experimental study into the confidence of the violation as a function of the number of individual single-bit frames used for the calculation, highlighting the benefit of the high frame rate enabled by the SPAD camera for fast characterisation of entanglement correlations. From our measurements in complementary position and momentum bases, we certify highdimensional entanglement in up to 48 dimensions. To this end, we used a robust method of entanglement certification developed by Bavaresco et al. 34 to compute the lower bound of the fidelity of the measured state with respect to a maximally entangled one in order to estimate its dimensionality. The short acquisition times required to reconstruct the bi-photon JPD and certify highdimensional entanglement pave the way for the implementation of quantum imaging and information processing schemes in real time. Indeed, we note that faster sensors, with frame rates up to 800 kiloframes s −1 , have been demonstrated 35 and that this work represents only the tip of the iceberg for the potential use of SPAD cameras within the field of quantum optics.

RESULTS
Experimental set-up Figure 1 shows the experimental apparatus used to measure spatial correlations between photon pairs. Spatially entangled photon pairs are produced via SPDC in a 0.5-mm-long β-barium borate (BBO) crystal, cut for type-I phase matching. The crystal is pumped by a 347-nm pulsed laser with a repetition rate of 100 MHz (pulse length of 10 ns), an average power of 50 mW and is spatially filtered and collimated (not shown) to a beam diameter of 0.7 mm (1/e 2 ). Spectral filters block the pump beam and select near-degenerate photon pairs at 694 ± 5 nm. Photon pairs are detected using an MPD-SPC3 SPAD camera with an array of 32 × 64 pixels, a fill factor of 80% and a 150 μm pixel pitch. The quantum efficiency at 694 nm is~9% and the dark count rate is 0.14 count per pixel s −1 . The SPAD camera has a nominal speed of 96 kiloframes s −1 and is operated at its minimum exposure time (10 ns) and dead-time (50 ns). Our photon flux over the full array was 22.5 mega-counts s −1 . To measure momentum correlations, we used a 4f-telescope (consisting of lenses f 2 and f 3 ) to image the Fourier plane of the BBO crystal surface (located 2f 1 away from the crystal) onto the SPAD camera (far-field (FF) configuration), as shown in Fig. 1a. Figure 1b shows the configuration used to measure position correlations of photons pairs at the surface of the BBO crystal (near-field (NF) configuration); the output surface of the crystal is imaged onto the SPAD camera using a 4ftelescope.
Measuring EPR entanglement We demonstrate the presence of spatial entanglement in SPDC light by violating the EPR criterion 33 where Δr = Δ(r 1 − r 2 ) and Δk = Δ(k 1 + k 2 ) are, respectively, measures of the correlation strength in position and momentum for a pair of photons labelled 1 and 2. Violation of a related EPR criterion has been achieved previously both in the case of scanning single-point detectors 33 and by full-field imaging using an electron-multiplying CCD (EMCCD) cameras 36,37 .
To determine the strength of the transverse position and momentum correlations, we measure the spatial JPD of photon pairs in the two configurations described in Fig. 1 using the SPAD camera. For a given pair of pixels located at positions r i and r j of the sensor, an element of the JPD denoted Γ(r i , r j ) represents the joint probability of detecting photon i of a pair at pixel r i and photon j at pixel r j . Γ(r i , r j ) is calculated from a set of N frames using the formula 38 where I l (r i ) ∈ {0, 1} is a binary value returned by the SPAD camera at pixel r i in the lth frame. The first term is an average value of the coincidence detection of photons belonging to either the same entangled pair (genuine coincidence) or different entangled pairs (accidental coincidence). Since multiple photon pairs can be detected during the time of an exposure, the contribution of accidental coincidences is generally greater than the genuine ones. The second term is an average value of accidental coincidences. Therefore, a subtraction between these two terms leaves only an average value of genuine coincidences that is Γ(r i , r j ) (see 'Methods'). Note that, while subtracting the accidental coincidences is important for our approach to work, this process also has the drawback of increasing the noise in the reconstructed JPD.
Reconstruction of bi-photon correlations Figure 2 shows results of JPD measurements performed in the FF configuration to study momentum correlations. Figure 2a shows the direct intensity image reconstructed from the sum of 10 7 frames, which corresponds to a total acquisition time of 140 s. This image represents the probability of detecting a photon with a given momentum k = (k x , k y ), with no information about the relative position of photons within the same pair (i.e. it is the marginal probability distribution). The JPD Γ(k 1 , k 2 ) is computed from this set of frames using Eq. (2) and is visualised using conditional projections. For example, Fig. 2b, c shows the conditional spatial distributions Γ(k|A) and Γ(k|B) of two photons measured by the sensor in coincidence and with high signal-tonoise ratio (SNR), with one photon detected at arbitrarily chosen positions A and B, respectively. Due to momentum conservation in the SPDC process, the signal and idler photons are measured at π radians in the transverse plane at the centre of the marginal distribution. The centrosymmetry of momentum correlations is characterised by the presence of an intense peak of coincidences in the sum-coordinate projection of the JPD shown in Fig. 2d. The height of the peak corresponds to the sum of coincidences measured in all pairs of symmetric pixels, while its width gives the strength of the correlation, that is, the momentum correlation width Δk. Accounting for the effective magnification of our optical system, we measured Δk = 1.0666(7) × 10 −2 rad μm −1 using a Gaussian fitting model. This model follows from the double-Gaussian approximation applied to the two-photon wavefunction of photon pairs produced by SPDC in a thin crystal 39,40 used in many experimental studies 31,33,36,37 We repeated the above analysis in the NF configuration shown in Fig. 1b to extract the position correlation width Δr. The photons in entangled pairs are position correlated, that is, they are born at the same position in the crystal and are expected to arrive at the SPAD sensor in the same position. Figure 3a, b shows, respectively, the direct intensity image and the projection of the JPD along the minus-coordinate reconstructed from a set of 10 7 frames. A coincidence peak is observed at the centre of the minuscoordinate projection that demonstrates strong position correlation between photon pairs. The central pixel in the minuscoordinate projection has been set to zero because the SPAD camera does not resolve the number of photons detected per pixel and therefore cannot measure photon coincidences at the same pixel. Accounting for the optical magnification, we measured a position correlation width Δr = 4.3(2) μm by fitting with a Gaussian model (see 'Methods').
The measured values of transverse position and momentum correlations width violate the EPR criterion in Eq. (1): Δr ⋅ Δk = 4.6 (2) × 10 −2 < 1/2, thus demonstrating the presence of spatial entanglement. This violation has a confidence level of C = 227 using the following definition: where σ = 10 −3 is the uncertainty on the product Δr ⋅ Δk.
In both position and momentum space, we observe that the correlation width is at most, as large as a single pixel on the SPAD camera. On the one hand, this is due to the relatively large 150 μm pixel pitch that is for comparison, over 10× larger than that of the EMCCD iXon 888 from Andor technologies. On the other hand, The magnification of our optical system was not large enough to spread the correlation width over more pixels. The implication is that the true correlation width as given by the crystal is necessarily smaller than what we measured (see 'Methods'). However, a more accurate and precise measurement of the correlation width will only lead to a tighter EPR violation. Note that this can now be realised with state-of-the-art SPAD sensors that boast megapixel arrays and 2.2 μm pixel pitch 41,42 . We measured an SNR of 320 (b) and 258 (c). d Projection of the joint probability distribution (JPD) along the sum coordinates k 1 + k 2 . A measured momentum correlation width of Δk = 1.0666(7) × 10 −3 rad μm −1 is obtained using a Gaussian fit (see 'Methods'). Spatial coordinates are in pixels and the analysis was performed on a total of 10 7 images. Confidence level analysis As detailed in the 'Methods' section, σ is calculated from the standard deviation of the noise surrounding the coincidence peaks in the sum-and minus-coordinate projections of the JPD. For a fixed exposure time and a constant source intensity, σ depends only on the number of frames acquired to compute the JPD 43 . Figure 4a shows that the measured values of C for different total number of frames N (black crosses) are found to scale as ffiffiffi ffi N p (blue dashed curve, see 'Methods'). In particular, confident EPR violation (C > 5) is not achieved for N < 2 × 10 4 . Figure 4b-g shows examples of sum-and minus-coordinate projections obtained for a total number of frames N = 2 × 10 3 (Fig. 4b, e), N = 10 4 (Figs. 4c, f), and N = 10 7 (Fig. 4d, g). We clearly observe a decrease of the noise in the projection images with increasing number of frames. Indeed, the average over a larger N in Eq. (2) enables to obtain a better estimation of the accidental term, such that its subtraction yields a more precise JPD. Consequently, after subtraction of the accidentals, the SNR of the JPD increases with the total number of frames, as does our confidence in the EPR violation. However, we observed that the position correlation signal appears noisier than its momentum counterpart. This is because we do not detect a significant part of our signal originating from two photons incident on the same pixel. Given that our SPAD camera cannot resolve photon number, we measure the weaker correlation signal from photons incident on adjacent pixels, owing to the relatively high 80% fill factor.
Certification of high-dimensional entanglement Measurements in two mutually unbiased bases (MUBs) allow the use of a recently developed entanglement witness for certifying high-dimensional entanglement 34 . In this respect, discrete position and momentum basis can be used as two MUBs and are accessible in our experimental set-up using the NF and FF configurations 44 . As shown in Fig. 5a, b, we then select d = 196 pixels uniformly distributed over a central region in a grid of 14 × 14 pixels in both configurations. Modes associated with the chosen pixels are denoted f m j ig m2½0;dÀ1 (discrete position basis) and fp j i p2½0;dÀ1 g (discrete momentum basis). To certify entanglement dimensionality of the generated state ρ, correlation measurements are performed in the two MUBs to compute a lower bound for the fidelity of the state ρ with respect to a maximally entangled target state Φ j i ¼ 1 ffiffi d p P dÀ1 m¼0 mm j i. From the correlation matrices measured in Fig. 5c, d, we obtained a lower bound value of the fidelity:Fðρ; ΦÞ ¼ 0:252ð9Þ (see 'Methods').
The entanglement dimensionality d ent that is certifiable with this method is the maximal r such that Fðρ; ΦÞ > r À 1 d ; allowing us to certify d ent = 48 dimensions. Finally, it is essential to note that the total acquisition time for performing all the measurements (i.e. in the two MUBs) and certifying the entanglement dimension was only 140 s. As a comparison, one day was required to certify a similar amount of dimensions of entanglement (55) using a state-of-the-art single-outcome projective measurement technique 45 , representing a speed improvement of~500× with our method. (blue dashed line). b-d Sum-coordinate projections of the JPD measured in the FF configuration using b 2 × 10 3 frames, c 10 4 frames and d 10 7 frames. e-g Minuscoordinate projection of the JPD measured in the NF configuration using e 2 × 10 3 frames, f 10 4 frames and g 10 7 frames. c, d show the measured correlations between all the pixels in the grid, with each pixel labelling the spatial coordinates of photon: m j i in the near-field (position) andp j i in the far-field (momentum). We calculated a lower boundF = 0.252(9) of the fidelity with respect to a d = 196 maximally entangled state, leading to a certified entanglement dimensionality of 48. Spatial coordinates are in pixels and the analysis was performed on a total of 10 7 images.

DISCUSSION
We have used a SPAD camera to characterise and quantify highdimensional entanglement between photon pairs. By measuring position and momentum correlations using 10 7 intensity images, we showed a violation of an EPR criterion by 227 sigmas for an acquisition time of 140 s. While EPR violation has been demonstrated by acquiring very few frames with a highly sensitive EMCCD camera 46 , quantum imaging approaches based on correlation measurements between spatially entangled photon pairs require the measurement of a large number of frames 26,28,29,31,47,48 , typically on the order of 10 6 -10 7 . This is to ensure high enough SNR on the conditional projections to reconstruct the image by exploiting photon-pair correlations. Using an EMCCD for example, the reconstruction of the spatial JPD is performed at a frame rate of the order of 100 frames s −1 , amounting to a total acquisition time of many hours. Our ability to reduce the JPD measurement time by a factor of 1000× will allow quantum imaging proof-of-principle experiments to evolve towards practical applications.
Furthermore, we have certified high-dimensional entanglement up to 48 dimensions in just 140 s, a significant development in quantum information processing with high-dimensional quantum states. While certifying high-dimensional entanglement is now routinely performed using single-outcome projective measurement techniques 34,45,49,50 , these approaches are tedious and prohibitively time-consuming. Our SPAD camera-based technique outperforms these approaches in term of speed and scalability, with a demonstrated speed-up by more than a factor 500× and the potential access to up to 4 × 10 6 two-photon modes ((32 × 64) 2 ) in parallel, indicating a key role to be played in future quantum information processing systems based on highdimensional entanglement. Finally, a natural combination with recently developed spatial mode sorting devices 51,52 holds promise for exploring high-dimensional entanglement in other types of spatial mode bases, including those carrying orbital angular momentum.

Details on Γ reconstruction
Equation (2) enables the reconstruction of the spatial JPD from a finite number of frames N acquired with the SPAD camera. This equation is derived from a theoretical model of photon-pair detection detailed in ref. 38 . In this work, a link is established between the JPD and the measured frames at the limit N → +∞: hIðr i ÞIðr j Þi À hIðr i ÞihIðr j Þi ð1 À hIðr i ÞiÞð1 À hIðr i ÞiÞ ; where A is a constant coefficient that depends on both the quantum efficiency of the sensor and the power of the pump laser, and Equation (5) is obtained under hypotheses 38 that are all verified in our work, including that (i) the quantum efficiency is the same for all pixels of the sensor and (ii) the number of pairs produced by SPDC during the time an exposure follows a Poisson distribution 53 . Moreover, in our experiment the probability of detecting a photon per pixel per frame is much lower than one (〈I(r)〉 ≪ 1), which allows us to express Eq. (5) as follows: Γðr i ; r j Þ % hIðr i ÞIðr j Þi À hIðr i ÞihIðr j Þi: In the practical case where only a finite number of frames N is measured, the first term on the right-hand side in Eq. (8) is estimated by multiplying pixel values within the same frame: The second term on the right-hand side in Eq. (8) is estimated by multiplying the averaged intensity values: I m ðr i ÞI n ðr j Þ: Combining Eqs. (8), (9) and (10) finally leads to Eq (2).

Δr and Δk measurements and uncertainties
Transverse position and momentum correlation widths, Δr and Δk, are estimated by fitting the sum-and minus-coordinate projections of the JPD measured in FF and NF configurations by a Gaussian model 36,37 of the form f(r) = a expðÀr 2 =2Δ 2 Þ, where a is a fitting parameter and Δ is the desired correlation width value (Δr or Δk). Note that, in the case of the position correlation width measurement, the central point of the minus-coordinate projection is excluded from the fitted data. The presence of noise in the sum-and minus-coordinate images induce uncertainties on values Δr and Δk returned by the fitting process. The standard deviation of the noise Σ is measured in an area composed of 40 × 40 pixels surrounding the central peak of coincidence. The link between the correlation width uncertainty δ Δ (δ Δ = δ Δr or δ Δk ) and Σ is given by calculating the value of grad[f] at the position r = Δ: and expanding it at the first order in r: In our case, the variations δf and δr identify to the uncertainty quantities Σ and δ Δ , respectively, which finally leads to: All correlation width values and uncertainties are expressed in the coordinate system of the crystal, after taking into consideration the magnifications introduced by the imaging systems detailed in Fig. 1. Variation of the confidence level C with the number of images N As defined in Eq. (3), the confidence levels C depends on both Δr ⋅ Δk and its uncertainty σ. For a given non-linear crystal and a stationary pump, Δr ⋅ Δk is constant, while σ depends on the quality of our measurement, including the total number of acquired frames N. To establish the theoretical link between C and N, we first relate σ to the uncertainties in position and momentum correlation widths (δ Δr and δ Δk ) by error propagation: Then, we replace δ Δk and δ Δr using Eq. (13) to show that σ ∝ Σ. Finally, we use the fact that Reichert et al. 43 have shown that Σ / 1= ffiffiffi ffi N p for a constant average pump power and a fixed exposure time. Thus, we conclude that C / ffiffiffi ffi N p : As shown in Fig. 4, this theoretical model fits successfully with the experimental data (R 2 = 0.998).

Effect of pixelation on the Gaussian fit
As shown in Figs. 2 and 3, the correlation widths of the photon pairs produced in our experiment are on the order of the pixel size. This pixelation effect introduces uncertainties on the Gaussian fits and the extracted values of the correlation widths. While this error is difficult to quantify, it is essential to note it can only lead to an overestimation of the correlation width values. In the limiting case where the correlation width is much smaller than the size of a pixel (i.e. the correlation image has only one non-zero pixel surrounded by zeros), the Gaussian fit will always return a value that is approximately equal to 20% of the pixel size. Therefore, the confidence value calculated from the correlation width values is a lower bound of the more accurate value that would be measured using a SPAD with more pixels. As a result, the pixelation effect cannot change our conclusion on the presence of entanglement. Simulations of the pixelation effect are discussed in Supplementary Note 1, and graphically illustrated in Supplementary Fig. 1.

Discrete position and momentum basis
In our experiment, we use discrete photonic transverse position and momentum bases given by a set pixel defined on the SPAD camera, which we refer to as the discrete position basis f m j ig m2½1;d and discrete momentum basis fp j ig p2½1;d . Our approach is based on the protocol proposed by Erker et al. 44 , where these bases are used as two MUBs to certify high-dimensional entanglement. More precisely, they are linked according to: where ω = e 2πi/d . Experimentally, these bases are accessed using simple optics, that is, lenses to image or Fourier image the output of the non-linear crystal producing the photon pairs. A subset of pixels is then selected in the illuminated areas of the sensor to optimise coincidence signals measured. In our case, we selected 196 pixels evenly separated from each other by 1 pixel and located on a square grid of 14 × 14 pixels at the centre of the sensor. The only difference between the scheme proposed by Eker et al. and our work is that only one image is produced on the camera in our case, against two in their proposal using a collinear beam splitter. This set-up prevents us from accessing the coincidence rate at the same pixel (e.g. same spatial modes) as already pointed out in Fig. 3. Rather, the coincidence rate of photon pairs in the same pixel is inferred by measuring the coincidence rate between each pixel and its neighbour. This inference leads to a lower value of the intra-pixel coincidence rate in the NF configuration. This approximation, therefore, underestimates the actual entanglement witness value.

Derivation of the dimensionality witness
To certify the presence of high-dimensional entanglement in the measured state ρ, we employ a recently developed witness that uses correlations in at two MUBs 34 . In our experiment, the two MUBs are the discrete position basis f m j ig m2½1;d and discrete momentum basis fp j ig p2½1;d (see previous 'Methods' section). Using only coincidence measurements in these two bases, one can determine a lower bound for the fidelity F(ρ, Φ) of the state ρ to a pure bipartite maximally entangled target state Φ j i. Since the fidelity to a target entangled state also provides information about the dimensionality of entanglement, we use this bound for certifying the dimension of entanglement of the state produced in our experiment.
We consider a maximally entangled target state written as: with d = 196. The fidelity F(ρ, Φ) of the state ρ to the target state Φ j i is defined as: where The entanglement dimensionality can be deduced from the fidelity taking into account that for any state ρ of Schmidt number r ≤ d, the fidelity of Eq. (20) is bound by: Hence, any state with F(ρ, Φ) > B r (Φ) must have an entanglement dimensionality of at least r + 1. Our goal is therefore to obtain a lower bound on the fidelity as large as possible for the target state whose Schmidt rank is as close as possible to the local dimension d. To achieve this experimentally, the method described in ref. 34 works the following way: Step 1: Matrix elements {〈mn|ρ|mn〉} m,n are calculated from the coincidence counts fN mn g mn measured in the discrete position basis via: These elements are shown in the matrix in Fig. 5c. They enable to calculate directly the term F 1 (ρ, Φ) = 0.00170(1) from the definition, Eq. (18).
Step 2: Matrix elements fhpṽjρjpṽig p;v are calculated from the coincidence counts fÑ pv g pv measured in the discrete momentum basis via: hpṽjρjpṽi ¼Ñ pv P k;lÑ kl : These elements are shown in the matrix in Fig. 5d. These matrix elements, together those of the discrete position basis, allow us to bound the fidelity term F 2 (ρ, Φ). This lower boundF 2 ðρ; ΦÞ ¼ 0:250ð9Þ is calculated via: where the prefactor γ mnm 0 n 0 is given by A derivation of Eq. (23) can be found in the 'Methods' section of ref. 34 .

Assumptions
Using this approach, we note that no assumptions are directly made about the underlying quantum state ρ. However, an assumption is made about our measurement process. Indeed, by using Eq. (2) to measure the JPD of photon pairs, we effectively perform a subtraction of accidental counts. Correcting for accidental coincidence is acceptable in our experiment since we trust our measurement devices and the final goal is only to assess the presence of entanglement and its dimension. However, such an assumption would not be acceptable in an adversarial scenario such as quantum key distribution as it is likely to compromise the security of the protocol. Furthermore, this assumption cannot be used if one wants to perform a loophole-free nonlocality tests, because it is likely to violate the fair-sampling assumption.

DATA AVAILABILITY
The experimental data that support the findings presented here are available at https://doi.org/10.5525/gla.researchdata.1071.

CODE AVAILABILITY
The codes used to process the data are available from the corresponding authors upon reasonable request.