# Phase Retrieval in Incoherent Diffractive Imaging using higher-order photon
correlation functions

M. Bojer  Friedrich-Alexander-Universität Erlangen-Nürnberg, Quantum Optics
and Quantum Information, Staudtstr. 1, 91058 Erlangen, Germany  J. Eckert
Friedrich-Alexander-Universität Erlangen-Nürnberg, Quantum Optics and Quantum
Information, Staudtstr. 1, 91058 Erlangen, Germany  S. Karl  Friedrich-
Alexander-Universität Erlangen-Nürnberg, Quantum Optics and Quantum
Information, Staudtstr. 1, 91058 Erlangen, Germany  S. Richter  Friedrich-
Alexander-Universität Erlangen-Nürnberg, Quantum Optics and Quantum
Information, Staudtstr. 1, 91058 Erlangen, Germany  J. von Zanthier
Friedrich-Alexander-Universität Erlangen-Nürnberg, Quantum Optics and Quantum
Information, Staudtstr. 1, 91058 Erlangen, Germany

###### Abstract

To obtain spatial information about an arbitrary object in x-ray structure
analysis, the standard method is to measure the intensity in the far field,
i.e., the first-order photon correlation function of the coherently scattered
x-ray photons (coherent diffractive imaging). Recently, it was suggested to
record alternatively the incoherently scattered photons and measure the
second-order photon correlation function to reconstruct the geometry of the
unknown object (incoherent diffractive imaging). Yet, besides various
advantages of the latter method, both techniques suffer from the so-called
phase retrieval problem. Lately, an ab-initio phase retrieval algorithm to
reconstruct the phase of the so-called structure factor of the scattering
objects based on the third-order photon correlation function was reported. The
algorithm makes use of the so-called closure phase, which contains important,
yet incomplete phase information, well-known from triple correlations and
their bispectrum in speckle masking and astronomy applications. Here, we
provide a detailed analysis of the underlying scheme and quantities in the
context of x-ray structure analysis. In particular, we explicitly calculate
the third-order photon correlation function in a full quantum mechanical
treatment and discuss the uniqueness of the closure phase equations
constructed from it. In this context, we recapitulate the sign problem of the
closure phase and how it can be lifted using redundant information. We further
show how the algorithm can be improved using even higher-order photon
correlation functions, e.g., the fourth-order correlation function, delivering
new phase relations appearing in the four-point correlations.

##  I Introduction

The goal of so-called coherent diffractive imaging (CDI) in x-ray structure
analysis is to reconstruct the geometry of an object via the detection of the
intensity of the coherently scattered x-ray photons in the far field. The
detected intensity distribution is thereby proportional to the modulus squared
of the so-called structure factor, the Fourier transform of the real space
distribution of the object. To reconstruct the object via an inverse Fourier
transform, the phase of the structure factor is thus needed. However, since
only the intensity of the scattered radiation is measured, the phase
information is missing. This is the so-called phase problem of CDI. Several
so-called phase retrieval algorithms have been developed over the last years
to recover the phase when only the intensity of the scattered light is known
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10].
Recently, it was suggested to measure the incoherently scattered light, which
for CDI is usually considered to be detrimental since it leads merely to a
constant background, in order to reconstruct the geometry of the object via
so-called incoherent diffractive imaging (IDI) [11, 12, 13]. The idea is to
use the second-order photon correlation function instead of the intensity,
i.e., the first-order photon correlation function. It can be shown that the
second-order correlation function of the incoherently scattered photons
equally gives access to the modulus squared of the structure factor. Thus,
IDI, equally to CDI, requires the phase of the structure factor to be
determined in an independent manner. It is well-known from speckle masking and
astronomy applications that triple correlations and their bispectrum contain
some useful yet incomplete phase information about the structure factor via
the so-called closure phase [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25].
However, since only the cosine of the closure phase is available, a sign
ambiguity arises. In that context, an ab-initio phase retrieval algorithm
based on the third-order photon correlation function was recently proposed
[26], which demonstrates how to get rid of the sign ambiguity using redundant
information. Here, we analyse the sign ambiguity problem in more detail and
show how even higher-order photon correlation functions, e.g., the fourth-
order photon correlation function, can be used to further reduce the
ambiguities.
The paper is organized as follows: In Sec. II, we introduce the so-called
structure factor, the Fourier transform of the real space distribution of the
object. In Sec. III, we recapitulate the basic working principle of coherent
diffractive imaging and explain the phase problem. In Sec. IV, we first
introduce the so-called complex degree of coherence and show its relation to
the structure factor. We then explain in detail how the second- and third-
order photon correlation functions can be used to reconstruct the missing
phase via the so-called closure phase. In this context, we analyse the sign
ambiguity of the closure phase and recapitulate the proposed algorithm of Ref.
[26]. We end this section by demonstrating how the new phase relations
occurring for four-point correlations in the fourth-order photon correlation
function can be used to further reduce the previously discussed ambiguity. We
summarise our results and finally conclude in Sec. VI.

##  II Structure factor

In order to model the structure of a complex object, we consider $N$ point-
like atoms, assumed to be two-level systems, located at positions
$\bm{R}_{i}$, $i\in\\{1,...,N\\}$. The real space distribution is thus given
by

$S(\bm{r})=\sum_{i=1}^{N}\delta(\bm{r}-\bm{R}_{i})\,.$ (1)

The Fourier transform of the real space distribution is the so-called
structure factor denoted by

$\displaystyle\tilde{S}(\bm{q})$ $\displaystyle=\int
d^{3}r\,S(\bm{r})e^{-i\bm{q}\bm{r}}$ $\displaystyle=\sum_{i=1}^{N}\int
d^{3}r\,\delta(\bm{r}-\bm{R}_{i})e^{-i\bm{q}\bm{r}}=\sum_{i=1}^{N}e^{-i\bm{q}\bm{R}_{i}}\,.$
(2)

We can write the structure factor in terms of its magnitude and phase as

$\tilde{S}(\bm{q})=|\tilde{S}(\bm{q})|e^{i\phi(\bm{q})}\,,$ (3)

which makes it evident that both the magnitude and phase are needed to
reconstruct the real space distribution via an inverse Fourier transform

$S(\bm{r})=\int
d^{3}q\,|\tilde{S}(\bm{q})|e^{i\phi(\bm{q})}e^{i\bm{q}\bm{r}}\,.$ (4)

Now, assuming a real-valued real space distribution $S(\bm{r})\in\mathbb{R}$,
we find

$\displaystyle\tilde{S}(-\bm{q})$
$\displaystyle=|\tilde{S}(-\bm{q})|e^{i\phi(-\bm{q})}=\int
d^{3}r\,S(\bm{r})e^{i\bm{q}\bm{r}}$
$\displaystyle=\tilde{S}^{*}(\bm{q})=|\tilde{S}(\bm{q})|e^{-i\phi(\bm{q})}$
(5)

and thus

$|\tilde{S}(-\bm{q})|=|\tilde{S}(\bm{q})|\,,\quad\phi(-\bm{q})=-\phi(\bm{q})\,,$
(6)

from which, in particular, $\phi(0)=-\phi(0)=0$ follows. By writing Eq. (1),
we implicitly defined an origin located at some point $\bm{R}$. Let us assume
that we would have chosen a different origin at $\bm{R}^{\prime}$ with
$\bm{R}^{\prime}=\bm{R}+\Delta\bm{R}$. With respect to the new origin, the
real space distribution reads

$S^{\prime}(\bm{r})=\sum_{i=1}^{N}\delta(\bm{r}-\bm{R}^{\prime}_{i})=\sum_{i=1}^{N}\delta(\bm{r}-\bm{R}_{i}+\Delta\bm{R})=S(\bm{r}+\Delta\bm{R})\,,$
(7)

where we used $\bm{R}^{\prime}_{i}=\bm{R}_{i}-\Delta\bm{R}$. That is, we
obtain a shifted real space distribution. Calculating the structure factor of
the shifted object leads to

$\displaystyle\tilde{S}^{\prime}(\bm{q})$
$\displaystyle=|\tilde{S}^{\prime}(\bm{q})|e^{i\phi^{\prime}(\bm{q})}=\int
d^{3}r\,S(\bm{r}+\Delta\bm{R})e^{-i\bm{q}\bm{r}}$
$\displaystyle=e^{i\bm{q}\Delta\bm{R}}\tilde{S}(\bm{q})=|\tilde{S}(\bm{q})|e^{i[\phi(\bm{q})+\bm{q}\Delta\bm{R}]}\,.$
(8)

Thus, the new magnitude and phase are given by
$|\tilde{S}^{\prime}(\bm{q})|=|\tilde{S}(\bm{q})|$ and
$\phi^{\prime}(\bm{q})=\phi(\bm{q})+\bm{q}\Delta\bm{R}$, i.e., the shift of
the origin leads to a linear shift in the phase. We will come to this outcome
later when we discuss the phase retrieval algorithm based on higher-order
photon correlation functions. We end this section by noting that the aim of
x-ray structure analysis is to determine the real space distribution
$S(\bm{r})$ via scattered x-ray photons.

##  III Coherent diffractive imaging

In this section, we recall some basic concepts of CDI, where it is assumed
that the incoming photons are coherently scattered by the atoms. The outgoing
electric field can be written as a superposition of the incoming field and a
scattered spherical wave field, which is proportional to the so-called
scattering factor $f(\bm{q})$, where $\bm{q}$ denotes the momentum transfer of
the scattering process. The scattering factor is thus proportional to the
structure factor

$f(\bm{q})\propto\tilde{S}(\bm{q})=\int
d^{3}r\,S(\bm{r})e^{-i\bm{q}\bm{r}}\,.$ (9)

Since the energy of the photons is unchanged in a coherent scattering process,
the momentum transfer vectors $\bm{q}$ are confined to the so-called Ewald
sphere. Therefore, to get full 3D information via detection of photons on a
pixelated detector, the object has to be rotated. Moreover, to determine the
real space distribution $S(\bm{r})$ via a Fourier transform of
$\tilde{S}(\bm{q})$, full information about the complex function
$\tilde{S}(\bm{q})$ is needed, i.e., the modulus and the phase have to be
known. However, the measured intensity of the coherently scattered photons is
merely proportional to $|f(\bm{q})|^{2}\propto|\tilde{S}(\bm{q})|^{2}$. That
is, one has access to the absolute value, but not to the phase of the
scattering factor. This is the so-called phase problem of CDI. Many phase
retrieval algorithms have been developed to reconstruct the scattering phases,
e.g., where reasonable constraints are imprinted by use of suitable
projections [1, 2, 3, 6], via Cauchy-Riemann approaches using the Cauchy-
Riemann equations [4], ptychographical methods [8, 9], or so-called greedy
methods [5]. In the following, we consider and extend a phase retrieval
algorithm based on the third-order photon correlation function using
incoherently scattered photons [26].

##  IV Incoherent diffractive imaging

In this section, we consider the case of IDI, where the incoherently emitted
fluorescence light of the atoms is measured in the far field and the second-
order photon correlation function is computed. We consider again $N$ two-level
atoms with excited state $\ket{e}$ and ground state $\ket{g}$. Neglecting the
dipole radiation pattern and a general constant prefactor, which would provide
a meaningful physical unit, the positive electric field operator at time $t=0$
is given by [27]

$\hat{E}^{(+)}(\bm{r})=\sum_{i=1}^{N}e^{-ik_{0}\hat{\bm{r}}\bm{R}_{i}}\hat{S}_{-}^{(i)}\,.$
(10)

Here, $k_{0}=2\pi/\lambda$ denotes the wave number, $\hat{\bm{r}}$ the unit
vector in direction of the observation, $\bm{R}_{i}$ the position of the $i$th
atom, and $\hat{S}_{-}^{(i)}$ the pseudo-spin lowering operator of the $i$th
atom. For convenience, we define $\bm{k}\coloneqq k_{0}\hat{\bm{r}}$. In this
way, we can write the positive electric field operator as a function of
$\bm{k}$ as

$\hat{E}^{(+)}(\bm{k})=\sum_{i=1}^{N}e^{-i\bm{k}\bm{R}_{i}}\hat{S}_{-}^{(i)}\,.$
(11)

We start by calculating the far field intensity produced by the atoms, which
can be written as [28]

$\displaystyle I(\bm{k})$
$\displaystyle=G^{(1)}(\bm{k},\bm{k})=\braket{\hat{E}^{(-)}(\bm{k})\hat{E}^{(+)}(\bm{k})}$
$\displaystyle=\sum_{i,j=1}^{N}e^{i\bm{k}\bm{R}_{i}}e^{-i\bm{k}\bm{R}_{j}}\braket{\hat{S}_{+}^{(i)}\hat{S}_{-}^{(j)}}\,.$
(12)

Assuming the atoms to be initially in the fully excited state end thus emit
the photons incoherently, the expectation value
$\braket{\hat{S}_{+}^{(i)}\hat{S}_{-}^{(j)}}$ is zero unless $i=j$, whereas
$\braket{\hat{S}_{+}^{(i)}\hat{S}_{-}^{(i)}}=1\,\forall i$. The intensity is
then given by

$I(\bm{k})=\sum_{i=1}^{N}\braket{\hat{S}_{+}^{(i)}\hat{S}_{-}^{(i)}}=N\,,$
(13)

which is a constant and thus gives no phase information about the real space
distribution. However, in contrast to the intensity, the complex degree of
coherence (CDC) $g^{(1)}(\bm{k}_{1},\bm{k}_{2})$, i.e., the cross-correlation
between the electric fields in the directions $\bm{k}_{1}$ and $\bm{k}_{2}$
does reveal structure information since

$\displaystyle g^{(1)}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\frac{G^{(1)}(\bm{k}_{1},\bm{k}_{2})}{\sqrt{G^{(1)}(\bm{k}_{1},\bm{k}_{1})}\sqrt{G^{(1)}(\bm{k}_{2},\bm{k}_{2})}}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}e^{-i\bm{q}\bm{R}_{i}}=\frac{1}{N}\tilde{S}(\bm{q})=\frac{\tilde{S}(\bm{q})}{\tilde{S}(0)}\,,$
(14)

where $\bm{q}=\bm{k}_{2}-\bm{k}_{1}$ is the momentum difference vector. As
found in Eq. (IV), measuring the CDC of the emitted fluorescence light of the
$N$ atoms would give full information about the real space distribution
$S(\bm{r})$. However, determining experimentally the CDC is usually very hard
[11]. Instead, we investigate the possibility of using higher-order photon
correlation functions to get access to the structure factor when measuring
incoherent fluorescence light in the following. As we show in App. A, the
second-order photon correlation function calculates to

$\displaystyle
g^{(2)}(\bm{k}_{1},\bm{k}_{2})=1-\frac{2}{N}+\frac{1}{N^{2}}|\tilde{S}(\bm{q})|^{2}\,.$
(15)

According to Eq. (15), $g^{(2)}(\bm{k}_{1},\bm{k}_{2})$ also gives access to
the modulus of the structure factor, however not to the phase.

##  V Phase retrieval using higher-order photon correlation functions

In order to get access to the phase of the structure factor
$\tilde{S}(\bm{q})$, one can consider the phase relations appearing in the
third-order photon correlation function. We start indeed by considering the
third-order correlation function, however, will subsequently also investigate
the fourth-order photon correlation function for this purpose.

###  V.1 Third-order photon correlation function

The third-order photon correlation function calculates to (see App. B)

$\displaystyle g^{(3)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=1-\frac{6}{N}+\frac{12}{N^{2}}+\frac{N-4}{N^{3}}[|\tilde{S}(\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}]$
$\displaystyle+\frac{2}{N^{3}}|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|\cos[\phi(\bm{q}_{1}+\bm{q}_{2})-\phi(\bm{q}_{1})-\phi(\bm{q}_{2})]\,,$
(16)

where we denote $\bm{q}_{1}=\bm{k}_{2}-\bm{k}_{1}$ and
$\bm{q}_{2}=\bm{k}_{3}-\bm{k}_{2}$. As can be seen from Eq. (16), the third-
order correlation function involves the so-called closure phase

$\Phi(\bm{q}_{1},\bm{q}_{2})\coloneqq\phi(\bm{q}_{1}+\bm{q}_{2})-\phi(\bm{q}_{1})-\phi(\bm{q}_{2})\,,$
(17)

which is well-known also in astronomy. Since both the number of atoms $N$ as
well as the magnitude of the structure factor $|\tilde{S}(\bm{q})|$ can be
determined by the second-order photon correlation function, we can focus on
the closure phase in the following. We see that Eq. (16) gives actually only
access to the cosine of the closure phase, leading to a sign (and a $2\pi$
shift) ambiguity. However, as recently shown, the sign ambiguity can be lifted
in certain cases, such that a complete phase retrieval using the third-order
photon correlation function is possible. We will present a detailed analysis
of the number of non-trivial inequivalent equations and the number of unknowns
in the following and review the possibilities to unlift the sign ambiguity. We
then show how the algorithm can be improved by additionally using the fourth-
order photon correlation function.

###  V.2 Uniqueness of closure phase equations and sign ambiguity lifting

In the linear regime, i.e., in the far-field, the vectors $\bm{q}_{1}$ and
$\bm{q}_{2}$ can be considered as difference vectors between different
positions of pixels on a camera. Therefore, we will adapt our notation in the
following and write

$\Phi(\bm{m},\bm{n})=\phi(\bm{m}+\bm{n})-\phi(\bm{m})-\phi(\bm{n})$ (18)

with discrete difference vectors $\bm{m}$ and $\bm{n}$ with natural number
entries. We now check the uniqueness of these equations for two different
scenarios, namely a pixelated camera in one dimension and two dimensions,
respectively.
Let us start by assuming that we have a camera of $M$ pixels along one
dimension. Without loss of generality, we can set the first pixel of the
triple correlations to be the first pixel of the camera. Now, we consider an
arbitrary pixel defined by the difference value $u\in\\{0,...,M-1\\}$ with
respect to the first pixel. The question now is, how many possibilities do we
have to achieve $m+n=u$ with $m,n\leq u$ and $m,n\in\mathbb{N}$. There are
exactly $u+1$ possibilities. This leads to a total number of equations given
by

$\sum_{u=0}^{M-1}(u+1)=\sum_{u=1}^{M}u=\frac{1}{2}M(M+1)\,.$ (19)

However, not all of these equations are independent. In fact, all equations
with $m$ and $n$ interchanged are equivalent since $\Phi(m,n)=\Phi(n,m)$
(redundant equations). In addition, $\phi(0)=0$, such that
$\Phi(m,0)=\Phi(0,m)=0$ does not give additional information either (trivial
equations). Let us first count the number of redundant equations, which for
$M$ even reads $\frac{M^{2}}{4}$, whereas for $M$ odd is given by
$\frac{M^{2}-1}{4}$. Further, the number of trivial equations is given by $M$.
Therefore, the number of non-trivial inequivalent equations is given by

$\displaystyle F^{(1)}_{\mathrm{even}}(M)$
$\displaystyle\coloneqq\frac{1}{2}M(M+1)-\frac{1}{4}M^{2}-M=\frac{1}{4}M(M-2)\,,$
(20) $\displaystyle F^{(1)}_{\mathrm{odd}}(M)$
$\displaystyle\coloneqq\frac{1}{2}M(M+1)-\frac{1}{4}(M^{2}-1)-M=\frac{1}{4}(M-1)^{2}\,.$
(21)

In contrast, the number of unknown variables is
$F^{(1)}_{\mathrm{even/odd}}(M)+M-2$, where $F^{(1)}_{\mathrm{even/odd}}(M)$
gives the number of closure phases $\Phi(m,n)$, i.e., essentially signs since
the value of $|\Phi(m,n)|$ is known, and $M-2$ phases $\phi(m)$, namely
$\phi(2),...,\phi(M-1)$ since $\phi(0)=0$ and $\phi(1)$ can be arbitrarily
chosen, as it corresponds to a global shift, i.e., it defines the origin of
the coordinate system [see Eq. (II)]. This means, in principle the number of
unknowns is greater than the number of equations. However, since many of the
unknowns are only signs, in certain cases it is possible to find a unique
solution. To exemplify this statement, let us consider the example of a one-
dimensional camera with 9 pixel. In this case, we have
$F^{(1)}_{\mathrm{odd}}(9)=16$ non-trivial inequivalent equations, namely

$\displaystyle\Phi(1,1)$ $\displaystyle=\phi(2)-2\phi(1)$ (22)
$\displaystyle\Phi(2,1)$ $\displaystyle=\phi(3)-\phi(2)-\phi(1)$ (23)
$\displaystyle\Phi(3,1)$ $\displaystyle=\phi(4)-\phi(3)-\phi(1)$ (24)
$\displaystyle\Phi(4,1)$ $\displaystyle=\phi(5)-\phi(4)-\phi(1)$ (25)
$\displaystyle\Phi(5,1)$ $\displaystyle=\phi(6)-\phi(5)-\phi(1)$ (26)
$\displaystyle\Phi(6,1)$ $\displaystyle=\phi(7)-\phi(6)-\phi(1)$ (27)
$\displaystyle\Phi(7,1)$ $\displaystyle=\phi(8)-\phi(7)-\phi(1)$ (28)
$\displaystyle\Phi(2,2)$ $\displaystyle=\phi(4)-2\phi(2)$ (29)
$\displaystyle\Phi(3,2)$ $\displaystyle=\phi(5)-\phi(3)-\phi(2)$ (30)
$\displaystyle\Phi(4,2)$ $\displaystyle=\phi(6)-\phi(4)-\phi(2)$ (31)
$\displaystyle\Phi(5,2)$ $\displaystyle=\phi(7)-\phi(5)-\phi(2)$ (32)
$\displaystyle\Phi(6,2)$ $\displaystyle=\phi(8)-\phi(6)-\phi(2)$ (33)
$\displaystyle\Phi(3,3)$ $\displaystyle=\phi(6)-2\phi(3)$ (34)
$\displaystyle\Phi(4,3)$ $\displaystyle=\phi(7)-\phi(4)-\phi(3)$ (35)
$\displaystyle\Phi(5,3)$ $\displaystyle=\phi(8)-\phi(5)-\phi(3)$ (36)
$\displaystyle\Phi(4,4)$ $\displaystyle=\phi(8)-2\phi(4)-\phi(1)\,,$ (37)

but also $16$ unknown signs of the closure phases $\Phi(m,n)$ and $7$ phases
$\phi(m)$, i.e., $\phi(2),...,\phi(8)$. We will show how to solve this problem
later.
Let us now assume we have a two-dimensional camera of $M^{2}$ pixels. Without
loss of generality, we can set the first pixel of the triple correlations to
be the pixel at the left lower corner of the camera. Now, we consider a pixel
defined by the difference vector $\bm{u}$ with respect to the first pixel. The
question is, how many possibilities do we have to achieve
$\bm{m}+\bm{n}=\bm{u}$ with $m_{x/y},n_{x/y}\leq u_{x/y}$ and
$\bm{m},\bm{n}\in\mathbb{N}^{2}$. There are exactly $(u_{x}+1)(u_{y}+1)$
possibilities. This leads to a total number of equations given by

$\sum_{u_{x},u_{y}=0}^{M-1}(u_{x}+1)(u_{y}+1)=\frac{1}{4}M^{2}(M+1)^{2}\,.$
(38)

However, not all of these equations are independent. In fact, all equations
with $\bm{m}$ and $\bm{n}$ interchanged are equivalent since
$\Phi(\bm{m},\bm{n})=\Phi(\bm{n},\bm{m})$ (redundant equations). In addition,
$\phi(0)=0$, such that $\Phi(\bm{m},0)=\Phi(0,\bm{m})=0$ does not give
additional information either (trivial equations). The number of redundant
equations for $M$ even reads $\frac{1}{8}M^{3}(M+2)$, whereas for $M$ odd is
given by $\frac{1}{8}(M-1)(M+1)^{3}$. Further, the number of trivial equations
is given by $M^{2}$. Therefore, the number of non-trivial inequivalent
equations is given by

$\displaystyle F^{(2)}_{\mathrm{even}}(M)$
$\displaystyle\coloneqq\frac{1}{4}M^{2}(M+1)^{2}-\frac{1}{8}M^{3}(M+2)-M^{2}=\frac{1}{8}M^{2}[M(M+2)-6]\,,$
(39) $\displaystyle F^{(2)}_{\mathrm{odd}}(M)$
$\displaystyle\coloneqq\frac{1}{4}M^{2}(M+1)^{2}-\frac{1}{8}(M-1)(M+1)^{3}-M^{2}=\frac{1}{8}(M-1)^{2}[M(M+4)+1]\,.$
(40)

In contrast, the number of unknown variables is
$F^{(2)}_{\mathrm{even/odd}}(M)+M^{2}-3$, where
$F^{(2)}_{\mathrm{even/odd}}(M)$ gives the number of closure phases
$\Phi(\bm{m},\bm{n})$, i.e., essentially signs since the value of
$|\Phi(\bm{m},\bm{n})|$ is known, and $M^{2}-3$ phases $\phi(\bm{m})$, since
$\phi(0)=0$ and $\phi((1,0)^{T})$, $\phi((0,1)^{T})$ can be arbitrarily
chosen, as they correspond to a global shift, i.e, they define the origin of
the coordinate system [see Eq. (II)]. Let us consider the example of a camera
with 9x9 pixels, i.e., $M=3$. In this case, we have
$F^{(2)}_{\mathrm{odd}}(3)=11$ non-trivial inequivalent equations, namely

$\displaystyle\Phi((1,0)^{T},(1,0)^{T})$
$\displaystyle=\phi((2,0)^{T})-2\phi((1,0)^{T})$ (41)
$\displaystyle\Phi((0,1)^{T},(0,1)^{T})$
$\displaystyle=\phi((0,2)^{T})-2\phi((0,1)^{T})$ (42)
$\displaystyle\Phi((1,0)^{T},(0,1)^{T})$
$\displaystyle=\phi((1,1)^{T})-\phi((1,0)^{T})-\phi((0,1)^{T})$ (43)
$\displaystyle\Phi((1,0)^{T},(1,1)^{T})$
$\displaystyle=\phi((2,1)^{T})-\phi((1,0)^{T})-\phi((1,1)^{T})$ (44)
$\displaystyle\Phi((2,0)^{T},(0,1)^{T})$
$\displaystyle=\phi((2,1)^{T})-\phi((2,0)^{T})-\phi((0,1)^{T})$ (45)
$\displaystyle\Phi((0,1)^{T},(1,1)^{T})$
$\displaystyle=\phi((1,2)^{T})-\phi((0,1)^{T})-\phi((1,1)^{T})$ (46)
$\displaystyle\Phi((0,2)^{T},(1,0)^{T})$
$\displaystyle=\phi((1,2)^{T})-\phi((0,2)^{T})-\phi((1,0)^{T})$ (47)
$\displaystyle\Phi((1,0)^{T},(1,2)^{T})$
$\displaystyle=\phi((2,2)^{T})-\phi((1,0)^{T})-\phi((1,2)^{T})$ (48)
$\displaystyle\Phi((2,0)^{T},(0,2)^{T})$
$\displaystyle=\phi((2,2)^{T})-\phi((2,0)^{T})-\phi((0,2)^{T})$ (49)
$\displaystyle\Phi((2,1)^{T},(0,1)^{T})$
$\displaystyle=\phi((2,2)^{T})-\phi((2,1)^{T})-\phi((0,1)^{T})$ (50)
$\displaystyle\Phi((1,1)^{T},(1,1)^{T})$
$\displaystyle=\phi((2,2)^{T})-2\phi((1,1)^{T})\,,$ (51)

but also $11$ unknown signs of the closure phases $\Phi(\bm{m},\bm{n})$ plus
six unknown phases $\phi(\bm{m})$.
We now review the idea how to lift the sign ambiguity referring to Ref. [26].
Let us consider again our 1D example of 9 pixels. Since $\phi(1)$ can be
arbitrarily chosen, Eq. (22) gives two possible values for $\phi(2)$ due to
the unknown sign of $\Phi(1,1)$. Further, Eq. (23) gives four possible values
for $\phi(3)$ since $\phi(2)$ has two possible values and due to the unknown
sign of $\Phi(2,1)$. Then, Eq. (24) gives eight possible values for $\phi(4)$.
However, this equation is not the only equation including previous phases.
Instead, Eq. (29) is an additional equation for $\phi(4)$ leading to four
possible values. Since both equations have to give the same value for
$\phi(4)$, the idea is to check for the intersection of the set containing the
eight possible values coming from Eq. (24) and the set containing the four
possible values coming from Eq. (29). If there is only one equal number in
both sets, we know that $\phi(4)$ has to be this number, which additionally
determines also $\phi(2)$ and $\phi(3)$. We can then go on with Eqs. (25) and
(30) for $\phi(5)$ and check again the intersection. To find the remaining
phases, we just go on with this procedure. In certain cases, this leads to a
unique solution for the phase function of the structure factor, such that full
information is obtained. However, there is an additional option to further
reduce the ambiguity in the retrieval of the phase and to improve this
algorithm, namely by considering in addition the fourth-order photon
correlation function giving an additional restrictions.

###  V.3 Fourth-order photon correlation function

As we show in App. C, the fourth-order photon correlation function leads to
different additional cosine terms of phases, i.e., it gives further
restrictions to the possible phase values. Besides the usual closure phases
as, for instance,
$\phi(\bm{q}_{1}+\bm{q}_{2})-\phi(\bm{q}_{1})-\phi(\bm{q}_{2})$, we find for
example also terms like
$\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1})-\phi(\bm{q}_{2})-\phi(\bm{q}_{3})$
or
$\phi(\bm{q}_{1})-\phi(\bm{q}_{3})+\phi(\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2})$,
which include four phase terms. To show how these additional terms can be
beneficial, let us consider again the 1D example. Here, we will construct by
hand a scenario where the fourth-order correlation function can resolve an
ambiguity, which is not resolvable by the third-order correlation function.
Even though this is a particular example, it shows how the fourth-order
correlation function can be helpful if no unique intersection value is found
for the reconstructed phases via the closure relation based on
$g^{(3)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$. Let us assume that the first four
closure phases are $\Phi(1,1)=\Phi(2,1)=\Phi(3,1)=\Phi(2,2)=\pm 1$. Without
loss of generality, we can choose $\phi(1)=0$, such that
$\phi(2)=\Phi(1,1)=\pm 1$. Let us neglect a possible point reflection of the
configuration, such that we can keep the plus sign for simplicity. This leads
to two possible values for $\phi(3)$, namely $\phi(3)\in\\{2,0\\}$ via Eq.
(23). Now, using Eq. (24) we find four possible values for $\phi(4)$, namely
$\phi(4)\in\\{3,1,1,-1\\}$ with different combinations of $\Phi(3,1)$ and
$\phi(3)$. However, via Eq. (29) we have an additional restriction for
$\phi(4)$, which leads to the set of possible values given by
$\phi(4)\in\\{3,1\\}$. As we find, the intersection of the two sets is not
unique. However, considering the first non-equivalent value of $g^{(4)}$ we
find (see App. C:

$\displaystyle g^{(4)}(1,1,1)$
$\displaystyle=1-\frac{12}{N}+\frac{1}{N^{2}}(64+3|\tilde{S}(1)|^{2}+2|\tilde{S}(2)|^{2}+|\tilde{S}(3)|^{2})+\frac{1}{N^{3}}\\{-144-34|\tilde{S}(1)|^{2}-24|\tilde{S}(2)|^{2}-10|\tilde{S}(3)|^{2}$
$\displaystyle+4|\tilde{S}(1)|^{2}|\tilde{S}(2)|\cos[\phi(2)-2\phi(1)]+4|\tilde{S}(1)||\tilde{S}(2)||\tilde{S}(3)|\cos[\phi(3)-\phi(2)-\phi(1)]\\}$
$\displaystyle+\frac{1}{N^{4}}\\{|\tilde{S}(1)|^{4}+|\tilde{S}(2)|^{4}+2|\tilde{S}(1)|^{2}|\tilde{S}(2)|^{2}|\tilde{S}(1)|^{2}|\tilde{S}(3)|^{2}+96|\tilde{S}(1)|^{2}+68|\tilde{S}(2)|^{2}+32|\tilde{S}(3)|^{2}+4|\tilde{S}(4)|^{2}$
$\displaystyle-28|\tilde{S}(1)|^{2}|\tilde{S}(2)|\cos[\phi(2)-2\phi(1)]-28|\tilde{S}(1)||\tilde{S}(2)||\tilde{S}(3)|\cos[\phi(3)-\phi(2)-\phi(1)]$
$\displaystyle-4|\tilde{S}(2)|^{2}|\tilde{S}(4)|\cos[\phi(4)-2\phi(2)]-4|\tilde{S}(1)||\tilde{S}(3)||\tilde{S}(4)|\cos[\phi(4)-\phi(3)-\phi(1)]$
$\displaystyle+2|\tilde{S}(1)|^{3}|\tilde{S}(3)|\cos[{\color[rgb]{1,0,1}\phi(3)-3\phi(1)}]+2|\tilde{S}(1)||\tilde{S}(2)|^{2}|\tilde{S}(3)|\cos[{\color[rgb]{1,0,1}\phi(3)-2\phi(2)+\phi(1)}]\\}$
(52)

Thus, in the last line we find two new phase relations involving $\phi(3)$,
namely $\phi(3)-3\phi(1)$ and $\phi(3)-2\phi(2)+\phi(1)$. If we come back to
our explicit example, we would have $\phi(3)-3\phi(1)=\phi(3)\in\\{2,0\\}$ and
$\phi(3)-2\phi(2)+\phi(1)=\phi(3)-2\phi(2)\in\\{0,-2\\}$. Consequently, since
the cosines give different values for the different possible phase relations,
it is possible to determine the correct value of $\phi(3)$ via a $g^{(4)}$
measurement, which then reduces the possible values for the other phases. Let
us, for instance, assume that the measurement revealed that $\phi(3)=0$ is the
correct value. Then, the possible values for $\phi(4)$ from Eq. (24) reduce to
$\phi(4)\in\\{1,-1\\}$. Now, the intersection with the set $\\{3,1\\}$ is
unique and we can immediately identify the value of $\phi(4)$ to be $1$. We
emphasize again that the additional constraints given by the fourth-order
correlation function $g^{(4)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})$ is
always helpful if no unique intersection value is found for the reconstructed
phases via the closure relation based on
$g^{(3)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$.

##  VI Conclusion

In conclusion, we have shown how higher-order photon correlation functions can
be useful and beneficial to retrieve the phase of the structure factor
$\tilde{S}(\bm{q})$, which is not accessible via CDI or IDI. To that aim, we
have provided a complete quantum mechanical derivation of the photon
correlation functions of order $3$ and $4$, displaying the closure phase and
additional phase relations. We have demonstrated in particular that by using
the fourth-order photon correlation function where new non-equivalent phase
relations appear due to four-point correlations the phase retrieval process
can be improved. Furthermore, we have shown by explicitly counting the number
of equations and unknowns that even though the sign ambiguity of the closure
phase equations can be lifted in certain cases, the closure phase equations
are not unique. We believe that our findings provide new input in the route of
using higher-order photon correlation functions for phase reconstruction.

## Acknowledgements

M.B. and J.v.Z. gratefully acknowledge funding and support by the
International Max Planck Research School - Physics of Light. S.R. and J.v.Z.
gratefully acknowledge funding and support by the Erlangen Graduate School in
Advanced Optical Technologies (SAOT). This work was funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID
429529648 – TRR 306 QuCoLiMa (”Quantum Cooperativity of Light and Matter”).

##  Appendix A Quantum mechanical calculation of $g^{(2)}$

We consider the so-called normalized second-order photon correlation function

$\displaystyle g^{(2)}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\frac{G^{(2)}(\bm{k}_{1},\bm{k}_{2})}{G^{(1)}(\bm{k}_{1},\bm{k}_{1})G^{(1)}(\bm{k}_{2},\bm{k}_{2})}=\frac{1}{N^{2}}\braket{\hat{E}^{(-)}(\bm{k}_{1})\hat{E}^{(-)}(\bm{k}_{2})\hat{E}^{(+)}(\bm{k}_{2})\hat{E}^{(+)}(\bm{k}_{1})}$
$\displaystyle=\frac{1}{N^{2}}\sum_{i,j,k,l=1}^{N}e^{i\bm{k}_{1}\bm{R}_{i}}e^{i\bm{k}_{2}\bm{R}_{j}}e^{-i\bm{k}_{2}\bm{R}_{k}}e^{-i\bm{k}_{1}\bm{R}_{l}}\braket{\hat{S}_{+}^{(i)}\hat{S}_{+}^{(j)}\hat{S}_{-}^{(k)}\hat{S}_{-}^{(l)}}\,.$
(53)

Since the atoms are assumed to be initially in the excited state and emit
incoherently, we only get a contribution if we have pairs of raising and
lowering operators with the same index. In addition, since the atoms can emit
at most one photon per pulsed excitation cycle, the case of all indices being
equal gives zero. Therefore, we can write

$\displaystyle\braket{\hat{S}_{+}^{(i)}\hat{S}_{+}^{(j)}\hat{S}_{-}^{(k)}\hat{S}_{-}^{(l)}}$
$\displaystyle=(1-\delta_{ij})(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\,.$
(54)

For the following calculations, it is convenient to describe the contributing
terms by permutations of the corresponding symmetric permutation groups. With
respect to the phase factors in Eq. (53) we can identify the index $i$ with
”1” (index of $\bm{k}_{1}$), $j$ with ”2”, $k$ with ”2”, and $l$ with ”1”.
Therefore, the Kronecker delta decomposition in Eq. (54) can be interpreted as
the two permutation group elements of the symmetric permutation group $S_{2}$,
where the Kronecker deltas correspond to the columns in Cauchy’s two-line
notation of a permutation. Explicitly, $\delta_{ik}$ corresponds to
$1\rightarrow 2$, $\delta_{jl}$ to $2\rightarrow 1$, $\delta_{il}$ to
$1\rightarrow 1$, and $\delta_{jk}$ to $2\rightarrow 2$. In Cauchy’s two-line
notation and cycle notation we can write

$\displaystyle\delta_{ik}\delta_{jl}\;$
$\displaystyle\hat{=}\left(\begin{array}[]{cc}1&2\\\
2&1\end{array}\right)=(1\;2)\,,$ (57) $\displaystyle\delta_{il}\delta_{jk}\;$
$\displaystyle\hat{=}\left(\begin{array}[]{cc}1&2\\\
1&2\end{array}\right)=()\,,$ (60)

where the first line denotes the transposition of the two elements, whereas
the second line is the identity. This notation will be especially useful when
we consider the third- and fourth-order photon correlation functions.
As mentioned earlier, the case, in which all indices are equal is zero. In Eq.
(54), this is accounted for by including the factor $(1-\delta_{ij})$. Thus,
we can write the second-order correlation function as

$\displaystyle g^{(2)}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\frac{1}{N^{2}}\sum_{i,j=1}^{N}(1-\delta_{ij})\left(1+e^{-i(\bm{k}_{2}-\bm{k}_{1})\bm{R}_{i}}e^{i(\bm{k}_{2}-\bm{k}_{1})\bm{R}_{j}}\right)\,,$
(61)

where we made use of Eq. (54). Simplifying $g^{(2)}(\bm{k}_{1},\bm{k}_{2})$
leads us to

$\displaystyle g^{(2)}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=\frac{1}{N^{2}}\left[N^{2}-2N+\sum_{i,j=1}^{N}e^{-i\bm{q}\bm{R}_{i}}e^{i\bm{q}\bm{R}_{j}}\right]$
$\displaystyle=1-\frac{2}{N}+\frac{1}{N^{2}}|\tilde{S}(\bm{q})|^{2}\,.$ (62)

Note that for $\bm{k}_{1}=\bm{k}_{2}=\bm{k}$, i.e., a difference vector
$\bm{q}=\bm{0}$, the correlation function reduces to

$g^{(2)}(\bm{k},\bm{k})=1-\frac{2}{N}+1=2\left(1-\frac{1}{N}\right)\,.$ (63)

Therefore, measuring the second-order photon correlation function for
$\bm{q}=\bm{0}$ gives access to the number of atoms $N$. In addition, we find
that measuring the second-order photon correlation function of the
fluorescence light, i.e., utilizing incoherent diffractive imaging, gives
access to the magnitude of the CDC.

##  Appendix B Quantum mechanical calculation of $g^{(3)}$

Next, we calculate the normalized third-order photon correlation function for
$N$ initially excited atoms, defined as

$\displaystyle g^{(3)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=\frac{G^{(3)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})}{G^{(1)}(\bm{k}_{1},\bm{k}_{1})G^{(1)}(\bm{k}_{2},\bm{k}_{2})G^{(1)}(\bm{k}_{3},\bm{k}_{3})}=\frac{1}{N^{3}}\braket{\hat{E}^{(-)}(\bm{k}_{1})\hat{E}^{(-)}(\bm{k}_{2})\hat{E}^{(-)}(\bm{k}_{3})\hat{E}^{(+)}(\bm{k}_{3})\hat{E}^{(+)}(\bm{k}_{2})\hat{E}^{(+)}(\bm{k}_{1})}$
$\displaystyle=\frac{1}{N^{3}}\sum_{i,j,k,l,m,n=1}^{N}e^{i\bm{k}_{1}\bm{R}_{i}}e^{i\bm{k}_{2}\bm{R}_{j}}e^{i\bm{k}_{3}\bm{R}_{k}}e^{-i\bm{k}_{3}\bm{R}_{l}}e^{-i\bm{k}_{2}\bm{R}_{m}}e^{-i\bm{k}_{1}\bm{R}_{n}}\braket{\hat{S}_{+}^{(i)}\hat{S}_{+}^{(j)}\hat{S}_{+}^{(k)}\hat{S}_{-}^{(l)}\hat{S}_{-}^{(m)}\hat{S}_{-}^{(n)}}\,.$
(64)

As before, we only get a contribution if we have pairs of raising and lowering
operators, which have the same index, whereas the expectation value vanishes
in the case that two pairs have the same index. In particular, we can write

$\braket{\hat{S}_{+}^{(i)}\hat{S}_{+}^{(j)}\hat{S}_{+}^{(k)}\hat{S}_{-}^{(l)}\hat{S}_{-}^{(m)}\hat{S}_{-}^{(n)}}=\delta_{il}\delta_{jm}\delta_{kn}+\delta_{il}\delta_{jn}\delta_{km}+\delta_{im}\delta_{jl}\delta_{kn}+\delta_{im}\delta_{jn}\delta_{kl}+\delta_{in}\delta_{jl}\delta_{km}+\delta_{in}\delta_{jm}\delta_{kl}$
(65)

for $i\neq j,i\neq k,j\neq k$ pairwise different. We can identify the
different delta contributions with permutation elements of the symmetric
permutation group $S_{3}$. In Cauchy’s two-line notation and in cycle notation
we find

$\displaystyle\delta_{il}\delta_{jm}\delta_{kn}\;$
$\displaystyle\hat{=}\left(\begin{array}[]{ccc}1&2&3\\\
3&2&1\end{array}\right)=(1\;3)\,,$ (68)
$\displaystyle\delta_{il}\delta_{jn}\delta_{km}\;$
$\displaystyle\hat{=}\left(\begin{array}[]{ccc}1&2&3\\\
3&1&2\end{array}\right)=(1\;3\;2)\,,$ (71)
$\displaystyle\delta_{im}\delta_{jl}\delta_{kn}\;$
$\displaystyle\hat{=}\left(\begin{array}[]{ccc}1&2&3\\\
2&3&1\end{array}\right)=(1\;2\;3)\,,$ (74)
$\displaystyle\delta_{im}\delta_{jn}\delta_{kl}\;$
$\displaystyle\hat{=}\left(\begin{array}[]{ccc}1&2&3\\\
2&1&3\end{array}\right)=(1\;2)\,,$ (77)
$\displaystyle\delta_{in}\delta_{jl}\delta_{km}\;$
$\displaystyle\hat{=}\left(\begin{array}[]{ccc}1&2&3\\\
1&3&2\end{array}\right)=(2\;3)\,,$ (80)
$\displaystyle\delta_{in}\delta_{jm}\delta_{kl}\;$
$\displaystyle\hat{=}\left(\begin{array}[]{ccc}1&2&3\\\
1&2&3\end{array}\right)=()\,.$ (83)

Therefore, we only have to calculate the sum for three of the summands
corresponding to the identity, one of the transpositions, and one of the
3-cycles. The remaining three summands are then easily deduced by index
replacements. We note again that we have to restrict the summation to all
$i,j,k$ being mutually different, which is accounted for by multiplication
with the factor $(1-\delta_{ij})(1-\delta_{ik})(1-\delta_{jk})$. Note that we
can write this product also as

$(1-\delta_{ij})(1-\delta_{ik})(1-\delta_{jk})=1+\delta_{ij}\delta_{ik}+\delta_{ij}\delta_{jk}-\delta_{ij}-\delta_{ik}-\delta_{jk}\,.$
(84)

We first consider the identity term

$\sum_{i,j,k=1}^{N}(1+\delta_{ij}\delta_{ik}+\delta_{ij}\delta_{jk}-\delta_{ij}-\delta_{ik}-\delta_{jk})=N^{3}+2N-3N^{2}\,.$
(85)

Second, we choose the transposition term $(1\;2)$, which gives

$\displaystyle\sum_{i,j,k=1}^{N}(1+\delta_{ij}\delta_{ik}+\delta_{ij}\delta_{jk}-\delta_{ij}-\delta_{ik}-\delta_{jk})e^{-i\bm{q}_{1}\bm{R}_{i}}e^{i\bm{q}_{1}\bm{R}_{j}}$
$\displaystyle=2N-N^{2}+(N-2)|\tilde{S}(\bm{q}_{1})|^{2}\,,$ (86)

where we defined $\bm{q}_{1}=\bm{k}_{2}-\bm{k}_{1}$. Last, we need one of the
3-cycle terms. We choose $(1\;3\;2)$, which gives

$\displaystyle\sum_{i,j,k=1}^{N}(1+\delta_{ij}\delta_{ik}+\delta_{ij}\delta_{jk}-\delta_{ij}-\delta_{ik}-\delta_{jk})e^{-i\bm{q}_{3}\bm{R}_{i}}e^{i\bm{q}_{2}\bm{R}_{k}}e^{i\bm{q}_{1}\bm{R}_{j}}$
$\displaystyle=\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}_{r}(\bm{q}_{3})-\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2})\tilde{S}_{r}(\bm{q}_{3})-\tilde{S}^{*}(\bm{q}_{2})\tilde{S}_{r}(\bm{q}_{3}-\bm{q}_{1})-\tilde{S}^{*}(\bm{q}_{1})\tilde{S}_{r}(\bm{q}_{3}-\bm{q}_{2})+2\tilde{S}_{r}(\bm{q}_{3}-\bm{q}_{1}-\bm{q}_{2})$
$\displaystyle=2N-|\tilde{S}(\bm{q}_{1})|^{2}-|\tilde{S}(\bm{q}_{2})|^{2}-|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}+\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}_{r}(\bm{q}_{1}+\bm{q}_{2})\,,$
(87)

where $\bm{q}_{2}=\bm{k}_{3}-\bm{k}_{2}$ and
$\bm{q}_{3}=\bm{k}_{3}-\bm{k}_{1}=\bm{k}_{3}-\bm{k}_{2}+\bm{k}_{2}-\bm{k}_{1}=\bm{q}_{1}+\bm{q}_{2}$.
Accounting for all six permutations ${\color[rgb]{0,0,0}()}$,
${\color[rgb]{0,0,0}(1\;2)}$, ${\color[rgb]{0,0,0}(2\;3)}$,
${\color[rgb]{0,0,0}(1\;3)}$, ${\color[rgb]{0,0,0}(1\;3\;2)}$, and
${\color[rgb]{0,0,0}(1\;2\;3)}$, the third-order photon correlation function
calculates to

$\displaystyle g^{(3)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=\frac{1}{N^{3}}\left[{\color[rgb]{0,0,0}N^{3}+2N-3N^{2}}+{\color[rgb]{0,0,0}2N-N^{2}+(N-2)|\tilde{S}(\bm{q}_{1})|^{2}}\right.$
$\displaystyle\left.+{\color[rgb]{0,0,0}2N-N^{2}+(N-2)|\tilde{S}(\bm{q}_{2})|^{2}}+{\color[rgb]{0,0,0}2N-N^{2}+(N-2)|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}}\right.$
$\displaystyle\left.+{\color[rgb]{0,0,0}2N-|\tilde{S}(\bm{q}_{1})|^{2}-|\tilde{S}(\bm{q}_{2})|^{2}-|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}+\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}_{r}(\bm{q}_{1}+\bm{q}_{2})}\right.$
$\displaystyle\left.+{\color[rgb]{0,0,0}2N-|\tilde{S}(\bm{q}_{1})|^{2}-|\tilde{S}(\bm{q}_{2})|^{2}-|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}+\tilde{S}_{r}(\bm{q}_{1})\tilde{S}_{r}(\bm{q}_{2})\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2})}\right]$
$\displaystyle=1-\frac{6}{N}+\frac{12}{N^{2}}+\frac{N-4}{N^{3}}[|\tilde{S}(\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}]$
$\displaystyle+\frac{2}{N^{3}}|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|\cos[\phi(\bm{q}_{1}+\bm{q}_{2})-\phi(\bm{q}_{1})-\phi(\bm{q}_{2})]\,,$
(88)

where we have written the structure factor in terms of its magnitude and phase
as

$\tilde{S}(\bm{q})=|\tilde{S}(\bm{q})|e^{i\phi(\bm{q})}\,,$ (89)

and also used

$|\tilde{S}(-\bm{q})|=|\tilde{S}(\bm{q})|\,,\quad\phi(-\bm{q})=-\phi(\bm{q})\,,$
(90)

which holds for an arbitrary real-valued atomic distribution $S(\bm{r})$.

##  Appendix C Quantum mechanical calculation of $g^{(4)}$

We consider the normalized fourth-order photon correlation function defined as

$\displaystyle
g^{(4)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})=\frac{G^{(4)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})}{G^{(1)}(\bm{k}_{1},\bm{k}_{1})G^{(1)}(\bm{k}_{2},\bm{k}_{2})G^{(1)}(\bm{k}_{3},\bm{k}_{3})G^{(1)}(\bm{k}_{4},\bm{k}_{4})}$
$\displaystyle=\frac{1}{N^{4}}\braket{\hat{E}^{(-)}(\bm{k}_{1})\hat{E}^{(-)}(\bm{k}_{2})\hat{E}^{(-)}(\bm{k}_{3})\hat{E}^{(-)}(\bm{k}_{4})\hat{E}^{(+)}(\bm{k}_{4})\hat{E}^{(+)}(\bm{k}_{3})\hat{E}^{(+)}(\bm{k}_{2})\hat{E}^{(+)}(\bm{k}_{1})}$
$\displaystyle=\frac{1}{N^{4}}\sum_{i,j,k,l,m,n,o,p=1}^{N}e^{i\bm{k}_{1}\bm{R}_{i}}e^{i\bm{k}_{2}\bm{R}_{j}}e^{i\bm{k}_{3}\bm{R}_{k}}e^{i\bm{k}_{4}\bm{R}_{l}}e^{-i\bm{k}_{4}\bm{R}_{m}}e^{-i\bm{k}_{3}\bm{R}_{n}}e^{-i\bm{k}_{2}\bm{R}_{o}}e^{-i\bm{k}_{1}\bm{R}_{p}}\braket{\hat{S}_{+}^{(i)}\hat{S}_{+}^{(j)}\hat{S}_{+}^{(k)}\hat{S}_{+}^{(l)}\hat{S}_{-}^{(m)}\hat{S}_{-}^{(n)}\hat{S}_{-}^{(o)}\hat{S}_{-}^{(p)}}\,.$
(91)

Again, the expectation value of the normally ordered product of raising and
lowering operators is only non-zero if we have pairs of raising and lowering
operators with the same index. This gives $4!=24$ different contributions,
which correspond to the $24$ permutation elements of the symmetric permutation
group $S_{4}$. The group consists of one identity permutation, six
transpositions, three double transpositions, eight 3-cycles and six 4-cycles,
which are listed in Tab. 1 in cycle notation.

permutation type | elements in cycle notation
---|---
identity | $()$
transposition |  $(1\;2)$, $(1\;3)$, $(1\;4)$, $(2\;3)$, $(2\;4)$, $(3\;4)$
double transposition |  $(1\;2)(3\;4)$, $(1\;3)(2\;4)$, $(1\;4)(2\;3)$
3-cycle |  $(1\;2\;3)$, $(1\;3\;2)$, $(2\;3\;4)$, $(2\;4\;3)$, $(3\;4\;1)$, $(3\;1\;4)$, $(4\;1\;2)$, $(4\;2\;1)$
4-cycle |  $(1\;2\;3\;4)$, $(1\;2\;4\;3)$, $(1\;3\;2\;4)$, $(1\;3\;4\;2)$, $(1\;4\;2\;3)$, $(1\;4\;3\;2)$
Table 1: Elements of the symmetric permutation group $S_{4}$.

Therefore, we find that we only have to calculate five different terms
corresponding to the five different permutation types. In addition, the
normally ordered expectation value gives zero if two creation and annihilation
operator pairs have the same index. Thus, the remaining indices $i,j,k,l$ need
to be mutually different. This is accounted for by multiplication with

$\displaystyle\text{del}(i,j,k,l)\coloneqq(1-\delta_{ij})(1-\delta_{ik})(1-\delta_{il})(1-\delta_{jk})(1-\delta_{jl})(1-\delta_{kl})$
$\displaystyle=1-\delta_{ij}\delta_{ik}\delta_{il}-\delta_{ij}\delta_{il}\delta_{jk}-\delta_{ij}\delta_{ik}\delta_{jl}-\delta_{ij}\delta_{jk}\delta_{jl}-\delta_{ij}\delta_{ik}\delta_{kl}-\delta_{ij}\delta_{jk}\delta_{kl}+\delta_{ij}\delta_{kl}+\delta_{il}\delta_{jk}+\delta_{ik}\delta_{jl}$
$\displaystyle+\delta_{ij}\delta_{ik}+\delta_{ij}\delta_{jk}+\delta_{ij}\delta_{il}+\delta_{ij}\delta_{jl}+\delta_{ik}\delta_{il}+\delta_{ik}\delta_{kl}+\delta_{jk}\delta_{jl}+\delta_{jk}\delta_{kl}-\delta_{ij}-\delta_{ik}-\delta_{il}-\delta_{jk}-\delta_{jl}-\delta_{kl}\,.$
(92)

Now, we evaluate the five different contribution types.
The identity type gives

$\displaystyle{\color[rgb]{0,0,0}p_{1}}=\sum_{i,j,k,l=1}^{N}\text{del}(i,j,k,l)=N^{4}-6N^{3}+11N^{2}-6N\,.$
(93)

The transposition type yields

$\displaystyle{\color[rgb]{0,0,0}p_{2}(\bm{q})}$
$\displaystyle=\sum_{i,j,k,l=1}^{N}\text{del}(i,j,k,l)e^{-i\bm{q}(\bm{R}_{k}-\bm{R}_{l})}=N^{2}\tilde{S}^{*}(\bm{q})\tilde{S}(\bm{q})-N^{3}+5N^{2}-5N\tilde{S}^{*}(\bm{q})\tilde{S}(\bm{q})-6N+6\tilde{S}^{*}(\bm{q})\tilde{S}(\bm{q})$
$\displaystyle=-N^{3}+5N^{2}-6N+(N^{2}-5N+6)|\tilde{S}(\bm{q})|^{2}\,.$ (94)

The double transposition type is given by

$\displaystyle{\color[rgb]{0,0,0}p_{3}(\bm{q}_{1},\bm{q}_{2})}$
$\displaystyle=\sum_{i,j,k,l=1}^{N}\text{del}(i,j,k,l)e^{-i\bm{q}_{1}(\bm{R}_{i}-\bm{R}_{j})-i\bm{q}_{2}(\bm{R}_{k}-\bm{R}_{l})}$
$\displaystyle=N^{2}-N\tilde{S}^{*}(\bm{q}_{1})\tilde{S}(\bm{q}_{1})-N\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{2})-6N+\tilde{S}^{*}(\bm{q}_{1})\tilde{S}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{2})-\tilde{S}^{*}(\bm{q}_{1})\tilde{S}(\bm{q}_{2})\tilde{S}^{*}(\bm{q}_{2}-\bm{q}_{1})$
$\displaystyle-\tilde{S}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{2}-\bm{q}_{1})-\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{1}+\bm{q}_{2})-\tilde{S}(\bm{q}_{1})\tilde{S}(\bm{q}_{2})\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2})+\tilde{S}^{*}(\bm{q}_{2}-\bm{q}_{1})\tilde{S}(\bm{q}_{2}-\bm{q}_{1})$
$\displaystyle+\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2})\tilde{S}(\bm{q}_{1}+\bm{q}_{2})+4\tilde{S}^{*}(\bm{q}_{1})\tilde{S}(\bm{q}_{1})+4\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{2})$
$\displaystyle=N^{2}-6N-(N-4)[|\tilde{S}(\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}]+|\tilde{S}(\bm{q}_{2}-\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}$
$\displaystyle+|\tilde{S}(\bm{q}_{1})|^{2}|\tilde{S}(\bm{q}_{2})|^{2}-2|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{2}-\bm{q}_{1})|\cos[\phi(\bm{q}_{1})-\phi(\bm{q}_{2})+\phi(\bm{q}_{2}-\bm{q}_{1})]$
$\displaystyle-2|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{2})-\phi(\bm{q}_{1}+\bm{q}_{2})]\,.$
(95)

The 3-cycle type reads

$\displaystyle{\color[rgb]{0,0,0}p_{4}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3})}$
$\displaystyle=\sum_{i,j,k,l=1}^{N}\text{del}(i,j,k,l)e^{-i\bm{q}_{3}\bm{R}_{j}+i\bm{q}_{1}\bm{R}_{k}+i\bm{q}_{2}\bm{R}_{l}}$
$\displaystyle=N\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{3})-N\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{3}-\bm{q}_{1})-N\tilde{S}^{*}(\bm{q}_{1})\tilde{S}(\bm{q}_{3}-\bm{q}_{2})-N\tilde{S}(\bm{q}_{3})\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2})$
$\displaystyle+2N\tilde{S}(\bm{q}_{3}-\bm{q}_{1}-\bm{q}_{2})-3\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{3})+3\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{3}-\bm{q}_{1})+3\tilde{S}^{*}(\bm{q}_{1})\tilde{S}(\bm{q}_{3}-\bm{q}_{2})$
$\displaystyle+3\tilde{S}(\bm{q}_{3})\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2})-6\tilde{S}(\bm{q}_{3}-\bm{q}_{1}-\bm{q}_{2})$
$\displaystyle=(N-3)[\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{3})-\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{3}-\bm{q}_{1})-\tilde{S}^{*}(\bm{q}_{1})\tilde{S}(\bm{q}_{3}-\bm{q}_{2})$
$\displaystyle-\tilde{S}(\bm{q}_{3})\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2})+2\tilde{S}(\bm{q}_{3}-\bm{q}_{1}-\bm{q}_{2})]\,.$
(96)

The 4-cycle type calculates to

$\displaystyle{\color[rgb]{0,0,0}p_{5}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3},\bm{q}_{4})}$
$\displaystyle=\sum_{i,j,k,l=1}^{N}\text{del}(i,j,k,l)e^{-i\bm{q}_{4}\bm{R}_{i}+i\bm{q}_{1}\bm{R}_{j}+i\bm{q}_{2}\bm{R}_{k}+i\bm{q}_{3}\bm{R}_{l}}$
$\displaystyle=\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}^{*}(\bm{q}_{3})\tilde{S}(\bm{q}_{4})-\tilde{S}^{*}(\bm{q}_{3})\tilde{S}(\bm{q}_{4})\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2})-\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{4})\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{3})$
$\displaystyle-\tilde{S}^{*}(\bm{q}_{1})\tilde{S}(\bm{q}_{4})\tilde{S}^{*}(\bm{q}_{2}+\bm{q}_{3})-\tilde{S}^{*}(\bm{q}_{2})\tilde{S}^{*}(\bm{q}_{3})\tilde{S}(\bm{q}_{4}-\bm{q}_{1})-\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{3})\tilde{S}(\bm{q}_{4}-\bm{q}_{2})$
$\displaystyle-\tilde{S}^{*}(\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{4}-\bm{q}_{3})+2\tilde{S}^{*}(\bm{q}_{3})\tilde{S}(\bm{q}_{4}-\bm{q}_{1}-\bm{q}_{2})+2\tilde{S}^{*}(\bm{q}_{2})\tilde{S}(\bm{q}_{4}-\bm{q}_{1}-\bm{q}_{3})$
$\displaystyle+2\tilde{S}^{*}(\bm{q}_{1})\tilde{S}(\bm{q}_{4}-\bm{q}_{2}-\bm{q}_{3})+2\tilde{S}(\bm{q}_{4})\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})+\tilde{S}(\bm{q}_{4}-\bm{q}_{1})\tilde{S}^{*}(\bm{q}_{2}+\bm{q}_{3})$
$\displaystyle+\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{3})\tilde{S}(\bm{q}_{4}-\bm{q}_{2})+\tilde{S}^{*}(\bm{q}_{1}+\bm{q}_{2})\tilde{S}(\bm{q}_{4}-\bm{q}_{3})-6\tilde{S}(\bm{q}_{4}-\bm{q}_{1}-\bm{q}_{2}-\bm{q}_{3})\,.$
(97)

Having the five different contribution types evaluated, we define the
difference vectors $\bm{q}_{1}=\bm{k}_{2}-\bm{k}_{1}$,
$\bm{q}_{2}=\bm{k}_{3}-\bm{k}_{2}$, $\bm{q}_{3}=\bm{k}_{4}-\bm{k}_{3}$,
$\bm{q}_{4}=\bm{k}_{3}-\bm{k}_{1}=\bm{q}_{1}+\bm{q}_{2}$,
$\bm{q}_{5}=\bm{k}_{4}-\bm{k}_{1}=\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3}$,
$\bm{q}_{6}=\bm{k}_{4}-\bm{k}_{2}=\bm{q}_{2}+\bm{q}_{3}$. Now, the normalized
fourth-order photon correlation function is simply the sum over the $24$ terms
corresponding to the different permutations of $S_{4}$, which can be
calculated via the functions $p_{1}$ to $p_{5}$ by properly plugging in the
correct difference vectors for the different permutations, i.e.,

$\displaystyle
N^{4}g^{(4)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})={\color[rgb]{0,0,0}p_{1}}{\color[rgb]{0,0,0}+p_{2}(\bm{q}_{1})+p_{2}(\bm{q}_{2})+p_{2}(\bm{q}_{3})+p_{2}(\bm{q}_{4})+p_{2}(\bm{q}_{5})+p_{2}(\bm{q}_{6})}$
$\displaystyle{\color[rgb]{0,0,0}+p_{3}(\bm{q}_{1},\bm{q}_{3})+p_{3}(\bm{q}_{2},\bm{q}_{5})+p_{3}(\bm{q}_{4},\bm{q}_{6})}{\color[rgb]{0,0,0}+p_{4}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{4})+p_{4}(-\bm{q}_{2},\bm{q}_{4},\bm{q}_{1})+p_{4}(\bm{q}_{1},\bm{q}_{6},\bm{q}_{5})}$
$\displaystyle{\color[rgb]{0,0,0}+p_{4}(-\bm{q}_{6},\bm{q}_{5},\bm{q}_{1})+p_{4}(\bm{q}_{2},\bm{q}_{3},\bm{q}_{6})+p_{4}(-\bm{q}_{3},\bm{q}_{6},\bm{q}_{2})+p_{4}(-\bm{q}_{3},\bm{q}_{5},\bm{q}_{4})+p_{4}(\bm{q}_{4},\bm{q}_{3},\bm{q}_{5})}$
$\displaystyle{\color[rgb]{0,0,0}+p_{5}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3},\bm{q}_{5})+p_{5}(-\bm{q}_{2},-\bm{q}_{3},\bm{q}_{5},\bm{q}_{1})+p_{5}(\bm{q}_{1},-\bm{q}_{3},\bm{q}_{6},\bm{q}_{4})}$
$\displaystyle{\color[rgb]{0,0,0}+p_{5}(-\bm{q}_{6},\bm{q}_{4},\bm{q}_{3},\bm{q}_{1})+p_{5}(-\bm{q}_{2},\bm{q}_{4},\bm{q}_{6},\bm{q}_{5})+p_{5}(-\bm{q}_{6},\bm{q}_{2},\bm{q}_{5},\bm{q}_{4})}\,.$
(98)

Thereby, the individual contributions are given by

$\displaystyle{\color[rgb]{0,0,0}p_{1}}=N^{4}-6N^{3}+11N^{2}-6N\,,$ (99)
$\displaystyle{\color[rgb]{0,0,0}p_{2}(\bm{q}_{1})+p_{2}(\bm{q}_{2})+p_{2}(\bm{q}_{3})+p_{2}(\bm{q}_{4})+p_{2}(\bm{q}_{5})+p_{2}(\bm{q}_{6})}=$
$\displaystyle-6N^{3}+30N^{2}-36N+(N^{2}-5N+6)[|\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|^{2}$
$\displaystyle+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{3})|^{2}]\,,$
(100)
$\displaystyle{\color[rgb]{0,0,0}p_{3}(\bm{q}_{1},\bm{q}_{3})+p_{3}(\bm{q}_{2},\bm{q}_{5})+p_{3}(\bm{q}_{4},\bm{q}_{6})}=3N^{2}-18N$
$\displaystyle+|\tilde{S}(\bm{q}_{1})|^{2}|\tilde{S}(\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}|\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}|\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|^{2}$
$\displaystyle-(N-4)[|\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{3})|^{2}]$
$\displaystyle+2[|\tilde{S}(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{3}-\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{3})|^{2}]$
$\displaystyle-2\\{|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{3}-\bm{q}_{1})|\cos[\phi(\bm{q}_{1})-\phi(\bm{q}_{3})+\phi(\bm{q}_{3}-\bm{q}_{1})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{3}-\bm{q}_{1})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1}+\bm{q}_{2})+\phi(\bm{q}_{3}-\bm{q}_{1})-\phi(\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{2})+\phi(\bm{q}_{1}+\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1}+\bm{q}_{2})+\phi(\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{2})+\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})]\\}\,,$
(101)
$\displaystyle{\color[rgb]{0,0,0}p_{4}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{4})+p_{4}(-\bm{q}_{2},\bm{q}_{4},\bm{q}_{1})+p_{4}(\bm{q}_{1},\bm{q}_{6},\bm{q}_{5})+p_{4}(-\bm{q}_{6},\bm{q}_{5},\bm{q}_{1})+p_{4}(\bm{q}_{2},\bm{q}_{3},\bm{q}_{6})}$
$\displaystyle{\color[rgb]{0,0,0}+p_{4}(-\bm{q}_{3},\bm{q}_{6},\bm{q}_{2})+p_{4}(-\bm{q}_{3},\bm{q}_{5},\bm{q}_{4})+p_{4}(\bm{q}_{4},\bm{q}_{3},\bm{q}_{5})}=16N^{2}-48N$
$\displaystyle-4(N-3)[|\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{3})|^{2}]$
$\displaystyle+2(N-3)\\{|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{2})+\phi(\bm{q}_{3})-\phi(\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{2})-\phi(\bm{q}_{1}+\bm{q}_{2})]$
$\displaystyle+|\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{3})+\phi(\bm{q}_{1}+\bm{q}_{2})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]\\}\,,$
(102)
$\displaystyle{\color[rgb]{0,0,0}p_{5}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3},\bm{q}_{5})+p_{5}(-\bm{q}_{2},-\bm{q}_{3},\bm{q}_{5},\bm{q}_{1})+p_{5}(\bm{q}_{1},-\bm{q}_{3},\bm{q}_{6},\bm{q}_{4})}$
$\displaystyle{\color[rgb]{0,0,0}+p_{5}(-\bm{q}_{6},\bm{q}_{4},\bm{q}_{3},\bm{q}_{1})+p_{5}(-\bm{q}_{2},\bm{q}_{4},\bm{q}_{6},\bm{q}_{5})+p_{5}(-\bm{q}_{6},\bm{q}_{2},\bm{q}_{5},\bm{q}_{4})}=$
$\displaystyle-36N+10[|\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{3})|^{2}]$
$\displaystyle+2[|\tilde{S}(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{3}-\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{3})|^{2}]$
$\displaystyle-6\\{|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{2})+\phi(\bm{q}_{3})-\phi(\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{3})+\phi(\bm{q}_{1}+\bm{q}_{2})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{2})-\phi(\bm{q}_{1}+\bm{q}_{2})]\\}$
$\displaystyle-2\\{|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{3}-\bm{q}_{1})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1}+\bm{q}_{2})+\phi(\bm{q}_{3}-\bm{q}_{1})-\phi(\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{2})+\phi(\bm{q}_{1}+\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1}+\bm{q}_{2})+\phi(\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{2})+\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{3}-\bm{q}_{1})|\cos[\phi(\bm{q}_{1})-\phi(\bm{q}_{3})+\phi(\bm{q}_{3}-\bm{q}_{1})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{3})]\\}$
$\displaystyle+2\\{|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})-\phi(\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2})+\phi(\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{2})+\phi(\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|$
$\displaystyle\times\cos[\phi(\bm{q}_{2})-\phi(\bm{q}_{1}+\bm{q}_{2})-\phi(\bm{q}_{2}+\bm{q}_{3})+\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]\\}\,.$
(103)

In total, we have

$\displaystyle
N^{4}g^{(4)}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3})=N^{4}-12N^{3}+60N^{2}-144N+(N^{2}-10N+32)\times$
$\displaystyle\times[|\tilde{S}(\bm{q}_{1})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}+|\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|^{2}]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})|^{2}|\tilde{S}(\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|^{2}|\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{2})|^{2}|\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|^{2}$
$\displaystyle+4[|\tilde{S}(\bm{q}_{1}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})|^{2}+|\tilde{S}(\bm{q}_{3}-\bm{q}_{1})|^{2}]$
$\displaystyle-4\\{|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{3}-\bm{q}_{1})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1}+\bm{q}_{2})+\phi(\bm{q}_{3}-\bm{q}_{1})-\phi(\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{2})+\phi(\bm{q}_{1}+\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1}+\bm{q}_{2})+\phi(\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{2})+\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+2\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{3}-\bm{q}_{1})|\cos[\phi(\bm{q}_{1})-\phi(\bm{q}_{3})+\phi(\bm{q}_{3}-\bm{q}_{1})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{3})]\\}$
$\displaystyle+2(N-6)\\{|\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{3})+\phi(\bm{q}_{1}+\bm{q}_{2})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{2})-\phi(\bm{q}_{1}+\bm{q}_{2})]$
$\displaystyle+|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{2})+\phi(\bm{q}_{3})-\phi(\bm{q}_{2}+\bm{q}_{3})]\\}$
$\displaystyle+2\\{|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})-\phi(\bm{q}_{3})+\phi(\bm{q}_{2}+\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2})]$
$\displaystyle+|\tilde{S}(\bm{q}_{1})||\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\cos[\phi(\bm{q}_{1})+\phi(\bm{q}_{2})+\phi(\bm{q}_{3})-\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]$
$\displaystyle+|\tilde{S}(\bm{q}_{2})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2})||\tilde{S}(\bm{q}_{2}+\bm{q}_{3})||\tilde{S}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})|\times$
$\displaystyle\times\cos[\phi(\bm{q}_{2})-\phi(\bm{q}_{1}+\bm{q}_{2})-\phi(\bm{q}_{2}+\bm{q}_{3})+\phi(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})]\\}\,,$
(104)

which can be used as an additional equation to check the validity of the
phases of the structure factor reconstructed from the normalized third-order
photon correlation function $g^{(3)}(\bm{q}_{1},\bm{q}_{2})$.

## References

  * Fienup [1978] J. R. Fienup, Reconstruction of an object from the modulus of its Fourier transform, Optics Letters 3, 27 (1978)
  * Fienup [1982] J. R. Fienup, Phase retrieval algorithms: a comparison, Applied Optics 21, 2758 (1982)
  * Fienup and Kowalczyk [1990] J. R. Fienup and A. M. Kowalczyk, Phase retrieval for a complex-valued object by using a low-resolution image, JOSA A 7, 450 (1990)
  * Belen’kii and Holmes [2004] M. S. Belen’kii and R. B. Holmes, Investigation of the Cauchy–Riemann equations for one-dimensional image recovery in intensity interferometry, JOSA A 21, 697 (2004)
  * Shechtman  _et al._ [2014] Y. Shechtman, A. Beck, and Y. C. Eldar, Gespar: Efficient phase retrieval of sparse signals, IEEE Transactions on Signal Processing 62, 928 (2014)
  * Strekalov  _et al._ [2014] D. V. Strekalov, I. Kulikov, and N. Yu, Imaging dark objects with intensity interferometry, Opt. Exp. 22, 12339 (2014)
  * Shechtman  _et al._ [2015] Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, Phase retrieval with application to optical imaging: A contemporary overview, IEEE Signal Processing Magazine 32, 87 (2015)
  * Cao  _et al._ [2016] S. Cao, P. Kok, P. Li, A. M. Maiden, and J. M. Rodenburg, Modal decomposition of a propagating matter wave via electron ptychography, Phys. Rev. A 94, 063621 (2016)
  * Wang  _et al._ [2018] W. Wang, H. Chen, Y. Yuan, Q. Han, G. Wang, H. Zheng, J. Liu, and Z. Xu, Ptychographical intensity interferometry imaging with incoherent light, Opt. Exp. 26, 20396 (2018)
  * Nakasako  _et al._ [2020] M. Nakasako, A. Kobayashi, Y. Takayama, K. Asakura, M. Oide, K. Okajima, T. Oroguchi, and M. Yamamoto, Methods and application of coherent x-ray diffraction imaging of noncrystalline particles, Biophysical Reviews 12, 541 (2020)
  * Classen  _et al._ [2017] A. Classen, K. Ayyer, H. N. Chapman, R. Röhlsberger, and J. von Zanthier, Incoherent diffractive imaging via intensity correlations of hard x rays, Phys. Rev. Lett. 119, 053401 (2017)
  * Schneider  _et al._ [2018] R. Schneider, T. Mehringer, G. Mercurio, L. Wenthaus, A. Classen, G. Brenner, O. Gorobtsov, A. Benz, D. Bhatti, L. Bocklage, _et al._ , Quantum imaging with incoherently scattered light from a free-electron laser, Nature Physics 14, 126 (2018)
  * Trost  _et al._ [2020] F. Trost, K. Ayyer, and H. N. Chapman, Photon statistics and signal to noise ratio for incoherent diffraction imaging, New Journal of Physics 22, 083070 (2020)
  * Gamo [1963] H. Gamo, Triple correlator of photoelectric fluctuations as a spectroscopic tool, Journal of Applied Physics 34, 875 (1963)
  * Lohmann  _et al._ [1983] A. W. Lohmann, G. Weigelt, and B. Wirnitzer, Speckle masking in astronomy: triple correlation theory and applications, Appl. Opt. 22, 4028 (1983)
  * Lohmann and Wirnitzer [1984] A. Lohmann and B. Wirnitzer, Triple correlations, Proceedings of the IEEE 72, 889 (1984)
  * Cornwell [1989] T. J. Cornwell, The applications of closure phase to astronomical imaging, Science 245, 263 (1989)
  * Kang  _et al._ [1991] M. G. Kang, K.-T. Lay, and A. K. Katsaggelos, Phase estimation using the bispectrum and its application to image restoration, Optical Engineering 30, 976 (1991)
  * Marathay  _et al._ [1994] A. S. Marathay, Y. Z. Hu, and L. Shao, Phase function of spatial coherence from second-, third-, and fourth order intensity correlations, Optical Engineering 33, 3265 (1994)
  * Ofir and Ribak [2006a] A. Ofir and E. N. Ribak, Offline, multidetector intensity interferometers – I. Theory, MNRAS 368, 1646 (2006a)
  * Ofir and Ribak [2006b] A. Ofir and E. N. Ribak, Offline, multidetector intensity interferometers – II. Implications and applications, MNRAS 368, 1652 (2006b)
  * Malvimat  _et al._ [2013] V. Malvimat, O. Wucknitz, and P. Saha, Intensity interferometry with more than two detectors?, MNRAS 437, 798 (2013)
  * Dravins  _et al._ [2013] D. Dravins, S. LeBohec, H. Jensen, and P. D. Nuñez, Optical intensity interferometry with the cherenkov telescope array, Astroparticle Physics 43, 331 (2013)
  * Nuñez and Domiciano de Souza [2015] P. D. Nuñez and A. Domiciano de Souza, Capabilities of future intensity interferometers for observing fast-rotating stars: imaging with two- and three-telescope correlations, MNRAS 453, 1999 (2015)
  * Shoulga and Ribak [2017] G. Shoulga and E. N. Ribak, Toward spectral intensity interferometry, Appl. Opt. 56, A23 (2017)
  * Peard  _et al._ [2023] N. Peard, K. Ayyer, and H. N. Chapman, Ab initio spatial phase retrieval via intensity triple correlations, Opt. Express 31, 25082 (2023)
  * Agarwal [1974] G. S. Agarwal,  _Quantum statistical theories of spontaneous emission and their relation to other approaches - Quantum Optics_ , edited by G. Höhler (Springer-Verlag, 1974)
  * Glauber [1963] R. J. Glauber, The quantum theory of optical coherence, Phys. Rev. 130, 2529 (1963)