# Pixel super-resolution using spatially-entangled photon pairs

Hugo Defienne${}^{\,1,^{*}}$ Patrick Cameron 1 Bienvenu Ndagano 1 Ashley Lyons
1 Matthew Reichert 3 Jiuxuan Zhao 2 Edoardo Charbon 2 Jason W. Fleischer 3
Daniele Faccio${}^{\,1,^{\dagger}}$
1School of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, UK
2Advanced Quantum Architecture Laboratory (AQUA), Ecole Polytechnique Federale
de Lausanne (EPFL), 2002 Neuchatel, Switzerland
3Department of Electrical and Computer Engineering, Princeton University,
Princeton, USA

(August 28, 2024)

###### Abstract

Pixelation occurs in many imaging systems and limits the spatial resolution of
the acquired images. This effect is notably present in quantum imaging
experiments with photon pairs, in which the number of pixels used to detect
coincidences is often limited by the sensor technology or the acquisition
speed. Here, we introduce a pixel super-resolution technique based on
measuring the full spatially-resolved joint probability distribution (JPD) of
spatially-entangled photon pairs. Without shifting optical elements or using
prior information, our technique doubles the pixel resolution of the imaging
system and enables retrieval of spatial information lost due to undersampling.
We demonstrate its use in various quantum imaging protocols, including quantum
illumination, entanglement-enabled quantum holography, and in a full-field
version of N00N-state quantum holography. Our JPD super-resolution technique
will benefit any full-field imaging system limited by the sensor spatial
resolution, including all already established and future photon-pairs-based
quantum imaging schemes, bringing these techniques closer to real-world
applications.

††preprint: APS/123-QED††Corresponding author:
∗<EMAIL_ADDRESS>
†<EMAIL_ADDRESS>
Acquiring high-resolution images over a large field of view is essential in
most optical imaging applications. In this respect, the widespread development
of digital cameras made of millions of pixels has strongly contributed to
create imaging systems with large space-bandwidth products. In classical
imaging, it is therefore mainly the imaging systems with strong spatial
constraints that suffer from pixelation and undersampling, such as lensless
on-chip microscopes Bishara  _et al._ (2010, 2011). However, this effect is
very present in quantum imaging schemes where it severely hinders the progress
of these techniques towards practical applications.

Quantum imaging systems harness quantum properties of light and their
interaction with the environment to go beyond the limit of classical imaging
or to implement unique imaging modalities Moreau  _et al._ (2019). Of the many
approaches, imaging schemes based on entangled photon pairs have turned out to
be among the most promising. Proof-of-principle demonstrations range from
improving optical resolution Boto  _et al._ (2000); Reichert  _et al._ (2018)
and imaging sensitivity Ono  _et al._ (2013); Israel  _et al._ (2014); Brida
_et al._ (2010) to the creation of new imaging protocols, such as ghost
imaging Pittman  _et al._ (1995); Aspden  _et al._ (2013), quantum
illumination Lopaeva  _et al._ (2013); Defienne  _et al._ (2019); Gregory  _et
al._ (2020) and quantum holography Devaux  _et al._ (2019); Defienne  _et al._
(2021a). Contrary to classical imaging, photon-pairs-based imaging systems
operate by measuring photon coincidences between many spatial positions of the
image plane in parallel (except from induced-coherence imaging approaches
Lemos  _et al._ (2014); Kviatkovsky  _et al._ (2020); Vanselow  _et al._
(2020)). In practice, this process is much more delicate than forming an
intensity image by photon accumulation and therefore requires specific
photodetection devices. Although such a task was originally performed with
raster-scanning single-photon detectors Pittman  _et al._ (1995), today most
implementations use single-photon sensitive cameras such as electron
multiplied charge coupled devices (EMCCD) Moreau  _et al._ (2012); Edgar  _et
al._ (2012), intensified complementary metal-oxide-semiconductor (iCMOS)
Chrapkiewicz  _et al._ (2014) and single-photon avalanche diode (SPAD) cameras
Unternährer  _et al._ (2016).

EMCCD cameras are the most widely used devices for imaging photon-pair
correlations thanks to their high quantum efficiency, low noise and high pixel
resolution. However, these cameras suffer from a very low frame rate
($\sim$100 fps) due to their electronic amplification process operating in
series and therefore require very long acquisition times ($\sim$10 hours) to
reconstruct correlation images Toninelli  _et al._ (2019); Defienne  _et al._
(2019); Gregory  _et al._ (2020); Defienne  _et al._ (2021a). Using a smaller
sensor area or a binning technique can reduce the acquisition time, but at the
cost of a loss in pixel resolution. By using an image intensifier, iCMOS
cameras do not use such slow electronic amplification processes and can
therefore reach higher frame rates ($\sim$1000 fps). However, so far these
cameras have only enabled correlation images with relatively small number of
pixels (and still several hours of acquisition), mostly because of their low
detection efficiency and higher noise level Chrapkiewicz  _et al._ (2014,
2016). Finally, SPAD cameras are an emerging technology that can detect single
photons at low noise while operating at very high frame rate (up to $800$
kfps). However, most devices available on the market are not yet suitable for
photon pair imaging because of their low quantum efficiency ($\sim$20$\%$),
low fill factor ($\sim$50$\%$) and small resolution ($\sim$1000 pixels). For
example, the commercial SPAD camera that currently has the highest pixel
resolution only has $32\times 64$ pixels (MicroPhotonDevice) Bronzi  _et al._
(2014).

Figure 1: Experimental results of JPD pixel super-resolution. (a) Experimental
setup. Light emitted by a diode laser (405 nm) illuminates $\beta$-barium
borate (BBO) crystal with a thickness of $0.5$ mm to produce spatially-
entangled pairs of photons by type I SPDC. Long-pass and band-pass filters at
$810\pm 5$nm (BPF) positioned after the crystal remove pump photons. A two-
lens system $f_{1}-f_{2}$ images the crystal surface onto an object $t$, that
is itself imaged onto the EMCCD camera by another two-lens system
$f_{3}-f_{4}$. $t$ is a grid-shaped amplitude object with a five-pixel period.
(b) Intensity image. (c) Sum-coordinate projection of the JPD
$P_{+}(\mathbf{r^{+}})$ after filtering. (d). Spatial frequency spectrum of
the intensity image (red curve) and sum-coordinate projection (blue curve)
obtained by a discrete Fourier transform and averaging aver the $x$-axis.
Spectral curves (blue and red) are normalised to the intensity of their
respective zero-order peak. $5.10^{6}$ frames were acquired to retrieve (b,c)
in approximately $13$ hours of acquisition.

It is clear from the above analysis that nearly all sensor technologies
currently used in photon-pairs-based imaging experiments suffer from a poor
pixel resolution, either directly when the technology does not have cameras
with enough pixels or indirectly when the experiment can only operate in a
reasonable time with a small number of pixels. In these systems, objects are
therefore often undersampled, resulting a loss of spatial information and the
creation of artefacts in the retrieved images. Here, we demonstrate a quantum
image processing technique that doubles the pixel resolution of any imaging
sensor used in photon-pairs-based imaging schemes. Our technique harnesses the
large dimensionality of the measured spatial joint probability distribution
(JPD) of photon pairs to reconstruct correlation images composed of twice as
many pixels as the camera used to measure it. We demonstrate the method in
three common photon-pairs-based imaging schemes: two quantum illumination
protocols using (i) a near-field illumination configuration with an EMCCD
camera Defienne  _et al._ (2019) and a (ii) far-field configuration with a
SPAD camera Defienne  _et al._ (2021b), and (iii) one entanglement-enabled
holography system Defienne  _et al._ (2021a). In addition, we use our JPD
pixel super-resolution method in (iv) a full-field version of N00N-state
quantum holography, a scheme that has only been demonstrated so far using a
scanning approach Ono  _et al._ (2013). Note that we refer to our technique as
’pixel super-resolution’ to avoid confusion with the term ’super-resolution’
describing imaging techniques capable of overcoming the classical diffraction
limit.

JPD pixel super-resolution technique. Figure 1.a describes the experimental
setup used to demonstrate the principle of our technique. In this scheme,
spatially entangled photon pairs produced by type-I SPDC in a thin
$\beta$-barium Borate (BBO) crystal illuminate an object $t$ using a near-
field illumination configuration (i.e. the output surface of the crystal is
imaged onto the object). The object is then imaged onto an EMCCD camera. In
our experiment, $t$ is a horizontal square-modulation amplitude grating. The
recorded intensity image is shown in Figure 1.b. In the zoomed area, it is
clear that the image is pixelated. This pixelation results in a loss of object
spatial information, confirmed when analysing the spatial spectrum of the
image (Figure 1.d, red curve). The spectrum shows only the first two harmonics
at $0.2$ cycles/pixel and $0.4$ cycles/pixel; all harmonics higher than the
Nyquist frequency of $0.5$ cycles/pixel are not resolved, and their associated
information is lost.

Figure 2: Spatial frequency analysis of images of a square-wave amplitude
grating with different periods, produced using a spatial light modulator (SLM)
operating in amplitude modulation (graings shown in insets). The spatial
frequency spectrum of the intensity (red curves) and sum-coordinate projection
(blue curve) images of gratings are shown, with periods of (a) $3$ pixels, (b)
$2$ pixels and (c) $1.6$ pixels. Each spectrum is obtained by performing a
discrete Fourier transform and averaging aver the $x$-axis. All spectral
curves (blue and red) are normalised to the intensity of their respective zero
order peak. $2.5.10^{7}$ frames were acquired to retrieve intensity image and
JPD in each case in approximately $20$ hours of acquisition. More details in
Methods.

Using an EMCCD camera, one can also measure the spatially-resolved JPD of
photon pairs $\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$ by identifying photon
coincidences between any arbitrary pair of pixels $\mathbf{r_{1}}$ and
$\mathbf{r_{2}}$ using a method described in Defienne  _et al._ (2018). In
practice, the information contained in the JPD can be used for various
purposes. For example, it was used in pioneering works with EMCCD cameras to
estimate position and momentum correlation widths of entangled photons pairs
in direct imaging of the Einstein-Podolsky-Rosen paradox Moreau  _et al._
(2012); Edgar  _et al._ (2012). In the experiment here, displaying the
diagonal component of the JPD, $\Gamma(\mathbf{r},\mathbf{r})$, reconstructs
an image of the object. Indeed, assuming an ideal imaging system and near-
perfect spatial correlations between photons Abouraddy  _et al._ (2002), the
JPD can be approximated as
$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})=\left|t(\mathbf{r_{1}})t(\mathbf{r_{2}})\delta(\mathbf{r_{1}}-\mathbf{r_{2}})\right|^{2}$
(where $\delta$ is a Dirac delta function), resulting in a diagonal image
proportional to the object shape i.e.
$\Gamma(\mathbf{r},\mathbf{r})=|t(\mathbf{r})|^{4}$. Such a diagonal image is
the quantity that is conventionally measured and used in all photon-pairs-
based imaging schemes using a near-field illumination configuration Reichert
_et al._ (2017, 2018); Toninelli  _et al._ (2019); Defienne  _et al._ (2019)
(see section 7 of supplementary document).

Figure 3:  Results in quantum illumination imaging. (a) Experimental setup for
quantum illumination in a near-field configuration. Photon pair source and the
EMCCD camera are the same as those used in Figure 1.a. A two-lens system
$f_{5}-f_{6}$ images the crystal surface onto an object $t_{1}$, that is in
turn imaged onto the EMCCD camera by a single-lens imaging system $f_{7}$. An
object $t_{2}$ is inserted and illuminated by a spatially filtered light-
emitting diode (LED) and spectrally filtered at $810\pm 5$ nm. Images of both
objects are superimposed on the camera using a beam splitter (BS). $t_{1}$ and
$t_{2}$ are negative and positive amplitude USAF resolution charts,
respectively. (b) Intensity image. (c) Sum-coordinate projection of the
filtered JPD $P_{+}(\mathbf{r^{+}})$. $6.10^{6}$ frames were acquired to
retrieve intensity image and JPD in approximately $20$ hours of acquisition.
(d) Experimental setup for quantum illumination in a far-field configuration.
Light emitted by a pulsed laser ($347$ nm) illuminates a BBO crystal ($1$ mm
thickness) to produce spatially entangled photon pairs ($694$ nm) by type-I
SPDC. After the crystal, pump photons are filtered out by a BPF at $694\pm 10$
nm. A lens $f_{8}$ is used to Fourier-image photon pairs onto a reflective
object $t_{3}$, that is then imaged onto the SPAD camera by a two-lens system
$f_{9}-f_{10}$. Another laser ($633$ nm) illuminates an object $t_{4}$ that is
also imaged onto the camera using a BS. $t_{3}$ is a cat-shape reflective
object and $t_{4}$ is a transmission object with the same cat-shape,
positioned so that the cat shapes overlap on the sensor. (e) Intensity image.
(f) Minus-coordinate projection of the filtered JPD $P_{-}(\mathbf{r^{-}})$.
$10^{7}$ frames were acquired to retrieve the intensity image and JPD in
approximately $7.5$ hours of acquisition.

Interestingly, there is another way to retrieve an image of the object without
using the diagonal component of the JPD. For that, one can project the JPD
along its sum-coordinate axis, achieved by integrating $\Gamma$:

$P_{+}(\mathbf{r^{+}})=\int\Gamma(\mathbf{r},\mathbf{r^{+}}-\mathbf{r})d\mathbf{r}=\left|t\left(\frac{\mathbf{r^{+}}}{{2}}\right)\right|^{4},$
(1)

where $P_{+}$ is defined as the sum-coordinate projection of the JPD and
$\mathbf{r^{+}}$ is the sum-coordinate variable. Each value
$P_{+}(\mathbf{r^{+}})$ is obtained by adding all the values
$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$ located on an anti-diagonal of the JPD
defined by $\mathbf{r_{1}}+\mathbf{r_{2}}=\mathbf{r^{+}}$. The essential point
in Eq. 1 is the presence of a scaling factor $1/2$ in the spatial variable: in
a discrete space, this indicates that $P_{+}$ is an image of the object
sampled over twice as many points as the number of pixels of the sensor. The
consequence is pixel super-resolution and the ability to retrieve lost spatial
information due to undersampling.

To understand the underlying principle of our approach, it is essential to
realize that the JPD is a 4D object containing much richer information than a
conventional 2D intensity image. In particular, the JPD contains correlation
information not only between photons detected at the same pixel on its main
diagonal $\Gamma(\mathbf{r},\mathbf{r})$, but also correlation information
about photons detected between nearest-neighbour pixels on the 8 other
diagonals directly above and bellow the main diagonal i.e.
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{x}})$,
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{y}})$ and
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{x}}\pm\mathbf{e_{y}})$ (where
$\mathbf{e_{x}}$ and $\mathbf{e_{y}}$ are the unit pixel vectors in transverse
directions $x$ and $y$, respectively). The $9$ total diagonals are equivalent
to $9$ low-resolution images acquired by shifting the sensor by half of a
pixel in the two transverse directions. Using an analogy with the shift-and-
add technique Farsiu  _et al._ (2004), a high-resolution image is then
retrieved from the diagonals by performing the integration in Equation 1 (see
section 2 of the supplementary document).

In practice, calculating directly $P_{+}$ from the the measured JPD results in
an image that is dominated by noise, in which the object is not visible (see
Figure 6.a in Methods). Indeed, experimentally measured JPD correlation values
are noisy, and all the noise from each diagonal adds up in the integration in
Equation 1, ultimately producing a noise of greater variance that dominates
the signal in the final image. To solve this issue, we apply a filtering
technique to the JPD to remove the weakest correlation values that have a
negligible contribution in the final image. By doing so, we also remove the
noise associated with these values and significantly improve the quality of
the final image. This filtering technique does not require prior information
of the object nor the source. After filtering, the calculated sum-coordinate
projection $P_{+}$ is normalised to remove the so-called ’motion error
artefact’ (see Methods). In the case of the experiment in Figure 1.a, by
filtering the JPD, we effectively set to zero all correlation values except
those in the main diagonal and the off-centre diagonals that are kept. Figure
1.c shows the resulting sum-coordinate projection of the measured JPD after
filtering. It reveals the grating amplitude with an increased pixel resolution
compared to the direct intensity image. Figure 1.d shows the spatial frequency
spectrum of the image projected along the x-axis (blue curve) and confirms
that our approach retrieves additional spatial information content (harmonics
at $0.6$ cycles/pixel and $0.8$ cycles/pixel) that cannot be resolved by the
conventional intensity measurement.

Experimental results with different grating periods, ranging from 3 to 1.6
pixels are shown in Figures 2.a-c. These results further confirm that our
technique retrieves spatial information content of an object that is not
visible in the intensity image, even in the case where the first harmonic has
a higher frequency than the Nyquist limit of $0.5$ cycles/pixel (Figure 2.c).
In this latter case, we also observe the presence of an aliased frequency peak
in the intensity spectrum at $0.39$ cycles/pixel (red curve), which is
expected when an image is undersampled. Note also that a similar peak appears
at the same position in the JPD projection image (blue curve), that is however
not predicted by Equation 1. This small residual frequency peak ($25\%$ of the
first harmonic at $0.61$ cycles/pixel) is an artefact originating from
imperfections in the JPD measured with an EMCCD camera Defienne  _et al._
(2018) (see Methods). The experimental results obtained in Figure 1 and 2 are
in very good accordance with simulations (see section 3 of supplementary
document).

Figure 4:  Results in entanglement-enabled quantum holography.(a) Experimental
setup. Light emitted by a laser diode at $405$ nm and polarized at 45∘
illuminates a stacked BBO crystal pair ($0.5$ mm thickness each) whose optical
axes are perpendicular to each other to produce pairs of photons entangled in
space and polarization by type-I SPDC. After the crystals, pump photons are
filtered out by a band pass filter at $810\pm 5$nm. A lens $f_{11}$ is used to
Fourier-image photon pairs onto an optical plane where an SLM with a phase
object $t_{6}$. A two-lens imaging system $f_{12}-f_{13}$ is then used to
image the SLM and the object onto two different parts of an EMCCD camera. A
polarizer (P) at 45∘ is positioned before the camera. (b) Photo of the
birefringent object used in the experiment (section of a bird feather). (c)
Spatial phase of the object reconstructed by combining four anti-diagonal
images $\Gamma(\mathbf{r},-\mathbf{r})$ measured for four phase shifts
$\\{0,\pi/2,\pi/2,3\pi/4\\}$ programmed by the SLM. (d) Spatial phase of the
same object reconstructed from four minus-coordinate projections
$P_{-}(\mathbf{r^{-}})$. $2.5.10^{6}$ frames were acquired to retrieve
intensity image and JPD in each case in approximately $17$ hours of
acquisition. White scale bar is $1$ mm.

Applications to quantum illumination. We demonstrate our technique on two
different experimental schemes based on quantum illumination protocols
described in Defienne  _et al._ (2019) and Defienne  _et al._ (2021b). In the
first implementation (Figure 3.a), an amplitude object $t_{1}$ is illuminated
by photon pairs and imaged onto an EMCCD camera using a near-field imaging
configuration. Simultaneously, another object $t_{2}$ is illuminated by a
classical source and also imaged on the camera. The intensity image, Figure
3.b, shows a superposition of both quantum and classical images. The goal of
this protocol is to segment the quantum object and therefore retrieve an image
showing only object $t_{1}$ illuminated by photon pairs, effectively removing
any classical objects or noise. Using our image processing method, we achieve
such a task and simultaneously retrieve a pixel super-resolved image of
$t_{1}$ (Figure 3.d).

In the second quantum illumination protocol, an object $t_{3}$ (a reflective
amplitude object) is illuminated using a far-field configuration and its image
is measured with a SPAD camera (Figure 3.d). Another object $t_{4}$ (a
transmission object with the same shape as $t_{3}$) is illuminated by a
classical source and also imaged on the SPAD camera. Objects are positioned so
that they overlap exactly on the camera, as shown in the intensity image in
Figure 3.e. In such a far-field configuration, it is important to understand
that information of the JPD is concentrated around the anti-diagonal because
photon pairs are spatially anti-correlated in the object plane Schneeloch and
Howell (2016). Consequently, an image of $t_{3}$ can then be retrieved by
displaying the main anti-diagonal component $\Gamma(\mathbf{r},-\mathbf{r})$.
Such an anti-diagonal image is the quantity that is conventionally measured
and used in all photon-pairs-based imaging schemes using a far-field
illumination configuration Gregory  _et al._ (2020); Defienne  _et al._
(2021a, b) (see section 7 of supplementary document). In this case, the JPD
pixel super-resolution technique must be adapted by using the minus-coordinate
projection $P_{-}$ of the filtered JPD in place of the sum-coordinate
projection $P_{+}$ to retrieve the high-resolution image (Figure 3.f). Even
though this image has a better pixel resolution, we also observe that it has a
lower signal-to-noise ratio (SNR) and the presence some horizontal dark stripe
artefacts. These effects are due to the poor quantum efficiency ($2.6\%$) and
fill-factor ($10.5\%$) of the SPAD camera. In addition, note that the object
appears twice in the image (the vertical cat-shape and its rotated symmetric
copy), which is a known consequence of using a far-field illumination
configuration Defienne  _et al._ (2021b, a) (see Methods).

Application to quantum holography. We now demonstrate our technique on two
different experimental quantum holography schemes. The first protocol under
consideration was demonstrated in Defienne  _et al._ (2021a) and its
experimental configuration is described in Figure 4.a. Pairs of photons
entangled in space and polarisation illuminate an SLM and a birefringent phase
object $t_{6}$ positioned in each half of an optical plane using a far-field
illumination configuration. The object and SLM are then both imaged onto two
different parts of an EMCCD camera. The object is a section of a bird feather
(Fig. 4.b). The SLM is used to compensate for optical aberrations and
implement a four phase-shifting holographic process by displaying uniform
phase patterns with values $\\{0,\pi/2,\pi,3\pi/2\\}$ (see Methods). In the
original protocol, four different images are obtained for each step of the
process by measuring the anti-diagonal component
$\Gamma(\mathbf{r},-\mathbf{r})$ of the JPD. These images are then combined
numerically to reconstruct the spatial phase of the birefringent object
(Figure 4.c). Using our new JPD pixel super-resolution approach, the four
images are now obtained from the minus-coordinate projection of the JPD and
recombined to retrieve a phase image with improved spatial resolution (Figure
4.d).

The second protocol is a full-field version of a N00N-state holographic scheme
based on photon pairs (N=2) Ono  _et al._ (2013). N00N states provide phase
measurements with N-times better sensitivity than classical holography Israel
_et al._ (2014). As shown in Figure 5, photon pairs entangled in space and
polarisation pass through an SLM and illuminate a birefringent object $t_{7}$.
In this case, we use a near-field configuration and the object, SLM and EMCCD
camera are all positioned in conjugate image planes. In our experiment,
$t_{7}$ is a checkerboard birefringent phase pattern composed of 9 square
areas with different phase values (Figure 5.b). The SLM is used to perform a
phase-shifting holographic process at four phase values
$\\{0,\pi/4,\pi/2,3\pi/4\\}$. Sum-coordinate projections of four filtered JPDs
are measured at each step of the process and then combined to reconstruct the
object phase image (see Methods) shown in Figure 5.c. This image has better
pixel resolution than the phase image of the same object retrieved using an
equivalent classical holographic system, shown in Figure 5.d. Furthermore, it
is worth noting that phase values in the image retrieved using photon pairs
are twice those of the object and the classical phase image, as expected for
N00N-state holography and also as demonstrated in previous scanning-based
imaging approach with N00N states Ono  _et al._ (2013); Israel  _et al._
(2014). In addition to preserving this sensitivity enhancement, our full-field
version of N00N-state holography then also provides a better pixel resolution,
an advantage that could only be matched with the scanning version by
increasing the acquisition time by at least a factor of four. The range of
applications can be extended to non-birefringent phase objects by using
birefringent plates, as done in differential interference contrast microscopy
Allen and David (1969); Terborg  _et al._ (2016).

Conclusion. We have introduced a pixel super-resolution technique based on the
measurement of a spatially-resolved JPD of spatially-entangled photon pairs.
This approach retrieves spatial information about an object that is lost due
to pixelation, without shifting optical elements or relying on prior
knowledge. We demonstrated that this JPD pixel super-resolution approach can
improve the spatial resolution in already established quantum illumination and
entanglement-enabled holography schemes, as well as in a full-field version of
N00N-state quantum holography, using different types of illumination (near-
field and far-field) and different single-photon-sensitive cameras (EMCCD and
SPAD). Our approach has the advantage that it can be used immediately in
quantum imaging schemes based on photon pairs, and even in some cases by only
reprocessing already acquired data. In addition, our approach can also be
implemented in any classical imaging system limited by pixelation, after
substituting the classical source by a source of entangled photon pairs with
similar properties. Indeed, contrary to classical pixel super-resolution
techniques, such as shift-and-add approaches Farsiu  _et al._ (2004),
wavelength scanning Luo  _et al._ (2016) and machine-learning-based algorithms
Shi  _et al._ (2016), our JPD super-resolution approach has the strong
advantage that it does not require displacing optical elements in the system
or having prior knowledge about the object being imaged. Our JPD pixel super-
resolution approach can therefore benefit all quantum and classical imaging
applications in which the spatial resolution is limited by the camera pixel
size.

Figure 5: Results in full-field N00N-state quantum holography. (a)
Experimental setup. Pairs of photons entangled in space and polarization are
produced by type I SPDC using the same source of Figure 4.a. Three consecutive
two-lens systems $f_{14}-f_{15}$, $f_{17}-f_{18}$ and $f_{19}-f_{20}$ are used
to successively image the crystal surface onto a SLM, an object $t_{7}$ and
onto an EMCCD camera. A polarizer (P) at 45∘ is positioned before the camera.
(b) The object $t_{7}$ is a checkerboard birefringent phase pattern composed
of 9 squared areas with different phase values uniformly distributed in
$[0,2\pi]$ created using another SLM. (c) Spatial phase image obtained by
combining four sum-coordinate projections $P_{+}(\mathbf{r^{+}})$ measured for
four different phase shifts $\\{0,\pi/4,\pi/2,3\pi/4\\}$ programmed by the
SLM. (d) Spatial phase of the same object reconstructed by classical
holography after replacing the source of photon pairs by a laser at $810$ nm
and $45^{\circ}$ polarisation and measuring four intensity images measured for
four different phase shifts $\\{0,\pi/2,\pi,3\pi/2\\}$. $2.10^{7}$ frames were
acquired to retrieve the JPD in each case in approximately $30$ hours of
acquisition.

### Methods

Experimental layouts. Experiment in Figure 1.a: BBO crystal has dimensions of
$0.5\times 5\times 5$ mm and is cut for type I SPDC at $405$ nm with a half
opening angle of $3$ degrees (Newlight Photonics). It is slightly rotated
around horizontal axis to ensure near-collinear phase matching of photons at
the output (i.e. ring collapsed into a disk). The pump is a continuous-wave
laser at $405$ nm (Coherent OBIS-LX) with an output power of approximately
$200$ mW and a beam diameter of $0.8\pm 0.1$ mm. A $650$ nm-cut-off long-pass
filter is used to block pump photons after the crystals, together with a band-
pass filter centred at $810\pm 5$ nm. The camera is an EMCCD (Andor Ixon Ultra
897) that operates at $-60^{\circ}$C, with a horizontal pixel shift readout
rate of $17$ MHz, a vertical pixel shift every $0.3$ $\mu$s, a vertical clock
amplitude voltage of $+4$V above the factory setting and an amplification gain
set to $1000$. It has a $16\mu$m pixel pitch. Exposure time is set to $2$ ms.
The camera speed is about $100$ frames per second (fps) using a region of
interest of $100\times 100$ pixels. The two-lens imaging system $f_{1}-f_{2}$
in Fig. 1.a is represented by two lenses for clarity, but is composed of a
series of $6$ lenses with focal lengths $45$ mm - $75$ mm - $50$ mm - $150$ mm
- $100$ mm - $150$ mm. The first and the last lens are positioned at focal
distances from the crystal and the object, respectively, and the distance
between two lenses in a row equals the sum of their focal lengths. Similarly,
the second two-lens imaging system $f_{3}-f_{4}$ in Figure 1.a is composed of
a series of $4$ consecutive lenses with focal lengths $150$mm - $50$ mm - $75$
mm - $100$ mm arranged as in the previous case.
Experiment in Figure 3.a: Laser, camera and photon-pairs source are the same
than in Figure 1.a. Exposure time is set to $6$ ms. The two-lens imaging
system $f_{5}-f_{6}$ is composed of two lenses with focal lengths $35$ mm and
$75$ mm. The single-lens imaging system $f_{7}$ is composed of one lens with
focal length $50$ mm positioned at a distance of $100$mm from both objects and
the camera. The same setup was used in Defienne  _et al._ (2019).
Experiment in Figure 3.d: BBO crystal has dimensions $5\times 5\times 1$mm cut
for type-I SPDC at $355$nm with a half opening angle of $3$ degrees (Newlight
Photonics). It is slightly rotated around horizontal axis to ensure near-
collinear phase matching of photons at the output (i.e. ring collapsed into a
disk). The pump is the third harmonic at $347$ nm of a femtosecond pulsed
laser with $100$ MHz repetition rate, $80$ mW average power and beam-diameter
of approximately $0.5$ mm (Chromacity). The SPAD camera is the SwissSPAD2 Ulku
_et al._ (2019). It has a total of $512\times 512$ pixels, a fill-factor of
$10.5\%$, a photon detection efficiency of $2.6\%$, a pixel pitch of $16.83$
$\mu$m and a frame rate of $370$ fps. The Fourier-imaging system $f_{8}$ is
composed of one lens with focal length $45$ mm. The two-lens imaging system
$f_{9}-f_{10}$ is composed of two lenses with focal lengths $100$ mm and $50$
mm. The same setup was used in Defienne  _et al._ (2021b).
Experiment in Figure 4.a: The paired set of stacked BBO crystals have
dimensions of $0.5\times 5\times 5$ mm each and are cut for type I SPDC at
$405$ nm. They are optically contacted with one crystal rotated by $90$
degrees about the axis normal to the incidence face. Both crystals are
slightly rotated around horizontal and vertical axis to ensure near-collinear
phase matching of photons at the output (i.e. rings collapsed into disks). The
pump laser and camera are the same than in Figure 1.a. A $650$ nm-cut-off
long-pass filter is used to block pump photons after the crystals, together
with a band-pass filter centred at $810\pm 5$ nm. The SLM is a phase only
modulator (Holoeye Pluto-2-NIR-015) with $1920\times 1080$ pixels and a $8$
$\mu$m pixel pitch. For clarity, it is represented in transmission in Figure
4.a, but was operated in reflection. Exposure time is set to $3$ ms. The
camera speed is about $40$ frame per second using a region of interest of
$200\times 200$ pixels. The single-lens Fourier imaging system $f_{11}$ is
composed of a series of three lenses of focal lengths $45$ mm - $125$ mm -
$150$ mm. The first and last lenses are positioned at focal distance from the
crystal and the object-SLM plane, respectively, and the distance between each
pair of lenses equals the sum of their focal lengths. The two-lens imaging
system $f_{12}-f_{13}$ is in reality composed by a series of $4$ lenses with
focal lengths $150$ mm - $75$ mm - $75$ mm - $100$ mm. The first and the last
lens are positioned at focal distances from respectively the SLM-object and
the camera, respectively, and the distance between two lenses in a row equals
the sum of their focal lengths. The same setup was used in Defienne  _et al._
(2021a) and in section 5 of the supplementary document.
Experiment in Figure 5.a: Laser, camera, SLM, filters and photon-pairs source
are the same than in Figure 4.a. The two-lens imaging system $f_{14}-f_{15}$
is composed of a series of $4$ lenses with focal lengths $45$ mm - $75$ mm -
$50$ mm - $150$ mm. They are arranged so that the first and last lenses are
positioned at focal length from the crystal and the SLM, respectively, and the
distance between each pair of lenses equals the sum of their focal lengths.
The two-lens imaging system $f_{16}-f_{17}$ is composed of two lenses with
focal lengths $100$ mm and $150$ mm. Similarly, the second two-lens imaging
system $f_{18}-f_{19}$ in Figure 1.a is in reality composed of a series of $4$
lenses with focal lengths $150$mm - $50$ mm - $75$ mm - $100$ mm, arranged as
for $f_{14}-f_{15}$. More details in section 6 of the supplementary document.

JPD measurement with a camera. $\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$, where
$\mathbf{r_{1}}$ and $\mathbf{r_{2}}$ are positions of camera pixels, is
measured by acquiring a set of $N+1$ frames
$\\{I_{l}\\}_{l\in[\\![1,N+1]\\!]}$ using a fixed exposure time and then
processing them using the formula:

$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})=\frac{1}{N}\sum_{l=1}^{N}\left[I_{l}(\mathbf{r_{1}})I_{l}(\mathbf{r_{2}})-I_{l}(\mathbf{r_{1}})I_{l+1}(\mathbf{r_{2}})\right]$
(2)

In all the results of our work, $N$ was on the order of $10^{6}-10^{7}$
frames. However, it is essential to note that not all the JPD values can be
directly measured using this process. When using an EMCCD camera, (i)
correlation values at the same pixel, i.e. $\Gamma(\mathbf{r},\mathbf{r})$,
cannot be directly measured because Equation 2 is only valid for different
pixels $\mathbf{r_{1}}\neq\mathbf{r_{2}}$ Defienne  _et al._ (2018), and (ii)
and correlation values between vertical pixels, i.e.
$\Gamma(\mathbf{r},\mathbf{r}\pm k_{y}\mathbf{e_{y}})$ (where $k_{y}$ is an
integer that defines the position of a pixel above or bellow $\mathbf{r}$),
cannot be measured because of the presence of charge smearing effects. To
circumvent this issue in practice, these values are interpolated from
neighbouring correlation values of the JPD i.e.
$[\Gamma(\mathbf{r},\mathbf{r}+\mathbf{e_{x}})+\Gamma(\mathbf{r},\mathbf{r}-\mathbf{e_{x}})]/2\rightarrow\Gamma(\mathbf{r},\mathbf{r})$
and $[\Gamma(\mathbf{r},\mathbf{r}\pm
k_{y}\mathbf{e_{y}}+\mathbf{e_{x}})+\Gamma(\mathbf{r},\mathbf{r}\pm
k_{y}\mathbf{e_{y}}-\mathbf{e_{x}})]/2\rightarrow\Gamma(\mathbf{r},\mathbf{r}\pm
k_{y}\mathbf{e_{y}})$, as detailed in Reichert  _et al._ (2018). As a result,
the gain in resolution along the $x$-axis in the experiments using near-field
imaging configurations (Figure 1.a, 3.a, 6.a) is not optimal. However, it is
important to note that the gain in pixel resolution along the $y$-axis is not
affected by this interpolation technique. Therefore, the spectral analyses
performed in Figures 1 and 2 are also not impacted because the grid-objects
used are horizontal (no spectral component on the $x$-axis) and all the
resulting spectral curves are obtained by summing along the $x$-axis. In
addition, this interpolation also does not affect experiments using far-field
illumination configurations in Figures 3.d and 5.a because the JPD diagonals
do not contain any relevant imaging information (that is in the JPD anti-
diagonals). More details are provided in Defienne  _et al._ (2018) and in
section 1 of the supplementary document.

JPD pixel super-resolution in far-field illumination configuration. When an
object is illuminated by photon-pairs using a far-field illumination
configuration (i.e. the crystal is Fourier-imaged on the object), the JPD
super-resolution technique must be adapted and the sum-coordinate projection
$P_{+}$ cannot be used to retrieve the high-resolution image. Indeed, assuming
an ideal imaging system and near-perfect spatial correlations between photon
pairs Abouraddy  _et al._ (2002), the JPD measured by the camera using a far-
field illumination configuration can be written as:

$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})=\left|t(\mathbf{r_{1}})t(\mathbf{r_{2}})\delta(\mathbf{r_{1}}+\mathbf{r_{2}})\right|^{2},$
(3)

where $\delta$ is a Dirac delta function. Information about the object can
then be retrieved by displaying the anti-diagonal component of the JPD i.e.
$\Gamma(\mathbf{r},-\mathbf{r})=|t(\mathbf{r})|^{2}|t(-\mathbf{r})|^{2}$.
$\Gamma(\mathbf{r},-\mathbf{r})$ is the quantity that is conventionally
measured and used in all photon-pairs-based imaging schemes using a far-field
illumination configuration Gregory  _et al._ (2020); Defienne  _et al._
(2021b, a). Interestingly, it is also possible to retrieve information about
the object by integrating the JPD along its minus-coordinate axis:

$P_{-}(\mathbf{r^{-}})=\int\Gamma(\mathbf{r},\mathbf{r}+\mathbf{r^{-}})d\mathbf{r}=\left|t\left(\frac{\mathbf{r^{-}}}{2}\right)t\left(\frac{-\mathbf{r^{-}}}{2}\right)\right|^{2},$
(4)

where $P_{-}$ is defined as the minus-coordinate projection of the JPD and
$\mathbf{r^{-}}$ is the minus-coordinate variable. Each value
$P_{-}(\mathbf{r^{-}})$ is obtained by adding all the values
$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$ located on a diagonal of the JPD
defined by $\mathbf{r_{1}}-\mathbf{r_{2}}=\mathbf{r^{-}}$. An essential point
in Eq. 4 is the presence of a scaling factor $1/2$ in the spatial variable. In
a discrete space, it means that $P_{-}$ is an image of the object sampled over
twice as many points as the number of pixels of the sensor. In theory,
calculating the minus-coordinate projection of the JPD can therefore achieve
pixel super-resolution and potentially retrieved lost spatial information of
undersampled objects. However, it is important to note that the JPD anti-
diagonal $\Gamma(\mathbf{r},-\mathbf{r})$ and minus-coordinate projection
$P_{-}$ images are not directly proportional to $|t(\mathbf{r})|$ as in the
near-field configuration, but to $|t(\mathbf{r})||t(-\mathbf{r})|$, which does
not always enable to retrieve $|t|$ without ambiguity. In works using a far-
field illumination configuration, this problem is solved by illuminating $t$
with only half of the photon pairs beam (i.e. $t(x,y\leq 0)=1$). The object
then appears in double in the retrieved image (object and its symmetric), but
no information is lost.

JPD filtering. Figure 6.a shows the sum-coordinate projection $P_{-}$
calculated using equation 1 from the unfiltered JPD measured in Figure 1.a.
This image is very noisy and does not reveal the object. To solve this issue,
a filtering method is applied and consists in setting to $0$ all values of the
JPD that have a negligible weight (i.e. values close to zero) so that their
associated noise does not contribute when performing the integration. To
perform this selection efficiently, correlation values with the lowest weights
are not selected individually but by groups. To distinguish these groups, we
must differentiate the two cases of an object illuminated by photon pairs
using (i) a near-field illumination configuration (Figures 1.a, 3.a and 5.a)
and (ii) a far-field configuration (Figures 3.d and 4.a).
Near-field configuration: Using a near-field illumination configuration, each
group corresponds to all the correlation values located in a specific JPD
diagonal $\Gamma(\mathbf{r},\mathbf{r}+\mathbf{r^{-}})$, where
$\mathbf{r^{-}}$ defines the position of the diagonal. To assess the weights
of the different JPD diagonals, the minus-coordinate projection of the JPD
(noted $P_{-}$) is calculated using the integration defined in Eq. 4. Each
value $P_{-}(\mathbf{r^{-}})$ is the sum of all the correlation values
$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}}+\mathbf{r^{-}})$ located on the diagonal
defined by $\mathbf{r_{1}}-\mathbf{r_{2}}=\mathbf{r^{-}}$ and therefore
provides a direct estimation of the corresponding diagonal weight. For
example, Fig. 6.b shows the minus-coordinate projection obtained using the JPD
measured in the experiment of Fig. 1.a. We observe that only a small number of
diagonals (values around the center in $P_{-}$) have significant weights,
while all the others are almost zero. By applying a threshold, only the
highest values of $P_{-}$ are kept, which in the case of the experiment in
Fig. 1.a corresponds to keeping values of the main JPD diagonal
$\Gamma(\mathbf{r},\mathbf{r})$ (i.e. center pixel of $P_{-}$) and those of
the $8$ other diagonals directly above and bellow i.e.
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{x}})$,
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{y}})$ and
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{x}}\pm\mathbf{e_{y}})$ (direct-
neighbour around the center pixel in $P_{-}$). After filtering the JPD, the
integration in Eq. 1 is performed and a normalisation process is used (see
next Methods section) to retrieve the low-noise high-resolution image shown in
Fig. 1.c.
Far-field configuration: Using a far-field illumination configuration, the
process is the same as in the near-field configuration described above after
inverting the roles played by the minus-coordinate $P_{-}$ and sum-coordinate
$P_{+}$ projections. In this configuration, the JPD sum-coordinate projection
$P_{+}$ enables to assess the weights of the different group of correlations
values associated with specific JPD anti-diagonals. For example, Figure 6.c
shows $P_{+}$ measured using the experiment in Figure 4.a, in which a peak of
correlation is well visible. By applying a threshold on this image, only the
highest values of $P_{+}$ are kept (i.e. those in the correlation peak), which
corresponds to keeping values of the main JPD anti-diagonal
$\Gamma(\mathbf{r},-\mathbf{r})$ (i.e. pixel at the center of the peak) and
those of the $8$ other anti-diagonals directly above and bellow i.e.
$\Gamma(\mathbf{r},-\mathbf{r}\pm\mathbf{e_{x}})$,
$\Gamma(\mathbf{r},-\mathbf{r}\pm\mathbf{e_{y}})$ and
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{x}}\pm\mathbf{e_{y}})$ (i.e. pixel
surrounding the center of the peak). After filtering the JPD, the integration
in Eq. 4 is performed and a normalisation process is used (see next Methods
section) to retrieve a low-noise high-resolution image.
Symmetric roles of $P_{+}$ and $P_{-}$: We underline that both the sum-
coordinate $P_{+}$ and minus-coordinate $P_{-}$ projections of the JPD are
used in both near-field and far-field configurations, but their roles are
different. In a near-field configuration, $P_{-}$ is first calculated from the
measured JPD to apply the filtering process, and then $P_{+}$ is calculated to
retrieve the high-resolution image from the filtered JPD. In a far-field
configuration, $P_{+}$ is first calculated from the measured JPD to apply the
filtering process, and then $P_{-}$ is calculated to retrieve the high-
resolution image from the filtered JPD. More details are provided section 4 of
the supplementary document.

Figure 6:  (a) Sum-coordinate projection of the JPD $P_{+}(\mathbf{r^{+}})$
measured in the experiment in Figure 1.a without filtering. (b) Minus-
coordinate projection image of the JPD $P_{-}(\mathbf{r^{-}})$ measured in the
experiment in Figure 1.a without filtering. (c) Sum-coordinate projection of
the JPD $P_{+}(\mathbf{r^{+}})$ measured in the experiment in Figure 4.a
without filtering. (d) Sum-coordinate projection of the JPD
$P_{+}(\mathbf{r^{+}})$ measured in the experiment in Figure 1.a after
filtering and before normalisation.

Normalisation. Figure 6.d shows the sum-coordinate projection $P_{+}$ directly
calculated from the filtered JPD measured in the experiment in Figure 1.a. We
observe that this image has an artefact taking the form of inhomogeneous
horizontal and vertical stripes. This artefact is very similar to this
commonly observed in shift-and-add super-resolution techniques Farsiu  _et
al._ (2004) and is often refereed as a ’motion error artefact’. In our case,
it originates from the approximation made when substituting the integration in
Equation 1 by a summation in a discrete space to calculate $P_{+}$, and also
from the JPD filtering. Using an analogy with the shift-and-add technique, our
problem is equivalent to a situation in which the different shifted low-
resolution images were measured using sources with different averaged
intensity during the first step of the process. Then, when these low-
resolution images are recombined into a high-resolution one (second step of
shift-and-add), the artefact appears in the resulting image because some
pixels placed next to each others were not acquired using the same intensity.
In practice, the most basic way to remove this artefact is to normalise each
low-resolution image by its total average intensity, so that neighbouring
pixels in the high-resolution image are at the same level. We use such a
normalisation approach in our work by dividing all values of the non-zero
diagonals in the filtered JPD by their spatial average value i.e.
$\Gamma(\mathbf{r},\mathbf{r}+\mathbf{r^{-}})\rightarrow\Gamma(\mathbf{r},\mathbf{r}+\mathbf{r^{-}})/\sum_{\mathbf{r}}\Gamma(\mathbf{r},\mathbf{r}+\mathbf{r^{-}})$,
where $\mathbf{r^{-}}$ identifies a specific JPD diagonal. After
normalisation, Figure 1.c is obtained. In the case of far-field illumination,
the same normalisation is applied to the values of the non-zero anti-diagonals
in the filtered JPD i.e.
$\Gamma(\mathbf{r},-\mathbf{r}+\mathbf{r^{+}})\rightarrow\Gamma(\mathbf{r},-\mathbf{r}+\mathbf{r^{+}})/\sum_{\mathbf{r}}\Gamma(\mathbf{r},-\mathbf{r}+\mathbf{r^{+}})$,
where $\mathbf{r^{+}}$ identifies a specific JPD anti-diagonal.
In some cases, such as in Fig. 4.f, the artefact is reduced but still visible
in the resulting image even after normalisation. The persistence of this
artefact is due the fact that the anti-diagonals are not only different in
terms of their average intensity (that can be corrected by normalisation), but
they are also different in terms of their spatial shape. By analogy to shift-
and-add, a similar artefact would be present if the different shifted low-
resolution images were measured using sources of light with different
distribution of spatial intensity. In the experiment shown in Figure 4.d, the
difference in spatial shape of the JPD anti-diagonals originates from the poor
quality of the SPAD camera sensor, in particular from its very low fill-factor
($10.5\%$), that can create spatial inhomogeneities in the coincidence values
measured for different pair of pixels depending on how the sensor is
positioned relative to the source and the object. To further reduce this
artefact, one could use more complex normalisation algorithms, such as $L_{1}$
or $L_{2}$ norms minimisation approaches Farsiu  _et al._ (2006) and kernel
regression Takeda  _et al._ (2007), that are commonly used in shift-and-add.

Full-field N00N-state quantum holography. In the holographic process used in
full-field N00N-state holography (Fig. 5.a), different JPDs are measured for
different uniform phase values programmed on the SLM. For each phase value
$\alpha$, the sum-coordinate projection of the JPD (noted $P^{\alpha}_{+}$) is
proportional to $1+\cos(2\Delta\theta-2\alpha)$, where
$\Delta\theta=\theta_{H}-\theta_{V}$ with $\theta_{H}$ (resp. $\theta_{V}$)
the horizontal (resp. vertical) phase component of $t_{7}$. To reconstruct the
phase image (Figure 5.c), JPD projections measured for four phases
$\alpha\in\\{0,\pi/4,\pi/2,3\pi/4\\}$ are recombined using the formula:
$\mbox{arg}[P^{0}_{+}-P^{\pi/2}_{+}+i(P^{\pi/4}_{+}-P^{3\pi/4}_{+})]=2\Delta\theta$.
To retrieve the classical phase image (Fig. 5.d), the photon pair source is
first substituted by a collimated laser beam at $810$nm and polarized at
$45^{\circ}$. Then, as conventionally done in quantitative phase imaging
Yamaguchi and Zhang (1997), four intensity images (noted $I^{\alpha}$) are
measured for four different constant phase value programmed on the SLM
$\alpha\in\\{0,\pi/2,\pi,3\pi/2\\}$ and combined using the formula
$\mbox{arg}[I^{0}-I^{\pi}+i(I^{\pi/2}-I^{3\pi/2})]=\Delta\theta$ to retrieve
the spatial phase of the object $\theta$. More details are provided section 6
of the supplementary document, including images of JPD projections measured
for each phase value and a description of the equivalent classical holography
experiment used to retrieved the phase image shown in Fig. 5.d

Residual aliased peak in the JPD projection. As detailed in Defienne  _et al._
(2018), the use of Eq. 2 to measure the JPD with an EMCCD camera is not
perfect and a residual background $\mathcal{B}$ is also retrieved i.e.
$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})+\mathcal{B}$. This background
originates from the fluctuations of the EMCCD camera gain during the
acquisition process. $\mathcal{B}$ is proportional to $\alpha
I(\mathbf{r_{1}})I(\mathbf{r_{2}})+\beta(I(\mathbf{r_{1}})+I(\mathbf{r_{2}}))+\gamma$,
where $\alpha$, $\beta$ and $\gamma$ are constant experimental parameters, and
$I(\mathbf{r})$ is the intensity value at pixel $\mathbf{r}$. In most
experiments, $\mathcal{B}$ is negligible compared to $\Gamma$ and is just
ignored. However, when performing the experiment with a high-frequency grid-
shaped object of period $1.6$ pixels (Fig. 2.c), most of the information in
$\Gamma$ is lost (all harmonics $>1$ cycles/pixel are cut in the blue curve)
while some low frequency content of the intensity image (red curve), including
in particular the aliased peak at $0.39$ cycles/pixel, are still retrieved.
These asymmetric losses level the ratio between $\Gamma$ and $\mathcal{B}$ in
favour of $\mathcal{B}$, making the aliased peak at $0.39$ cycles/pixel in the
blue curve (due to $\mathcal{B}$) a bit more visible relative to this at
$0.61$ cycles/pixel (due to $\Gamma$). More details are provided section 1 of
the supplementary document.

Funding. D.F. acknowledges financial support from the Royal Academy of
Engineering Chair in Emerging Technology, UK Engineering and Physical Sciences
Research Council (grants EP/T00097X/1 and EP/R030081/1) and from the European
Union’s Horizon 2020 research and innovation programme under grant agreement
No 801060. H.D. acknowledges funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No. 840958.

Contributions. H.D conceived the original idea, designed and performed
experiments, analysed the data and prepared the manuscript. P.C. and H.D.
performed the simulations. M.R. and J.F. contributed to developing the quantum
illumination protocol with the EMCCD camera. J.Z and E.C. contributed to
developing the quantum illumination protocol with the SPAD camera. A.L. and
B.N. contributed to developing the entanglement-enabled quantum holography
protocol. All authors contributed to the manuscript. D.F. supervised the
project.

Competing interests. The authors declare no competing interests.

Data availability. All data that support the plots within this paper and other
findings of this study are available from the corresponding author upon
reasonable request.

## Supplementary document

###  I Details on the joint probability distribution (JPD) measurement

This section provides more details about the measurement of the spatially
resolved JPD of photon pairs using a Electron Multiplied Charge Coupled Device
(EMCCD) or single-photon avalanche diode (SPAD) camera. Further theoretical
details can be found in Defienne  _et al._ (2018), Reichert  _et al._ (2018)
and Defienne  _et al._ (2021b).

####  I.1 JPD measurement.

EMCCD ans SPAD cameras can be used to reconstruct the (a) intensity
distribution $I(\mathbf{r})$ and (b) the spatial joint probability
distribution $\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$ of photon pairs, where
$\mathbf{r}$, $\mathbf{r_{1}}$ and $\mathbf{r_{2}}$ are positions of pixels.
To do that, the camera first acquires a set of $N+1$ frames
$\\{I_{l}\\}_{l\in[\\![1,N+1]\\!]}$ using a fixed exposure time and without
using a threshold in the case of EMCCD camera. Then:

  1. (a)

The intensity image is reconstructed by averaging over all the frames:

$I(\mathbf{r})=\frac{1}{N}\sum_{l=1}^{N}I_{l}(\mathbf{r})$ (5)

  2. (b)

The JPD is reconstructed by performing the following substation:

$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})=\frac{1}{N}\sum_{l=1}^{N}I_{l}(\mathbf{r_{1}})I_{l}(\mathbf{r_{2}})-\frac{1}{N}\sum_{l=1}^{N}I_{l}(\mathbf{r_{1}})I_{l+1}(\mathbf{r_{2}})$
(6)

Under illumination by the photon pairs, intensity correlations in the left
term of Equation 6 originate from detections of both real coincidences (two
photons from the same entangled pair) and accidental coincidences (two photons
from two different entangled pairs), while intensity correlations in the
second term originate only from photons from different entangled pairs
(accidental coincidence) because there is zero probability for two photons
from the same entangled pair to be detected in two successive images. A
subtraction between these two terms leaves only genuine coincidences, that are
proportional to the spatial JPD of the pairs.

####  I.2 Interpolation of some JPD values measured with an EMCCD camera.

Using an EMCCD camera, not all correlation values of the JPD can be measured
using Equation 6, in particular those (a) at the same pixel i.e.
$\Gamma(\mathbf{r},\mathbf{r})$ and (b) and those between vertical
neighbouring pixels i.e. $\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{y}})$
(where $\mathbf{e_{y}}$ is the unit pixel vector in transverse direction $y$):

  1. (a)

As demonstrated in theory in Defienne  _et al._ (2018), Equation 6 is not
valid for $\mathbf{r_{1}}=\mathbf{r_{2}}$. For example, Figure 7.a shows a
conditional image $\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$ for a reference
pixel $\mathbf{r_{1}}=(-5,3)$ obtained from the JPD measured in the experiment
described in Figure 1.a of the manuscript using Equation 6. In this image, we
observe that the value at $\mathbf{r_{2}}=\mathbf{r_{1}}$ is more than
$10^{3}$ times brighter than the other images values, confirming that this
value is wrong. In practice, this issue is solved by interpolating all the
values in the diagonal $\Gamma(\mathbf{r},\mathbf{r})$ using neighbouring
correlation values on the same row:
$[\Gamma(\mathbf{r},\mathbf{r}+\mathbf{e_{x}})+\Gamma(\mathbf{r},\mathbf{r}-\mathbf{e_{x}})]/2\rightarrow\Gamma(\mathbf{r},\mathbf{r})$
($\mathbf{e_{x}}$ is the unit pixel vector in transverse direction $x$).
Figure 7.b shows the conditional image $\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$
for $\mathbf{r_{1}}=(-5,3)$ obtained after interpolating the pixel value at
$\mathbf{r_{2}}=\mathbf{r_{1}}$. More details in Reichert  _et al._ (2018).

  2. (b)

However, we still observe in Figure 7.b the presence of artefacts in the
values at positions $\mathbf{r_{2}}=\mathbf{r_{1}}\pm k_{y}\mathbf{e_{y}}$,
where $k_{y}$ is an integer that defines the position of a pixel above or
bellow $\mathbf{r_{1}}$. These artefacts are due to cross-talk between pixels
located on the same column of the EMCCD camera, originating from the so-called
charge smearing effect Reichert  _et al._ (2018). In practice, these values
are interpolated using neighbouring correlation values on the same row:
$[\Gamma(\mathbf{r},\mathbf{r}+\mathbf{e_{x}})+\Gamma(\mathbf{r},\mathbf{r}-\mathbf{e_{x}})]/2\rightarrow\Gamma(\mathbf{r},\mathbf{r})$
and $[\Gamma(\mathbf{r},\mathbf{r}\pm
k_{y}\mathbf{e_{y}}+\mathbf{e_{x}})+\Gamma(\mathbf{r},\mathbf{r}\pm
k_{y}\mathbf{e_{y}}-\mathbf{e_{x}})]/2\rightarrow\Gamma(\mathbf{r},\mathbf{r}\pm
k_{y}\mathbf{e_{y}})$. Figure 7.c shows the resulting conditional image
$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$ for $\mathbf{r_{1}}=(-5,3)$ after
interpolating the pixel values.

Figure 7:  Interpolation of JPD values when using EMCCD camera. (a)
Conditional image $\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$ for a reference
pixel $\mathbf{r_{1}}=(-5,3)$ directly obtained from the JPD measured in the
experiment described in Figure 1.a without interpolating the pixel values. (b)
Same conditional image obtained after interpolating the diagonal values
$\Gamma(\mathbf{r},\mathbf{r})$ of the measured JPD. (c) Same conditional
image obtained after interpolating all the values
$\Gamma(\mathbf{r},\mathbf{r}\pm k_{y}\mathbf{e_{y}})$ of the measured JPD,
where $k_{y}$ is an integer that correspond to an arbitrary number of vertical
pixel above or bellow $\mathbf{r}$, and $\mathbf{e_{y}}$ is the unit vector
pixel in the $y$-axis.

####  I.3 Interpolation of some JPD values measured with an EMCCD camera.

Similar issues are encountered when measuring the JPD with a SPAD camera, as
shown in Figure 3.d of the manuscript. (a) Equation 6 is not valid for
$\mathbf{r_{1}}=\mathbf{r_{2}}$ and therefore cannot be used to measure the
main diagonal $\Gamma(\mathbf{r},\mathbf{r})$. (b) The presence of cross-talk
between neighbouring pixels of the SPAD camera also corrupt the values in the
$8$ other diagonals
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{x}}\pm\mathbf{e_{y}})$. In
practice, this problem can be solved by setting all these values to $0$ and
using another configuration than a near-field illumination configuration to
perform imaging. In Figure 3.d of the manuscript, we did that by using a far-
field illumination so that the image information content is located in the
anti-diagonal components of the JPD and not on its diagonals. More details in
Defienne  _et al._ (2021b).

####  I.4 Residual term in the measured JPD using EMCCD camera.

The use of Equation 6 to measure the JPD with an EMCCD camera is not perfect
and a residual background $B$ is also retrieved:
$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})+\mathcal{B}$. This problem is studied
in Defienne  _et al._ (2018). This background term originates from the
fluctuations of the EMCCD camera gain during the acquisition process, for
example due to temperature drifts and variations in the high voltage clock
amplitude. $\mathcal{B}$ is a complex term proportional to $\alpha
I(\mathbf{r_{1}})I(\mathbf{r_{2}})+\beta(I(\mathbf{r_{1}})+I(\mathbf{r_{2}}))+\gamma$,
where $I(\mathbf{r})$ is the intensity value at pixel $\mathbf{r}$, and
$\alpha$, $\beta$ and $\gamma$ are constant experimental parameters depending
on the mean and variance values of the camera gain and electronic noise. In
practice, $\mathcal{B}$ is often negligible compared to $\Gamma$ and can be
ignored in most experiments. However, because $\Gamma$ and $\mathcal{B}$ do
not have the same spatial frequency content, spatial filtering effects can
create asymmetric losses and level these two terms in favor of $\mathcal{B}$.
In particular, it is the case in Figure 2.c of the manuscript where spatial
frequency contents are filtered asymmetrically by pixellation.

###  II JPD projections in a pixelated space.

In the experiments shown in our work, the measured JPDs $\Gamma$ takes the
form of a 4-dimensional matrix containing $(M_{Y}\times M_{X})^{4}$ elements,
where $M_{Y}\times M_{X}$ corresponds the size of the illuminated region of
the camera sensor. An element of the matrix is written
$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$, where
$\mathbf{r_{1}}=x_{1}\mathbf{e_{x}}+y_{1}\mathbf{e_{y}}$ and
$\mathbf{r_{2}}=x_{2}\mathbf{e_{x}}+y_{2}\mathbf{e_{y}}$ are the spatial
coordinates of two pixels. These coordinates are discrete:
${x_{i}}\in[1,M_{X}]$ and ${y_{i}}\in[1,M_{Y}]$ (where $i=1$ or $2$). In our
work, the information content of $\Gamma$ is analysed using two types of
projections:

  1. 1.

In theory, the sum-coordinate projection is defined as:

$P_{+}(\mathbf{r^{+}})=\int\Gamma(\mathbf{r},\mathbf{r^{+}}-\mathbf{r})d\mathbf{r}.$
(7)

In practice, it is calculated using the discrete space formula:

$P_{+}(\mathbf{r^{+}})=\sum_{x=1}^{M_{X}}\sum_{y=1}^{M_{Y}}\Gamma(\mathbf{r},\mathbf{r^{+}}-\mathbf{r}).$
(8)

where $\mathbf{r}=x\mathbf{e_{x}}+y\mathbf{e_{y}}$. It represents the
probability of detecting pairs of photons generated in all symmetric
directions relative to the position
$\mathbf{r^{+}}=x^{+}\mathbf{e_{x}}+y^{+}\mathbf{e_{y}}$.

  2. 2.

In theory, the sum-coordinate projection is defined as:

$P_{-}(\mathbf{r^{-}})=\int\Gamma(\mathbf{r},\mathbf{r^{+}}+\mathbf{r})d\mathbf{r}.$
(9)

In practice, it is calculated using the discrete space formula:

$P_{-}(\mathbf{r^{-}})=\sum_{x=1}^{M_{X}}\sum_{y=1}^{M_{Y}}\Gamma(\mathbf{r},\mathbf{r^{-}}+\mathbf{r}).$
(10)

where $\mathbf{r}=x\mathbf{e_{x}}+y\mathbf{e_{y}}$. This represents the
probability for two photons of a pair to be detected in coincidence between
pairs of pixels separated by an oriented distance
$\mathbf{r^{-}}=x^{-}\mathbf{e_{x}}+y^{-}\mathbf{e_{y}}$.

###  III Simulations of the JPD pixel super-resolution technique

####  III.1 Near-field illumination configuration.

Figure 8 shows a 1D simulation illustrating the principle of our technique
when using a near-field illumination configuration. In such a simple imaging
configuration, described in Figure 8.a, a thin type-I BBO crystal is
illuminated by a collimated blue pump laser to produce spatially-entangled
photon pairs. Crystal surface, object (grid-shape) and camera are positioned
in the same optical plane. In this simulation, photon-pairs are considered to
be strongly correlated in space Abouraddy  _et al._ (2002):
$\Gamma({x_{1}},{x_{2}})\approx\left|t(x_{1})t(x_{2})\delta(x_{1}-x_{2})\right|^{2}$,
where $\delta$ is a Dirac delta function. Figure 8.b shows the JPD of photon-
pairs measured without the object. The JPD is 2D-matrix because the imaging
system is 1D. The intense diagonal visible in the JPD is a clear signature of
the near-perfect position correlations between photon pairs. After inserting
the object, the diagonal of the JPD becomes modulated by the object shape
(Figure 8.c). Figure 8.d shows the JPD after pixellation. Figure 8.e shows the
JPD sum-coordinate projection, obtained by summing correlation values along
the anti-diagonals. This projection enables to retrieve a pixel super-resolved
image of the object i.e. $|t(x/2)|^{2}$. Figure 8.f shows the JPD minus-
coordinate projection, obtained by summing correlation values along the
diagonals. This projection shows an intensity correlation peak and is used in
our work to select diagonals with the highest weights in the JPD filtering
process (see section IV).

Figure 8: Simple illustration of JPD pixel super-resolution using a near-filed
illumination configuration. (a) 1D scheme for imaging with spatially-entangled
photon pairs using a near-field illumination configuration. Crystal surface,
object and camera are positioned in the same optical plane. (b) JPD without
object. (c) JPD with object. (d) JPD with object after pixellation. (e) Sum-
coordinate projection of the JPD. (f) Minus-coordinate projection of the JPD.

####  III.2 Far-field illumination configuration.

Figure 9 shows a 1D simulation illustrating the principle of our technique
when using a far-field illumination configuration. In this simple
configuration, described in Figure 9.a, a thin type-I BBO crystal is
illuminated by a collimated blue pump laser to produce spatially-entangled
photon pairs. The crystal surface is Fourier-imaged onto the object, that is
imaged onto the camera. In this simulation, photon-pairs are considered to be
strongly correlated Abouraddy  _et al._ (2002):
$\Gamma({x_{1}},{x_{2}})\approx\left|t(x_{1})t(x_{2})\delta(x_{1}+x_{2})\right|^{2}$
(where $\delta$ is a Dirac delta function). The object has a grid shape is
positioned in only half of the optical plane (the other half is empty). Figure
9.b shows the JPD of photon-pairs measured without the object. The JPD is a
2D-matrix because the imaging system is 1D. The intense anti-diagonal visible
in the JPD is a clear signature of the strong momentum correlations between
photon-pairs. After inserting the object, we observe that each half of the JPD
anti-diagonal become modulated by the object shape (Figure 9.d). Figure 8.d
shows the JPD after pixellation. Figure 8.e shows the JPD sum-coordinate
projection, obtained by summing correlation values along the anti-diagonals.
This projection shows a intensity correlation peak and is used in our work to
select anti-diagonals with the highest weights in the JPD filtering process
(see section IV). Figure 8.f shows the JPD sum-coordinate projection, obtained
by summing correlation values along the diagonals. This projection enables to
retrieve a pixel super-resolved image proportional to $|t(x/2)||t(-x/2)|^{2}$.
Because the object is positioned in only one half of the beam, the resulting
image contains two copies of the pixel-super-resolved image.

Figure 9: Simple illustration of JPD pixel super-resolution using a far-filed
illumination configuration. (a) 1D scheme for imaging with spatially-entangled
photon pairs using a far-field illumination configuration. Crystal surface is
Fourier-imaged onto the object positioned in one half of the optical plane.
Object is imaged on the camera. (b) JPD without object. (c) JPD with object.
(d) JPD with object after pixellation. (e) Sum-coordinate projection of the
JPD. (f) Minus-coordinate projection of the JPD.

####  III.3 1D simulations of JPD pixel super-resolution.

1D simulations were performed by adapting the integral formalism of entangled
photon pairs propagation provided in Abouraddy  _et al._ (2002) to a matrix
formalism computed with Matlab. Figure 10 shows simulated images of a grid-
shape object (5 pixels period) obtained using a pixelated sensor and a near-
field illumination configuration. Simulation parameters have been chosen to
match those used in the experiment shown in Figure 1.a of the manuscript.
Pixellation effect is well visible by eye when comparing the ground-truth
object (Figure 10.a) to the intensity image (Figure 10.b). Figure 10.c shows
the JPD sum-coordinate projection providing a better pixel resolution than the
intensity image. Such a better resolution is confirmed when comparing the
spatial frequency spectrum of both images in Figure 10.d. Indeed, the JPD
projection spectrum (blue curve) contains harmonics at $0.6$ cycles/pixel and
$0.8$ cycles/pixel, that are absent in the intensity image spectrum (red
curve). These simulation results are in very good accordance with the
experimental results obtained in Figure 1 of the manuscript.

Figure 10: 1D simulation of JPD pixel-super-resolution. (a) Ground-truth image
of a grid-shape object with $5$ pixels period (no pixellation). (b) Intensity
image of the object acquired using a pixelated sensor and a near-field
illumination configuration. (c) Sum-coordinate projection of the JPD acquired
using a pixelated sensor and a near-field illumination configuration. (d)
Spatial spectrum of the intensity (red curve) and sum-coordinate projection
(blue curve) calculated by discrete Fourier transform. Simulation parameters
have been chosen to match those used in the experiment shown in Figure 1.a of
the manuscript.

Figure 11 shows 1D simulations of spatial spectra calculated from images
obtained with grid-shape objects of different periods. Simulation parameters
have been chosen to be close to those used in the experimental results shown
in Figure 2 of the manuscript. These spectra further confirm that the JPD sum-
coordinate projection enables to retrieve spatial information that are not
present in the intensity image, even in the case where the first harmonics has
a higher frequency than the Nyquist frequency (Figure 11.c). These simulation
results are in very good accordance with the experimental results obtained in
Figure 2 of the manuscript. We however observe two small differences in Figure
11.c compared to Figure 2.c of the manuscript: (i) the absence of the residual
first order aliased peak in the JPD projection curve and (ii) the presence of
a second order aliased peak at $0.7$ cycles/pixel. (i) is explained by the
fact that the simulation do not use Equation 6 to produce the JPD and
therefore does not contain the residual background term $\mathcal{B}$ at the
origin of the first order aliased peak in Figure 2.c of the manuscript (see
section I). (ii) is explained by the fact that no optical frequency cut-off is
considered in the simulation, while it exists in the experiment. Such an
optical cut-off causes a more rapid decrease of the higher harmonics in the
experiment compared to the simulation, and therefore of their corresponding
aliased peaks. The second order aliased peak is then too weak to be visible in
the experimental spectrum shown in Figure 2.c, but is visible in the simulated
spectrum shown in Figure 11.c.

Figure 11: 1D simulations of spatial frequency spectra obtained from grid-
shape objects with different periods. Spatial frequency spectrum of intensity
image (red curve) and sum-coordinate projection (blue curve) obtained by
imaging grid-shape object with period of (a) $1/0.35\approx 2.9$ pixels, (b)
$1/0.5=2$ pixels and (c) $1/0.65\approx 1.5$ pixels. Simulation parameters
have been chosen to be close to those used in the experimental results shown
in Figure 2 of the manuscript.

###  IV Details on JPD filtering

The JPD filtering method consists in setting to $0$ values of the JPD that
have a negligible weight (i.e. values close to zero) so that their associated
noise do not contribute when performing the integration to calculate the JPD
projection. To perform this selection efficiently, correlations values with
the lowest weights are selected by groups. As explained in Methods, selecting
these groups of correlation values require to differentiate the two cases of
an object illuminated by photon pairs using (i) a near-field illumination
configuration (Figures 1.a, 3.a and 5.a of the manuscript) and (ii) a far-
field configuration (Figures 3.d and 4.a of the manuscript). In this section,
we will describe step-by-step the filtering process in the case of the near-
field illumination configuration. The same process can then be used in the
far-field case by inverting the roles of $P_{-}$ and $P_{+}$.

####  IV.1 Image of a JPD diagonal average values.

First, the minus-coordinate projection of the JPD (noted $P_{-}$) is
calculated using Equation 10. Each value $P_{-}(\mathbf{r^{-}})$ is the sum of
all the correlation values
$\Gamma(\mathbf{r_{1}},\mathbf{r_{2}}+\mathbf{r^{-}})$ located on the diagonal
defined by $\mathbf{r_{1}}-\mathbf{r_{2}}=\mathbf{r^{-}}$ and therefore
provides a direct estimation of the corresponding diagonal average weight.
Figure 12.a shows the minus-coordinate projection obtained using the JPD
measured in the experiment of Figure 1.a of the manuscript (only the center
area of $10\times 10$ pixels is shown). We observe that only a small number of
diagonals (pixels around the centrer) have significant weights, while all the
others are almost zero.

####  IV.2 Thresholding.

Second, a threshold $T$ is applied to the minus-coordinate projection to
select only the strongest values of $P_{-}$ and their associated diagonals in
the JPD. If $T=0$ (no threshold), all JPD values are kept and the resulting
JPD sum-coordinate projection is very noisy (Figure 12.e, same as Figure 6.a
of the manuscript). If $T=2\%$ of the image maximum value, only the diagonals
associated with the values located in the area delimited by the blue line in
Figure 12.d are kept, resulting into a sum-coordinate projection with less
noise (Figure 12.f). For the higher threshold values of $T=10\%$ (Figure 12.c)
and $T=50\%$ (Figure 12.d), the noise in the corresponding JPD sum-coordinate
projection shown in Figure 12.g and h is further reduced. In this case, the
best quality image is obtained using $T=50\%$ (Figure 12.g), which corresponds
to keeping only the correlations values of main JPD diagonal
$\Gamma(\mathbf{r},\mathbf{r})$ (i.e. center pixel of $P_{-}$) and those of
the $8$ other diagonals directly above and bellow i.e.
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{x}})$,
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{y}})$ and
$\Gamma(\mathbf{r},\mathbf{r}\pm\mathbf{e_{x}}\pm\mathbf{e_{y}})$ (direct-
neighbour around the center pixel in $P_{-}$). A similar filtering thresholds
was used also in all the other experiments shown in our work, resulting in
keeping only the $9$ diagonals (or anti-diagonals in far-field configuration)
with highest weights. However, note that using the highest threshold level to
keep only $9$ diagonals/anti-diagonals is not necessarily the best solution.
For example, if the correlation peak in the projection is asymmetric (e.g. due
to optical aberrations in the imaging system), using a too high threshold can
remove some diagonals in an asymmetric way and therefore create artefacts in
the retrieved high-resolution image.

Figure 12: JPD filtering. (a) Minus-coordinate projection of the JPD measured
in the experiment described in Figure 1.a of the manuscript before filtering
(centred zoom of $30\times 30$ pixels). (b) JPD minus-coordinate projection
obtained after using a threshold $T=2\%$ of the maximum. Values of the JPD
diagonals associated with pixel values outside the blue line are set to zero.
(c) JPD minus-coordinate projection obtained after using a threshold $T=10\%$
of the maximum. (d) JPD minus-coordinate projection obtained after using a
threshold $T=10\%$ of the maximum. (e) JPD sum-coordinate projection obtained
without filtering. (f) JPD sum-coordinate projection obtained after filtering
with threshold $T=2\%$ of the maximum. (g) JPD sum-coordinate projection
obtained after filtering with a threshold $T=10\%$ of the maximum. (h) JPD
sum-coordinate projection obtained after filtering with a threshold $T=50\%$
of the maximum.

###  V Details on entanglement-enabled quantum holography

This section provides more details about the entangled-enabled quantum
holography protocol used in Figure 4.a of the manuscript. This protocol has
been introduced in Defienne  _et al._ (2021a).

####  V.1 Theory.

Assuming that photon-pairs are near-perfectly correlated in space Abouraddy
_et al._ (2002), the vectorial two-photon wavefunction $\mathbf{\Psi}$ in the
object/SLM plane simplifies as:

$\displaystyle\mathbf{\Psi}(\mathbf{r_{1}},\mathbf{r_{2}})$ $\displaystyle=$
$\displaystyle\delta(\mathbf{r_{1}}+\mathbf{r_{2}})\big{[}t_{6}^{H}(\mathbf{r_{1}})t_{6}^{H}(\mathbf{r_{2}})\mathbf{e_{H}}$
$\displaystyle+$ $\displaystyle
e^{i(\alpha(\mathbf{r_{1}})+\alpha(\mathbf{r_{2}}))}t_{6}^{V}(\mathbf{r_{1}})t_{6}^{V}(\mathbf{r_{2}})\mathbf{e_{V}}\big{]}$

where $t_{6}^{H}$ (resp. $t_{6}^{V}$) is the horizontal (resp. vertical)
polarization component of the polarisation sensitive object $t_{6}$,
$\mathbf{e_{H}}$ and $\mathbf{e_{V}}$ are the unit horizontal and vertical
polarisation vectors, $\delta$ is the Dirac delta function and
$\alpha(\mathbf{r})$ is the phase programmed at position $\mathbf{r}$ of the
SLM. In the results shown in Figure 4 of the manuscript, $t_{6}$ is a
birefringent object i.e. $t_{6}^{H}=|t_{6}|e^{\theta_{H}}$ and
$t_{6}^{V}=|t_{6}|e^{\theta_{V}}$, where $\theta_{H}$ and $\theta_{V}$ are the
horizontal and vertical polarization components of the object phase,
respectively. As a result, after passing through the polariser at
$45^{\circ}$, the JPD measured by the camera takes the form of :

$\displaystyle\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$ $\displaystyle=$
$\displaystyle|t_{6}(\mathbf{r_{1}})t_{6}(\mathbf{r_{2}})\delta(\mathbf{r_{1}}+\mathbf{r_{2}})\big{[}1$
$\displaystyle+$
$\displaystyle\cos\big{(}\Delta\theta(\mathbf{r_{1}})+\Delta\theta(\mathbf{r_{2}})-\alpha(\mathbf{r_{1}})-\alpha(\mathbf{r_{2}})\big{)}\big{]}|^{2}$

where $\Delta\theta(\mathbf{r})=\theta_{H}(\mathbf{r})-\theta_{V}(\mathbf{r})$
is the object spatial phase that we are trying to retrieve.

####  V.2 Phase image retrieved using JPD anti-diagonals
$\Gamma(\mathbf{r},-\mathbf{r})$ (no pixel super-resolution).

To achieve quantum holography using the conventional protocol described in
Defienne  _et al._ (2021a), only the anti-diagonal component of the JPD
$\Gamma(\mathbf{r},-\mathbf{r})$ needs to be considered. Taking into account
that in our experiment (Figure 4.a of the manuscript) the object $t_{6}$ and
the SLM are positioned on each halves of the optical plane (i.e.
$t_{6}(x<0,y)=1$ and $\alpha(x>0,y)=0$), the JPD diagonal measured by the
camera then simplifies as:

$\displaystyle\Gamma(\mathbf{r},-\mathbf{r})$ $\displaystyle=$
$\displaystyle|t_{6}(\mathbf{r})|^{2}\left[1+\cos(\Delta\theta(\mathbf{-r})-\alpha(\mathbf{r}))\right]^{2}$
$\displaystyle\mbox{ , for all }\mathbf{r}=(x>0,y)$

In the holographic reconstruction process, four JPD anti-diagonals (noted
$\Gamma^{\alpha}(\mathbf{r},-\mathbf{r})$) are measured for four different
uniform phase shifts $\alpha(\mathbf{r})\in\\{0,\pi/2,\pi,3\pi/2\\}$
programmed on the SLM. Figures 13.c-f show the four anti-diagonal images
measured at each phase value. The spatial phase $\Delta\theta(\mathbf{r})$
shown in Figure 13.k (and also in Figure 4.b of the manuscript) is then
reconstructed from these four JPD anti-diagonals using the formula:
$\mbox{arg}[\Gamma^{0}(\mathbf{r},-\mathbf{r})-\Gamma^{\pi}(\mathbf{r},-\mathbf{r})+i(\Gamma^{\pi/2}(\mathbf{r},-\mathbf{r})-\Gamma^{3\pi/2}(\mathbf{r},-\mathbf{r}))]=\Delta\theta(\mathbf{r})$.
In practice, note also that a static and non-uniform phase distortion term
(noted $\Phi_{0}(\mathbf{r})$) is also present in Equation V.2, originating
from the photon generation process Hegazy and Obayya (2015). This term is not
shown here for clarify because it is directly compensated in the experiment by
superimposing a correction phase pattern $-\Phi_{0}(\mathbf{r})$ to all the
uniform phase $\alpha$ programmed on the SLM (see Defienne  _et al._ (2021a)
for more details about the phase correction).

Figure 13: Pixel super-resolution in entanglement-enabled quantum holography.
(a) and (b), Intensity images measured on each half of the camera. c-f, JPD
anti-diagonals $\Gamma^{\alpha}(\mathbf{r},-\mathbf{r})$ measured for
$\alpha=0$ (c), $\alpha=\pi/2$ (d), $\alpha=\pi$ (e) and $\alpha=3\pi/2$ (f).
(g-j), JPD minus-coordinate projections $P_{-}^{\alpha}(\mathbf{r^{-}})$
measured for $\alpha=0$ (g), $\alpha=\pi/2$ (h), $\alpha=\pi$ (i) and
$\alpha=3\pi/2$ (j). (k) Phase image retrieved using the anti-diagonal images
(same as Figure 4.c of the manuscript). (l) Phase image retrieved using the
minus-coordinate projections (same as Figure 4.d of the manuscript).

####  V.3 Phase image retrieved using JPD minus-coordinate projections
$P_{-}(\mathbf{r^{-}})$ (with pixel super-resolution).

To use our pixel super-resolution approach, the minus-coordinate projections
(noted $P_{-}^{\alpha}$) of the JPD is now considered instead of the anti-
diagonal. Taking into account that in our experiment (Figure 4.a of the
manuscript) the object $t_{6}$ and the SLM are positioned on each halves of
the optical plane (i.e. $t_{6}(x<0,y)=1$ and $\alpha(x>0,y)=0$), the JPD
minus-coordinate is calculated from Equation V.1 and simplifies as:

$\displaystyle P_{-}(\mathbf{r^{-}})$ $\displaystyle=$
$\displaystyle|t_{6}(\mathbf{r^{-}}/2)|^{2}\big{[}1+\cos(\Delta\theta(\mathbf{-r^{-}}/2)$
$\displaystyle-$ $\displaystyle\alpha(\mathbf{r^{-}}/2))\big{]}^{2}\mbox{, for
all }\mathbf{r}=(x>0,y)$

In the holographic reconstruction process, four JPD minus-coordinate
projections (noted $P_{-}^{\alpha}(\mathbf{r^{-}})$) are measured for four
different uniform phase shifts $\alpha(\mathbf{r})\in\\{0,\pi/2,\pi,3\pi/2\\}$
programmed on the SLM. Figures 13.g-j show the four minus-coordinate
projections measured at each phase value. The spatial phase $\Delta\theta$
shown in Figure 13.l (and also in Figure 4.c of the manuscript) is finally
reconstructed by combining the four JPD minus-coordination projections using
the formula:
$\mbox{arg}[P_{-}^{0}(\mathbf{r^{-}})-P_{-}^{\pi}(\mathbf{r^{-}})+i(P_{-}^{\pi/2}(\mathbf{r^{-}})-P_{-}^{3\pi/2}(\mathbf{r^{-}}))]=\Delta\theta(\mathbf{r^{-}}/2)$.
The resulting phase image has a better spatial resolution than the one shown
in Figure 13.k measured using the conventional entanglement-enabled quantum
holographic process.

###  VI Details on full-field NOON state quantum holography

The use of optical NOON-state (N=2) for holographic imaging has been
demonstrated by Ono et al Ono  _et al._ (2013) using a raster-scanning
approach. In our work, we extended it to full-field imaging using the
experiment described in Figure 5.a.

####  VI.1 Theory.

Assuming that photon-pairs are near-perfectly correlated in space Abouraddy
_et al._ (2002), the vectorial two-photon wavefunction $\mathbf{\Psi}$ in the
object plane of the experiment described in Figure 5.a can be written as:

$\displaystyle\mathbf{\Psi}(\mathbf{r_{1}},\mathbf{r_{2}})$ $\displaystyle=$
$\displaystyle\delta(\mathbf{r_{1}}-\mathbf{r_{2}})\big{[}t_{7}^{H}(\mathbf{r_{1}})t_{7}^{H}(\mathbf{r_{2}})\mathbf{e_{H}}$
$\displaystyle+$ $\displaystyle
e^{i(\alpha(\mathbf{r_{1}})+\alpha(\mathbf{r_{2}}))}t_{7}^{V}(\mathbf{r_{1}})t_{7}^{V}(\mathbf{r_{2}})\mathbf{e_{V}}\big{]}$

where $t_{7}^{H}$ (resp. $t_{7}^{V}$) is the horizontal (resp. vertical)
polarization component of the polarisation sensitive object $t_{7}$,
$\mathbf{e_{H}}$ and $\mathbf{e_{V}}$ are the unit horizontal and vertical
polarisation vectors, $\delta$ is the Dirac delta function and
$\alpha(\mathbf{r})$ is the phase programmed at position $\mathbf{r}$ of the
SLM. In the results shown in Figure 5 of the manuscript, $t_{7}$ is a
birefringent phase object produced by using a second SLM i.e.
$t_{7}^{H}=e^{\theta_{H}}$ and $t_{7}^{V}=e^{\theta_{V}}$, where $\theta_{H}$
and $\theta_{V}$ are the horizontal and vertical polarization components of
the phase. As a result, after passing through the polariser at $45^{\circ}$,
the JPD measured by the camera takes the form of :

$\displaystyle\Gamma(\mathbf{r_{1}},\mathbf{r_{2}})$ $\displaystyle=$
$\displaystyle\delta(\mathbf{r_{1}}-\mathbf{r_{2}})\big{[}1+\cos\big{(}\Delta\theta(\mathbf{r_{1}})$
$\displaystyle+$
$\displaystyle\Delta\theta(\mathbf{r_{2}})-\alpha(\mathbf{r_{1}})-\alpha(\mathbf{r_{2}})\big{)}\big{]}^{2}$

where $\Delta\theta(\mathbf{r})=\theta_{H}(\mathbf{r})-\theta_{V}(\mathbf{r})$
is the object spatial phase that we are trying to retrieve.

####  VI.2 Double-phase effect in NOON-state holography (N=2).

In the experimental configuration shown in Figure 5.a of the manuscript,
photon-pairs accumulate twice the phase of the object when passing through it
at a given spatial position. In practice, this double-phase effect is
precisely what provides this protocol a better phase sensitivity compared to
classical holography, as already demonstrated in previous works using raster-
scanning configurations Ono  _et al._ (2013); Israel  _et al._ (2014). A
simple way to visualise it is to measure JPD minus-coordinate projections
(noted $P_{-}^{\alpha}$) for different uniform phase pattern $\alpha$
programmed on the SLM with no object in the setup. In this case,
$P_{-}^{\alpha}$ can be calculated from Equation VI.1 and takes the following
theoretical form:

$P_{-}(\mathbf{r^{-}})=\left|1+\cos\left(2\alpha)\right)\right|^{2}\delta(\mathbf{r^{-}})$
(17)

where the double-phase effect is captured by the presence of the factor of 2
next to the SLM phase coefficient $\alpha$. Figures 14.a-c show JPD minus-
coordinate projections measured for different values of phases (a) $\alpha=0$,
(b) $\alpha=\pi/4$ and (c) $\alpha=\pi/2$ using the experimental setup in
Figure 5.a of the manuscript with no object. As predicted by Equation 17,
correlation value at the center peak (i.e. peak at
$\mathbf{r^{-}}=\mathbf{0}$) varies with the value of $\alpha$. Figure 14.d
shows correlation peak values for $29$ values of $\alpha\in[0,2\pi]$ (red
curve). This curve confirms the sine shape with a period of $1/2$ predicted in
Equation 17. As a comparison, average intensity values are measured using an
equivalent classical holography setup (see section VI.5) for
$\alpha\in[0,2\pi]$ are measured (blue curve). This curve shows the classical
sine shape curve with a period of $1$ (blue curve). Note that the contrast of
these sine shape curves are not unity because the horizontal and vertical
polarization components of entangled photons and classical light beams were
not perfectly balanced in our experiments.

Figure 14: Double-phase effect. (a - c) JPD minus-coordinate projections
measured using the experiment described in Figure 5.a of the manuscript with
no object for different uniform phase values programmed on the SLM (a)
$\alpha=0$, (b) $\alpha=\pi/4$ and (c) $\alpha=\pi/2$ (only central areas of
$30\times 30$ pixels are shown). (d) JPD minus-coordinate correlation peak
values (red curve) and average intensity values measured with an equivalent
classical holography setup (blue curve) are represented for different
$\alpha\in[0,2\pi]$. The equivalent classical holography setup is described in
section VI.5.

####  VI.3 NOON-state quantum holography using the JPD diagonal
$\Gamma(\mathbf{r},\mathbf{r})$ (no pixel super-resolution).

In the presence of a birefringent phase object $t_{7}$ in the experimental
setup, the JPD minus-coordinate projection cannot be used to retrieve its
phase component $\Delta\theta(\mathbf{r})$. Instead, the JPD diagonal
component $\Gamma(\mathbf{r},\mathbf{r})$ is considered:

$\Gamma(\mathbf{r},\mathbf{r})=\left|1+\cos(2\Delta\theta(\mathbf{r})-2\alpha(\mathbf{r}))\right|^{2}$
(18)

In the holographic reconstruction process, four JPD diagonals (noted
$\Gamma^{\alpha}(\mathbf{r},-\mathbf{r})$) are then measured for four
different uniform phase shifts
$\alpha(\mathbf{r})\in\\{0,\pi/4,\pi/2,3\pi/4\\}$ programmed on the SLM.
Figures 15.b-e show the four diagonal images measured at each phase value. The
spatial phase $2\Delta\theta(\mathbf{r})$ shown in Figure 15.j is finally
reconstructed by combining these four diagonal images using the formula:
$\mbox{arg}[\Gamma^{0}(\mathbf{r},\mathbf{r})-\Gamma^{\pi/2}(\mathbf{r},\mathbf{r})+i(\Gamma^{\pi/4}(\mathbf{r},\mathbf{r})-\Gamma^{3\pi/4}(\mathbf{r},\mathbf{r}))]=2\Delta\theta(\mathbf{r})$.
Note that, as predicted by the theory, the reconstructed phase is twice this
of the object (i.e. $2\Delta\theta(\mathbf{r})$) because of the double-phase
effect.

Figure 15: NOON-state quantum holography using JPD diagonals and sum-
coordinate projections. (a) Constant intensity image acquired by the camera
under photon-pairs illumination (not sensitive to the SLM phase). b-e JPD
diagonals measured for four different phase values programmed on the SLM (b)
$\alpha=0$, (c) $\alpha=\pi/4$, (d) $\alpha=\pi/2$ and (e) $\alpha=3\pi/4$.
f-i JPD sum-coordinate projections measured for four different phase values
programmed on the SLM (f) $\alpha=0$, (g) $\alpha=\pi/4$, (h) $\alpha=\pi/2$
and (i) $\alpha=3\pi/4$. (j) Phase image retrieved using the JPD diagonals.
(k) Phase image retrieved using the JPD sum-coordinate projections (same as in
Figure 5.c of the manuscript).

####  VI.4 NOON-state quantum holography using the JPD sum-coordinate
projection $P_{+}$ (with pixel super-resolution).

To use our pixel super-resolution approach, the sum-coordinate projections
(noted $P_{+}^{\alpha}$) of the JPD are now considered instead of the JPD
diagonals. The JPD sum-coordinate projection is calculated using Equation
VI.1:

$P_{+}(\mathbf{r^{-}})=\left|1+\cos(\Delta\theta(2\mathbf{r^{+}}/2)-2\alpha(\mathbf{r^{+}}/2))\right|^{2}$
(19)

In the holographic reconstruction process, four JPD sum-coordinate projections
(noted $P_{+}^{\alpha}(\mathbf{r^{-}})$) are measured for four different
uniform phase shifts $\alpha(\mathbf{r})\in\\{0,\pi/4,\pi/2,3\pi/4\\}$
programmed on the SLM. Figures 15.f-i show the four sum-coordinate projections
measured at each phase value. The spatial phase $2\Delta\theta$ shown in
Figure 15.l (and also in Figure 5.c of the manuscript) is finally
reconstructed by combining these four JPD sum-coordinate projection using the
formula:
$\mbox{arg}[P_{+}^{0}(\mathbf{r^{+}})-P_{+}^{\pi/2}(\mathbf{r^{+}})+i(P_{+}^{\pi/4}(\mathbf{r^{+}})-P_{+}^{3\pi/4}(\mathbf{r^{+}}))]=2\Delta\theta(\mathbf{r^{+}}/2)$.
The resulting phase image has a better spatial resolution than the one shown
in Figure 15.k measured using the JPD diagonal. Note that the reconstructed
phase is also twice this of the object (i.e. $2\Delta\theta(\mathbf{r})$)
because of the double-phase effect.

####  VI.5 Equivalent classical holography setup.

The classical holography setup used to retrieve the phase in Figure 5.d of the
manuscript is shown in Figure 16.a. It is the same setup than this described
in Figure 5.a of the manuscript after substituting the photon-pairs source by
a collimated classical laser beam at $810$nm and polarised at $45^{\circ}$. In
the holographic reconstruction process, four intensity images (noted
$I^{\alpha}$) are measured for four different uniform phase shifts
$\alpha(\mathbf{r})\in\\{0,\pi/2,\pi,3\pi/2\\}$ programmed on the SLM. Figures
16.c-f show the four intensity images measured at each phase value. The
spatial phase of the object $\Delta\theta$ (Figure 16.b and Figure 5.d of the
manuscript) is finally reconstructed by combining these four intensity images
using the formula:
$\mbox{arg}[I^{0}(\mathbf{r})-I^{\pi}(\mathbf{r})+i(I^{\pi/2}(\mathbf{r})-I^{3\pi/2}(\mathbf{r}))]=\Delta\theta(\mathbf{r}/2)$.

Figure 16: Equivalent classical holography setup. (a) Experimental setup used
to perform classical holography. It is the same than this described in Figure
5.a of the manuscript after substituting the photon-pairs source by a
collimated classical laser beam at $810$nm and polarised at $45^{\circ}$. (b)
Phase image retrieved by classical holography (same as Figure 5.d of the
manuscript). (c-f) Intensity images measured for four different uniform phase
values programmed on the SLM (c) $\alpha=0$, (d) $\alpha=\pi/2$, (e)
$\alpha=\pi$ and (f) $\alpha=3\pi/2$.

####  VI.6 Other artificial and real birefringent phase objects measured using
full-field NOON-state holography (N=2).

Figure 17 shows phase images of other artificial and real birefringent objects
retreived by NOON-state quantum holography and classical holography. Note that
the range of application of this full-field NOON-state holography technique
can be extended to non-birefringent phase objects by using birefringent
plates, as done in differential interference contrast microscopy Allen and
David (1969).

Figure 17: Phase images of other artificial and real birefringent objects.
(a-c) Phase images of a (a) smiley face and (b) the number 15 created using
another SLM retrieved by full-filed NOON-state quantum holography using the
sum-coordinate projections. Phase image of (c) pieces of scotch tape forming a
cross retrieved by full-filed NOON-state quantum holography (N=2). (d-f) Phase
images of the same (d) smiley face, (e) the number 15 and (f) pieces of scotch
tape retrieved by classical holography.

###  VII Diagonal and anti-diagonal images

Figure 18 shows diagonal and anti-diagonal images associated with the
experiments shown in Figure 1.a, Figure 3.a and d, and Figure 5.a of the
manuscript.

Figure 18: Diagonal and anti-diagonal images. (a) Diagonal image
$\Gamma(\mathbf{r},\mathbf{r})$ measured using the experiment shown in Figure
1.a of the manuscript. (b) Diagonal image $\Gamma(\mathbf{r},\mathbf{r})$
measured using the experiment shown in Figure 3.a of the manuscript. (c) Anti-
diagonal image $\Gamma(\mathbf{r},-\mathbf{r})$ measured using the experiment
shown in Figure 3.d of the manuscript. (d) Diagonal image
$\Gamma(\mathbf{r},\mathbf{r})$ measured using the experiment shown in Figure
5.a of the manuscript.

### References

  * Bishara  _et al._ (2010) W. Bishara, T.-W. Su, A. F. Coskun, and A. Ozcan, Optics Express 18, 11181 (2010), publisher: Optical Society of America.
  * Bishara  _et al._ (2011) W. Bishara, U. Sikora, O. Mudanyali, T.-W. Su, O. Yaglidere, S. Luckhart, and A. Ozcan, Lab on a Chip 11, 1276 (2011).
  * Moreau  _et al._ (2019) P.-A. Moreau, E. Toninelli, T. Gregory, and M. J. Padgett, Nature Reviews Physics 1, 367 (2019).
  * Boto  _et al._ (2000) A. N. Boto, P. Kok, D. S. Abrams, S. L. Braunstein, C. P. Williams, and J. P. Dowling, Physical Review Letters 85, 2733 (2000).
  * Reichert  _et al._ (2018) M. Reichert, H. Defienne, and J. W. Fleischer, Scientific reports 8, 7925 (2018).
  * Ono  _et al._ (2013) T. Ono, R. Okamoto, and S. Takeuchi, Nature communications 4, 2426 (2013).
  * Israel  _et al._ (2014) Y. Israel, S. Rosen, and Y. Silberberg, Physical Review Letters 112, 103604 (2014).
  * Brida  _et al._ (2010) G. Brida, M. Genovese, and I. R. Berchera, Nature Photonics 4, 227 (2010).
  * Pittman  _et al._ (1995) T. B. Pittman, Y. H. Shih, D. V. Strekalov, and A. V. Sergienko, Physical Review A 52, R3429 (1995).
  * Aspden  _et al._ (2013) R. S. Aspden, D. S. Tasca, R. W. Boyd, and M. J. Padgett, New Journal of Physics 15, 073032 (2013).
  * Lopaeva  _et al._ (2013) E. D. Lopaeva, I. Ruo Berchera, I. P. Degiovanni, S. Olivares, G. Brida, and M. Genovese, Physical Review Letters 110, 153603 (2013).
  * Defienne  _et al._ (2019) H. Defienne, M. Reichert, J. W. Fleischer, and D. Faccio, Science Advances 5, eaax0307 (2019).
  * Gregory  _et al._ (2020) T. Gregory, P.-A. Moreau, E. Toninelli, and M. J. Padgett, Science Advances 6, eaay2652 (2020).
  * Devaux  _et al._ (2019) F. Devaux, A. Mosset, F. Bassignot, and E. Lantz, Physical Review A 99, 033854 (2019).
  * Defienne  _et al._ (2021a) H. Defienne, B. Ndagano, A. Lyons, and D. Faccio, Nature Physics 17, 591 (2021a).
  * Lemos  _et al._ (2014) G. B. Lemos, V. Borish, G. D. Cole, S. Ramelow, R. Lapkiewicz, and A. Zeilinger, Nature 512, 409 (2014).
  * Kviatkovsky  _et al._ (2020) I. Kviatkovsky, H. M. Chrzanowski, E. G. Avery, H. Bartolomaeus, and S. Ramelow, Science Advances 6, eabd0264 (2020), publisher: American Association for the Advancement of Science Section: Research Article.
  * Vanselow  _et al._ (2020) A. Vanselow, A. Vanselow, P. Kaufmann, I. Zorin, B. Heise, H. M. Chrzanowski, S. Ramelow, and S. Ramelow, Optica 7, 1729 (2020), publisher: Optical Society of America.
  * Moreau  _et al._ (2012) P.-A. Moreau, J. Mougin-Sisini, F. Devaux, and E. Lantz, Physical Review A 86, 010101 (2012).
  * Edgar  _et al._ (2012) M. P. Edgar, D. S. Tasca, F. Izdebski, R. E. Warburton, J. Leach, M. Agnew, G. S. Buller, R. W. Boyd, and M. J. Padgett, Nature Communications 3, 984 (2012).
  * Chrapkiewicz  _et al._ (2014) R. Chrapkiewicz, W. Wasilewski, and K. Banaszek, Optics Letters 39, 5090 (2014).
  * Unternährer  _et al._ (2016) M. Unternährer, B. Bessire, L. Gasparini, D. Stoppa, and A. Stefanov, Optics Express 24, 28829 (2016).
  * Toninelli  _et al._ (2019) E. Toninelli, P.-A. Moreau, T. Gregory, A. Mihalyi, M. Edgar, N. Radwell, and M. Padgett, Optica 6, 347 (2019).
  * Chrapkiewicz  _et al._ (2016) R. Chrapkiewicz, M. Jachura, K. Banaszek, and W. Wasilewski, Nature Photonics 10, 576 (2016).
  * Bronzi  _et al._ (2014) D. Bronzi, F. Villa, S. Tisa, A. Tosi, F. Zappa, D. Durini, S. Weyers, and W. Brockherde, IEEE Journal of Selected Topics in Quantum Electronics 20, 354 (2014).
  * Defienne  _et al._ (2021b) H. Defienne, J. Zhao, E. Charbon, and D. Faccio, Physical Review A 103, 042608 (2021b).
  * Defienne  _et al._ (2018) H. Defienne, M. Reichert, and J. W. Fleischer, Physical review letters 120, 203604 (2018).
  * Abouraddy  _et al._ (2002) A. F. Abouraddy, B. E. Saleh, A. V. Sergienko, and M. C. Teich, JOSA B 19, 1174 (2002).
  * Reichert  _et al._ (2017) M. Reichert, H. Defienne, X. Sun, and J. W. Fleischer, Journal of Optics 19, 044004 (2017).
  * Farsiu  _et al._ (2004) S. Farsiu, M. Robinson, M. Elad, and P. Milanfar, IEEE Transactions on Image Processing 13, 1327 (2004), conference Name: IEEE Transactions on Image Processing.
  * Schneeloch and Howell (2016) J. Schneeloch and J. C. Howell, Journal of Optics 18, 053501 (2016).
  * Allen and David (1969) R. D. Allen and G. B. David, Zeitschrift fur wissenschaftliche Mikroskopie und mikroskopische Technik 69, 193 (1969).
  * Terborg  _et al._ (2016) R. A. Terborg, J. Pello, I. Mannelli, J. P. Torres, and V. Pruneri, Science Advances 2, e1600077 (2016).
  * Luo  _et al._ (2016) W. Luo, Y. Zhang, A. Feizi, Z. Göröcs, and A. Ozcan, Light: Science & Applications 5, e16060 (2016), number: 4 Publisher: Nature Publishing Group.
  * Shi  _et al._ (2016) W. Shi, J. Caballero, F. Huszar, J. Totz, A. P. Aitken, R. Bishop, D. Rueckert, and Z. Wang (2016) pp. 1874–1883.
  * Ulku  _et al._ (2019) A. C. Ulku, C. Bruschini, I. M. Antolović, Y. Kuo, R. Ankri, S. Weiss, X. Michalet, and E. Charbon, IEEE Journal of Selected Topics in Quantum Electronics 25, 1 (2019), conference Name: IEEE Journal of Selected Topics in Quantum Electronics.
  * Farsiu  _et al._ (2006) S. Farsiu, M. Elad, and P. Milanfar, in  _Visual Communications and Image Processing 2006_, Vol. 6077 (International Society for Optics and Photonics, 2006) p. 607703.
  * Takeda  _et al._ (2007) H. Takeda, S. Farsiu, and P. Milanfar, IEEE Transactions on Image Processing 16, 349 (2007), conference Name: IEEE Transactions on Image Processing.
  * Yamaguchi and Zhang (1997) I. Yamaguchi and T. Zhang, Optics Letters 22, 1268 (1997).
  * Hegazy and Obayya (2015) S. F. Hegazy and S. S. A. Obayya, JOSA B 32, 445 (2015).