Ab Initio Spatial Phase Retrieval via Intensity Triple Correlations

Second-order intensity correlations from incoherent emitters can reveal the Fourier transform modulus of their spatial distribution, but retrieving the phase to enable completely general Fourier inversion to real space remains challenging. Phase retrieval via the third-order intensity correlations has relied on special emitter configurations which simplified an unaddressed sign problem in the computation. Without a complete treatment of this sign problem, the general case of retrieving the Fourier phase from a truly arbitrary configuration of emitters is not possible. In this paper, a general method for ab initio phase retrieval via the intensity triple correlations is described. Simulations demonstrate accurate phase retrieval for clusters of incoherent emitters which could be applied to imaging stars or fluorescent atoms and molecules. With this work, it is now finally tractable to perform Fourier inversion directly and reconstruct images of arbitrary arrays of independent emitters via far-field intensity correlations alone.


Introduction
Coherent diffractive imaging uses the stationary far-field interference of elastically-scattered light to infer the geometry of a scattering potential via Fourier analysis. Since most photodetectors perform an intensity measurement, information about the relative phases ( ì ) of the scattered waves at pixels ì is lost and Fourier inversion to real space is incomplete [1]. This "phase problem" is shared across a variety of imaging modalities, including x-ray crystallography and optical microscopy, and research in each field has arrived at a variety of techniques to obtain the phase information.
Non-stationary or incoherent scattering processes are known to provide more information, as much as twice the information cut-off in an optical microscope utilising incoherent illumination or fluorescence as compared with plane-wave illumination [2]. The far-field intensity distribution of such a process is featureless, but the measurement of intensity-intensity correlations can nevertheless be used to extract the Fourier amplitude of the object's structure as first demonstrated by Hanbury Brown and Twiss on the radio emission of bright stars [3]. This approach is attractive in situations where lenses of high enough angular or spatial resolution do not exist. This is certainly the case in the X-ray regime where recent work has examined the possibility of using photon pair correlations to retrieve the Fourier spectrum of x-ray fluorescence emission [4][5][6][7][8][9]. One is still left with the phase problem, which can be solved using iterative phase retrieval [9] when the correlations are adequately and extensively sampled by detectors with large numbers of pixels. However, it has been known since the 1960s that intensity triple correlations can reveal partial phase information directly in the form of the so-called closure phase [10,11]. This has been heavily investigated in the field of radio astronomy [12][13][14][15][16] with the aim to develop the means to reconstruct an image of an arbitrary arrangement of emitters without the use of additional constraints.
Retrieval of the Fourier phase, ( ì ), begins by first computing the absolute value of the closure phase, |Φ( ì , ì)| = | ( ì + ì) − ( ì ) − ( ì)|, from the triple correlations. Unfortunately, we need the signed value of Φ to recover completely ( ì ) for the following reason: At the start of phase extraction, an estimate value is chosen for the first pixel. But for the next pixel, the sign ambiguity of |Φ| returns two possible values of ( ì + 1) and an additional two possible values for every subsequent pixel (except for special arrays where sgn(Φ) may be assumed constant). To avoid this exponential expansion of the solution space, we show how redundant information contained in |Φ| may be used to constrain the possible values of sgn(Φ). Multiple publications have described the concept but, to the best of our knowledge, no one has yet provided complete or useful details on calculating sgn(Φ( ì , ì)) [17][18][19][20][21][22][23][24][25][26][27][28][29][30]. Ab-initio phase retrieval from the third-order intensity correlations has thus remained incomplete for decades. With our method, it is now possible to solve for the phase of an arbitrary array of incoherent emitters from the third-order intensity correlations alone. Combined with the second-order intensity correlations, we have a completely general method for reconstructing images of arrays of incoherent emitters.
In this paper, we describe our solution to the sign problem of the closure phase, with the help of a simple 1D example using round numbers in section 3.1.1, and show a numerical implementation of our method with simulated data from classical independent light sources. This same method may be used to reconstruct images of star clusters or, with some corrections to account for the use of a quantum light source, arrays of fluorescent molecules or atoms.

Theory
The diagram in Fig. 1 depicts and contrasts structure determination via coherent scattering to that obtained from incoherent emission. When illuminated with a plane wave with a wave-vector ì , the elastically scattered field has stationary intensity given by for a number of point scatterers with scattering factors and positions ì relative to an arbitrary real-space origin. The photon momentum transfer ì is equal to the difference ì − ì . The phase of each scattered wave, ì · ì , is derived from the difference in the optical path along the directions of the incoming and outgoing waves as compared to a scatterer at the origin. The intensity pattern is thus proportional to the square modulus of the Fourier transform ( ì) of the distribution of scatterers in terms of spatial frequencies equated with ì. The origin of the pattern, ì = 0, is located in the direction of the incident beam. In this forward direction all scattered waves are in phase and there is strong constructive interference, with intensity generally falling with scattering angle. The pattern consists of speckles whose width is inversely proportional to the extent of the object. The recovery of the object's scattering potential is obtained by an inverse Fourier transform of ( ì), but only after the corresponding phases are obtained. If, instead, the object consists of a collection of incoherent point emitters, then there is no dependence on any incident beam and the phase of the emission, relative to that of an emitter at an arbitrary real-space origin, is ì · ì + . We assume that the emission phases are random and uncorrelated on timescales greater than the relevant system coherence time, , due to independent, spontaneous emission at random times. The total light field in this scenario is Fig. 1. Schematic sketch of coherent diffraction in the forward detection plane, intersecting with the incident beam ì . Fluorescence speckle is emitted by the atoms isotropically and the position of the second detector is not dependent on the incident beam. Coherent diffraction data is collected as a function of the scattering vector ì = ì − ì . Correlations between triples of fluorescence photons (intensities) at pixels separated by ì 1 , ì 2 , and ì 3 reveal the spatial phase information lost in the coherent diffraction experiment. often referred to as pseudo-thermal or chaotic and has intensity with the amplitude of electric field emission of the th emitter. The intensity pattern depends on the orientation of the object, and, given the complete independence of emission, at an instant of time this pattern has a uniform intensity modulated by speckles of the same size as the case for coherent scattering. When rapid exposures are measured with a photodetector, we can consider the phases are reset shot-to-shot, changing the instantaneous speckle pattern. From the right-hand side of Eqn. 2, we observe that the structure (sum of ì for all ) would be difficult to discern by averaging intensities over many shots-the random phase resets would drive the interference speckle visibility to zero. However, it remains possible to obtain structural information via intensity correlations [4].
In the following, we use the word atom to refer to any member of a collection of point fluorescent (atoms and molecules) or thermal (stars) light sources. We assume these atoms to undergo spontaneous emission independently, i.e. that each atom emits a field with a phase or time delay that is uncorrelated to the fields emitted by the other atoms.

Intensity Correlations
We consider photon emission vectors in reciprocal space, ì 1 , ì 2 , and ì 3 , and their vector differences as depicted in Fig. 1. The ensemble average of third-order intensity correlations of the light field, Eqn. 2, over all shots can be expressed as and is called the bispectrum. Similarly, the mean second-order intensity correlation function may be written as (2) This equation is often referred to as the Siegert Relation in quantum optics [31]. For a full derivation of Equations 7 and 9 please review the Supplement. In Eq. 7, we have an expression for (3) in terms of constants, the square modulus of (1) , and the real part of a product of complex-valued (1) . Since (1) may be acquired from (2) in Eq. 9, it is possible to extract the last term 7b alone. This term is referred to as the closure in the astronomy literature. We can rewrite the closure as where we have expressed (1) in polar coordinates in the complex plane. As the radial component ( (1) ) is easily obtained from (2) , the phase information, ( ì), can be isolated as follows.
Suppose we set ì 1 = ì and ì 2 = ì where ì , ì map to discrete pixels on a detector. The symmetry of (1) ( ì) and anti-symmetry of ( ì) allow us to rearrange the closure into The inverse cosine of this expression is known as the closure phase, which we represent via the symbol Just as in the Siegert Relation, the third-order correlation function encodes the phase ( ì) at pixels in ì -space beyond the physical spatial extent of the detector (| ì max | = 2 ì max ) as depicted in Fig. 1. Together, the double and triple correlations allow retrieval of the equivalent of a coherent diffraction pattern and its phase across an area of ì -space four times larger than the area of detector coverage [4,6].

Phase Retrieval
Equation 12 for the closure phase Φ( ì , ì) is a difference equation which can be used like a discrete derivative to estimate the slope of ( ì ) between pixels separated by ì. The antisymmetry of the phase pins ( ì = ì 0) = 0. Since overall translation in real-space results in phase ramps in reciprocal space, we can estimate the value of the phase at a nearest-neighbor pixel of ( ì = ì 0) = 0 without loss of generality. This estimate can be refined later by seeking to minimize the total error (see Section 3.2) of all pixels and treating the initial value as a parameter. The difference equation and calculated Φ( ì , ì) from experimental data reveals the value of the phase at the next-nearest-neighbor pixels and so forth until the phase on the entire pixel array has been calculated.
Once the second pixel of (next-nearest neighbor) in any direction is calculated, the interval of the difference equation (ì) may be increased to find the slope between every other (instead of every) pixel in the same direction. Essentially, the phase values calculated for pixels near the origin constrain the possible phase values of pixels far from the origin.

Determining the sign of the closure phase sgn(Φ)
Due to the sign ambiguity of the inverse cosine, every datum from Φ( ì , ì) points to two possible values of the phase ( ì + ì) ( ì + ì) = +Φ( ì , ì) + ( ì ) + ( ì) ≡ + (13) or for any ì and ì. Assuming a global sign often leads to an incorrect slope for , as shown in Fig. 2. The fact that multiple values of Φ( ì , ì) relate to the value of the phase at a single pixel allows us to determine the proper sign of Φ( ì , ì) for each ì and ì. Suppose for a given pixel at ì there exist sets ( ì , ì) in Φ( ì , ì) for which ì + ì = ì. Each set offers a pair of possible values for ( ì), giving 2 possible values for ( ì) altogether. We know that each and every one of the pairs contains the correct value, so comparing the pairs should reveal it. Ideally, the correct value is included times between the pairs and is found simply by taking the intersection of all pairs. Next, we show a simple 1D example to illustrate the principle.
The correct contour is recovered with the wrong slope in the sections highlighted in red. The blue highlight section has both the correct shape and slope, but is offset due to incorporation of incorrect early in the solution.
For purposes of this demonstration, let us assume we have correctly estimated (1) = 1. Then, we can try to determine (2) via which gives the following options Since this is a 1D example, there is no other point ( , ) such that we can deduce the correct value using Φ( , ). The same is true for the calculation of (3) since we only have the information Φ(1, 2) = 1 available to us. In a 2D phase retrieval we would have more values of Φ available to us to solve these pixels close to the origin, but for this 1D example, we will follow the branch we know to be correct with (2) = −3 and (3) = −1. Now the core principle of our algorithm will be demonstrated. We want to find (4) using (1), (2), (3), Φ(2, 2), and Φ(1, 3). This is the first pixel for which we have two constraints on (4) via Φ(2, 2) and Φ(1, 3) which gives us four possible values of (3) in two pairs labeled by the points on the Φ matrix that are associated with (Φ(2, 2) → {4, −16} and Φ(1, 3) → {4, −4}). If we assume that our calculation of up to (3) is accurate, then we should expect that the correct value of ( = 4) is contained in both pairs. Examining the intersection of (4, −4) and (4, −16), we may infer that Similarly, we may determine the correct value and sign of (5) via two constraints from |Φ|.
Again, the intersection of the pairs of possible solutions gives the correct answer (5) = 2.
The last entry requires finding the intersection among three pairs since there are three non-redundant pieces of data in |Φ| we may use.
The intersection of the three pairs leads us to conclude that (6) = 7.

Numerical Algorithm
In practice, the intersection of solution pairs is never exact and a numerical estimation subject to input noise is required. We devised an algorithm to accurately find the intersection of all pairs of possible solutions. For each ì = ì + ì there are, in general, multiple sets of ( ì , ì) positions with previously calculated phases. Each of these sets generates a pair of solutions for the two signs of Φ( ì , ì), say +, and −, . As seen in the example above, for every , one of these two values is approximately the same (mod 2 ). A visually instructive way to determine this common value is to plot each of these pairs as points in a two-dimensional plane, both as ( + , − ) and ( − , + ). Thus, either the horizontal and vertical components of each point are common with all the others. Graphically, this means that the points (approximately) form a vertical and horizontal line intersecting at the true solution ( ( ì), ( ì)). In order to stay in a bounded domain, it is simpler to consider ordered pairs (cos( + ), sin( − )) and (cos( − ), sin( + )) to constrain the search to ∈ [− , ] and remove 2 offsets of the value of ± (see Fig. 3A). Finding the intersection given some noise in the data ± then amounts to the minimization of the error function where the sum is over the pairs of possible solutions for our chosen ì. The optimal value is the desired value of at ì, ( ì) = opt . The landscape of the error function ( ) presents challenges for conjugate gradient optimization because it contains multiple local minima separated by large barriers. An example for a random pixel in a 2D phase array is shown in Fig. 3B. Since the value of the phase needs to be optimized for each pixel on the detector, a rapid and accurate method of determining the absolute minimum is desired. This is most straightforwardly accomplished by supplying an optimization algorithm with the minimum value of log ( ) on a grid in [− , ]; the logarithm increases the contrast of the absolute minima significantly, allowing accurate calculation of an initial guess for the conjugate gradient optimizer. The optimizer polishes the brute-force search to a precise final value. Since the value of ( ì + ì) depends on previous values of ( ì ) and ( ì), it is especially important that values calculated early in the retrieval are accurate. Depending on the quality of the data in Φ( ì , ì), the error function may present multiple deep local minima which cause the algorithm to, initially, choose an incorrect value of . In this case, log[ ( )] for subsequent pixels in the retrieval sharply increases, indicating that at least one previous pixel has assigned incorrectly. Plotting the total fit error for all pixels indicates the location of problematic pixels which require resolving by toggling candidate phase values until the total error of all pixels is minimized. The consequences of toggling alternate values of near the origin to minimize the error across all pixels are illustrated nicely in Fig. 4 by comparing the boxed and unboxed figures. (E-H) shows the same plots for the same set of atoms and number of shots but without alternate phase value toggling. When alternate phase values at pixels adjacent to spikes in the error function are toggled correctly as in (B), the error is reduced and the difference between the retrieved and true values is small. When alternates are not toggled correctly as in (F), the error is large and faithful structure retrieval is less likely. Phase information past the physical edge of the detector ( ì max = 2 in this example) is retrieved via the triple correlations, enhancing real space resolution that would be measured in a typical diffraction experiment.

Results
Using Equation 11 and the algorithm described in Section 3, we can calculate the Fourier phase from triple correlation data and compare the result to the true value. Fig. 4 shows the results of phase retrieval in a simulation with a 2D pixel detector. Note that in regions where the error Fig. 4B is small, the solved phase matches the true phase quite well. In this example, sufficient phase information is retrieved to fully resolve the seven simulated atoms (blue circles in Fig. 4D) with only 10 4 shots. The Fourier inversion was performed using the phase retrieved via our algorithm from the third-order correlation function and (1) calculated via the second-order correlation function. No use of coherent diffraction data was required.
Acquiring additional shots significantly improves the fidelity of phase retrieval. The primary practical limit on phase retrieval via the triple correlations is the computation of the bispectrum which, for a 2D detector, requires storage of a ( pix × pix ) 2 floating point array. For large detectors, the bispectrum can rapidly consume all available memory on small workstations. The 11 × 11 detector simulated in Fig. 4 with 10 4 shots was chosen as a reasonable compromise between memory usage, execution time, and visual impact for this demonstration on a 2015 MacBook Pro running a 2.8GHz Intel i7 processor.

Conclusions
We have described a mathematical solution to the sign problem in phase retrieval from triple correlations of fluorescent or thermal light. We provided a numerical demonstration of our method with simulated data which accurately retrieved the reciprocal space Fourier phase of an atom array. We envision this method as another potential solution to the phase problem in crystallography via photon correlations of fluorescence radiation from atoms pumped by x-ray free-electron lasers. The method presented in this paper may also prove interesting for pulse metrology [32][33][34][35], observing many-body correlations in ultracold atomic gases [36][37][38], imaging in turbid media [39][40][41], and for imaging with radio telescope arrays [12][13][14].

S1.2. Detection of a Photon
The reader may wish to refer to [42][43][44][45] in the next few sections where we derive the expression for intensity triple correlations used in the main text.
Consider some initial state of the photon field | with some detection device in the ground state. | is the final state of the field after our detection system has been excited by a photon. The transition amplitude for this event is then = |ˆ| (S1) and the transition probability is where is the photon field annihilation operator. The final state of the field is usually not observed directly, so we trace over the possible | This short calculation shows how photon detection is analogous to a measurement of the degree of first-order coherence.

S1.3. First-Order Coherence in Diffraction
In coherent diffraction experiments, we measure the degree of first-order coherence as a function of reciprocal space coordinates.
In this picture of elastic scattering, a photon in mode ì 0 is destroyed at site ì and a new photon of the same energy is created in the mode ì at site ì = ì . If we define the momentum transfer vector ì = ì − ì 0 , then For elastic scattering, we can generally regard the state of the light field to be classical. Furthermore, we consider a system with only one atomic species such that = = , and It follows that ( ì 0) = | | 2 so the normalized correlation function is (1) Photodetectors measure only the fringe visibility (the light intensity), eliminating the phase information of the first-order coherence (1)

S1.4. Pair Correlation Function
Light with no first-order coherence (e.g., fluorescence) can still exhibit second-order coherence. The pair-wise correlation function is written in second-quantization as where the angle brackets indicate the expectation value of a quantum state. Expanding singleparticle states in a single mode of momentum-spacê with ∈ [0, 2 ) a random phase that varies slower than the coherence time. In this picture, each atom at ì emits a photon into the mode ì with phase .
The pair-wise correlation averaged over an ensemble of measurements (shots) is The outer braces notate the ensemble mean. There are two cases where the ensemble mean pair-wise correlation is non-zero: We contract the fields accordingly and rewrite the correlation function as By writing out the first two terms with , indices over the number of emitters we double count the case that = . In Equation S15, we have subtracted this case once in the last term. The quantum-classical correspondence between field annihilation/creation operators and classical waves (ˆ( +) ↔ andˆ( −) ↔ * ) with normal ordered correlation functions is implemented with the use of the coherent states |{ } [45]. Fock states have no corresponding classical state and must be treated separately.

S1.4.1. Classical Light
Consider the case that each atom emits a single coherent field so that the initial state of the field is where | is the state contributed by the th atom. Note that since the coherent states are eigenstates ofˆ. Because the correlation function operator is normal ordered, we can directly resume from Equation S15 with Now, the ensemble cancellation of random phases in the E ( ) gives the following approximation Supposing identical emission from each atom, = = , and recalling |E ( )| 2 = 1, we can reduce the sums over constant terms (2) We can rewrite the last term where ì = ì 1 − ì 2 . Normalizing the entire correlation function by dividing out the square of (S23) we can express the normalized second-order correlation function (2) Observe that the last term is identical to the modulus square first-order correlation function from Equation S9. Apparently, the coherent diffraction pattern is recovered from photon pair correlations of fluorescence speckle intensities.
The imperfect cancellation of random phases causes the "phase noise" of Hanbury Brown and Twiss. The more atoms included, the longer random walk in the complex plane requires additional shots to reach the same level of phase noise as with fewer atoms.

S1.4.2. Quantum Light
Consider the case that each atom emits a single photon so that the initial state of the field is defined where | is the state contributed by the th atom. Note that requires an exact match of creation and annihilation operator indices to be non-zero, since Fock states form an orthonormal basis and are not eigenstates ofˆorˆ †. For coherent states, the normal-ordered correlation operator is already diagonalized. For Fock states, this is not so. The commutation relation [ˆ,ˆ †] = allows us to diagonalize the cases S13 and S14 as well as the double-counting correction term in S15 Writing the expectation value ˆ †ˆ = as the occupation number of the mode emitted by the th atom, the diagonalized second order correlation function (Eq. S15) reads Supposing single photon emission from identical atoms, = = 1 and writing E ( 2 ) 2 = 1, we get Normalizing by ˆ † ( ì 0 )ˆ( ì 0 ) 2 = 2 as in the previous section, we obtain which is quite similar to the case where we assume classical light. Note here, however, that the expression is an exact equality -this is due to the orthogonality of Fock states contracting only those random phases which cancel perfectly. There are no − terms where ≠ , so the ensemble average over shots is not subject to the Hanbury Brown and Twiss phase noise.

S1.5. Triple Correlation Function
The correlation function between triples of photons is written in second-quantization as where the field operators are defined as before. The triple correlation averaged over an ensemble of measurements is and again we must consider the field contractions for which the expression is non-zero. We find that there are six = , = , = (S38) The first case (S38) contracts fields with the same ì , and will thus be a constant after summation and averaging. For the remaining cases, we must carefully account for double-and triple-counting. We can quickly see that cases S41, S42, and S43 each contract one set of fields with the same ì , leaving a pair of free indices which do not contract to a constant. Each of these pairs over count S38 when the indices are equal, so we must subtract this case from each of S41, S42, and S43. The cases S39 and S40, however, have three free indices which do not produce field contractions to a constant. Firstly, we can contract any pair of the three indices in each of S39 and S40 to get the result of any one of S41, S42, and S43 again -this over counts each of S41, S42, and S43 once for each of S39 and S40. Additionally, this overcounting correction must itself be corrected just as S41, S42, and S43 were above. Finally, S39 and S40 each overcount the case where all three indices are equal. All cases (S38 through S43) over-counting corrections can be expressed as when contracted with the six-fold sum S37.

S1.5.1. Classical Light
Consider once more our scenario where classical light is emitted As with the pair correlation operator above, the triple correlation operator is diagonal in the coherent state basis. Proceeding analogously and defining we soon reach the expression remembering the over-counting rules from the previous section. Observe that the last two terms are products of independent sums and complex conjugates. Moreover, each of the second, third, and fourth terms is a first-order correlation function S7. After normalizing by the cube of | | 2 , we write The last term here is the subject of the main text and allows retrieval of the exact structural phase ( ì) sought in coherent diffraction experiments.

S1.5.2. Quantum Light
As with the pair correlations of quantum light, the triple correlator acting on non-classical states of the field picks up additional corrections during diagonalization. There are six terms (neglecting the over-counting correction terms, which must also be diagonalized) that must be diagonalized in the Fock basis. As an example, diagonalizing the first term of S44 produces the operator ˆ †ˆ †ˆ †ˆˆˆ = ˆ †ˆˆ †ˆˆ †ˆ−ˆ †ˆˆ †ˆ−ˆ †ˆˆ †ˆ+ˆ †−ˆ †ˆˆ †ˆ+ˆ †ˆ (S51) Clearly, obtaining the exact form of the triple correlations for chaotic quantum light is a tedious task. The result is not important for this paper, but it is important to realize that the classical and quantum cases lead to different expressions. We will skip this derivation for now. S1.6. Friedel's Law and Symmetries of Φ( ì , ì) Given a real function ( ), its Fourier Transform ( ) = F ( ( )) has the following properties: leads to some useful symmetries of the correlation functions and phase map (1)

S1.7. Coherent Diffraction Pattern Retrieval
Fig . S2 shows the result of calculating the pairwise correlation function (2) ( ì) from first-order incoherent fluorescence of atoms. The autocorrelation enables accurate retrieval of the coherent diffraction pattern out to 2| ì max | (past the physical edge of the detector at | ì max |) [6]. in the boxed region of (C) plus extra coherent diffraction data beyond the maximum scattering angle covered by the detector. The augmented ì -space window produces a real space (object) autocorrelation function with enhanced spatial resolution (E) compared to the real space autocorrelation (F) obtained from inverse Fourier inversion of (C) and (D).

S1.8. Extended Momentum Space Sampling Enhances Real Space Resolution
Fig . S3 shows that extended sampling of momentum space enhances the resolution of the real space autocorrelation and, if the phase information is available, the resolution of inverse Fourier Transform (the object).

S1.9. The Cross-Correlation Theorem
For correlations of one-dimensional or small two-dimensional pixel arrays, computationally, it is usually faster to work directly with the outer product of the arrays, sum the resulting matrix along the diagonal, and divide out the underlying support (which represents correlations of unbunched photons). However, for larger arrays, this computation rapidly becomes memorylimited; for example, the second-order correlation of a typical 1000 by 1000 pixel array involves the calculation of a four-dimensional array represented by 1000 4 floating point numbers. If each floating point number were 32 bits, we would need about 4 terabytes to store this array.
The autocorrelation (2) ( 1 , 2 ) of intensities ( ) on a detector is then where the last expression is the convolution of intensities. The triple correlation is most easily calculated from the bispectrum as where = 1 − 2 and = 2 − 3 .

S1.10. Coarse Binning and Interpolation of Phase Information
A major difficulty of performing a triple correlation experiment is the computational complexity of the method. For a square, two-dimensional detector with side length of pixels, the third-order correlation function (correlations between triples of pixels) is a six-dimensional function. This must either be stored in memory or rapidly generated on-the-fly during phase retrieval, neither of which is an attractive option for working with large detectors.
We can consider a case where the sampling rate of the triple correlation may be reduced, the phase retrieved and interpolated, and combined with densely sampled diffraction data. The upscaled phase information breaks the symmetry of the real space structure autocorrelation, selecting peaks at the true atomic positions. For example, Fig. S4 shows that phase information sampled at 30 pixels combined with the modulus from coherent diffraction data sampled at 200 pixels across the same region of ì -space retrieves a noisy, but discernible structure.

S1.11. Linear Phase Ramps and Sign Flips Shift and Invert the Fourier Transform
Fig . S5 demonstrates the effect of linear ramps and sign flips of the phase in real space. A linear phase ramp in momentum space produces a translation in real space. A sign flip of the entire phase inverts the real space structure. Both exact numerical solving and differential evolution structure fitting obtain the phase accurately up to a linear phase ramp or sign flip (or both) in most cases, maintaining the correct relative positions but producing incorrect absolute positions and orientations.

S1.12. Harmonic Inversion Resolves Real Space Positions from Fourier Data
In Fig. S5 we used harmonic inversion to retrieve precise real space positions of the atoms from the Fourier modulus and phase. The inverse Fourier transform gives peaks with a finite width and, thus, uncertainty in the position of the atom. Harmonic inversion is not limited by the Fourier uncertainty principle and superresolves the positions of the atoms directly from the momentum space data rather than the Fourier transform [46][47][48]. S1.13. 1D Phase Retrieval Example S1.13.1. Exact Coherent Phase Retrieval Using Equation 11 and the algorithm described in 3, we can calculate the Fourier phase from triple correlation data and compare the result to the true value. Figure S6 shows the results of phase retrieval in a simulation with a 1D pixel detector.

S1.13.2. Coherent Phase Retrieval via Fitting
Retrieval of the Fourier phase is highly sensitive to statistical and systematic noise. The examples of exact phase retrieval shown in the previous section were calculated for ideal experimental conditions with no sources of noise. Noise is especially problematic in the division step of Equation 11. Moreover, it is expected that the phase noise will become problematic for sources with large numbers of incoherent atoms [6]. Therefore, phase retrieval via an optimization algorithm which compares the closure phase (Equation 12) of a trial atom array to a measured closure phase may be more tractable for real experimental data.
We used a differential evolution algorithm to iterate on the spatial arrangement of the atoms until satisfactory convergence between the simulated closure phase and the closure phase of the trial structure was achieved. It can be seen in Fig. S6 that the exact method based on the algorithm and the structure optimization method achieve similar results.  )) and the same modulus used to produce (C). (D) suffers from additional noise compared to (C), but each of the five atoms remains well-resolved. S1.14. Alternate Phase Toggling Figure S7 illustrates how phase toggling is implemented in the algorithm from the main text.  . Phase retrieval results compared to truth values (blue) for a 1D array of three atoms using the exact numerical solution (green) of Section 3 and differential evolution (orange). (A) shows the retrieved (orange and green) and truth (blue) phases used, together with the modulus obtained from Equation 9, to execute the inverse Fourier transform in (B), (C), and (D). In (B), the truth object from the inverse Fourier transform and the known positions (vertical dots) of the atoms match, unsurprisingly. (C) shows that differential evolution correctly found atom positions (orange vertical dots) with the correct separation and ordering, but the structure appears flipped and shifted compared to the truth positions (blue vertical dots). This is due to the linear phase ramp in (A) compared to the truth phase. The effect of linear phase ramps and sign flips is detailed in Supplement S1.11. Similarly, the exact numerical solution in (D) also retrieved the correct interatomic separations but the phase (green) in (A) is inverted to the truth phase, causing a flip of the real space structure. Since the exact numerical algorithm finds the phase without directly determining atom positions, we used harmonic inversion (see Supplement S1.12) to determine a sub-linewidth position for the three atoms (green vertical dots) to compare to the truth atom positions (blue vertical dots). The Fourier transforms (solid curves) are each normalized to the height of the tallest peak. = − + , ∈ Z is determined beginning at the origin and ending in the upper right corner. In the first example (left), where toggling is not used, the solved value of in (A) has a large difference (F) from the true value of in (B). We observe that the value of log [ ( )] (E) is seen to spike (red box) along the third diagonal (green and red boxes). The value of the error function increases and propagates to surrounding pixels as additional diagonals are solved. This indicates that a pixel in the second diagonal (blue box) has been assigned incorrectly. Resolving the second diagonal by substituting alternate minima of log [ ( )] eventually arrives at a ( ì ) (C) which has uniformly low values of the error function (G) such that the difference (H) between the solution (C) and truth (D) is quite small. Note that the error function, (E) and (G), is defined to be zero at the origin and the origin nearest neighbors as these values of ( ì ) are free parameters of the solution process.