Interferometric imaging using shared quantum entanglement

Quantum entanglement-based imaging promises significantly increased resolution by extending the spatial separation of optical collection apertures used in very-long-baseline interferometry for astronomy and geodesy. We report a table-top entanglement-based interferometric imaging technique that utilizes two entangled field modes serving as a phase reference between two apertures. The spatial distribution of a simulated thermal light source is determined by interfering light collected at each aperture with one of the entangled fields and performing joint measurements. This experiment demonstrates the ability of entanglement to implement interferometric imaging.

Coherent measurement of light entering separate collection apertures of an imaging system, an approach called aperture synthesis that forms the basis of interferometric imaging, enables increased angular resolution beyond the single-aperture diffraction limit, as first demonstrated by Michelson and Pease [1].The resolution of interferometric imaging is limited, in-principle, only by the aperture separation (the 'baseline').Unlike interferometric imaging in the radio frequency (RF) band, where Earth-sized telescope arrays have been implemented based on locally recording the RF fields at each telescope [2], in the optical band (visible and nearinfrared) the maximum baseline length is limited, in part, by the difficulty in performing local coherent detection with low noise at the few-photon level.The desire to increase angular resolution in the optical band by increasing the baseline is motivated by a number of applications currently limited by resolution: the search for exoplanets and the study of their atmospheres, resolved imaging of black hole event horizons in the near-infrared to complement the mm-wave imaging of the Event Horizon Telescope (EHT), geodesy, imaging of planet-forming disks, and stellar surfaces beyond the Sun [3][4][5][6][7][8].This has prompted the quantum information science community to search for new tools such as shared optical entanglement between the receivers to extend optical interferometric imaging baselines [9][10][11][12][13][14][15][16][17].
There are three common approaches to achieving high angular resolution with synthetic apertures: 1) transporting the received fields to a central location and interfering them (direct interference) as depicted in Fig. 1(a); 2) distributing and interfering strong mutually coherent fields, known as local oscillators (LO), with two received fields, as shown in Fig. 1(b); or, 3) measuring intensity correlations between the collected fields as in the Hanbury Brown-Twiss (HBT) experiment [18][19][20][21].In the first approach, the challenge to increasing the baseline lies in constructing low-loss optical channels between telescopes, with the current state-of-the-art baseline be- a) Direct interference combines the collected fields at a beam splitter with a known variable phase shift, δ, and performs intensity measurements at the output.b) Indirect interference combines the collected fields with reference fields (laser light or path-entangled single-photon state) followed by local measurements at each aperture.Correlations between measurement outcomes at each telescope yield the mutual coherence function of the fields collected by the two telescope stations.
ing 330 m long [22].The second approach is limited by the photon-number uncertainty of the LO (also called shot noise) as noted by Townes [23].In the optical regime, direct interference is typically preferred due to its higher signal-to-noise ratio (SNR) compared to interference with a classical LO [24,25].In contrast, radiofrequency interferometric imaging (e.g., EHT) uses LOs since the effect of shot noise is significantly reduced: for the same average energy as one photon in the optical field, the RF field contains approximately 10 6 photons [23].The third (HBT) approach is limited to relatively bright sources since it requires detection of at least two photons from the astronomical source, whereas the other methods require only one photon.
In this Letter, we present an experimental realization of an interferometric imaging technique that utilizes quantum entanglement as proposed by Gottesman, Jennewein, and Croke (GJC) [9].To increase the baseline while delivering high SNR at optical wavelengths, GJC proposed the use of a distributed, path-entangled singlephoton state as a shared phase reference in conjunction with a quantum repeater network to alleviate transmission losses [3].A path-entangled reference state (PERS) can be generated by, e.g., splitting a single photon occupying a single traveling mode into two paths of (un)equal length using a beam splitter (BS) [26][27][28].The entangled state of the two modes (1 and 2) in the Fock basis is In the present application, the state provides a phase reference, δ, between two locations while minimizing the shot noise associated with the reference field.
The two paths are interfered with light collected by the apertures, as depicted in Fig. 1(b).A coincidence detection event between photon-counting detectors at each aperture, labeled ±1 and ±2 in Fig. 1(b), realizes a quantum joint measurement allowing coherent detection.For weak thermal-like sources, schemes using joint quantum measurements of the received fields (e.g., direct interference or the GJC protocol) can yield more information per received photon than any scheme using independent measurements at each receiver (e.g., using classical LOs) [15].Here independent measurements are defined generally to include local operations and classical communication.
Interferometric imaging is based on measurement of the complex degree of coherence of an electromagnetic field after propagation from a source [29].In the detection plane, the complex degree of coherence is given by j , where the electric field, E, is evaluated at two points labeled by transverse position vectors, ⃗ y 1 and ⃗ y 2 , and angular brackets imply ensemble averaging.The complex visibility, v(∆⃗ y)e iϕ(∆⃗ y) , arising from the interference of the fields collected at ⃗ y 1 and ⃗ y 2 is equal to j(⃗ y 1 , ⃗ y 2 ).The van Cittert-Zernike theorem for a spatially incoherent source parallel to the detection plane shows that v(∆⃗ y)e iϕ(∆⃗ y) is proportional to the Fourier transform of the normalized source intensity distribution- , where ⃗ x is the transverse position vector in the source plane-and depends on the displacement ∆⃗ y = ⃗ y 2 − ⃗ y 1 , between ⃗ y 1 and ⃗ y 2 , and the difference of their distance from the optical axis, Here z is the distance from the source plane to the observation plane and λ is the wavelength [30][31][32][33].In interferometric imaging, measuring the complex visibility enables reconstruction of the source intensity distribution.Our experimental setup demonstrating the GJC scheme is shown in Fig. 2. It consists of a pulsed heralded single-photon source and a pseudo-thermal light source in a single spectral-temporal mode mimicking filtered light from a star, which enables mode matching of the two sources.Both sources are derived from a titanium-sapphire (Ti:Sa) laser with 80 MHz repetition rate, 830 nm central wavelength and approximately 10 nm full-width half-max (FWHM) bandwidth.To simulate starlight, i.e., a spatially incoherent light source, we direct a 5 mW laser beam onto a double-slit mask placed 5 mm in front of a rotating, nearly-Lambertian diffuser.The scattered light propagates 1.02 m to the detection plane resulting in a transverse coherence area on the order of a few square millimeters.A balanced, free-space BS separates the field into two spatial paths, enabling simultaneous coupling into two bare, polarization-maintaining, single-mode fibers (PM-SMF; Thorlabs PM780-HP) with a core diameter of 5 µm, acting as small-aperture telescopes.One of the fibers, T1, is scanned in transverse position (y 1 ), to sample the field over a range of 6 mm in 10 µm steps.The other fiber, T2, is statically positioned in the center of the scattered field.Scattered light from the rotating diffuser creates a timevarying speckle pattern in the far-field so that a single spatial mode exhibits thermal (Bose-Einstein) photonnumber statistics (defined as p(n) = nn /(n + 1) n+1 ) with average photon number n ≈ 0.008 and second-order coherence g (2) (0) = 2.00 ± 0.05 [34].Because the intensity distribution in the far field is spread over a region greater than 10 cm, each fiber collects the same mean photon number, n, as long as the fibers are not displaced more than a few millimeters from one another.The joint state of the collected fields is described by a density matrix, ρ kl ij ; with i and k labeling creation and annihilation operators for the field at T1, and j and l playing similar roles at T2.Further details throughout the paper are given in the Supplemental Materials [35].
The main beam from the Ti:Sa (1.24 W) is frequencydoubled to a wavelength of 415 nm in a 700 µm-thick beta barium borate (BBO) crystal to produce a 100 mW beam that pumps an 8 mm-thick potassium dihydrogen phosphate (KDP) crystal phase-matched for degenerate, collinear, type-II spontaneous parametric downconversion (SPDC), which is designed to produce photons in the same spectral-temporal mode as the pseudothermal source.
The source produces a two-mode squeezed-vacuum state that is nearly spectrally separable with a central wavelength of 830 nm where the herald mode has ∼5 nm of bandwidth FWHM (horizontally-H-polarized), while the heralded mode has a bandwidth of ∼12 nm FWHM (vertically-V-polarized) [36].FIG. 2. A titanium-sapphire (Ti:Sa) laser is separated into two paths by a beam splitter (BS).One path creates a single spectral-temporal-mode thermal-like state by passing through a double slit (D) and then a time-varying scatterer (S).The scattered light is collected by fibers T1 and T2 after a BS.The other path undergoes second-harmonic generation in a BBO crystal followed by parametric down conversion in a KDP crystal.The generated photon pair is split at a polarizing BS (PBS); one is used for heralding.The heralded photon is sent to a BS to generate the path-entangled reference state along two paths, with one path having a time-dependent length controlled by a PZT to introduce a phase shift, δ(t), driven by a signal generator (SG).Each path is coupled into single-mode fibers (S1 and S2), and combined with the fields collected by T1 and T2 in two fiber BS (FBS).The outputs of the FBSs are monitored by superconducting nanowire single-photon counting detectors (SNSPDs) and time taggers.
A polarizing BS separates the H and V fields, where the V field is sent to a single-mode fiber connected to a superconducting nanowire detector (SNSPD).We measure an unheralded second-order coherence g (2) (0) = 1.66 ± 0.06 (which implies a Schmidt number of 1.51 and near separability [37]) and an average photon number of approximately 0.01, after correcting for detection efficiency.A detection event at the heralding detector indicates the presence of a single photon (or more) in the H-polarized beam.The density operator of the heralded field in the photon-number basis is approximately ρ heralded ≈ η p(1)|1⟩⟨1|+(2−η)η p(2)|2⟩⟨2| where η ≈ 0.28 is the measured heralding efficiency and p(n) are thermal photonnumber statistics, indicating the possibility of more than one photon in the pulse.
The heralded H-polarized field is directed to a 50:50 BS realizing the PERS in Eq. 1.A mirror bonded to a piezo-electric translator (PZT) in one path introduces a controllable, relative phase difference between the two paths.The PZT is controlled by a signal generator (SG) that applies a 600 Hz triangular waveform, giving a timedependent phase, δ(t), which varies between 0 and approximately 4π.The heralded PERS distributed to the telescopes is approximately given by Eq. 1, which neglects the vacuum and higher-order photon numbers (a proof that the full state is entangled is given in the Supplemental Materials [35]).The two paths are coupled into PM-SMFs, indicated by S1 and S2 in Fig. 2.
The complex visibility of the simulated starlight, as defined in Eq. ( 2), is measured by interfering the two modes of the PERS, S1 and S2, with the fields collected by the two fibers acting as telescopes, T1 and T2 in Fig. 2, at fiber BSs (FBS).The FBS outputs are sent to SNSPDs.All detection events are time-tagged using a time-todigital converter.We keep coincidence events only when they are registered between different telescopes.The coincidence events, occurring probabilistically at times t j , can be associated with a known PZT phase, δ(t j ), by comparison to a phase-locked square voltage pulse used as a clock indicating the change in direction of the PZT.The set of event phases are used to estimate the complex visibility (see Supplemental Materials, [35] and [38] for details).
Theory predicts the estimated visibility measured by our experiment to be (see Supplemental Materials [35]) where ± is the phase difference acquired from the FBSs determined by which detectors register events, ρ is the density matrix of the collected star photons that depends on Eq. 2, p(n) is the Bose distribution of the two-mode squeezed state from the SPDC, η is the heralding efficiency, Γ ≈ 0.5 denotes the coherent overlap of the idler field mode from the SPDC and the field collected from the pseudo-thermal source.Given the experimental pa- rameters, the predicted maximum value of |ṽ| from Eq. 3 is 0.24.
In this experiment, the magnitude of the complex visibility versus baseline is measured.The relative phase of the interferometer paths fluctuates on a 10 s time scale, which is typically longer than necessary to acquire sufficient coincidence events to determine the amplitude and phase of the complex visibility at a given baseline, but does not allow a stationary phase reference between the measurements at different separations.However, the Fourier transform of the phase-independent, modulussquared complex visibility as a function of baseline, |v(∆⃗ y)| 2 , yields the autocorrelation of the normalized source intensity distribution, C(δx) = I(x)I(x + δx)dx (see Supplemental Materials [35]).Such a reconstruction without interferometer stability is sufficient to verify a source distribution using a model fit.
For two known source distributions created from 0.5 mm wide vertically-oriented double slits separated by 1 mm or 2 mm, illuminated by a Gaussian-profile beam, Fig. 3(a) and (b) shows the magnitude of the complex visibility as a function of baseline in the horizontal direction, ∆⃗ y.Each data point in the measured visibility results from processing around 10,000 collected events-the PZT phases at the time of a coincidence, δ(t j )-over 10 s (see Supplemental Details [35] and [38]).and w) for the double slits with 0.5 mm slit width, and slit separations of 1 mm (d = 0.5 mm) and 2 mm (d = 1 mm), which are illuminated by the same Gaussian beam.
The Fourier transform of the squared visibility is proportional to the autocorrelation of the source intensity distribution, shown in Fig. 3(c) and (d).The estimated errors, which are equal in size to the circular points, account for statistical errors and phase fluctuations, but do not account for systematic errors.We also show a best fit of the reconstructed autocorrelation data to a model based on the autocorrelation of an intensity distribution describing the experiment (see Supplemental Materials [35]).Parameters in the model are the distance from the center of a slit to the midpoint between the slits, d; the width of the slits, w; and the radius of the common Gaussian beam illuminating the slits, σ.
The best-fit parameters and errors determined using maximum likelihood estimation for the sources are shown in Table I and found to lie within fabrication tolerances of the slits (±20 µm).We also find that, as expected, the fitted Gaussian beam radius parameter in the two cases are equal within error.Overall, the fit to the expected intensity distributions for the two slits indicates that the GJC scheme is capable of source reconstruction.
In conclusion, we have shown that an approximately single-photon state distributed across two paths can act as an entangled-state phase and amplitude reference for use in distributed VLBI, as proposed by GJC [9].The complex degree of coherence of a distant incoherent source was determined and served to reconstruct its intensity distribution.To achieve extended baselines, future work will study true thermal sources, reduce higher-order photon number contributions, and extend the scheme to detection of multiple spectral-temporal modes with the purpose of increasing detection rates.Lengthening baselines will require a quantum repeater chain to distribute the path-entangled light with lower loss than using optical fiber alone [39,40].Correcting for turbulence-induced phase distortion by the Earth's atmosphere is possible in arrays with three or more telescopes using closure phase [41], and progress toward such an implementation with PERSs can be found in [15].Further, it was proposed that a quantum network in which single-photon states could be stored in quantum memories and used on demand could provide a further increase in multiplexing ability [11,12,16,17].While functioning quantum networks are years in the future, if successful the proposed scheme could become one of the first practical uses of a quantum network and would open new

Supplementary Material
Interferometric imaging using shared quantum entanglement Matthew R. Brown, Markus Allgaier, Valérian Thiel, John D. Monnier, Michael G. Raymer, and Brian J. Smith, Here we provide several background and technical details that explain and support the methods used in the main body of the paper.

A: Model of Coincidence Count Probability
To predict and model the coincidence count probability in the VLBI scheme, we need two important matrices: one describing the 'faux' star electromagnetic field collected by the single-mode-fiber 'telescopes,' and the other describing the signal field from the down-converted pair, which is used as a path-entangled reference field (PRF).To model the density matrix describing the field collected by the telescopes, we first realize that the act of rotating the diffuser creates (by the central-limit theorem) Gaussian statistics in the complex field variables, as is well known in statistical optics [33].Further, we know from the van Cittert-Zernike Theorem (Eq. 1 in the main text) that the fields collected by the telescopes are correlated [33].We can then use the optical equivalence theorem (relating semi-classical field states to quantum states [30]) to write the density operator for the mutual state of light arriving at the telescopes as, ).
The density matrix elements are: , where, in the experiment,  ̅ 1 ≈ 0.008 and  ̅ 2 ≈ 0.008 are the average photon numbers collected by the telescopes T1 and T2 respectively, and (Δ)  ϕ(Δy) is the complex interference visibility determined by the van Cittert-Zernike theorem (see Eq. ( 1) in the main text).
The intermediate state of the signal field post-selected on detecting a photon in the herald mode is approximated by where  is the heralding efficiency (the probability that a photon in the herald mode heralds a photon in the signal arm) and () is the Bose distribution.The signal is then propagated to a beam splitter and each output mode is then propagated one to each telescope to be interfered with the collected "faux" star field.As depicted in Fig. 1(b) in the main text, one mode is chosen to be have a phase change given by .Using this state we can calculate the resulting probability to observe a coincidence event between pairs of detectors at different telescopes to be ), which is similar to the classical expression for interference, where the first term is the total probability of a coincidence event, and the second indicates the coherent variation of probability of a coincidence event as either the visibility is changed by moving a telescope or the phase, , is changed by varying one path length in the interferometer.The ± is determined by which output port of the beam splitter yielded the measured event.

B: Path Entanglement of the Heralded Field
To demonstrate the entanglement of the heralded field, we must first notice that the general state (density operator) of the heralded field including loss and non-unit efficiency of the heralding detector is given by thermal difference states [42], represented by the form: where , , , and  are known functions of the average photon number of the initial down conversion, the quantum efficiency of the heralding detector, and the loss the heralded field experiences prior to the beam splitter (which we assume to be upper bounded by 10%) and are defined in [42].To verify entanglement between the fields in the two paths after the heralded field is split by a beam splitter, we show that the projection to the qubit subspace of the two Hilbert spaces defining the bipartite state is entangled, which is sufficient to show the entire state is entangled [43].To show the entanglement of the qubit subspace, it is sufficient to show that the partial transpose of the density matrix is not positive by, for instance, determining that the determinant is negative (implying that at least one of the eigenvalues is negative) [44,45].The partial transpose of the qubit density matrix, derived from the state above, is ) .
Hence, the determinant is ).
For our measured experimental parameters and an estimated upper bound on the loss (10%) from the down conversion crystal to the beam splitter ( ≈ 351.38,  ≈ 1,  ≈ 0.72,  ≈ 0.009), we find the determinant to be ≈-0.038thus showing that the is entangled.

C: Aligning Telescope Fibers and Measuring Scattered Light
One important aspect of the experiment is alignment of the telescope fibers so that they collect, approximately, the same transverse part of the field.To first order, we first try to find the peak of the intensity distribution (usually Gaussian) without the diffuser since each fiber is much smaller than the spatial distribution in the far field.However, this method becomes less precise as the intensity distribution at the diffuser becomes larger.
The only way to truly center the fibers, then, is to raster scan and search for interference fringes in the coincidences as described by the van Cittert-Zernike theorem.
Further, we also need to ensure that the light scattered by our diffuser has properties consistent with a pseudothermal source and has no ballistic throughput.With our fibers aligned, we measure a  (2) (0) = 2.00 ± 0.05, consistent with thermal-field statistics.However, a  (2) (0) value alone does not ensure that we truly have a thermal state; therefore, we image, on a camera (Thorlabs DC1545M), the spatial structure 48.For near-singlemode thermal-like light, we expect this form to be a negative exponential for large intensities with a peak at low intensities due to the nonzero size of the individual pixels [46,47].Finally, we want the source coherence area to be as small as possible.To test for this, we first realize that there is a generalized van Cittert-Zernike theorem that deals with a nonzero coherence area at the source [33].If the coherence area were too large, we would notice the intensity distribution in the far field to have a small angular spread, that is, to not be Lambertian.We measured that the intensity angular spread indicates a source coherence roughly less than 6 microns, small enough to be considered point-like when examining the far-field correlations only near the beam axis.

D: Timing and Synchronization
In order for the non-local oscillator field and the collected thermal-like field to interfere at each of the fiber beam splitters, the temporal envelopes of the pulses must overlap.To accomplish maximum overlap, Hong-Ou-Mandel (HOM) interference is observed between the unheralded single photon and the effective thermal field collected by the telescopes.There are two delay stages to control the relative timing between the 100 fs, thermal field pulses and the unheralded SPDC 100 fs pulses.The stages are then set to the position of maximum two-photon interference [48].Once the delays are set to the position maximizing the HOM interference, the experiment is ready to record data.Below are some technical details covering the synchronization processing of time tags collected during measurement.
The data comes in the form of time-tags with a timing jitter of 13 ps from two ID900 time-to-digital converters (TDC).The "start measuring" command from the controlling computer had a different time of arrival between the two TDCs, hence, they must be synchronized for coincidences to be detected.Specifically, this task is to find the timing offset between the two TDCs and add this time to the time tags from the TDC that received the start message first.To measure the absolute timing delay between the TDCs we use a 600 Hz square-pulse voltage source, phase-locked to the 600 Hz triangle-wave signal that drives the PZT, which makes a timevarying phase between the different paths of the single photon.To make the square pulse common between the two TDCs, it is routed through a delay generator (SRS DG535) where the signal is split into two channels, one sent to each TDC.Any differing delays accrued by the cable lengths were accounted for by measuring the timing difference on a single TDC.In our case, the timing delay was much less than the 6 ns coincidencewindow (<100 ps) and is negligible.Given that each TDC has a common signal, we note that during the internal phase-locking procedure of the square-pulse generator there was a slight timing change between subsequent square pulses as the square pulse is phase-locked to the triangle wave.Since each TDC will see this sudden shift in timing, we used this shift to calculate the absolute timing between the TDCs by the following procedure: From one TDC, using the stream of time tags resulting from the square wave, calculate the inverse of the difference of subsequent time tags (hereafter called the frequency).The list of frequency values can be searched for any frequencies that are different (by usually more than a few Hz; a threshold can be applied by measuring the standard deviation of these arrival frequencies of the square wave and is in the sub-Hz regime).
The last position where the frequency deviates from 600 Hz was selected.Since the frequencies are made of two consecutive time tags ordered by their time-of-arrival (first and last), we took the time tag directly after the last time tag of the pair to ensure that there are no further drifts in frequency due to the phase-locking procedure.The resulting time tag (called the synchronization time) at each TDC are compared and the difference of these time tags yields the absolute timing difference between the two TDCs to the nanosecond level or better.To confirm synchronization, coincidences between the herald channel on one of the TDCs and one of the detection channels (which contains a single photon from the other half of the PDC pair) connected to the other TDC are measured.Finally, all time tags prior to the synchronization time were dropped since the square pulse was not yet phase locked to the triangle wave driving the PZT.

E: Determining Complex Visibility from Time Tags
The time tags produced by measurement are processed in software.In our case, the information is contained in the probability of a coincidence as found in section A. We must then take each of the measured time-tags and compute the coincidences, with a 6 ns window, between the detectors of interest including the herald detector.There are four sets of time-tags corresponding to the four possible detector combinations that give information (the other two possible groupings, coincidences at the same telescope, cause the loss of a factor of two in data rate, as explained in [10]).Furthermore, two of the four possible combinations exhibit a π-phase shift compared to the other two.To make sure that the phasors from all four arrangements are in phase, an extra πphase is added to the phases measured from one of the arbitrarily chosen, but constant, pairs of detector combinations.The heralded data is then created by finding coincidences (again with a 6 ns window) with the herald signal from PDC and each of the four coincidence streams found previously.It is important to note that all data currently is still in time-stamp format.
Each coincidence event is associated with a specific PZT phase found by reference to the known triangle waveform.This phase (δ  ), also called a phase tag, can be determined by the known function applied by the signal generator that drives the PZT (600 Hz triangle pulse with ~4π amplitude) given by the time difference of the coincidence time-tag to the time-tag of the nearest square pulse, which is phase-locked so that its rising edge starts at the beginning of a triangle wave (see section B).Effectively, these square pulses also act as a clock.The amplitude of the triangle pulse is measured by aligning a classical laser beam into the same paths (and fibers) as the distributed single photon.By then interfering the outputs of the single-photon collection fibers using a fiber beam splitter and measuring the interference pattern of the output as the PZT is swept over time (by e.g. a photodiode and analog-to-digital converter), one can estimate the approximate phase applied by the PZT through curve fitting.The peak-to-peak voltage applied to the PZT is adjusted so that a relative phase the two paths of roughly 4π is applied to the optical fields.
We estimate the measured complex visibility in a single trial recording many events by averaging over the resulting phasors to yield the mean phasor, The brackets indicate averaging over N samples, where N is of the order ten thousand.This technique is called phase-tagged photon counting [38].This estimator for the complex visibility can be shown to be unbiased.

F: Phase Fluctuations, Autocorrelations, and Error Propagation
Environmental phase fluctuations in our system occur on an estimated times scale of 10 s.Given that the phase fluctuations may be of order of 2 , we cannot currently achieve full imaging as defined in the van Cittert-Zernike theorem.Rather, we use the fact that the autocorrelation, ( ⃑), of the source intensity is related to the squared modulus of the measured visibility [33] through Thus, we are interested in estimating quantities related to the magnitude of the visibility.
In theory, each trial should yield a complex value for the mean phasor that is normally distributed with equal variance in the real and imaginary axes, with no correlation between them, and centered around the complex visibility by the central limit theorem.However, due to the uncontrolled environmental phase fluctuations between trials, the complex distribution may be more complicated.Thus, a method of determining the variance in the complex direction of the mean value (our estimate for the visibility) is needed.This is achieved by rotating the measured distribution (of approximately 10,000 complex samples) such that the new mean phasor points along the positive real axis.One can then compute the variance in the real direction, which corresponds to the variance of the original distribution in the complex direction of the visibility [49].With this definition, one can use standard error propagation from |(Δ ⃑)| 2 to ( ⃑); however, propagation from (Δ ⃑) (in the direction of (Δ ⃑)) to |(Δ ⃑)| 2 requires the use of a second-order delta method for cases where the mean value is small [50].Given that our error-propagation function, () =  2 , is a polynomial of degree two, this error propagation is exact.In particular, this means that Here,  is the number of samples, (,  2 ) is the normal distribution centered at  with variance  2 ,  1 2 is the chi-squared distribution with one degree of freedom, and  → symbolizes convergence in distribution.So, the new mean and variance are given by  2 +  2 and 4  2  2 + 2  2 respectively, where  is the mean of (Δ ⃑) and  2 is the asymptotic variance (i.e.variance due to the central limit theorem) of (Δ ⃑) in the direction of the mean.Both  and  2 can be determined directly from the sampled distributions.This procedure yields error estimates that are smaller than the circles representing data points in the reconstructed autocorrelation functions in the main text, Fig. 3.

G: Model Fitting
In this paper, we fit the inferred source intensity autocorrelation function () to a theoretical model that represents well the known source distribution.Here we detail the model and the fitting parameters we assume, along with the fitting parameters for the 1 mm and 2 mm slits.We define the model for the normalized intensity distribution as (; , , ) = where  is the radius of the Gaussian beam that impinges on the double  is the width of an individual slit, 2 is the distance between the centers of the two slits, and the UnitBox() function is defined as 1 when || ≤ 1 2 and 0 otherwise.We can then calculate the autocorrelation of (; , , ), which is omitted here due to its complexity.However, there are a few complications that must be accounted for: the visibility is not unity in our experiment, the discrete nature of the data, and the non-zero nature of the resulting data.The first two of these issues require us to use a multiplicative scaling parameter, , while the second forces an offset parameter, .Given that we need a multiplicative scaling factor, we do not need to worry about normalization.Hence, we will use the unnormalized intensity  ̃(; , , ) =  ).Thus, () = ∬  ̃(; , , ) ̃( + ; , , ) 2  and the final model is given as  () + , where we fit to all parameters listed above using Mathematica's NonlinearModelFit function.The full parameter list of the fits found for the data is shown below.The beam radii inferred are comparable; and, the larger error is attributed to the slight differences in the curvature of the side peaks and the necessity to fit a scaling parameter.
To fit the measured visibility, we use these parameters in the Fourier transform of  ̃(; , , ) and allow for scaling and offset parameters as usual.Below, we show the functional form of the modeled absolute value of the visibility since the phase is unimportant in this experiment: These expressions allow us to propagate the model fit of the autocorrelation back to our measured visibility.

FIG. 1 .
FIG.1.Schematics of two approaches to interferometric imaging.Light from a distant incoherent source with intensity distribution I(x) is collected by apertures with baseline ∆y.a) Direct interference combines the collected fields at a beam splitter with a known variable phase shift, δ, and performs intensity measurements at the output.b) Indirect interference combines the collected fields with reference fields (laser light or path-entangled single-photon state) followed by local measurements at each aperture.Correlations between measurement outcomes at each telescope yield the mutual coherence function of the fields collected by the two telescope stations.

FIG. 3 .
FIG. 3. (a), (b) show the measured magnitude of the complex visibility, |ṽ(∆y)|, versus telescope separation, ∆y, in mm and, u, normalized by wavelength for double slits with 1 and 2 mm separation, respectively.The points are linearly interpolated to guide the eye.The black line is a model fit that takes parameters found from the fits shown in (c), (d) and adjusts scaling and offset.(c), (d) show the reconstructed autocorrelation of the source distribution for the 1 and 2 mm slit separations versus position in mm and arcmin.Experimental error bars are subsumed by the circular points used (see Supplemental Materials [35] for details).The dashed black lines are the best model fits.Insets show the unnormalized source intensity distributions estimated from the model fit parameters.
,,=0, where  ̂ and  ̂ † are the annihilation and creation operators for mode  and the density matrix,    , 26 cm away from the diffuser and near the center of the transverse structure shown in Figure B.1a).From this image, we gain two major things: first, we have good spatial scattering without a ballistic component; second, we can use the image to make a histogram of the instantaneous intensity values, as shown in Figure B.1b).

Figure B. 1 .
Figure B.1.Scattering of the light through the diffuser.In a) we show a far-field camera image of the scattered light.Note the lack of a ballistic component (bright center spot) in the scattering.In b), we show the histogram of the number of pixels sharing a specific intensity captured by the camera.This shows the nearly Gaussian and Lambertian character of the scattering medium.

TABLE I .
Autocorrelation fit parameters (σ, d, where  1 ,  2 ,  are the detector efficiencies at telescope 1, telescope 2, and for the herald, and Γ is the total mode overlap between the PRF and the thermal-like field collected by the telescopes.This form can be rewritten as