Spatial Wavefunction Characterization of Femtosecond Pulses at Single-Photon Level

Reading quantum information of single photons is commonly realized by quantum tomography or the direct (weak) measurement approach. However, these methods are time-consuming and face enormous challenges in characterizing single photons from an ultrafast light source due to the stringent temporal mode matching requirements. Here, we retrieve the spatial wavefunction of indistinguishable single photons from both a continuous wave source and a femtosecond light source using a self-referencing interferometer. Our method only requires nine ensemble-averaged measurements. This technique simplifies the measurement procedure of single-photon wavefunction and automatically mode matches each self-interfering single photon temporally, which enables the measurement of the spatial wavefunction of single photons from an ultrafast light source.

intensified charge-coupled device (ICCD) camera (model iStar DH734). The gated pulse that activates the intensifier tube of the ICCD camera is set to the period of the ultrashort pulses.
The probability of two photons landing on the same pixel is negligible. The acquired interferograms are accumulations of tens of thousands of frames. The interferograms with shearing amount of s x = 0, s x = 300µm, s y = 0, and s y = 300µm at phase-shift of 0, π/2, π and 3π/2 are obtained by rotating and translating the CDSI. The phases of these interferograms are manipulated numerically to extract the wavefront shape of the incident laser beam, which is at a single photon level unless stated otherwise. For a subset of the full wavefront information, in particular the odd order of x of the wavefront, only a subset of the mentioned interferograms are needed. After the data acquisition, they are processed using Matlab, a numerical computing software. The details of the formulas and the calculations involved for the wavefront extraction are provided in the 'Wavefront extraction' subsection.

Photon source
The photon source is a mode-locked oscillator of the Model MTS Mini Ti-Sapphire Laser Kit that generates ultrashort pulses with an 800 nm central wavelength, a 56 nm bandwidth, and 106 fs pulse duration at a 89 MHz repetition rate. The beam is attenuated by 9 orders of magnitude with ND-filters to the photon rate of ∼1 photon per pulse.

Imaging with ICCD camera
The ICCD camera used is of the model iStar DH734. Even though the minimum exposure time of the camera is 1 ms, the time resolution of this ICCD camera is determined by the 25 ps time resolution of the gated pulse that controls the intensifier tube. The gated pulse is set to be the pulse to pulse separation of the laser, which is 11.225 ps ≈ 1/89 MHz. The laser is attenuated such that the expected photon number in a pulse is approximately 1. The ICCD is thermoelectrically cooled for noise reduction. Single photons can be resolved by increasing the gain sufficiently high to surpass the CCD readout noise while setting a threshold to be above the CCD readout noise in the photon counting mode. The image readout mode is used with a sub-region selected just large enough to capture the output beams. The vertical and horizontal binnings are set to be 8 to generate 'superpixel' array that sums data on-chip prior to readout and significantly reduce the chance for a single photon event to register in more than one pixel.
This, in combination with sub-region-selection, increases the frame rate of the ICCD detector.
The data are acquired in accumulation mode with the mentioned settings. A video showing the typical shearing interferogram made up of accumulation of 20, 40, 60, 80, 100, 200, 500, 1000, 5000, and 15000 frames is found in Visualization 2. The camera frame rate in this experiment is about 20 frames per second, so each set of measurement with 15000 frames takes about 750 seconds. The measurement time required for our method scales with the dimension d of the state, or the number of pixels.

Probability of undercounting
Undercounting occurs when two photons land on the same pixel in a single camera frame. This can be estimated using conditional probability of Prob(≥ 2 photons in 1 frame) × Prob(2 photons lands on same pixel|2 photons in 1 frame). The pulse energy of the oscillator is attenuated to 0.72 photons per camera frame for ultrashort pulses and 1.04 photons per camera frame for the CW beam. For a source that follows the Poisson distribution at such intensities, the probability of having 2 or more photons in a single frame is 1 − Prob(0 photon) − Prob(1 photon) = 0.1628 for ultrashort pulses and 0.2790 for the CW beam. The gate pulse of the ICCD camera is set to the pulse-to-pulse separation of the oscillator. Therefore, each frame captures exactly one pulse. We now need to estimate the probability of two photons landing on the same pixel. According to the probability distribution of the beam profile measured in Fig. 2 of the ultrashort pulse, the pixel with the highest probability has a probability of 0.004, which is an overestimation of two photons landing on the same pixel given two photons in 1 camera frame. Similarly, this value is approximately 0.002 for the CW beam as shown in Fig. S7.
Therefore, the probability of having two or more photons in a camera frame and two photons landing on the same pixel must be less than 0.1628 * 0.004 = 0.00065 for ultrashort pulses and less than0.2790 * 0.002 = 0.00056 for the CW beam.

Operation principle of CDSI
As shown in figure 1 in the main manuscript, the final wavefunction at exit face 3 is the superposition of the left reflected portion and the right transmitted portion of the initial wavefunction.
The optical path difference between the two portions can be written as where n and n air are the refractive index of the BSC and air respectively, t(y) = 2πy sin α/λ + t(0) ≈ 2πyα/λ + t(0) is the thickness difference measured in radians of the two entrance faces resulted from the y-wedge angle α of face 2, and s is the shearing amount. The function W (x, y) is the wavefront of the incident beam, which can be separated into odd and even order terms and then the terms containing the wavefront become W e (x, y) − W e (x + s, y) To simplify the terms, we taylor expand the sheared wavefront as follows: The shearing amount s is chosen to be small for wavefront retrieval, so we keep only the two highest order terms of s as such: Phase-shifting of the CDSI

Michelson vs. CDSI on vibration effects
As mentioned previously, in comparison to the Michelson interferometer, our CDSI with this BSC is 3 orders of magnitude less vulnerable to vibration in the direction that causes phase shift.
Visualization 3 shows the direct comparison of the vibration effects on Michelson interferometer vs. CDSI. Both interferometers were mounted on speakers that run in series with a function generator as a source. The interferogram of the Michelson interferometer immediately vanishes when vibration is introduced while that of the CDSI had no visible difference. Vibration on the orthogonal directions has no effect for Michelson interferometers. For CDSI, vibration in the beam propagation direction has hardly any effect on the interferogram. Vibration in the shearing direction will change the shearing amount s and affect the interferogram that depends on the wavefront shape according to Eq. (2) in the main manuscript.

Avoiding unwanted diffraction effect
A very sharp beam cut off can be observed in the interferogram as in Fig. 2(a) -2(h) of the main manuscript if the entrance plane of the BSC is imaged. Otherwise, diffraction will become pronounce and stripes parallel to the cut off can be seen. In addition, the beam cut off will smear out into tilted fringes as well. These diffraction effects are shown in Fig. S5. The intensity of the beam is adjusted beyond the saturation of the camera to better illustrate the diffraction effects. The CCD camera used to take this diffraction measurement is of the model Firefly MV FMVU-03MTM.

Wavefront extraction
The normalized 4-bin phase-shifting phase retrieval is based on four normalized interferograms described by: is the average of the probability density of the two portions of the incident wavefunction on the CDSI, | is the geometric mean of the probability density of the two portions of the incident wavefunction, and φ(x, y) is the phase difference between the interfering beam. Each interferogram is normalized by dividing the total photon number, which is measured by capturing all outputs of an interferometer. Solving for the phase φ(x, y) yields φ(x, y) = tan −1 P r(x, y; 3π/2) − P r(x, y; π/2) P r(x, y; 0) − P r(x, y; π) (S9) which is Eq. (3) in the main manuscript.
After the phase retrieval, it's important to extract the wavefront. However, one should avoid using the technique of curve fitting in this step and manipulate the phase instead. The actual wavefront is maintained this way and further analysis using curving fitting of Zernike polynomials is available. For a typical shearing interferometer, wavefront extraction can easily be done by numerical integrating. The CDSI, however, requires a more sophisticated manipulation of the phase to extract the wavefront. The full procedure is as follows.
Recalls that the phase of the CDSI with the shearing direction ofx is where s is the shearing amount in the +x direction, the subscript ox denotes odd order of x and ex denotes even order of x. The first two terms are due to the geometry of the setup where one of the entrance faces of the beam splitter cube (BSC) has a y-wedge angle. By measuring the phase of an interferogram with s = 0, we will obtain By Eq. (S11) and the property of odd function (W ox (0, y) = 0), we can solve for the odd component of the wavefront. The solution is Now we isolate the derivative of the wavefront by subtracting Eq. (S10) by (S11): Numerically integrating the above equation and dividing by s x yields The even component of the wavefront can be computed by subtracting Eq. (S12) by Eq. (S14): The total wavefront is where f (y) = W oy (0, y) + W ey (0, y). These wavefront components can be calculated using the interferogram with the shearing direction ofŷ following the same derivation above. Their expressions are where s y is the shearing amount in theŷ direction and should be equal to s x . Adding up all of the odd and even components in Eq. (S12, S15, S17, S18) will give us the total wavefront by Eq. (S16). However, because the interferograms of our CDSI only consists of half of the beam, we only obtain the total wavefront in the first quadrant. To obtain the wavefront in all four quadrants, we utilize the property of odd and even function to get where x > 0 and y > 0. As a summary, the wavefront in the four quadrants are extracted using the four phase distributions φ x (x, y; s = 0), φ x (x, y; s = s x ), φ y (x, y; s = 0), φ y (x, y; s = s y ).
If only the odd order wavefront are of interest, not all four phase distributions are needed. In particular, extracting W ox (x, y) (W oy (x, y)) only requires one phase distribution φ x (x, y, s = 0) (φ y (x, y, s = 0)).

Elimination of noise from intensity fluctuation by normalization
In the main manuscript, we discussed how measuring over all the outputs in an interferometer and normalizing the measurements enhance the accuracy of measurement by eliminating the uncertainty due to photon number fluctuation. Here, we will show that the normalization results in a phase sensitivity that scales as N −1/2 independent of photon number fluctuation. Let N ∼ P r(N = n) be the photon number distribution incident on the interferometer for a particular set of measurement conditions (integration time, quantum efficiency, incident beam) and X ∼ B(n, p i,j )/n be the normalized binomial distribution with n trials and p i,j being the probability that the photon lands on a particular pixel (x i , y j ). From now on, the subscript i, j will be dropped. The measured probability p will be a compound binomial distribution: Here we do not consider n = 0 as a measurement. The variance of the compound distribution is given by  Z|N = n)).

Accuracy enhancement by normalization for other types of interferometer
In general, this type of accuracy enhancement by normalization mentioned in the last section can be applied to most interferometric phase retrieval methods.In practice, certain interferometers require some minor modification in order to apply this accuracy enhancement. The setups must be modified to measure all the outputs of the interferometer. This is because the accuracy enhancement relies on the fact that the probability of the photons landing on the detector sums up to one. In particular, for polarization phase-shifting techniques that use polarizer (36), the polarizer must be replaced with a polarizing beam splitter and all of the outputs must be mea-

sured. For an interferometer with a back-propagating beam such as Michelson interferometer
and Signac interferometer (17), a beam splitter (BS) must be placed before the interferometer to measure the back-propagating output. However, this also leads to an attenuation of the backpropagating beam so the other output must be attenuated by the same factor using an identical BS or simply by multiplying this factor in the result. This results in a probability distribution that adds up to a simple fraction that is independent of the phase distribution. For interferometers that utilize diffraction (37), the diffraction efficiency must be known. These spatially uniform losses are fundamentally identical to the effect of non-unitary quantum efficiency η at photodetectors. The pixelated probability distribution at the detector is still distributed identically, however, with values multiplied by η. In addition to the triangular wavefront shape imprinted on the spatial light modulator (SLM) that is retrieved in Fig. 3 of the main manuscript, we have performed the same exact experiment with a sinusoidal wavefront shape. The related results are shown in Fig. S6. This experiment is also performed with a continuous wave (CW) laser beam at a single-photon level by breaking the mode-locking operation of the oscillator. The results are shown in Fig. S7. The experiment performed using a triangular wavefront shape and a CW laser beam at a single-photon level is shown in Fig. S8. The results of the experiments conducted using pulsed and CW operation are nearly identical except the beam size of the CW mode is slightly larger and the beam quality is not as Gaussian as the pulsed operation.

Simulations of complex spatial structures
A spatially uniform beam with an orbital angular momentum of 1 and its shearing interferogram with shear s x = 0, s x = 1 pixel, s y = 0, and s y = 1 pixel at four different phase-shifts of 0, π/2, π, and 3π/2, are simulated. The intensity profile is chosen to be uniform to make the fringes easy to recognize. The wavefront is then extracted using the method described in the section "wavefront extraction." The results are shown in Fig. S9. The retrieved wavefront is near identical to the simulated wavefront. However, the abrupt 2π phase change and the discontinuity at the singularity causes artifact along thex direction. Excluding this artifact, the error is very close to zero.

Matlab data files and codes
All of the data in this paper are processed using Matlab. The Matlab data files and codes may be requested from the authors. The function unwrap_phase that is used in most of the codes is developed by Kasim, et al. (38, 39).

Statistical Analysis
The gated pulse of the intensifier tube that intensifies and collects optical signal is set to be the period of the ultrashort laser pulses. Thus, the attenuated pulses follow a Poisson distribution with a photon rate of about η photons per frame where η is the quantum efficiency of the ICCD camera. The probability of landing on a particular pixel in the detector is described by the probability density in Eq.
(1) of the main manuscript after normalization: where K is the normalization constant such that the probability sums up to 1 over all the pixels.
The error in measured value of p of each pixel is p(1 − p)/N where N is the total photon number measured over both output of the BSC for each measurement. The phase is retrieved from four interferograms using: φ(x, y) = tan −1 P r(x, y; 3π/2) − P r(x, y; π/2) P r(x, y; 0) − P r(x, y; π) (S25) We use the variance formula of error propagation to find the absolute error as follows [29]: ∆φ(x, y) = 3 i=0 dφ(x, y) dP r(x, y, iπ/2) 2 [∆P r(x, y, iπ/2)] 2 (S26) The derivative can be easily calculated using chain rule, and the trigonometric identities d tan −1 (θ)/dθ = 1/(1 + θ 2 ). The absolute error was calculated in the main manuscript which is just the standard deviation of a normalized binomial distribution. The expression becomes where N k is the total number of photons over all outputs used for the measurement of P r(x, y, iπ/2). It is of great interest to determine the noise limit. In fact, the phase noise in Eq. (S27) scales as 1/ √ N k , the standard quantum limit.

Reaching the standard quantum limit
For any type of phase retrieval that involves multiple measurement of intensity (probability) distribution, normalization results in a phase error of [29]: as derived in Eq. (5) of the main manuscript. The absolute error of ∆P r k (x, y) is P r k (x, y)(1 − P r k (x, y))/N k where N k is the total number of photons over all outputs used for the measurement of P r k (x, y). The derivative term is independent of N k because both φ(x, y) and the probability distribution P r k (x, y) are independent of N k . Given that the total number of photons used of each intensity measurement are the same, the phase noise dφ(x, y) scales as 1/ √ N k , the standard quantum noise limit.
There is a special case in which the standard quantum noise limit is achieved. This is when the visibility is one and there is only one pixel on each of the detectors over two complement outputs of an interferometer. Because there is only a single pixel over an output, the phase needs to be spatially uniform. For such a configuration, the photon count over an output can be described by the equation N 1 = N * p 1 = (N + N cos φ)/2, where N 1 , N is the total number of photons registered over detector 1, both detectors, respectively. The phase uncertainty can be calculated as follows: Among the series of N photons, the N 1 photons that landed on an output can be modeled by a binomial distribution. This allows us to calculate the right hand side of Eq. (S29), resulting in Substituting p 1 = (1 + cos φ)/2, and sin φ = √ 1 − cos 2 yields exactly ∆φ = 1/ √ N .

Number-phase uncertainty relation
The number-phase uncertainty relation (27)