Biphoton State Reconstruction via Phase Retrieval Methods

The complete measurement of the quantum state of two correlated photons requires reconstructing the amplitude and phase of the biphoton wavefunction. We show how, by means of spatially resolved single photon detection, one can infer the spatial structure of bi-photons generated by spontaneous parametric down conversion. In particular, a spatially resolved analysis of the second-order correlations allows us to isolate the moduli of the pump and phasematching contributions to the two-photon states. When carrying this analysis on different propagation planes, the free space propagation of pump and phasematching is observed. This result allows, in principle, to gain enough information to reconstruct also the phase of pump and phasematching, and thus the full biphoton wavefunction. We show this in different examples where the pump is shaped as a superposition of orbital angular momentum modes or as a smooth amplitude with a phase structure with no singularities. The corresponding phase structure is retrieved employing maximum likelihood or genetic algorithms. These findings have potential applications in fast, efficient quantum state characterisation that does not require any control over the source.


I. INTRODUCTION
The tomography of quantum states plays a fundamental role in modern quantum technologies [1][2][3][4][5].At the same time, such a task can be particularly challenging when dealing with systems of many particles and/or many degrees of freedom [6,7].If a projective measurement approach is adopted, with no prior information, the number of required measurements scales quadratically with the dimensionality of the Hilbert space [1].The case of the reconstruction of high dimensional two-photon states is of particular interest given their applications in fundamental quantum mechanics [8][9][10][11], highdimensional quantum communications [12,13], and quantum imaging [14].As an alternative to projective approaches based on mode sorting, recent works introduced the possibility of exploiting variations of classical interferometric techniques for directly reconstructing the spatial structure of phase and amplitude of an unknown biphoton state [15].This approach is practically faster and more reliable than traditional methods thanks to the use of time-stamping cameras [16][17][18], which is proving to be a promising resource for quantum optics experiments [15,[19][20][21][22][23][24][25][26][27].However, interferometric approaches will require the phase locking of the reference biphoton state with the unknown one, a task which can be harder to achieve in cases where the unknown source is not easily accessible.It is thus desirable to have the possibility of obtaining the phase structure of the unknown biphoton from measurements that do not require any control of the source.Here, we propose a method based on measuring the spatial coincidence distribution in different propagation planes and reconstructing the phase of the biphoton employing phase retrieval algorithms.
We demonstrate how the analysis of coincidence images allows us to extract the square moduli, hereafter referred to as "intensities", of the two complex functions that contribute to the biphoton state: the pump and the phasematching.Intriguingly, this separation can be performed in any propagation * aderrico@uottawa.caplane (assuming free space propagation and eventual lenses that can be present in the setup).The intensity of the extracted pump or phasematching function in a given plane is related to the respective fields in another plane at a distance d by paraxial propagation through a distance d/2, or assuming the pump and phasematching functions to have half the wavelength of the biphoton state.Retrieving pump and phasematching intensities at different planes and knowing their relationship opens the possibility of reconstructing the phase of these fields and, thus, the full biphoton wavefunction.We demonstrate theoretically and experimentally how to isolate the pump and phasematching intensities from coincidence images, and we give some examples of their phase reconstruction based on the exploitation of maximum likelihood approaches -in the case a decomposition in orthogonal mode sets can be conveniently used-or genetic algorithms -in the case of aberrated pumps with smooth amplitudes-.

II. THEORY
Consider a two-photon state with fixed wavelengths and polarisations, in a paraxial propagation regime.In this case, assuming the z-axis corresponding to the mean propagation direction, the system can be described in terms of the probabilities of photons passing through the transverse position X = (x, y) in a given plane.The quantum state is thus described in terms of the creation operators â † X .Defining the transverse position eigenstate as |X⟩ := â † X |0⟩, with |0⟩ the vacuum state, a general biphoton state is given by where Ψ(X i , X s ) is called the biphoton wavefunction.The indexes i and s denote, conventionally, idler and signal photons, respcetively.Equation ( 1) is a representation of the state in a basis of transverse coordinates eigenstates for a given value of z.The free space propagation of the biphoton wavefunction  from z to z ′ is where G ∆z (X ′ , X) is a free space propagator through the distance ∆z = z ′ − z and X ′ = (x ′ , y ′ ) are the transverse coordinates in the plane z ′ .For paraxial fields, the free space propagation is described in terms of the Fresnel propagator A common and interesting case is correlated biphoton states with the following wavefunction, In particular, in the case of biphotons generated via spontaneous parametric down-conversion (SPDC) in Type-I crystals, E p is the spatial, slowly varying amplitude of the pump laser, and ϕ is the phasematching function [28,29].This structure simplifies the propagation of the biphoton wavefunction significantly.Performing the change of variables R = (X i + X s )/2, ∆ = (X i − X s )/2 [24], and expanding the squares (X ′ − X) 2 in the Fresnel propagator, it is straightforward to verify that: This allows one to separate the 4-dimensional (4D) integral describing the SPDC propagation in free space into the product of two bi-dimensional integrals: where N is a normalisation constant.Thus, the propagated biphoton wavefunction has the structure , where E ′ p and ϕ ′ are obtained applying a single Fresnel propagator on the functions E p and ϕ.It is important to emphasize that the propagation of E p and ϕ must be evaluated for ∆z/2 if ∆z is the propagation distance considered for the biphoton state.
Experimentally, from a spatially resolved coincidence detection, one can extract the pump and phasematching contributions in each measurement plane.Coincidence detection allows to retrieve the 4D function C From this data set, the pump and phasematching squared moduli can be extracted when post-selecting on spatially correlated and anticorrelated states, respectively.This is done by evaluating the quantities and Apart from a constant shift c, the post-selected coincidences give patterns that are spatially distributed as the pump and the phasematching intensity.Here, we considered a generic case in which correlations and anti-correlations are chosen with a constant offset c.This allows for increasing the statistics without the need for long exposure times in the experimentsee the following section and Supplementary Figure S1.
Using the technique described, we can extract the individual functions of the pump E p and phasematching ϕ that contribute to the biphoton state.We can then observe how these functions propagate through free space.By examining the distributions of the pump and phasematching functions in different planes, it is possible to obtain information on their phase structure without the need for interferometric or tomographic measurements.

III. EXPERIMENTAL RESULTS
We observed the free space propagation of the pump and phasematching contribution of an SPDC state generated in a 0.5 mm thick Type-I BBO crystal pumped by a pulsed 405 nm laser beam (pulse duration 150 fs, repetition rate 80 MHz).A Liquid Crystal Spatial Light Modulator (LC-SLM) was used to shape the spatial structure of the pump laser.A conceptual scheme of the experiment is sketched in Figure 1 a, and a detailed setup is described in the Methods.Idler and signal photons can be spatially separated -with 50% probability-by a non-polarising beamsplitter and imaged on two regions of a time-stamping camera (Timepix3D) [16,18,30].If the spatial distributions of the idler and signal do not overlap on the camera sensor, one can choose two regions of interest (ROI), defined as the groups of pixels hit by the idler and signal photons, respectively.The time-stamping resolution of the Timepix3D is less than 10 ns and allows to extract coincidences between pixels contained in the ROIs, C(X i , X s ).From the 4-dimensional coincidence distribution is possible to extract marginals, e.g. the spatial correlations along x: C x (x i , x s ) = yi,ys C(x i , y i , x s , y s ).Alternatively, 2D sections of the 4D coincidence pattern can be extracted.As shown by Eqs. ( 7) and (8), this allows to obtain the intensity distributions , where the constant c can be chosen at will.To fully exploit the accumulated data, one can sum up the intensities obtained for each value of c after appropriately shifting each distribution by ±c.In this way, it is possible to achieve smooth reconstructions from data obtained with a few minutes of exposure (see Supplementary Materials, Figure S1).A first example of the reconstruction is shown in Fig. 1 b and c, for the case of a Gaussian pump.The analysis was carried out in two different propagation planes at distances z 1 = 7.5 cm and z 2 = 34.5 cm from the image plane of the crystal.The spatial correlations show the evolution from a spatially correlated state (signal and idler photon are localised in roughly the same transverse position) to a spatially anti-correlated state (signal and idler photons are localised in opposite transverse positions with respect to the biphoton state propagation axis).The correlations do not appear particularly sharp in these two planes.However, they show width of the order of a single pixel in both the crystal image plane and the corresponding Fourier plane.The extracted pump intensities display the expected Gaussian distribution (slightly narrower in z 2 due to a wavefront curvature of the pump laser on the crystal plane).The extracted phasematching intensity shows the formation of the characteristic SPDC cone in a collinear phasematching configuration.
It is well known that measuring the intensity distribution of a coherent beam in two different propagation planes can be used to retrieve the phase.This non-interferometric phase retrieval approach was first proposed by Gerchberg and Saxton [31] for the case where the two intensity distributions are measured in conjugate planes (i.e. the corresponding fields are one the Fourier transform of the other).Unfortunately, doing this for an SPDC state implies that, for instance, the reconstructed phasematching intensity will have a width of one or a few pixels in the first recording plane (and a similar effect will be observed for the pump distribution in the far field).This is due to the fact that the spatial extension of pump and phasematching functions are extremely different on these two planes.One can then think of looking for a compromise and choose two intermediate planes sufficiently far away from each other but where the spatial extension of the pump and phasematching are of the same order.However, in these cases, the Gerhberg-Saxton (GS) -as well as variations like the Fienup algorithm [32,33]-tends to have a worse convergence and be sensitive to noise.This is acceptable for applications where a phase pattern is generated to achieve target intensities but is less practical for phase retrieval since the GS can easily Due to these difficulties, we focused on retrieving the phase structure of phasematching and pump using optimization algorithms.However, these algorithms relied on some assumptions about the modal structure of the measured quantities.This is not a big issue since the SPDC physics for thin crystals has been extensively studied, and SPDC spatial mode structures are expected to follow the field continuity.
As a first step, we consider the phase retrieval of the phasematching function ϕ.In Fig. 2-a and b, we report the measured intensity |ϕ(X, z)| 2 in two different planes.ϕ(X, z) propagates essentially as a strongly diverging beam.One can expect -as also predicted by the thin-crystal theory of SPDC-that the phase of ϕ in a given plane is quadratic: arg ϕ(X, z) = πX 2 /(λz).When applying this phase structure on the amplitude measured in z 2 and numerically propagating back to z 1 , one obtains the intensity displayed in Fig. 2-c, which is in good agreement with the experimental result.The reconstructed phasematching function is shown in Fig. 2-c.The same phase structure was obtained via an optimization algorithm based on the decomposition of the phasematching function in orthogonal modes.We describe this process of reconstruction of the pump phase in detail below.
Pump lasers prepared in spatial structures different from a single Gaussian have been considered, e.g. to shape correlations in the Orbital Angular Momentum (OAM) degree of freedom [7] or to control multimode Hong-Ou-Mandel bunching or anti-bunching [34].In the latter case, one considers a pump field that is an asymmetric function of one transverse coordinate.This is obtained by applying a π-phase jump on an input Gaussian beam. Figure 3-a shows the pump intensity contribution to the SPDC state in two different propagation planes.To retrieve the phase structure of the corresponding field, we consider its approximation to a finite superposition of paraxial modes f κ (X, z), Here, κ is, in general, a set of discrete indices.If the pump intensity in two different planes is I i (X) = |E p (X, z i )| 2 , with i = 1, 2, then the optimal modal decomposition can be found minimizing the functional The choice of the modes set and the range of values for κ varies case-by-case.Typical choices can be the Hermite-Gauss or Laguerre-Gauss sets [35].However, spatially shaped beams are often generated by applying a phase pattern with line or point singularities on an input Gaussian beam.This typically requires a large number of coefficients in the decomposition equation (9).In this situation, another convenient choice is the over-complete set of Hypergeometric-Gauss modes HyGG −|ℓ|,ℓ (X) [36] -the detailed expression used in this work is given in the Methods.We recall that the index ℓ refers to the Orbital Angular Momentum (in units of ℏ) carried by the corresponding mode.The coefficients c ℓ will thus give the OAM spectrum of the analyzed field.shows the field intensity of the decomposition (9) obtained minimising Eq. ( 10) where f κ → HyGG −|ℓ|,ℓ .The excellent agreement with the experimental data is also highlighted in Fig. 3-c, where the experimental amplitude distribution along y, for fixed x = 1.9 mm, is compared with the reconstructed one.Figure 3-c shows the reconstructed field in z = 0.A phase variation of ∼ π between the two intensity maxima can be observed, together with an overall smooth quadratic phase due to the imperfect collimation of the pump on the crystal.In Fig. 3-e, the OAM power spectrum is reported, showing how the main contribution comes from two OAM modes with ℓ = ±1.The smaller ℓ = 0 contribution (associated with a Gaussian mode) is due to a residual misalignment of the SLM phase mask.In Fig. 4, we consider a pump created in an unbalanced superposition of OAM modes.Again, a finite decomposition in 15 HyGG modes (with ℓ = −7, . . ., 7) yields a reconstruction that matches well the experimental data.
The use of a finite set of orthogonal modes can be less efficient for fields with no phase singularities.In most applications SPDC is generated by a pump in the fundamental mode of the laser cavity, however, the actual phase and amplitude can be altered by imperfections in the experimental setup.Thus, one can expect that the pump contribution has a smooth phase factor that can be expanded in Zernike polynomials.We provide proof of principle for these applications by introducing specific optical aberrations on the Gaussian pump with a UV-SLM.Figure 5 a-d shows the extracted pump shapes in the crystal image plane z 1 = 0 and at z 2 = 19 cm for cases in which a coma and a second order astigmatism phase was introduced on the pump.We assumed the phase ξ p (X, z = 0) to be a superposition ξ p (X, z = 0) = n,m γ n,m Z m n (X), where Z m n (X) are Zernike polynomials [37] and γ n,m are real coefficients, while the amplitude is given by the square root of the experimentally retrieved intensity.Due to the lack of analytical expression for the propagation of aberrated modes, we relied on the use of a genetic algorithm instead of the maximum likelihood approach used in the previous examples.Random choices of γ n,m were used to define the individuals that initialize the genetic algorithm.A numerical Fresnel propagation approach was used to calculate the fields in z 2 resulting from the different individuals and compare the intensity with the experimental one.The function in Eq. ( 10) was here used as the Fitness function of the genetic algorithm (details are given in the Methods section and the Supplementary materials).Panels f and h of Fig. 5 show the best-reconstructed pump field at z 2 -which is in good agreement with the experiment-, and the corresponding amplitude and phase in z 1 are shown in panels e and g.

IV. DISCUSSION AND CONCLUSIONS
In conclusion, we have shown a powerful application of coincidence imaging of biphoton states.Spatially resolved second-order correlations allow the extraction of information about the two main physical contributions to the SPDC biphoton states, the spatial structure of the pump beam and the phasematching function, which is determined by the physical properties of the nonlinear crystal used for the biphoton state generation.The intensity of these two functions can be extracted at any distance from the crystal, and the relationship of the obtained intensities at different distances is given by an appropriate paraxial propagation.We exploited these results to extract, employing optimization methods, the phase of the two investigated functions and, thus, the full biphoton state.It must be stressed that this high-dimensional state reconstruction requires only two spatially resolved coincidence measurements, which, thanks to modern time-stamping cameras, can be performed in a few minutes without any control of the biphoton source.However, the renunciation of interferometric methods, and thus of direct phase measurements, comes at the expense of designing the proper algorithm to find the best phase structure that describes the experimental results.For smooth phases, if the field intensities are retrieved in propagation planes too close to each other, there will be a higher uncertainty in the reconstructed phase.Moreover, the separation of pump and phasematching can be rigorously achieved in the thin crystal limit, when the transverse walk-off between pump and downconverted photons is negligible.The more general scenario will require more robust approaches to analyze the 4D secondorder correlations in different propagation planes and extract the full phase patterns.We expect that this problem could be tackled by utilizing properly trained Neural Networks.Lastly, the ability to separate pump and phasematching contributions opens new opportunities for quantum imaging applications, some of which will be explored in future works.
V. METHODS

A. Experimental Setup
The experimental setup is described in detail in Fig. 6.A 405 nm pump laser is generated as the second harmonic of a Ti:Sa pulsed laser (Chameleon Vision II).Then, a magnifying system of lenses with a pinhole in the beam's focus is used to generate the desired Gaussian beam that is then sent to the Spatial Light Modulator (SLM).The beam's phase and amplitude are structured using the amplitude-phase masking technique [38], which requires appropriate phase masks displayed on the SLM and selecting the first diffraction order.The latter is achieved by placing an iris in a de-magnifying system of lenses after the SLM.Sending the resulting structured pump beam through the BBO-Type I crystal, the SPDC is generated and collimated by a 75 mm lens.After the crystal, the pump beam is filtered by a low-pass filter.Signal and idler pho-tons are separated with 50 % probability, by a half wave-plate and a polarizing beamsplitter (PBS), instead of using a nonpolarising beamsplitter.This approach guarantees that the two photons are now orthogonally polarised.Using two mirrors and a PBS, one can simultaneously match the optical paths of the two photons and propagate them along parallel paths, with a relative displacement smaller than the camera sensor.A bandpass filter (centred at 810 nm and with 10 nm bandwidth) is mounted on the camera intensifier to select frequency degenerate photons.The Timepix camera collects spatially resolved time stamps with ∼ 1 ns resolution [18] from which the 4D coincidence distribution is extracted.

B. Hypergeometric-Gaussian modes
Hypergeometric-Gaussian modes HyGG p,ℓ with p = −|ℓ| are a set of paraxial modes whose expression in the waist plane is given by HyGG −|ℓ|,ℓ (r, ϕ, z = 0) ∝ exp −r 2 /w 2 exp(iℓϕ), where (r, ϕ, z) are cylindrical coordinates, ℓ an integer number, and w a real parameter corresponding to the waist radius.Here, we consider the more general case where a quadratic phase is added on the waist plane: where J |ℓ| (.) is a cylindrical Bessel Function of order |ℓ| and 1/σ := 1/ξ − 1/ζ.The solution to the last integral (see Ref. [39]) can be expressed in terms of modified Bessel func-tions:

C. Details of genetic algorithm
Genetic algorithms evolve a population of candidate solutions toward optimal solutions to the problem of minimizing a cost function.In our case, the cost function is given by Eq. ( 10).The analytical expressions for the propagation of families of optical modes at a finite distance z, such as Hermite-Gauss or Laguerre-Gauss modes, enable the adoption of standard minimization routines.This is not possible in the case of aberrated beams, where there is no analytical formula for the propagated field.Genetic algorithms offer a suitable alternative, as these iterate from random guesses that evolve to the physical solutions of the optimization problem [40].In the following, we detail the sequence of operators used in our genetic algorithm.An initial population of N individuals is randomly generated within the range [−20 , 20], where each individual is a set of real coefficients γ (genes) for the chosen Zernike polynomials.In groups of k, these individuals compete for the possibility to reproduce.Only the individual better minimizing the cost function within each pool is given access to the reproduction stage.This is the so-called tournament mechanism, working as a selection operator [41].The reproduction occurs in the form of blend crossover, largely employed to mate real-valued individuals [42].When two individuals reproduce, two newborns originate as a weighted mixture of the parents.To emulate the mutation of individuals in a natural environment, genetic mutations are included in the workflow as Gaussian noise with mean µ and standard deviation σ, potentially affecting each gene of newborn individuals [43].Blend crossover and mutation are non-deterministic operators, occurring with probability p c and p m , respectively.To push the algorithm to a faster convergence, our algorithm is equipped with elitism, i.e., the best individual from the parent generation is guaranteed a place in the next one, replacing the worst individual of the offspring.The algorithm ends when a certain condition is verified [44].In our implementation, the maximum number of generations N gen is adopted as termination criterion.The user defines by hand all the parameters determining the evolutionary sequence (also called hyperparameters).The complete set of hyperparameters used in our algorithm is listed in Table I.

Figure 1 .
Figure 1.Experimental layout and coincidence analysis.a-Simplified sketch of the experimental setup.A pump laser (central wavelength 405 nm) is incident on a 0.5 mm thick BBO Type-I crystal that generates photon pairs with the same polarisation.The photons propagate through an imaging system (depicted here as a single lens L; see Methods for the detailed experimental layout).Signal and idler photons are spatially separated with a 50/50 beamsplitter (BS) and impinge on a time-stamping camera (Tpx3D), which allows for retrieving spatially resolved coincidence counts.The camera was moved along the propagation direction to measure spatial correlations in different planes and extract the corresponding pump and phasematching function intensities.An example with a Gaussian pump is shown in panels b and c.

Figure 2 .
Figure 2. Full reconstruction of phasematching function.a-b.Experimentally extracted intensity of the phasematching function of an SPDC state in two different propagation planes.c.Reconstructed phasematching intensity in z1 assuming the quadratic phase shown in panel d and applied on the experimental amplitude in z2.The intensity in z1 is obtained by Fresnel propagation.The excellent agreement with the experimental intensity in z1 indicates the correctness of the assumed phase structure.

Figure 3 .
Figure 3. Phase reconstruction of the anti-symmetric state.a. Experimentally extracted pump intensities for an anti-symmetric SPDC state.b,c.Reconstructed intensities from an optimal superposition of HyGG modes.Panel c shows the comparison with the experiment (gold-coloured triangles) along the line x = 1.6 mm.d.Reconstructed phase and amplitude of the pump contribution in z1.e.The OAM power spectrum of the pump, which was obtained from the reconstruction.

Figure 4 .
Figure 4. Phase reconstruction of OAM superposition states.a. Experimentally retrieved pump intensity in a superposition of, nominally, OAM=4 and OAM=-2.b.Intensities of the reconstructed field obtained by superimposing 15 HyGG modes.c.Retrieved phase and amplitude of the pump contribution at the crystal image plane.d.The OAM power spectrum of the reconstructed field.

Figure 5 .
Figure 5. Phase reconstruction of aberrated Gaussian pumps.Panels a to d show the extracted pump shapes in two propagation planes, z = 0 and 19 cm.For a-b, a coma aberration was imposed through an SLM on the pump laser.For c-d, instead, second-order astigmatism was applied on the pump laser.Panels e and g show the reconstructed phase and amplitude distributions in z = 0.In contrast, f and h show the intensity obtained by numerically propagating the fields in e and g to z = 19 cm (to be compared with panels b and d, respectively).The phase distributions in e and g are shown with the tip and tilt contributions removed since these are due only to the imperfect centring of the intensity patterns in the two planes.

4 Crossover probability pc = 0. 9 Figure S1 .
Figure S1.Illustration of pump and phasematching function extraction.The process of pump and phasematching intensity extraction is illustrated.The first two rows show the spatial correlations with a 1-pixel wide green band which indicates the region of coincidence postselection.The third row shows the coincidence distribution after the considered postselection.It is possible to observe how the spatial shift of this distribution is affected by the position of the postselection region.After compensating for this shift the results can be accumulated to obtain a less noisy distribution of the pump and phasematching intensity, as shown in the last row.

Figure S2 .
Figure S2.Detailed results of Genetic algorithm.Intensities of the best individual at each generation for a, coma, and b, second-order astigmatism.c and d show, respectively, the experimental results used as a target of the genetic algorithm.In panels e and f we show the evolution of the similarity S between the best individual and the target intensity.Panels g and h show the Zernike decomposition corresponding to the best individual of the last generation.