Introduction

The van Cittert-Zernike theorem constitutes one of the pillars of optical physics1,2. As such, this fundamental theorem provides the formalism to describe the modification of the coherence properties of optical fields upon propagation1,2,3,4. Over the last decades, extensive investigations have been conducted to explore the evolution of spatial, temporal, spectral, and polarization coherence of diverse families of optical beams5,6,7,8. In the context of classical optics, the investigation of the van Cittert-Zernike theorem led to the development of schemes for optical sensing, metrology, and astronomical interferometry9,10,11. Nowadays, there has been interest in exploring the implications of the van Cittert-Zernike theorem for quantum mechanical systems12,13,14,15,16,17. Recent efforts have been devoted to study the evolution of the properties of spatial coherence of biphoton systems12,13,15,17,18,19. Specifically, the van Cittert-Zernike theorem has been extended to analyze the spatial entanglement between a pair of photons generated by parametric down conversion12,13,20,21. The description of the evolution of spatial coherence and entanglement of propagating photons turned out essential for quantum metrology, spectroscopy, imaging, and lithography12,13,15,17,22,23,24,25,26. Nevertheless, previous research has not explored the evolution of the excitation mode of the field that establishes the quantum statistical properties of light fields12,13,14,15,16,17,18,19,20,21,22,23,24,25,26.

There has been important progress on the preparation of multiphoton systems with quantum mechanical properties27,28,29,30,31. The interest in these systems resides in the complex interference and scattering effects that they can host25,30,32,33. Remarkably, these fundamental processes define the statistical fluctuations of photons that establish the nature of light sources27,28,29,34,35,36. Furthermore, these quantum fluctuations are associated to distinct excitation modes of the electromagnetic field that determine the quantum coherence of a light field34,36. In the context of quantum information processing, the interference and scattering among photons have enormous potential to perform operations that are intractable on classical systems25,32. However, the manipulation of multiphoton systems requires complex light-matter interactions that are hard to achieve at the few-photon level28,37. Indeed, it has been assumed that light-matter interactions are needed to modify the excitation mode of an optical field38,39. These challenges have motivated interest in linear optical circuits for random walks, boson sampling, and quantum computing32,40,41. Moreover, the interaction of multiphoton quantum systems with the environment leads to the degradation of their nonclassical properties upon propagation33,42. Indeed, quantum-to-classical transitions are unavoidable in nonclassical systems interacting with realistic environments42. These vulnerabilities have prevented the use of nonclassical states of light for the sensing of small physical parameters with sensitivities that surpass the shot-noise limit24,43. The possibility of using nonclassical multiphoton states to demonstrate scalable quantum sensing has constituted one of the main goals of quantum optics for many decades43.

In contrast to well-established paradigms in the field of quantum optics12,13,14,15,16,17,18,19,20,21,22,23,24,25,26, our work demonstrates that the quantum statistical fluctuations of multiphoton light fields can be modified upon propagation in the absence of conventional light-matter interactions. These understood as optical processes taking place in the absence of photon absorption and emission44. We introduce the quantum van Cittert-Zernike theorem to describe the underlying scattering effects that give rise to photon statistics modifications in multiphoton systems. Remarkably, our work unveils the conditions under which sub-shot-noise quantum fluctuations could potentially be extracted from thermal light sources. These effects remained elusive for multiple decades due to the limitations of the classical theory of optical coherence to describe multiphoton scattering12,13,14,15,16,17,18,19,20,21,22,23,24,25,26. This stimulated the idea that the quantum statistical fluctuations of light fields were not affected upon propagation. Our work provides an all-optical alternative to prepare multiphoton systems with sub-Poissonian-like statistics. Previously, similar functionalities have been demonstrated in nonlinear optical systems, photonic lattices, plasmonic systems, and Bose-Einstein condensates27,29,38,45.

Results and discussion

Description of the quantum van Cittert-Zernike theorem

We demonstrate the quantum van-Cittert-Zernike theorem by analyzing the propagation of multiphoton two-mode correlations in the setup depicted in Fig. 1. In general, each mode can host a multiphoton system with an arbitrary number of photons. We consider a thermal, spatially incoherent, unpolarized beam that interacts with a polarization grating. This grating modifies the polarization of the thermal beam at different transverse spatial locations x according to πx/L. Here, L represents the length of the grating. The thermal beam propagates to the far-field, where it is measured by two point detectors46,47,48. We then post-select on the intensity measurements made by these detectors to quantify the correlations between different modes of the beam.

Fig. 1: The proposed setup for investigating the multiphoton quantum van Cittert-Zernike theorem.
figure 1

We consider an incoherent, unpolarized beam interacting with a polarization grating of length L at z = 0. After interacting with the grating, the beam propagates a distance of z onto the measurement plane, where two point detectors are placed ΔX apart.

Full size image

The multiphoton quantum van Cittert-Zernike theorem can be demonstrated for any incoherent, unpolarized state, the simplest of which is an unpolarized two-mode state49. The two-mode state can be produced by a source emitting a series of spatially independent photons with either horizontal (H) or vertical (V) polarization, giving an initial state27,50

$$\begin{array}{lll}\hat{\rho }={\hat{\rho }}_{1}\otimes {\hat{\rho }}_{2}\\\;\;\;=\frac{1}{4}\left({\left\vert H\right\rangle }_{1}{\left\vert H\right\rangle }_{2}{\left\langle H\right\vert }_{1}{\left\langle H\right\vert }_{2}\right.+{\left\vert H\right\rangle }_{1}{\left\vert V\right\rangle }_{2}{\left\langle H\right\vert }_{1}{\left\langle V\right\vert }_{2}\\ \qquad+\,{\left\vert V\right\rangle }_{1}{\left\vert H\right\rangle }_{2}{\left\langle V\right\vert }_{1}{\left\langle H\right\vert }_{2}\left.\,+\,{\left\vert V\right\rangle }_{1}{\left\vert V\right\rangle }_{2}{\left\langle V\right\vert }_{1}{\left\langle V\right\vert }_{2}\right),\end{array}$$
(1)

where the subscripts denote the mode. For simplicity, we begin by considering the case where a single photon is emitted in each mode.

We can find the state immediately after the polarization grating shown in Fig. 1 to be

$$\begin{array}{r}{\hat{\rho }}_{{{{\rm{pol}}}}}=\hat{P}\left({x}_{1}\right){\hat{\rho }}_{1}\hat{P}\left({x}_{2}\right)\otimes \hat{P}\left({x}_{3}\right){\hat{\rho }}_{2}\hat{P}\left({x}_{4}\right),\end{array}$$
(2)

where \(\hat{P}(x)\) is the projective measurement given by

$$\begin{array}{r}\hat{P}\left(x\right)=\left[\begin{array}{ccc}{\cos }^{2}\left(\frac{\pi x}{L}\right)&&\cos \left(\frac{\pi x}{L}\right)\sin \left(\frac{\pi x}{L}\right)\\ \cos \left(\frac{\pi x}{L}\right)\sin \left(\frac{\pi x}{L}\right)&&{\sin }^{2}\left(\frac{\pi x}{L}\right)\end{array}\right].\end{array}$$
(3)

For ease of calculation, we utilize the Heisenberg picture, back-propagating the detector operators to the polarization grating. The point detector is modeled by \({\hat{O}}_{j,k,z}(X)={\hat{a}}_{j,z}^{{\dagger} }(X){\hat{a}}_{k,z}(X)\), where z is the distance between the grating and the measurement plane. The ladder operator \({\hat{a}}_{j,z}(X)\) is defined as

$$\begin{array}{r}{\hat{a}}_{j,z}(X)=\int\nolimits_{-\frac{L}{2}}^{\frac{L}{2}}dx{\hat{a}}_{j,0}\left(x\right)\exp[-\frac{2\pi i}{z\lambda }xX],\end{array}$$
(4)

where X is the position of the detector on the measurement plane, λ is the wavelength of the beam and j, k is the polarization of the operator. Eq. (4) describes the contribution of each point on the polarization grating plane to the detection measurement. Since we wish to keep the information of each interaction on the screen, we choose to calculate the four-point auto covariance by51,52

$$\begin{array}{rl}{G}_{jklm}^{\left(2\right)}\left({{{\bf{X}}}},z\right)={{{\rm{Tr}}}}[{\hat{\rho }}_{{{{\rm{pol}}}}}{\hat{a}}_{j,z}^{{\dagger} }\left({X}_{1}\right){\hat{a}}_{k,z}\left({X}_{2}\right){\hat{a}}_{l,z}^{{\dagger} }\left({X}_{3}\right){\hat{a}}_{m,z}\left({X}_{4}\right)],\end{array}$$
(5)

where \({{{\bf{X}}}}=\left[{X}_{1},{X}_{2},{X}_{3},{X}_{4}\right]\), allowing for the measurement of a post-selected coherence. We then set X2 = X1 and X4 = X3, since we are working with two point detectors. We allow the operators of the two detectors to commute, recovering the well-known expression for second-order coherence53. The second-order coherence of any post-selected measurement is then found to be

$$\begin{array}{l}{G}_{jklm}^{(2)}({{{\bf{X}}}},z)=\int\,d{x}_{1}\int\,d{x}_{2}\int\,d{x}_{3}\int\,d{x}_{4}{C}_{jklm}({{{\boldsymbol{x}}}})F({{{\boldsymbol{x}}}},{{{\bf{X}}}},z)\\\qquad\qquad\qquad\; \times\,[\delta \left({x}_{1}-{x}_{2}\right)\delta \left({x}_{3}-{x}_{4}\right)+\delta \left({x}_{1}-{x}_{4}\right)\delta \left({x}_{3}-{x}_{2}\right)],\end{array}$$
(6)

where the limits of integration for each integral is − L/2 to L/2, x = [x1, x2, x3, x4], \({C}_{jklm}\left({{{\boldsymbol{x}}}}\right)\) is the coefficient of the \({\left\vert j\right\rangle }_{1}{\left\vert k\right\rangle }_{2}{\left\langle l\right\vert }_{1}{\left\langle m\right\vert }_{2}\) element of the density matrix \({\hat{\rho }}_{{{{\rm{pol}}}}}\) in Eq. (2), and j, k, l, m ∈ {H, V}. Furthermore, \(F\left({{{\boldsymbol{x}}}},{{{\bf{X}}}},z\right)\) is given as

$$\begin{array}{r}F\left({{{\boldsymbol{x}}}},{{{\bf{X}}}},z\right)=\exp[\frac{2\pi i}{\lambda z}\left({X}_{4}{x}_{4}-{X}_{3}{x}_{3}+{X}_{2}{x}_{2}-{X}_{1}{x}_{1}\right)].\end{array}$$
(7)

We set X2 = X1 and X4 = X3, which properly describes the two point detectors allowing Eq. (6) to become a 2D Fourier transform6. By observing Eq. (6), it is important to note that there are two spatial correlations that contribute to the coherence at the measurement plane. One is the correlation of a photon with itself which existed prior to interacting with the polarizer, while the other is the spatial correlation gained by two photon scattering. It is important to note that the key difference between our formalism and the classical formalism is the existence of the scattering term as discussed in the Supplementary Notes 1 and 3. Due to the nature of projective measurements in Eq. (2), the density matrix \({\hat{\rho }}_{{{{\rm{pol}}}}}\) will no longer be diagonal in the horizontal-vertical basis, allowing for the beam to temporarily gain and lose polarization coherence6. The self-coherence of a photon results in the minimum coherence throughout all measurements in the far-field. The correlations from two photon scattering sets the maximum coherence and determines how it changes with the distance between the detectors.

To extend the description of a two-mode system comprising two photons to a multiphoton picture capable of handling any state, we need to propagate a value other than the photon statistics. Attempting to propagate a multiphoton field under the Schrodinger and Heisenberg pictures becomes computationally hard, scaling on the order of O(2nn!) where n represents the number of photons54. As a result, we propagate multiphoton states using a quantum version of the beam coherence-polarization (BCP) matrix6,55. This formalism allows us to estimate the evolution of the four-point-correlation matrix, reducing the total elements of interest. Consequently, the two-photon calculation represents the simplest case that the BCP matrix can handle and is in agreement with the general multiphoton picture. We begin by defining the BCP matrix as

$$\left\langle \hat{J}\left({X}_{1},{X}_{2},z\right)\right\rangle =\left[\begin{array}{ll}\left\langle {\hat{E}}_{H}^{{\dagger} }\left({X}_{1},z\right){\hat{E}}_{H}\left({X}_{2},z\right)\right\rangle &\left\langle {\hat{E}}_{H}^{{\dagger} }\left({X}_{1},z\right){\hat{E}}_{V}\left({X}_{2},z\right)\right\rangle \\ \left\langle {\hat{E}}_{V}^{{\dagger} }\left({X}_{1},z\right){\hat{E}}_{H}\left({X}_{2},z\right)\right\rangle &\left\langle {\hat{E}}_{V}^{{\dagger} }\left({X}_{1},z\right){\hat{E}}_{V}\left({X}_{2},z\right)\right\rangle \end{array}\right].$$
(8)

Here, the angle brackets denote the ensemble average, whereas the quantities \({\hat{E}}_{\alpha }^{{\dagger} }\left(X,z\right)\) and \({\hat{E}}_{\alpha }\left(X,z\right)\) represent the negative- and positive-frequency components of the α-polarized (with α = H, V) field-operator at the space-time point (X, z; t), respectively. We can then propagate the BCP matrix through the grating and to the measurement plane, by considering an initial BCP matrix of the form: \({I}_{2}\delta \left({X}_{1}-{X}_{2}\right)\)55. The details of the calculation can be found in the Supplementary Notes 1, 2 and 4. Upon reaching the measurement plane, we can find the second-order coherence matrix given by56

$$\begin{array}{r}{{{{\bf{G}}}}}^{(2)}({{{\boldsymbol{X}}}},z)=\langle \hat{J}\left({X}_{1},{X}_{2},z\right)\otimes \hat{J}\left({X}_{3},{X}_{4},z\right)\rangle .\end{array}$$
(9)

Each element of the G(2) matrix is a post-selected coherence matching each combination of polarizations shown in Eqs. (5)-(6). As shown in the Supplementary Note 2, the result obtained is equivalent to the approach described in Eqs. (1)-(7).

In order to demonstrate the results of our calculation, we first look at the second-order coherence of the horizontal mode in the far-field. By normalizing either Eq. (5) or the matrix element of Eq. (9), we find the coherence of the horizontal mode to be

$$\begin{array}{lll}{g}_{{{{\rm{HHHH}}}}}^{(2)}\,(\nu )=1+\frac{1}{16}{{{{\rm{sinc}}}}}^{2}\left(2-\nu \right)+\frac{5}{8}{{{{\rm{sinc}}}}}^{2}\left(\nu \right)+\frac{1}{16}{{{{\rm{sinc}}}}}^{2}\left(2+\nu \right)\\\qquad\qquad\;\;+\;\frac{1}{4}{{{\rm{sinc}}}}\left(2-\nu \right){{{\rm{sinc}}}}\left(1-\nu \right)+\frac{3}{8}{{{{\rm{sinc}}}}}^{2}\left(1-\nu \right)\\\qquad\qquad\;\;+\;\,\frac{1}{4}{{{\rm{sinc}}}}\left(1+\nu \right)\left({{{\rm{sinc}}}}\left(2+\nu \right)+{{{\rm{sinc}}}}\left(1-\nu \right)\right)\\\qquad\qquad\;\;+\;\frac{3}{8}{{{{\rm{sinc}}}}}^{2}\left(1+\nu \right)+\frac{1}{8}{{{\rm{sinc}}}}\left(\nu \right)\left(\right.{{{\rm{sinc}}}}\left(2-\nu \right)\\\qquad\qquad\;\;+\;{{{\rm{sinc}}}}\left(2+\nu \right)+6{{{\rm{sinc}}}}\left(1-\nu \right)+6{{{\rm{sinc}}}}\left(1+\nu \right)\left)\right.,\end{array}$$
(10)

where \({g}_{jklm}^{(2)}\) is the normalized second-order coherence. Here \({{{\rm{sinc}}}}(\nu )=\sin (\pi \nu )/\left(\pi \nu \right)\) and \(\nu =L{{\Delta }}X/\left(\lambda z\right)\). Therefore, \({g}_{jklm}^{(2)}\) depends on the distance between the detectors ΔX = X1 − X2, the length of the polarization grating L, the wavelength λ, and the distance in the far field z. The same holds true for all other \({g}_{jklm}^{(2)}\), where each expression can be found in the Supplementary Notes 4 and 5. Since Eq. (10) applies to all incoherent unpolarized states, we will perform the analysis for two mode thermal states, as the statistical properties are well studied35,36,38,53.

Quantum Coherence of a Multiphoton System upon Propagation

As shown in Fig. 2, increasing the separation ΔX of the detectors causes the correlations to gradually decrease. Once ν ≈ 2.7, the detectors become uncorrelated. We note that g(2)(ν) = 1, ν ≠ 0 represents an uncorrelated measurement since this can only be true when the two spatial modes become separable. In addition, when one of the two measured modes is no longer contributing to the measurement we get a g(2)(ν) = 0. Interestingly, by fixing the distance ΔX between the two detectors, we can increase the correlations by moving the measurement plane further into the far-field. This is equivalent to decreasing ν, causing correlations to increase to a possible maximum value of g(2)(0) = 1.62. By measuring \({g}_{{{{\rm{grating}}}}}^{(2)}(0)\) immediately after the polarization grating at x = 0, a horizontally polarized beam is measured with a \({g}_{{{{\rm{grating}}}}}^{(2)}(0)=2\). We note the theory we presented only applies to the far-field, therefore these two values do not contradict each other. While the exact transition between the near and the far-fields are beyond the scope of the paper, we note that the horizontal mode along the central axis becomes more coherent as it propagates to the far-field, as predicted by the van Cittert-Zernike theorem14,51.

Fig. 2: The second-order coherence for various post-selected measurements in the far field.
figure 2

The x-axis is how the g(2) changes as a function of \(\nu =L{{\Delta }}X/\left(\lambda z\right)\) while keeping L, λ and z fixed. As the detectors move further apart, the spatial correlations created by the polarization grating decrease until they diminish entirely at ν ≈ 2.7. In addition, certain post-selected measurements allow us to quantify the coherence between two fields that possess sub-Poissonian-like statistics, suggesting the possibility of sub-shot-noise measurements. Note, these measurements can be also performed using quantum state tomography57.

Full size image

Setting one detector to measure the vertical mode and the other detector the horizontal mode, given by \({g}_{{{{\rm{HHVV}}}}}^{(2)}\) in Fig. 2, we can measure the coherence between the horizontal and vertical mode. This post-selective measurement results in a different effect from when we only measured only the horizontal mode. Placing the detectors immediately after the polarization grating at x = 0, we measure \({g}_{{{{\rm{grating}}}}}^{(2)}(0)=0\) since there is no vertically polarized mode. However, we measure g(2)(0) ≈ 1.1 when the beam is propagated into the far-field. In this case, the polarization grating leads to the thermalization of the beam38. This can be verified by removing the polarization grating and repeating the measurement giving \({g}_{{{{\rm{initial}}}}}^{(2)}(0)=1\) since the two modes are completely uncorrelated.

The measurements of \({g}_{{{{\rm{HHHH}}}}}^{(2)}\) and \({g}_{{{{\rm{HHVV}}}}}^{(2)}\) can be performed using point detectors, however we predict more interesting effects that can be observed through the full characterization of the field. This information can be obtained through quantum state tomography57,58. We find that the second-order coherence \({g}_{{{{\rm{VHHV}}}}}^{(2)}\), \({g}_{{{{\rm{HHVH}}}}}^{(2)}\) and \({g}_{{{{\rm{HVHV}}}}}^{(2)}\) is below one suggesting sub-Poissonian-like statistics—i.e. a photon distribution narrower than the characteristic Poissonian distribution of coherent light-sources—which potentially may allow for sub-shot-noise measurement59. This feature is found using the BCP matrix approach, therefore it is true for all unpolarized incoherent fields. The sub-Poissonian-like statistics were achieved only with the use of post-selection without nonlinear interactions24,25,27. Our theoretical formalism unveils the possibility of modifying the photon statistics of the electromagnetic field in free space without resorting to complex light-matter interactions28,38,45,60. Another interesting feature is that the \({g}_{{{{\rm{VHHV}}}}}^{(2)}\) decays and resurrects at ν ≈ 1.6 before decaying again. The fact that these elements are nonzero contrasts classical analyses where off-diagonal correlations of the system are zero since the system is unpolarized in the far-field6,49. In fact the classical analysis fails to describe any emergent phenomena shown since the final density matrix is the identity matrix (see Supplementary Note 3). This new quantum degree of polarization is likely caused by the two photon scattering induced by the polarization grating. Furthermore, the creation of nonzero off-diagonal elements suggests that we induced correlations between orthogonally-polarized fields in our system.

The sub-Possonian statistics are exclusive to unpolarized systems. Returning to Eq. (6) we have that there are two correlations contributing to the final coherence, one from the photons self-coherence that existed prior to interaction with the screen and another coherence term that comes from photon scattering. For unpolarized states there is no initial correlations in the off-diagonal elements of the density matrix, since by definition the off-diagonal elements are zero49. This results in the first term of Eq. (6) to be zero for all off-diagonal elements. As noted above, this term sets the minimum value of the coherence measurement to zero, allowing for the measurement of sub-Poissonian-like statistics.

Finally, we would like to highlight the fact that the quantum statistical properties of multiphoton systems can change upon propagation without complex light-matter interactions due to the scattering of their constituent single-mode photons carrying different polarizations. This effect is quantified through the second-order quantum coherence g(2)(τ = 0) defined as \({g}^{(2)}(\tau =0)=1+(\langle {({{\Delta }}\hat{n})}^{2}\rangle -\langle \hat{n}\rangle )/{\langle \hat{n}\rangle }^{2}\)34,36. In this case, the averaged quantities in g(2)(τ = 0) are obtained through the density matrix of the system’s state, as described in Eq. (1), at different spatial coordinates (X, z). In Fig. 3, we report the photon-number distribution of the combined vertical-horizontal multiphoton field. In this case, a single photon-number-resolving detector was placed at X = 0.4L35. Note that by selecting the proper propagation distance z, one could, in principle, generate on-demand multiphoton systems with sub-Poissonian-like or Poissonian statistics27,38, see the Supplementary Note 6 for details on the combined-field photon-number distribution calculation. As indicated in Figs. 2 and 3, the evolution of quantum coherence upon propagation lies at the heart of the quantum van Cittert-Zernike theorem for multiphoton systems. Finally, we note that a generalized form of the Hanbury Brown and Twiss effect suggests possible connections with the quantum van Cittert-Zernike theorem61. We believe that it will be interesting to pursue further investigations along this research direction.

Fig. 3: The modification of the photon-number distribution and quantum coherence of a thermal multiphoton system upon propagation.
figure 3

In this case, the multiphoton system comprises a mixture of single-mode photons with either vertical or horizontal polarization. We assumed a single photon-number-resolving detector placed at different propagation distances: a z = 0, b z = 50L, c z = 100L, d z = 150L, e z = 200L, f z = 250L, g z = 300L, h z = 350L. In the transverse plane, the photon-number-resolving detector is placed at X = 0.4L.

Full size image

In conclusion, we have investigated new mechanisms to control the photon statistics of multiphoton systems. We describe these interactions using a quantum version of the van Cittert-Zernike theorem. Specifically, by considering a polarization grating together with conditional measurements, we show that it is possible to control the quantum coherence of multiphoton systems. Moreover, we unveiled the possibility of producing multiphoton systems with sub-Poissonian-like statistics without complex light-matter interactions28,38,45,60. This possibility cannot be explained through the classical theory of optical coherence4. Thus, our work demonstrates that the multiphoton quantum van Cittert-Zernike theorem will have important implications for describing the evolution of the properties of quantum coherence of many-body bosonic systems28,33. As such, our findings could offer alternatives to creating novel states of light by controlling the collective interactions of many single-photon emitters30.

Methods

Calculation of the quantum BCP matrix

We now present the detailed calculation where two-photon correlations are built up during propagation. Let us start from a spatially incoherent and unpolarized source, whose four-point correlation matrix is written, in the \(\{\left\vert HH\right\rangle ,\left\vert HV\right\rangle ,\left\vert VH\right\rangle ,\left\vert VV\right\rangle \}\) basis, as

$$\begin{array}{ll}{G}_{{{{\rm{ini}}}}}\left({x}_{1},{x}_{2};{x}_{3},{x}_{4},0\right)={j}_{{{{\rm{ini}}}}}\left({x}_{1},{x}_{2},0\right)\otimes {j}_{{{{\rm{ini}}}}}\left({x}_{3},{x}_{4},0\right)\\\qquad\qquad\qquad\qquad\quad\, ={\lambda }^{4}{I}_{0}^{2}\left[\delta \left({x}_{2}-{x}_{1}\right)\delta \left({x}_{3}-{x}_{4}\right)+\delta \left({x}_{2}-{x}_{3}\right)\delta \left({x}_{1}-{x}_{4}\right)\right]\\\qquad\qquad\qquad\qquad\quad\;\;\;\; \times \left[\begin{array}{cccc}1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1\end{array}\right],\end{array}$$
(11)

where

$${j}_{{{{\rm{ini}}}}}\left({x}_{1},{x}_{2},0\right)={\lambda }^{2}{I}_{0}\delta \left({x}_{2}-{x}_{1}\right)\left[\begin{array}{cc}1&0\\ 0&1\end{array}\right]$$
(12)

stands for the two-point correlation matrix for a spatially incoherent and unpolarized photon source, and I0 describes a constant-intensity factor. Note that, for the sake of simplicity, we have restricted ourselves to a one-dimensional case, i.e., we have taken only one element of the transversal vector r = (x, y).

To polarize the source, we make use of a linear polarizer. Specifically, we cover the source with a linear polarization grating whose angle between its transmission axis and the x-axis is a linear function of the form θ = πx/L, with L being the length of the grating. The four-point correlation matrix after the polarization grating can thus be written as

$$\begin{array}{lll}{G}_{{{{\rm{out}}}}}\left({x}_{1},{x}_{2};{x}_{3},{x}_{4},0\right)=\left[{P}^{{\dagger} }\left({x}_{1}\right){j}_{{{{\rm{ini}}}}}\left({x}_{1},{x}_{2},0\right)P\left({x}_{2}\right)\right]\otimes \left[{P}^{{\dagger} }\left({x}_{3}\right){j}_{{{{\rm{ini}}}}}\left({x}_{3},{x}_{4},0\right)P\left({x}_{4}\right)\right]\\\qquad\qquad\qquad\qquad\qquad\;\, \times {{{\rm{rect}}}}\left({x}_{1}/L\right){{{\rm{rect}}}}\left({x}_{2}/L\right){{{\rm{rect}}}}\left({x}_{3}/L\right){{{\rm{rect}}}}\left({x}_{4}/L\right)\\\qquad\qquad\qquad\qquad\qquad\;\, \times {\lambda }^{4}{I}_{0}^{2}\left[\delta \left({x}_{2}-{x}_{1}\right)\delta \left({x}_{3}-{x}_{4}\right)+\delta \left({x}_{2}-{x}_{3}\right)\delta \left({x}_{1}-{x}_{4}\right)\right]\end{array}$$
(13)

where the product of \({{{\rm{rect}}}}\left(\cdots \right)\) functions describe the finite size of the source, and the action of the polarization grating is given by Eq. (3).

Furthermore, the elements defined by the previous equation follow a propagation formula of the form62

$$\begin{array}{l}{G}_{jklm}\left({{{{\boldsymbol{r}}}}}_{1},{{{{\boldsymbol{r}}}}}_{2};{{{{\boldsymbol{r}}}}}_{3},{{{{\boldsymbol{r}}}}}_{4},z\right)=\int\int\int\int\,{G}_{jklm}\left({{{{\boldsymbol{\rho }}}}}_{1},{{{{\boldsymbol{\rho }}}}}_{2};{{{{\boldsymbol{\rho }}}}}_{3},{{{{\boldsymbol{\rho }}}}}_{4},0\right){K}^{* }\left({{{{\boldsymbol{r}}}}}_{1},{{{{\boldsymbol{\rho }}}}}_{1},z\right)K\left({{{{\boldsymbol{r}}}}}_{2},{{{{\boldsymbol{\rho }}}}}_{2},z\right)\\\qquad\qquad\qquad\qquad\qquad\;\; \times\, {K}^{* }\left({{{{\boldsymbol{r}}}}}_{3},{{{{\boldsymbol{\rho }}}}}_{3},z\right)K\left({{{{\boldsymbol{r}}}}}_{4},{{{{\boldsymbol{\rho }}}}}_{4},z\right){d}^{2}{{{{\boldsymbol{\rho }}}}}_{1}{d}^{2}{{{{\boldsymbol{\rho }}}}}_{2}{d}^{2}{{{{\boldsymbol{\rho }}}}}_{3}{d}^{2}{{{{\boldsymbol{\rho }}}}}_{4},\end{array}$$
(14)

with the Fresnel propagation kernel defined by

$$K({{{\boldsymbol{r}}}},{{{\boldsymbol{\rho }}}},z)=\frac{-i\exp (ikz)}{\lambda z}\exp \left[\frac{ik}{2z}{({{{\boldsymbol{r}}}}-{{{\boldsymbol{\rho }}}})}^{2}\right],$$
(15)

where k = 2π/λ. Interestingly, in the context of the scalar theory, a spatially incoherent source is characterized by means of a delta-correlated intensity function, which indicates that subfields—making up for the whole source—at any two distinct points across the source plane are uncorrelated36.

By substituting Eq. (13) into the one dimensional version of Eq. (14), we can obtain the explicit form of the polarized, four-point correlation matrix elements. As an example, we can find that, in the far-field, the normalized four-point correlation function for H-polarized photons reads as

$$\begin{array}{lll}{G}_{{{{\rm{HHHH}}}}}\left({\nu }_{1},{\nu }_{2},{\nu }_{3},{\nu }_{4};z\right)=\,\left[{{\mathrm{sinc}}}\,\left({\nu }_{1}\right)+\frac{1}{2}{{\mathrm{sinc}}}\,\left({\nu }_{1}-1\right)+\frac{1}{2}{{\mathrm{sinc}}}\,\left({\nu }_{1}+1\right)\right]\\ \quad\qquad\qquad\qquad\qquad\qquad\times \left[{{\mathrm{sinc}}}\,\left({\nu }_{2}\right)+\frac{1}{2}{{\mathrm{sinc}}}\,\left({\nu }_{2}-1\right)+\frac{1}{2}{{\mathrm{sinc}}}\,\left({\nu }_{2}+1\right)\right]\\ \quad\qquad\qquad\qquad\qquad\qquad+\frac{1}{16}\left[\right.{{\mathrm{sinc}}}\,\left(2+{\nu }_{3}\right)\left({{\mathrm{sinc}}}\,\left({\nu }_{4}\right)+2{{\mathrm{sinc}}}\,\left(1-{\nu }_{4}\right)+{{\mathrm{sinc}}}\,\left(2-{\nu }_{4}\right)\right)\\ \quad\qquad\qquad\qquad\qquad\qquad+2{{\mathrm{sinc}}}\,\left(1+{\nu }_{3}\right)\left(\right.3{{\mathrm{sinc}}}\,\left({\nu }_{4}\right)+3{{\mathrm{sinc}}}\,\left(1-{\nu }_{4}\right) \\ \quad\qquad\qquad\qquad\qquad\qquad+{{\mathrm{sinc}}}\,\left(2-{\nu }_{4}\right)+{{\mathrm{sinc}}}\,\left(1+{\nu }_{4}\right)\left)\right.\\ \quad\qquad\qquad\qquad\qquad\qquad+{{\mathrm{sinc}}}\,\left(2-{\nu }_{3}\right)\left({{\mathrm{sinc}}}\,\left({\nu }_{4}\right)+2{{\mathrm{sinc}}}\,\left(1+{\nu }_{4}\right)+{{\mathrm{sinc}}}\,\left(2+{\nu }_{4}\right)\right)\\ \quad\qquad\qquad\qquad\qquad\qquad+2{{\mathrm{sinc}}}\,\left(1-{\nu }_{3}\right)\left(\right.3{{\mathrm{sinc}}}\,\left({\nu }_{4}\right)+{{\mathrm{sinc}}}\,\left(1-{\nu }_{4}\right)\\ \quad\qquad\qquad\qquad\qquad\qquad+3{{\mathrm{sinc}}}\,\left(1+{\nu }_{4}\right)+{{\mathrm{sinc}}}\,\left(2+{\nu }_{4}\right)\left)\right.\\ \quad\qquad\qquad\qquad\qquad\qquad+{{\mathrm{sinc}}}\,\left({\nu }_{3}\right)\left(\right.10{{\mathrm{sinc}}}\,\left({\nu }_{4}\right)+6{{\mathrm{sinc}}}\,\left(1-{\nu }_{4}\right)+{{\mathrm{sinc}}}\,\left(2-{\nu }_{4}\right) \\ \quad\qquad\qquad\qquad\qquad\qquad+6{{\mathrm{sinc}}}\,\left(1+{\nu }_{4}\right)+{{\mathrm{sinc}}}\,\left(2+{\nu }_{4}\right)\left)\right.\left]\right.\end{array}$$
(16)

with

$${\nu }_{1}=L\frac{{x}_{2}-{x}_{3}}{\lambda z};{\nu }_{2}=L\frac{{x}_{4}-{x}_{1}}{\lambda z};{\nu }_{3}=L\frac{{x}_{2}-{x}_{1}}{\lambda z};{\nu }_{4}=L\frac{{x}_{4}-{x}_{3}}{\lambda z}.$$
(17)

Finally, by realizing that when monitoring the two-photon correlation function with two detectors, at the observation plane in z, we must set: x2 = x3 and x1 = x4, we find that

$${G}_{{{{\rm{HHHH}}}}}\left({\nu }_{1},{\nu }_{2},{\nu }_{3},{\nu }_{4};z\right)={G}_{{{{\rm{HHHH}}}}}\left(0,0,{\nu }_{1},-{\nu }_{1};z\right).$$
(18)

We can follow a similar procedure as above to obtain the remaining terms of the four-point correlation matrix.