Shaping the spatial correlations of entangled photon pairs

Quantum imaging enhances imaging systems performance, potentially surpassing fundamental limits such as noise and resolution. However, these schemes have limitations and are still a long way from replacing classical techniques. Therefore, there is a strong focus on improving the practicality of quantum imaging methods, with the goal of finding real-world applications. With this in mind, in this tutorial we describe how the concepts of classical light shaping can be applied to imaging schemes based on entangled photon pairs. We detail two basic experimental configurations in which a spatial light modulator is used to shape the spatial correlations of a photon pair state and highlight the key differences between this and classical shaping. We then showcase two recent examples that expand on these concepts to perform aberration and scattering correction with photon pairs. We include specific details on the key steps of these experiments, with the goal that this can be used as a guide for building photon-pair-based imaging and shaping experiments.


I. INTRODUCTION
Quantum imaging is the field of optics that aims to exploit the quantum properties of light to improve the performance of imaging systems [1].Among the existing non-classical sources of light, those that emit pairs of entangled photons have proven to be particularly fruitful.For example, imaging with sensitivity below the shot-noise limit [2] and with resolution beating classical diffraction [3,4] were both achieved with photon pair states.In addition, these sources also enabled the development of new imaging modalities.Examples include ghost imaging [5], imaging with undetected photons [6], entanglement-enabled holography [7] and Hong-Ou-Mandel microscopy [8].Thanks to recent technological advancements, photon-pair-based imaging systems have become relatively simple to implement.Firstly, entangled photon pairs can be produced using spontaneous parametric down conversion (SPDC) [9].This is a non-linear optical process that converts a single, highenergy photon, called the pump, into two lower energy photons, traditionally called the signal and idler.It does not require complex laser systems or highly specialised non-linear crystals.For example, many experiments -like the one presented here -use standard room-temperature β−Barium Borate (BBO) crystals and a blue diode laser with a power of a few tens of milliwatts.Then, another factor that adds to the ease of implementation is the recent advent of commercial single-photon-sensitive cameras.Paired with advanced image processing techniques, these cameras enable the detection of photon coincidences across many transverse spatial positions simultaneously -a key measurement when working with entangled photon pairs.For exam-ple, electron-multiplying charge-coupled devices (EMC-CDs) and single-photon avalanche diodes (SPADs) cameras have enabled the measurement of intensity correlations functions i.e. the G (2) for photon pairs across thousands of spatial modes [10][11][12][13].Despite this recent progress, the low brightness of SPDC sources (i.e.tens of picowatts) and long acquisition times (e.g.sometimes up to several hours) still restrict their application to the laboratory environment and specific niche areas.
Besides generation and detection, it is also important to control how photon pairs and their characteristics (e.g.correlations) propagate within the optical system.For example, spatial light modulators (SLMs) can be placed in the pump beam or the photon pairs path to manipulate their spatial properties.They enable tailoring of the spatial entanglement between photons propagating in free space by manipulating their orbital angular spectrum [14][15][16][17][18][19][20], by shaping momentum-position correlations [21][22][23] or by tuning the pump spatial coherence [24,25].In addition, they can also be used to control photon pair propagation through complex media such as multimode fibers [26][27][28] and scattering layers [29][30][31][32].These efforts are mainly directed towards applications in communication and information processing.
Though less explored, such control can also be of interest for quantum imaging.Indeed, if we consider the classical domain, light structuring has many applications in imaging, particularly in microscopy [33].Examples range from structured illumination microscopy [34], which aims to improve imaging resolution by illuminating the sample with specifically tailored pattern of light, to adaptive optics [35] and wavefront shaping [36], which use an SLM to correct for optical aberrations and scattering, with applications also for contrast-enhanced and quantitative phase imaging [37,38].It is therefore interesting to see how we could apply these concepts to quantum imaging and it is within this context that this tutorial is written.Here, we describe two basic experiments that both use an SLM to shape the spatial correlations between entan-gled photon pairs.We detail their practical implementation and highlight their differences from their classical counterparts.Additionally, we show examples of applications that illustrate the practical uses of this two-photon shaping, including a case in imaging.We also include detailed descriptions of many of the key steps in this work in the Supplementary Material (S.M.), and provide example code in Ref. [39].

II. EXPERIMENTAL SETUP
When shaping spatial correlations between pairs of photons, two experimental configurations can be considered.These are called the far-field (FF) and near-field (NF) shaping configurations, named for the position of the SLM relative to the nonlinear crystal.Note that in other works the configurations may instead be named for the position of the detector relative to the crystal.Fig- ures 1a,b shows schematics of these two configurations.In both cases a nonlinear crystal, a thin β-Barium Borate (BBO) crystal, is illuminated by a continuous-wave (CW) pump laser at 405nm.This produces pairs of photons at 810nm via Type I SPDC.Immediately after the crystal, the pump photons are filtered out using a 650nm-cut-off long-pass filter, allowing only the photons of higher wavelengths generated by SPDC to propagate through the system.In practice, it is common to also include a shortpass filter in the pump beam's path before the crystal (not shown).This filter removes residual low-frequency light emitted by the laser, typically found at the output of blue diode lasers.After propagating through the system, the pairs are detected with an EMCCD camera, and the spatial correlations are measured following the method described in Ref. [12] (see also Methods).In both configurations, the pairs are controlled by a spatial light modulator that is positioned in a Fourier plane of the camera.A bandpass filter at 810 ± 10nm is used to remove ambient light and most non-degenerate pairs.Due to the phase matching conditions in the SPDC process, the photons at the plane of the crystal are strongly correlated in transverse position (r) and strongly anticorrelated in transverse momentum (k) [40].Both of these types of correlations can be exploited, hence the two experimental configurations.
The FF shaping configuration is shown in Figure 1a.As suggested by the name, the SLM is positioned in a Fourier plane of the crystal, and its surface is imaged on the EMCCD.Since the pairs are correlated in position at the crystal plane, they are also correlated in position at the camera, so they arrive at approximately the same pixel.Figure 1c shows an example of the intensity at the camera.The shape of this is mostly dependent on the intensity profile of the pump.Figure 1d shows the so-called minus-coordinate projection of the measured G (2) of the photon pairs.The way this is measured is detailed in the Methods section.It represents the probability of detecting both photons of a pair simultaneously on two pixels of the camera separated by a distance (x 2 − x 1 ,y 2 − y 1 ), expressed in pixels.This type of measurement is essential when working with entangled photon pairs because it provides information about the spatial correlations between the pairs, something that an intensity image alone does not provide.A bright, narrow peak, as seen in the figure, indicates that the pairs are strongly correlated in position.This means that, when one photon in a pair is detected at a position (x 1 , y 1 ), there is a high probability to detect its twin within a very small area around this pixel (x 2 , y 2 ) ≈ (x 1 , y 1 ).For the measurement of the data shown in Figure 1d, the SLM did not display any phase mask.
The NF shaping configuration is shown in Figure 1b.Here, the SLM is positioned in an conjugate plane to the surface of the crystal, and the Fourier plane of the crystal is imaged onto the camera.Now, due to momentum conservation in the SPDC process, the pairs arrive at the camera at diametrically opposite pixels, relative to the centre of the beam.Figure 1e shows an example of the intensity at the camera.Here, the typical SPDC ring can be seen.The thickness of this ring is proportional to the bandwidth of the pairs, and the radius is dependent on the angle between the pump optical axis and the normal to the crystal surface.Figure 1f shows the so-called sum-coordinate projection of the measured G (2) .The method of measuring this projection is also explained in the Methods section.It represents the probability of detecting two entangled photons at any pair of pixels of the camera (x 1 , y 1 ) and (x 2 , y 2 ) with a given barycenter value (x 1 + x 2 , y 1 + y 2 ).For example, the central value of the sum-coordinate projection corresponds to the sum of all coincidence rates detected between pairs of pixels with a barycenter at (0, 0), which means all the pairs of pixels that are exactly anti-symmetric i.e. (x 2 , y 2 ) = (−x 1 , −y 1 ).In our experiment, the presence of an intense peak at the center indicates strong spatial anti-correlation.This means that when a photon from a pair is detected at a position (x 1 , y 1 ), there is a high probability that its twin will be detected within a very small area around the symmetric pixel (x 2 , y 2 ) ≈ (−x 1 , −y 1 ).No phase mask was displayed on the SLM to perform the measurement shown in Figure 1f.

III. INTENSITY AND CORRELATION SHAPING THEORY
In this section we will describe the theory behind twophoton correlation shaping and compare it to classical shaping.When comparing classical and quantum imaging, it is useful to work with the spatial intensity correlation functions, specifically, the first and second order functions G (1) and G (2)  [41].The quantities we measure in experiment can be written in terms of these, and they can also be used to derive expressions for propagation through an optical system.
First we introduce the relevant measured quantities, and how they are written in terms of G (1) and G (2) .When imaging and shaping with classical coherent light, the quantity we measure with the camera is the intensity of the electric field, I(r) = |E(r)| 2 .This can be written in terms of the first-order spatial correlation function of the field as: where Ê(+) and Ê(−) are the positive and negative frequency component of the quantum operator associated with the electric field, respectively.In practice, the intensity is conventionally measured using a camera by accumulating photons on each pixel.
When working with photon pairs, the quantity that interests us the most is the second-order spatial correlation function of the intensity, G (2) (r 1 , r 2 ).This is also often referred as the joint probability distribution of the photons.It is written in terms of the second-order spatial correlation function of the field as: (3) The G (1) describes the field correlations and the second-order correlation function can be seen as the intensity-correlation analogue of this.In practice, it is measured by detecting photon coincidence between pairs of spatial positions r 1 and r 2 .As detailed in Ref. [12], it can also be reconstructed by measuring the intensity covariance between pairs of pixels on a camera.Assuming we have pure two-photon states, G (2) (r 1 , r 2 ) gives the probability of simultaneously detecting two photons from a pair at positions r 1 and r 2 .This quantity is related to the spatial two-photon wave function ϕ, as: In this way, ϕ is to G (2) in the two-photon case as the field E is to I in the case of coherent light.
Since we want to see how we can shape the distributions I and G (2) , it is useful to know how they propagate.Any linear system can be described by its coherent pointspread function (PSF).Given the PSF, written h(r ′ , r), the correlation functions can be propagated through any linear system as: and For example, let's consider the two experimental configurations described in Figures 1a and b (NF and FF shaping).In both cases, the camera is positioned in a Fourier plane of the SLM.Since the action of the SLM is to impart a phase profile, θ(r), onto the beam, the PSF associated with the propagation between the SLM (shaping plane) and the camera (detection plane) can be written as where f is the focal length of the lens immediately after the SLM, λ is the wavelength of the light being used, r and r ′ are the transverse coordinates in the SLM and camera plane, respectively.Note that in each configuration there are in reality several lenses performing a magnification in addition to the Fourier transform, which can be taken into account by changing the effective value of f in the previous equation.

A. Shaping coherent light
For a perfectly spatially coherent source, e.g. a laser, the first-order correlation function is given by where E(r) is the electric field at position r in the SLM plane and * denotes the complex conjugate.Thus, from equations 2 and 4, the intensity after propagation through the system is given by Using equation 6 and simplifying, we obtain the intensity at the camera where F[...] is the 2-dimensional Fourier transform, and * denotes the 2-dimensional convolution operation.This is the expected result from Fourier-optics and says that the field measured at the camera is simply the (scaled) Fourier transform of the phase mask on the SLM.The intensity correlation function for a coherent source is simply For coherent light, G (2) can thus be fully computed from the intensity measurements and therefore does not contain any additional information.

B. Shaping entangled photon pairs
Now let's consider a two-photon state.In the input plane i.e. the SLM plane, we can write our state in the position basis as where ϕ(r 1 , r 2 ) is the two-photon wavefunction expressed in the position basis, r 1 and r 2 are the transverse positions in the SLM plane for photon 1 and photon 2 respectively.In our case, both photon have the same polarisation and spectral frequency.This results from the pair generation process that we have chosen to use in our experiments, which is Type I SPDC.In this simplified notation, |r 1 , r 2 ⟩ = |r 1 ⟩ 1 ⊗ |r 2 ⟩ 2 denotes the state in which photon 1 is at position r 1 and photon 2 is at position r 2 .The G (1) for such a state is given by and the direct intensity is then Additionally, one can compute the second-order field correlations as: giving the intensity correlations: Using equations 4 and 5, I pairs and G pairs can be expressed in the camera plane: and where r 1 ′ and r 2 ′ are the transverse positions in the camera plane, ϕ(r 1 , r 2 ) is the two-photon wavefunction in the SLM plane and h is the PSF.Now we consider the two shaping configurations shown in Figures 1a and b.In both cases, the imaging system from the SLM to the camera, h, is the same (up to a different magnification factor), but the input states are different.We start with the configuration of Figure 1a.Here, the SLM is positioned in the Fourier plane of the crystal.In our experimental conditions the collimated pump beam diameter is much larger than the crystal thickness.Therefore, we can assume that photon pairs are near-perfectly anti-correlated in the plane immediately before the SLM [42].That is, the wavefunction , where ϕ 0 is the amplitude envelope of the two-photon wavefunction in the crystal Fourier plane.It is linked to the intensity measured in the SLM plane as I(r) = |ϕ 0 (r)| 2 .In practice, it takes the shape of a disk or a ring, as shown in Figure 1.e, and its spatial phase is assumed to be uniform.After performing the change of variables r + = (r 1 + r 2 )/2 and r − = (r 1 − r 2 )/2, the intensity correlation in the camera plane can be expressed as where ψ(k) = θ(k) + θ(−k) and all global constants are omitted for clarity.Now we consider the configuration in Figure 1b.Here we are imaging the crystal plane onto the SLM.Since we have thin crystal and a large pump diameter, we can assume that the photon pairs are perfectly correlated, i.e. [42], where ϕ ′ 0 is the amplitude envelope of the two-photon wavefunction in the crystal plane.It is linked to the intensity measured in the SLM plane as: In practice, it takes the shape of the pump beam, as shown in Figure 1.c, and its spatial phase is assumed to be uniform.Then, the intensity correlations at the camera are given by

) By performing a similar propagation calculation with G
(1) pairs , we find that the intensities on the camera plane are spatially uniform in both configurations and do not depend on the phase displayed on the SLM i.e.I pairs (r ′ ) is a constant.Thus, as would be the case for a perfectly spatially incoherent source, the SLM does not modulate the intensity measured in the Fourier plane.This nearperfect spatial incoherence is indeed observed experimentally and is a consequence of our experimental conditions, specifically the use of a collimated pump with a diameter much larger than the thickness of the crystal.By changing the crystal and the illumination conditions [43], it would be possible to work in an intermediate regime with partially spatially coherent light, allowing modulation of both the intensity and the intensity correlations.
Comparing equations 9, 18 and 19 we see that the spatial intensity correlations can be shaped in a manner that is almost equivalent to spatial intensity shaping in the classical case.Under our experimental conditions, F[E(r)], F[ϕ 0 (r)] and F[ϕ ′ 0 (r)] are very sharply peaked, and can therefore be approximated as Dirac Delta functions.This allows for the simplification of the equations to focus on the role played by the SLM: Interestingly, with entangled photon pairs we find that depending on the configuration, the shaping behaves slightly differently.In the case of FF shaping, we see that the phase on the SLM affects the projection in the form of the function ψ(r) = θ(r) + θ(−r).This suggests that what shapes the spatial correlations is effectively the SLM mask plus a spatially inverted version of the mask.In the case of NF shaping, we see that the projection is affected by a 2θ(r) phase term, i.e. the pairs 'see' a phase that is double that of the actual mask on the SLM.

IV. RESULTS
We now show how to experimentally measure and modulate the photon-pair G (2) , following equations 18 and 19.For this, we perform shaping experiments using the FF and NF configurations shown in Figures 1a and b, respectively.
From a practical standpoint, it is important to first clarify a key aspect regarding the measurement of the spatial G (2) that we perform.The G (2) is measured using the approach described in Ref. [12] and in Methods.This quantity takes the form of a four-dimensional matrix.To visualize it, we project it into two-dimensional images as seen in examples in Figures 1d and f.Depending on whether we use the FF or NF configuration, we perform projections onto the minus-or sumcoordinates, respectively.These projections are formally defined as spatial averages of the measured G (2) function: (2) (r − , r + )dr − for the minus-coordinate projection, and C + (r + ) = S G (2) (r − , r + )dr − for the sum-coordinate projection, where G (2) is expressed here using the variables (r + , r − ) and S is the integration area.As seen in equations 18 and 19, these projections follow the symmetries of the measured G (2) that we expect in each configuration: These projections effectively reveal the spatial modulation of G (2) , which is what we aim to observe, while max-imising the signal-to-noise ratio in the measured signal.In practice, this spatial averaging reduces the acquisition time needed to resolve the correlations from several hours to just a few minutes.
Figure 2 shows the results for the FF shaping configuration, where we are imaging the near-field of the crystal on the camera (Fig. 1a).As a simple demonstration, we use a π/2-modulated phase grating pattern, shown in Figure 2a.According to equation 18, the photon pairs in this system 'see' the actual grating plus a spatially inverted version.Therefore, if θ(r) + θ(−r) = const, then we expect to see no modulation of the minus-coordinate projection.To demonstrate this, the grating is displayed on the SLM and translated laterally, measuring C − at each lateral shift.Figure 2b shows the values of each diffraction order as a function of lateral shift.If the grating is positioned such that one of the steps is at the centre of the SPDC beam, then we have ψ(r) = π/2 = constant and we expect no modulation in the ideal case.Figure 2c shows the minus-coordinate projection measured in such a case.Even though we do not observe a complete extinction of higher orders, mainly due to the fact that the phase pattern is not perfectly asymmetric, in practice the measurements in Figure 2a confirm that these higher orders are minimized while the zero order is maximized.Instead, if the grating is positioned 1/4 of a grating period away from this, then we get ψ(r) = 2θ(r), and the correlations should be maximised in the first-order diffraction peaks.This can be seen in Figure 2d.At the intermediate grating positions, ψ(r) contains higher frequency components, and we observe second-order diffraction peaks.Furthermore, it's interesting to note that replicating this experiment with classical coherent light would not yield any change in the intensity-measured diffraction pattern.Indeed, a lateral shift in the phase pattern in the Fourier plane would only influence the spatial phase component of the diffraction pattern, which the camera is not sensitive to.
According to equation 19, we expect a different phenomenon in the case of the NF shaping configuration (Fig. 1b), where we are now imaging the Fourier plane of the crystal on the camera.As before, we display a phase grating on the SLM (Fig 3a) but, since we instead have the 2θ(r) phase term, lateral translation will have no effect on C + .Instead, we vary the amplitude of the phase grating α from 0 to 2π radians, and measure C + at each step.Figure 3b shows the values of the diffraction orders at each grating amplitude.As a comparison, we also performed the same experiment using classical coherent light and recorded the first-order diffraction peak intensity in function of α (pink curve).As expected, it follows a sinusoidal pattern reaching its maximum at α = π.In the case of entangled photon pairs, we see the oscillation but, since C + (δr + ) is modulated by twice the phase mask, the frequency of this oscillation is doubled.Figures 3c-e show examples of the sum-coordinate projections measured for α equals to π/6, π/2 and 2π/9, respectively.Finally, it is important to note that in se two experiments described, the intensity images measured on the camera do not change regardless of the phase pattern displayed on the SLM: they remain the same as that shown in Figures 1c and e.
In light of these results, it is interesting to discuss the physical interpretation of such a G (2) shaping demonstration, in particular using the corpuscular picture.By modulating G (2) , we effectively change the joint probability of detecting both photons of a pair at two specific camera pixels.However, we achieve this without changing the probability of detecting a photon at a given camera pixel (i.e. the marginal probability), since the intensity remains uniform.For instance, in the case of Figure 2d, the probability has been modified in such a way that the photons no longer arrive on the camera strongly correlated in position; instead, they are now forced to be separated from each other by a distance defined by the period of the phase grating programmed on the SLM.In essence, this means that we can control the collective behaviour of photon pairs without changing their individual behaviour.

V. APPLICATIONS
In this section we show two examples of experiments that use these concepts of two-photon shaping.The first is the method of quantum-assisted adaptive optics (QAO), as reported in Ref. [44].In this work the spatial correlations of entangled photons are exploited to gain information about the optical aberrations within an imaging system.Effectively, the central value of the sumcoordinate projection is used as a guidestar to implement an adaptive optics (AO) scheme.The key result underpinning this is effectively an extension of Equation 19.In the low-aberration (i.e.non-scattering) regime, the sumcoordinate projection directly encodes the point-spread function of the system via where * denotes the convolution operation.This means that if we can find the phase mask on an SLM that optimises C + , it will also optimise h and restore imaging performance.Figure 4a-c shows results of an AO experiment using this QAO method.The experimental setup used to obtain these results is the same as that described in Ref. [44], where the SLM is located in an image plane of the crystal, and the camera is positioned in a Fourier plane (i.e.NF shaping configuration).This approach is interesting because, unlike many classical AO methods, it does not rely on choosing a good image metric and the performance is independent of the sample structure.
Additionally, it should be noted that the quantum correlations only need to be used as an optimisation target and the actual imaging can be fully classical, avoiding the long acquisition times typically involved in quantum imaging experiments.A detailed experimental setup can be found in Ref. [44].a, Reference intensity and sum-coordinate projection of G (2) images with no aberrations present.b, Intensity and projection images in the presence of aberrations.c, Intensity and projection images after the projection peak has been optimised.d Example phase mask retreived by optimising the projection peak.e-h, Example 2: Entanglement transmission through a scattering medium.A thin scattering medium is placed in the Fourier-plane of a BBO crystal.The transmission matrix of this medium is measured using a classical beam, and is used to compute a correction mask displayed on an SLM in a conjugate plane to the scatterer and the crystal surface (i.e.NF shaping configuration).Without the scatterer, the spatial entanglement of the pairs can be verified.With the scatterer, the entanglement can no longer verified.Applying the correction on the SLM allows measuring entanglement again after the scattering medium.d, Sum-coordinate projection without the scatter.e, Sum-coordinate projection of G (2) in the presence of the scatterer.f, Sum-coordinate projection with the scatterer and after correction.h, Example phase mask retreived by inverting the classically measured transmission matrix.A detailed experimental setup can be found in Ref. [30] In some imaging systems, aberrations become so complex that we enter in the regime of light scattering.Correcting optical scattering using classical light is a complex task [45,46] and it becomes even more challenging when using photon pairs.Indeed, each photon of the pairs is scattered in random directions, causing them to lose their correlations.Since these correlations are at the core of photon pair-based quantum imaging, preserving them is essential for these systems to function.This situation is explored in Ref. [30], which we show here as an example.In this work, the authors propagate a two-photon stateidentical to that generated in the experiments shown in Figures 1a and b -through a scattering medium (parafilm layer).The SLM is located in an image plane of the crystal, and the measurement is performed in a Fourier plane (i.e.NF shaping configuration).The experimental setup is described in full detail in Ref. [30].As we can see in the sum-coordinate projections shown in Figures 4d.e., the strong spatial correlations between the pairs measured in free-space (Fig. 4d) are lost in the presence of the scattering medium (Fig. 4e).
To restore the correlations, the authors leverage the concepts of transmission matrix based wavefront shaping, introduced in Ref. [46].The transmission matrix, T , is the discrete version of the point-spread function.In this formalism the two-photon state Ψ in is propagated with matrix multiplications.For a system with an SLM followed by the scattering medium, we have: where D is a diagonal matrix representing propagation through the SLM and Ψ out is the matrix associated with the two-photon state at the output.In this method, the transmission matrix is first measured classically using a coherent light source at 810nm, the same wavelength as the photon pairs.Then, once the transmission matrix is known, it can be used to compute a phase mask for the SLM that negates the effects of the scattering layer.This, in turn, restores the spatial correlations of the pairs.Figures 4e show the sum-coordinate projection measured after applying the correction mask on the SLM.
Both of these methods aim to recover spatial correlations by correcting for optical aberrations.In the first example, the correction is found by directly optimising the quantum correlations.This is possible because it targets the case when the aberrations are low-order (i.e.smooth), meaning the correlation signal (C + ) can be resolved in a relatively short time (∼ 2 minutes).If the aberrations become sufficiently high-order so that there is scattering, as is the case with the second example, then the acquisition times become prohibitive.Therefore, the correction is found via a classical transmission matrix measurement instead.

CONCLUSION
In summary, we have shown how the spatial correlations of an entangled two-photon state can be shaped with both a theoretical description and an experimental examples.Structuring such correlations is very similar to shaping intensity with classical light, but with some important differences.Depending on the position of the SLM relative to the crystal (i.e. the spatial basis in which we modulate the phase of the photons, position or momentum), the photon pairs 'see' a different version of the phase mask displayed on the SLM.For an SLM located in the position basis (i.e.conjugate plane to the crystal surface), correlations are modulated according to twice the phase mask.For an SLM positioned in the momentum basis (Fourier plane of the crystal), they are modulated by the combination of the phase mask and its spatial inverse.We demonstrate this experimentally, showing that there are in fact cases where we see no modulation of the correlations even when a non-flat phase mask is displayed on the SLM.
The aim of this article is to introduce spatial correlation-shaping between entangled photons with a simple example that has a well-known classical analogue.We provides the necessary technical details for readers to replicate the described experiments.We believe that these simple experiments can be expanded to inspire new methods of quantum imaging, much like what has been done in the examples shown in Figure 4.
The experimental setup is depicted in Figure 1.The pump is a collimated continuous-wave laser at 405nm (Coherent OBIS-LX) with an output power of 200mW and a beam diameter of 0.8±0.1 mm.The BBO crystal has dimensions 0.5×5×5 mm and is cut for Type I SPDC at 405nm with a half opening angle of 3 degrees (Newlight Photonics).The crystal is slightly rotated around the horizontal axis to ensure near-colinear phase matching of photons at the output.A 650 nm cut-off long pass filter is used to remove the pump photons after the crystal.This is important, as the pump can cause fluorescence in the subsequent lenses.A band-pass filter at 810 ± 10nm selects near-degenerate photon pairs.The SLM is a PLUTO-NIR-15 from Holoeye.It is a liquidcrystal-on-silicon device with a resolution of 1920×1080 pixels, and a pixel pitch of 8µm.The beam radius on the SLM is approximately 2.4mm, corresponding to 300 pixels.The EMCCD is the model iXon Ultra 897 from Andor, which has a resolution of 512x512 pixels with a pixel pitch of 16µm.The camera is operated with a region of interest (ROI) of 100 × 100 pixels.Far-field shaping configuration.In Figure 1a, for clarity, lens f 1 is depicted as a single lens; however, in reality it represents three lenses arranged in the confocal configuration 50 mm -150 mm -100 mm.The first lens is positioned 50 mm after the crystal, and the last lens is situated 100 mm before the SLM.Lenses f 2 − f 4 each have a focal length of 100 mm, with f 2 placed 100 mm after the SLM and f 4 placed 100 mm before the camera.
Near-field shaping configuration.
In Figure 1b, for clarity, lenses f 1 and f 2 are depicted as two lenses; however, in reality it represents four lenses arranged in the confocal configuration 45 mm -75 mm -50 mm -150 mm.The first lens is positioned 45 mm after the crystal, and the last lens is situated 150mm before the SLM.Lenses f 3 , f 4 and f 5 have a focal length of 100 mm, 50 mm and 75 mm, respectively, with f 3 placed 100mm after the SLM and f 5 placed 75 mm before the camera.
B. Measurement G (2) , C − and C + In each experiment of our work, we measure the spatially-resolved second-order intensity correlation function G (2) , and then use it to compute the sum and minus coordinate projections.G (2) takes the form of a 4dimensional matrix containing (N x × N y ) 2 pixels, where N x × N y corresponds to the region of the sensor used to capture data.An element of the matrix is written ijkl , where (i, j) and (k, l) are pixel labels corresponding to spatial positions (x i , y j ) and (x k , y l ).Assuming that the two-photon state detected by the camera is pure, G (2) ijkl can be obtained directly by measuring the covariance matrix, as detailed in Ref. [12].Such a matrix is measured by acquiring a set of M + 1 frames, denoted {I (l) } l∈[[1,M +1]] , using a fixed exposure time of 0.002s and then processing them using the formula: The first term in the sum is the covariance matrix between a frame and itself and corresponds to an estimate of the real coincidences from photon pairs.The second term is the covariance matrix between a frame and a subsequent frame.Since we know the true pairs will never be detected across two different frames, this corresponds to an estimate of accidental coincidences coming from uncorrelated photons being detected in the same frame.Since the EMCCD camera cannot resolve the number of photons incident on a single pixel, the photon coincidences at the same pixel cannot be measured, and so the corresponding values G ijkl is a discrete version of the continuous second-order intensity correlation function G (2) (r 1 , r 2 ) = |ϕ(r 1 , r 2 )| 2 where ϕ is the spatial two-photon wave-function associated with the photon pairs.Such a formalism has been employed in many studies describing the propagation of entangled photon pairs [42,47,48].Then, all the coordinate projections can be computed from G (2) .C − (r − ) is defined as where G ± is G (2) expressed using the variables (r + , r − ), such that G . Using the variables (r 1 , r 2 ), the equation becomes: Similarly, C + (r + ) is defined as: = G (2) (r, 2r + − r 1 )dr 1 .
In the manuscript, we use indifferently G ± and G (2) for clarity.In practice, C − is calculated using the discretevariable formula: and C + using: C. Example of G (2) Measurement in Practice The previous section gives a mathematical description of a G (2) measurement, but it is also useful to have a more practical description.
For simplicity, let us assume that we have N x = N y = N so each frame contains N 2 pixels.A typical acquisition requires around 10 million frames.It is not practical -and normally not even possible -to store this many frames for each acquisition, so instead we break the processing down into smaller blocks and continuously sum the results of each block.The real and accidental coincidence terms from Equation 24 are computed separately and the subtraction is done after the acquisition.The m th block results in two N 2 × N 2 arrays corresponding to these terms, called R (m) and A (m) respectively.These are summed for each block to get the final arrays so that R = m R (m) and A = m A (m) , and finally The following is a step by step breakdown of the acquisition procedure: 1. Initialise empty N 2 × N 2 arrays corresponding to R and A. R (m) and A (m) will be added to these arrays after each block, then discarded.
2. Fill the camera's internal buffer with frames and download them as a block to the processing computer.
3. When the buffer has been emptied begin the next acquisition to fill the buffer.At the same time, begin the processing on the most recent block of frames.See following for details on processing.
4. Add R (m) and A (m) to R and A and delete the current block of frames from memory.

Return to
Step 2 and repeat until the required number of frames have been processed.
The processing for each block is done as follows.For the m th block of M frames, labelled I (m) , will be in the form of an N 2 × M array where each column is a frame that has been unwrapped into a 1-dimensional vector.From this, R (m) can be computed via an outer product.In practise, this is simply a matrix multiplication with the array of frames: where X T denotes the transpose of X.A (m) is computed slightly differently.We define as I (m) with the last frame (i.e.column) removed.Similarly, we define as I (m) with the first frame (column) removed.I are both N 2 × (M − 1) arrays.Now, the accidental coincidence matrix is computed as .
The two terms are the covariance matrix of a given frame with the next frame, and a given frame with the previous frame.For N × N pixel frames, G (2) is actually an N × N × N × N element 4-d array.However, this processing gives G (2) in an N 2 × N 2 form, where each column/row corresponds to a 2-d conditional probability distribution that has been unwrapped into a 1-d column/row vector.This processing is done in parallel to the acquisition of the next block of frames.For an EMCCD camera, the processing is faster than the acquisition, so the camera speed is still the limiting factor.See Ref. [39] for a template MATLAB script to perform the acquisition and processing.At the end of the processing we typically save R, A, and the relevant projection.The intensity images taken from the camera are all discarded as they are processed.Therefore, we also typically save a total intensity image that is simply the sum of all of the frames that have been acquired.

D. Alignment
Aligning a photon pair-based imaging experiment can be difficult since the power of the down-converted light is on the order of tens of picowatts, making it invisible to the eye, and most cameras.In addition to this, it is essential to optimise the correlation width of the pairs, otherwise signal-to-noise will be prohibitively low.This is very sensitive to lens misalignment, meaning rougher alignment methods that may work in classical imaging systems such as simply positioning lenses with a ruler, will not be sufficient.
In this section we will briefly describe the alignment process.Due to the low flux, it is impractical to align the system directly with the SPDC light.Instead, it is easier to use a bright, spatially coherent source at the same wavelength as the pairs, e.g. a laser.Here, we use a superluminescent diode (SLED) that is spatially coherent, and has a spectrum of ∼ 800 − 820 nm.The same bandpass filter as the one used to filter photon pairs is used to select the desired wavelength at 810 nm from this range.
If this alignment beam is co-aligned with the SPDC pump laser, then co-propagating photon pairs will follow the same general trajectory, and so it can be used to align the optics after the nonlinear crystal.The lens directly after the crystal should be relatively short -i.e. have a sufficiently high NA -to ensure all of the highwavevector photons are collected.The position of this lens is critical, so mounting it on a translation stage for fine control is recommended.The position of the crystal itself is also key, so it is recommended that the crystal also be mounted on a stage.In general, the lenses should be placed starting from the camera, and working backwards, finishing with the short lens after the crystal.For a full description of the alignment process, see SM.

E. Calibrating the SLM
An SLM is an array of liquid crystal pixels whose birefringence is controlled by an applied voltage.Effectively, it is a display on which grayscale images can be show.The grayscale value at each pixel determines the voltage applied across the pixel, and therefore controls its birefringence.This variation in birefringence across the SLM imparts a spatially dependent phase to light that is incident on it.Ideally, the imparted phase will be a linear function of pixel value, with grayscale value of 0 giving 0 phase shift, and 255 giving 2π.In practise, this is not generally the case.The imparted phase is wavelengthdependent, so a grayscale range of 0-255 could correspond to a phase range of more or less than 0-2π.The response may also be non-linear.To account for nonlinearity, we calibrate the pixel response via the decorrelation of a reference speckle as a function of pixel grayscale value.Note that there are plenty of other methods to calibrate an SLM, such as the more common one described in Ref. [49].The one we describe here has the advantage of being less sensitive to alignment errors.
Two things are needed for this method: a spatially coherent source of illumination at the desired wavelength (e.g. a SLED), and a scattering medium.The coherent light must be (roughly) collimated at the SLM, and illu-minate a sufficiently large region.The scattering medium can be anything that produces a speckle pattern.A thin piece of semi-transparent plastic e.g. from a plastic sleeve, or a layer of parafilm stretched across a microcope slide works well.By placing the scatterer in the beam path, we measure a reference speckle at the camera.Then, by randomly selecting approximately half of the SLM pixels and scanning their grayscale value from 0 to 255, the resulting speckle images will be changed relative to the reference as a function of the actual phase imparted.One can demonstrate that the value of spatial correlation (corr2 function of Matlab) between this speckle and the reference speckle varies like the cosine of the global phase between the two speckles.With this information, a calibration curve can be found to determine the correct relationship between SLM gray scale values and optical phase values.For a full description of the calibration process, see SM. Example code and data can be found in [39].
(f) Finish by placing the lens f 1 .Generally, this lens has a short focal length, typically between 25 and 50mm.As stated in the Methods, this lens should be mounted on a translation stage for fine control of its alignment.More precise alignment will be done in the following step, but you should still aim for the best alignment possible by hand first.
FIG. 6. Example of a rotatable crystal mount.a, Image of full mount and translation stage.The crystal is mounted on a circular lens mount and post.This itself is mounted horizontally on a rotation mount to enable simple control over the crystal angle.Finally, everything is mounted on a (manually actuated) translation stage for precise control over the crystal position.b, Top-down image of the same crystal mount to show the crystal followed by a short-pass spectral filter and short focal length lens (f1).The filter is tilted to direct the reflection of the pump to a desired location.

4.
Precisely align lens f 1 and crystal.
The photon-pair correlations are particularly sensitive to the positions of lens f 1 and the crystal.The reason is that it is generally necessary to use a lens with a short focal length just after the crystal.Indeed, this allows for a sufficiently large numerical aperture to collect all the k-vectors emitted by the crystal.It is therefore important to align these well.
(a) First, choose between the lens with focal length f n or 0.5f n in order to image the front focal plane of lens f 1 onto the camera.This configuration is called near-field imaging configuration.Now we are imaging the plane in which we want to put the crystal.The crystal is transparent to the SLED, but it will likely have dust/imperfections on its surface.It can be put in roughly the correct plane by getting these imperfections in-focus.If there are no visible imperfections, another option is use a cross target or similar object that can be easily swapped with the crystal without moving the entire crystal mount.For an example of the mount we have used, see Figure 6.
(b) Now that all of the lenses and the crystal are positioned, the SPDC light should be visible on the camera.
In the following steps we will measure and use the photon-pairs spatial correlations directly to precisely adjust the positions of the lens and crystal.From here, remove as much background light as possible by covering the setup and turning off all other light sources.
(c) Now, put the setup in the far-field imaging configuration i.e. choose between the lens with focal length f n or 0.5f n in order to image the back focal plane of lens f 1 onto the camera.We will align the lens f 1 first since, in the far-field configuration, the correlation width does not depend on the distance between the lens and crystal.
(d) Tilt the crystal around its horizontal axis until a ring (or rings if the crystal is a Type II or paired Type I) is visible on the camera.Slowly tilt the crystal until this ring is collapsed almost to a disk, and fits into an area of at most 200 × 200 pixel.The exact size of this disk will depend on the camera being used but, if it is bigger than 200 × 200 pixels, it is likely too large and will slow down the correlation processing.In this case, it is recommended to reduce the magnification of your imaging system.
(e) Do a short G (2) acquisition (i.e. up to few minutes) and compute the sum-coordinate projection.There should be a peak in the centre.If there is a peak, simply find the position of the lens translation stage that maximises it.If there is no peak, move the lens and try again.If, after using the full travel of the stage, you still see no peak, then it is likely a problem with the alignment of other lenses.Check their alignment with the SLED and, if necessary, start again from step 3.
(f) Once the peak of the sum-coordinate projection in the far-field imaging configuration is optimised, move to the near-field imaging configuration.Now, repeat the process above for the position of the crystal instead of the lens, and optimise the peak in the minus-coordinate projection of G (2) instead.
5. Aligning the SLM.Ideally, the SLM should be exactly in the Fourier-plane of the camera.However, since it acts in reflection, this is not always easy to accomplish this.After aligning all lenses between the SLM and camera, choose between the lens with focal length f n or 0.5f n in such a way as to image a plane close to that of the SLM.Display a high-contrast mask on the SLM with many sharp edges (e.g. a 0/π grating).These sharp edges are visible due to diffraction effects, and can be used to position the SLM in the correct plane.Moving the SLM will change the distance between the lenses before and after it, introducing defocus error to the system.Therefore, it is best to use longer lenses directly before and after the SLM so that the depth of focus is long and defocus errors are reduced.
If you have completed all of these steps successfully, your system should be well aligned.The peaks in the near-field and far-field imaging configurations should be similar in size to those in Figures 1d,f  where A f inal and B are two constants.Indeed, the intensity at a position r on the camera results from the interference between a speckle s P generated by the passive part of the SLM and a speckle s A generated by the active part.When we phase-shift the active part with respect to the passive part by a global phase θ, then the intensity measured on the camera can be written as follows: I θ (r) = |s P | 2 +|s A | 2 +2|s P s A |cos (α A (r) − α P (r) + θ), where α P (r) and α A (r) are the phase components of s A (r) and S P (r), respectively.Calculating the correlation M between I θ and I 0 using Matlab's corr2 function involves spatially averaging the products I θ (r) and I 0 (r ′ ) for all pairs of positions r and r ′ .Assuming that each speckle is well developed, then phases α P and α A are randomly distributed between 0 and 2π across all the camera pixels, leading to the following results: + 2⟨[|s A (r ′ )| 2 +|s P (r ′ )| 2 ]|s A (r)||s P (r)cos (α A (r) − α P (r))⟩ r,r ′ + 2⟨|s A (r)s P (r)s A (r ′ )s P (r ′ )|cos (α A (r) − α P (r − α A (r ′ ) + α P (r ′ ))⟩ r,r ′ + 2⟨|s A (r)s P (r)s A (r ′ )s P (r ′ )|cos (θ)⟩ r,r ′ = A + Bcos(θ). (D2) If the SLM is perfectly calibrated, then θ = 255 2π G, leading to equation D1.However, in practice this is never the case, and the relationship between θ and G i.e. θ = f (G), is not so simple.It is precisely this function f that we are aiming to experimentally measure and determine here.
To achieve this, we then start form the correlation curve measured in Figure 8a.Firstly, it is necessary to ensure that the SLM implements sufficient phase shifting, i.e. if the cosine is cut off before reaching a maximum, then the pixels are not modulating all the way to 2π.Generally, an SLM will come with software to control the voltage that is applied across the pixels.If possible, use this to adjust the maximum applied voltage so you get one full oscillation, erring towards more than a full oscillation.Once a calibration curve with more than one oscillation is obtained, the curve is unlikely to be a perfect cosine because f is generally not linear.f is determined using the following procedure: 1. Rescale the data so that M ranges between -1 and 1. Figure 8b shows such a re-scaled curve (blue).
For the right-side model, if the fit was done with shifted data, then the coefficients can be redefined: with Y 0 = max(Y lef t ) ensuring the fits can be merged at G π .We use these models to create the function f (e.g. in the form of a Matlab function) that transforms a grayscale level G into its corresponding phase value.
7. Similarly, we need to create a function corresponding to f −1 .For this, the models must be inverted.

FIG. 1 .
FIG.1.Experimental Setup.Spatially entangled photon pairs centred at 810nm are produced via Type I spontaneous parametric down conversion (SPDC) using a collimated (0.8mm diameter), continuous-wave laser at 405nm and a thin β−Barium Borate nonlinear crystal (NLC).Pump photons are filtered out by a long-pass filter (LP) at 650nm.A bandpass filter at 810 ± 5nm before the single-photon sensitive camera filters out any non-degenerate pairs.a, Diagram of a far-field shaping configuration (FF) .The surface of the NLC is imaged onto the camera using two 4-f relays (f1 − f2 and f3 − f4).The spatial light modulator (SLM) is placed in the Fourier plane, or far-field, of the crystal.In this configuration, the photons are correlated at the camera, and anti-correlated at the SLM.b, Diagram of a near-field shaping configuration (NF).The Fourier plane (FP) of the NLC, obtained via the lens f1, is imaged into the camera with two 4f relays (f2 − f3 and f4 − f5).The SLM is placed in the conjugate plane to the surface of the NLC.In this configuration, the photons are anti-correlated at the camera, and correlated at the SLM.c,e, Direct intensity images from the camera in the FF shaping (c) and NF shaping (e) configurations.d, Minus-coordinate projection of the measured G(2) in the FF configuration.f, Sum-coordinate projection of the measured G(2) in the NF configuration.c,d are from an acquisition of 2.5 × 10 5 frames.e,f are from an acquisition of 6 × 10 6 frames.M -mirror.

FIG. 2 .
FIG. 2. Shaping the two-photon correlations in a far-field shaping (FF) configuration.A phase grating (a) is displayed on the illuminated region of the SLM.Line plot shows a single row of the grating.Red dashed circles are illustration of approximate position of SPDC ring on the SLM.Grating width is exaggerated for illustrative purposes.The actual period used in acquisition was 75 pixels.b, Magnitude of correlation peaks as a function of grating lateral offset β.Crosses are data points, solid lines are sine fits.Data is normalised for each order individually to allow better visual comparison.c-e, Minus-coordinate projections of G (2) for β = 0 (c), β = 35 pixels (d), and β = 70 pixels (e) peaks (cropped).Correlations for each grating offset are from an acquisition of ∼ 1.8 × 10 6 frames.

FIG. 3 .
FIG. 3. Shaping the two-photon correlations in a near-field shaping (NF) configuration.a, Phase grating displayed on the illuminated region of the SLM.Line plot shows a single row of the grating.b, Magnitude of correlation peaks as a function of grating amplitude α.Crosses are data points, solid lines are sinusoidal fits.Data is normalised for each order individually to allow better visual comparison.Classical data (purple) has scale of intensity, rather than correlations.c,d Sum-coordinate projections of G (2) for grating amplitude α = π/6 (c), α = π/2 (d) and α = 2π/9 (e).Correlations for each grating amplitude are from an acquisition of ∼ 5 × 10 5 frames.For details on the acquisition of the classical data, see S.M.

FIG. 4 .
FIG. 4.Examples of applications of two-photon shaping.a-dExample 1: Adaptive optics with entangled photons.The sample, a mosquito pupa, is illuminated with photon pairs and imaged onto an EMCCD camera.The direct intensity (large grayscale images) and sum-coordinate projection (inset images) are measured.A layer of transparent polymer (polydimethylsiloxane) placed after the sample is used to introduce aberrations, which distorts the intensity images and projection.An optimisation algorithm is used to find the phase mask on an SLM that maximises the projection central value, which also restores imaging performance.The SLM is positioned in the conjugate image plane of the crystal (i.e.NF shaping configuration).A detailed experimental setup can be found in Ref.[44].a, Reference intensity and sum-coordinate projection of G(2) images with no aberrations present.b, Intensity and projection images in the presence of aberrations.c, Intensity and projection images after the projection peak has been optimised.d Example phase mask retreived by optimising the projection peak.e-h, Example 2: Entanglement transmission through a scattering medium.A thin scattering medium is placed in the Fourier-plane of a BBO crystal.The transmission matrix of this medium is measured using a classical beam, and is used to compute a correction mask displayed on an SLM in a conjugate plane to the scatterer and the crystal surface (i.e.NF shaping configuration).Without the scatterer, the spatial entanglement of the pairs can be verified.With the scatterer, the entanglement can no longer verified.Applying the correction on the SLM allows measuring entanglement again after the scattering medium.d, Sum-coordinate projection without the scatter.e, Sum-coordinate projection of G(2) in the presence of the scatterer.f, Sum-coordinate projection with the scatterer and after correction.h, Example phase mask retreived by inverting the classically measured transmission matrix.A detailed experimental setup can be found in Ref.[30]

( 2 )
ijij are set to the mean of the neighbouring values.G .

FIG. 7 .
FIG. 7. Speckle patterns.a, Reference speckle pattern measured with a flat phase programmed on the SLM.b, Speckle pattern measured with random subset of SLM pixels set to the gray value G = 117.Images were acquired with with Thorlabs Zelux 1.6 MP Monochrome CMOS camera.The SLM is an Holoeye Pluto NIR-II.

2 . 3 .
Record the gray values G of the first maximum, labelled G 0 ; the minimum, labelled G π ; and the second maximum, labelled G 2π .Record the G values where M = 0 (or closest).Label the lower G π/2 and higher G 3π/2 .See Figure8b.Split the data in half at G π , so we haveM lef t = M (G lef t ) for G lef t ∈ [G 0 , G π ] and M right = M (G right ) for G right ∈ [G π + 1, G 2π ].

4 .
Compute Y lef t = arccos(M lef t ) and Y right = −arccos(M right ) + 2π, and plot Y vs G, as shown in Figure 8c.This is the pixel response of the SLM i.e. the function f .

5 .
Fit Y lef t and Y right as functions of G lef t and G right , respectively.It is usually sufficient to use quadratic polynomials for the models.In the example code, to improve the fit consistency, we shift the right data so that the first point is at the origin.That is, we fit Y ′ right = Y right − π as a function of G ′ right = G right − G π .The results can be re-shifted after the fit.Figures 9ab show examples of fitted models.6.Now you have two models describing the pixel response from 0 to G π and G π to G 2π .If they are quadratics, they are of the form If they are quadratics, thenG model = −a 2 + a 2 2 − 4a 1 (a 3 − Y ) 2a 1 .(D5)Now you can define a function that computes G model−lef t for Y ∈ [0, π], using the coefficients for the left side fit, and G model−right for Y ∈ [π, 2π] for the right-side fit.See Figure9cfor the final combined model.

FIG. 9 .
FIG. 9. Quadratic fits of the measured SLM pixel response.a,b, Data points and fitted functions for the left and right data, respectively.c, Combined inverse function that maps the target phase to the grayscale pixel value that gives this phase.