Emergence of multiphoton quantum coherence by light propagation

The description of the evolution of spatial, temporal, and polarization properties of the light field gave birth to the classical theory of optical coherence [1, 2, 3, 4, 5, 6, 7]. Naturally, these properties of light can be fully characterized through the classical electromagnetic theory [3]. Furthermore, there has been interest in describing the evolution of propagating quantum optical fields endowed with these classical properties [8, 9]. This has been accomplished by virtue of the Wolf equation and the van Cittert-Zernike theorem [1, 8, 10, 11]. Nevertheless, there is a long-sought goal in describing the evolution of the inherent quantum mechanical properties of the light field that define its nature and kind [12, 13]. Such formalism would enable modeling the evolution of the excitation mode of propagating electromagnetic fields in the Fock number basis. Given the large number of scattering and interference processes that can take place in a quantum optical system with many photons, this ambitious description remains elusive [14, 15, 16, 17, 18]. Although, it is essential to describe the evolution of propagating multiphoton wavepackets in diverse quantum photonic devices [18, 19, 20, 21].

The quantum theory of optical coherence developed by Glauber and Sudarshan provides a description of the excitation mode of the electromagnetic field [13, 12, 22]. This elegant formalism led to the identification of different kinds of light that are characterized by distinct quantum statistical fluctuations and noise levels [15, 23, 12, 22, 24]. As such, a particular quantum state of light is typically associated with a specific emission process and a light source [22]. Moreover, the quantum theory of electromagnetic radiation enables describing light-matter interactions [25]. These consist of absorption and emission processes that can lead to the modification of the excitation mode of the light field and consequently to different kinds of light [25, 22]. This possibility has triggered interest in achieving optical non-linearities at the single-photon level to engineer and control quantum states of light [26, 27, 28, 29]. Thus, it is believed that the excitation mode of the light field, and its quantum properties of coherence, remain unchanged upon propagation in free space [22, 25].

We demonstrate that the statistical fluctuations of thermal light fields, and their quantum properties of coherence, can be modified upon propagation in the absence of light-matter interactions [9, 33]. This effect results from the scattering of multiphoton wavepackets that propagate in free space. The large number of interference effects hosted by multiphoton systems with up to twenty particles leads to a modified light field with evolving quantum statistical properties [16, 14]. Further, we show that the evolution of multiphoton quantum coherence can be described by the nonclassical formulation of the van Cittert-Zernike theorem [9]. Interestingly, our description of propagating multiphoton quantum coherence unveils conditions under which multiphoton systems with sub-shot-noise quantum properties are formed [34]. Remarkably, these quantum multiphoton systems are produced upon propagation in the absence of optical nonlinearities [26, 27, 28, 29]. As such, we believe that our findings provide an all-optical alternative for the preparation of multiphoton systems with nonclassical statistics. Given the relevance of photonic quantum control for multiple quantum technologies, similar functionalities have been explored in nonlinear optical systems, photonic lattices, plasmonic systems, and Bose-Einstein condensates [35, 17, 36, 26, 27, 28, 29].

The optical system under consideration is depicted in Figure 1a. In this case, an unpolarized thermal field is scattered by an optical grating to produce multiphoton wavepackets with distinct polarization properties at different transverse spatial positions [5]. The scattered photons contained in the thermal beam interfere upon propagation to change the statistical fluctuations of the field [9]. Interestingly, these effects enable the modification of the quantum properties of coherence of the initial multiphoton thermal system in free space. As discussed in the Supplementary Information (SI), we describe our initial thermal system as an incoherent superposition of coherent states [13, 12, 37]

$$\hat{\rho}=\int d\Sigma\bigotimes_{\bm{s}}\Big{(}|\alpha\rangle\langle\alpha|_{\Sigma,H,\bm{s}}+|\alpha\rangle\langle\alpha|_{\Sigma,V,\bm{s}}\Big{)},$$

(1)

where $|\alpha\rangle_{\Sigma,B,\bm{s}}$ represents a coherent state of amplitude $\alpha$ with random mode-structure $\Sigma$, where $\hat{a}_{\Sigma,B,\bm{s}}=\int d\bm{\rho}\text{Rect}[(\bm{s}-\bm{\rho})/d]\Sigma(\bm{\rho})\hat{a}_{B}(\bm{\rho})$ and polarization $B\in\{H,V\}$ (see SI). Furthermore, the tensor product over positions $\bm{s}$ represents the pixelated transverse beam profile where each pixel has sidelength $d$.

After the polarization grating, the resulting state is obtained via the transformation

$\displaystyle\hat{a}_{B}(\bm{\rho})\rightarrow P_{HB}(\bm{\rho})\hat{a}_{H}(\bm{\rho})+P_{VB}(\bm{\rho})\hat{a}_{V}(\bm{\rho})+P_{\emptyset B}(\bm{\rho})\hat{a}_{\emptyset}(\bm{\rho}),$

(2)

where $\hat{a}_{\emptyset}(\bm{\rho})$ is the mode for photon loss and $P_{AB}$ are the components of the transformation

$$\bm{P}(x)=\begin{bmatrix}\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)\\ \sin\left(\frac{\pi x}{L}\right)&\cos\left(\frac{\pi x}{L}\right)\end{bmatrix}$$

(3)

where $A\in\{H,V,\emptyset\}$ and $L$ is the length of the polarization grating. The beam described by Eq. (1) is then propagated by a distance of $z$ before being measured by two pairs of photon-number-resolving (PNR) detectors [23, 38]. This propagation can be modeled through the Fresnel approximation on the mode structure of the initial beam [39]. We can then compute the second-order correlation function
$G^{(2)}_{ijkl}(\bm{r}_{1},\bm{r}_{2},z)$ for the post-selected polarization measurements in the detection plane as [2]

		$\displaystyle G^{(2)}_{ijkl}(\bm{r}_{1},\bm{r}_{2},z)=\langle\hat{a}^{\dagger}_{i}(\bm{r}_{1})\hat{a}^{\dagger}_{j}(\bm{r}_{2})\hat{a}_{k}(\bm{r}_{1})\hat{a}_{l}(\bm{r}_{2})\rangle$		(4)
		$\displaystyle=I_{0}\int d\bm{\rho}_{1}d\bm{\rho}_{2}d\bm{\rho}_{3}d\bm{\rho}_{4}F(\bm{r}_{1},\bm{r}_{2},\bm{\rho}_{1},\bm{\rho}_{2},\bm{\rho}_{3},\bm{\rho}_{4},z)$
		$\displaystyle\times\Big{[}\delta(\bm{\rho}_{1}-\bm{\rho}_{3})\delta(\bm{\rho}_{2}-\bm{\rho}_{4})+\delta(\bm{\rho}_{1}-\bm{\rho}_{4})\delta(\bm{\rho}_{2}-\bm{\rho}_{3})\Big{]}.$

Remarkably, the Dirac-delta functions in Eq. (4) demonstrate the presence of nontrivial correlations. Given the complexity of $I_{0}$ and $F(\bm{r}_{1},\bm{r}_{2},\bm{\rho}_{1},\bm{\rho}_{2},\bm{\rho}_{3},\bm{\rho}_{4},z)$, their explicit expressions are given in the SI. These describe the coherence of a photon with itself, which existed prior to interacting with the grating, and the spatial coherence gained by multiphoton scattering. These terms unveil the possibility of modifying quantum coherence of multiphoton systems upon propagation [9]. We use this approach to describe the correlation properties of the multiphoton wavepackets that constitute our light beam. This allows us to use an equivalent density matrix for the system $\hat{\rho}_{ijkl}(z)$ (see SI) at the detection plane to compute its corresponding joint photon-number distribution $p_{ijkl}(n_{1},n_{2},z)$ as

$\displaystyle p(n_{1},n_{2},z)=\text{Tr}\left[\hat{\rho}_{ijkl}(z)|n_{1},n_{2}\rangle\langle n_{1},n_{2}|\right].$

(5)

As we shall see in the next section, these formulae allow for the prediction of very interesting correlation effects. Specifically, they predict that the statistical make-up of the light field is changing upon propagation in free space. The classical analogue to this behavior is explained by the van Cittert-Zernike theorem [10, 11], which predicts that the classical coherence properties of a light source can change upon free-space propagation. Therefore, we interpret our results in Eq. (4) as those of a quantum van Cittert-Zernike theorem. This is because they predict the modification of quantum coherence upon free-space propagation, and that is directly analogous to the classical theorem’s predictions. Specifically, Eq. (4) predicts this free-space quantum modification through the nontrivial scattering effects induced by the Dirac-delta functions. Interestingly, these delta functions arise from the unique coherence properties of thermal light (see SI for further details). Furthermore, Eq. (5) allows us to study multiparticle quantum coherence, which is also changing upon free-space propagation. These quantum van Cittert-Zernike effects, therefore, are not only arising in polarization subsystems, but also in multiphoton subsystems. This showcases the fundamental and intrinsic quantum impacts of free-space propagation on our state.

We explore the modification of the quantum coherence properties of propagating multiphoton systems using the experimental setup in Figure 1b. We use a combination of waveplates and a spatial light modulator (SLM) to rotate the polarization properties of our multiphoton system at any transverse position [30]. In addition, this arrangement enables us to characterize the polarization and photon-number distribution of multiphoton systems at different propagation planes. Specifically, we perform measurements at different propagation planes associated with the propagation parameter $\nu=L\Delta X/(\lambda z)$. Here, the transverse distance between detectors is described by $\Delta X$ and the wavelength of the beam by $\lambda$. As demonstrated in Figure 2, the many interference effects hosted by the propagating multiphoton system modify the photon-number distribution of the polarized components of the initial beam [15, 17, 20]. These processes lead to multiphoton systems with different quantum fluctuations and quantum properties of coherence [15, 40]. Each multiphoton system is characterized through the degree of second-order self coherence

$$g^{(2)}_{\nu}(0)=\frac{G^{(2)}_{\text{HHHH}}(\bm{r},\bm{r},z)}{G^{(1)}_{\text{HH}}(\bm{r},z)^{2}},$$

(6)

where $G^{(1)}_{i,j}(\bm{r},z)=\langle\hat{a}^{\dagger}_{i}(\bm{r})\hat{a}_{j}(\bm{r})\rangle=\sqrt{I_{0}}L/(2z^{2}\lambda^{2})$. Interestingly, the multiphoton system in Figure 2a is nearly thermal [22]. However, propagation leads to different kinds of multiphoton wavepackets. We show these from Figure 2a to f. The coherence properties of the multiphoton system in Figure 2f approach those observed in coherent light beams [34]. Remarkably, the conversion processes of the multiphoton system, described by Eq. (5), take place in free space in the absence of light-matter interactions [25, 26, 27, 28, 35, 17, 36].

The polarization and photon-number properties of the propagating light beam at different transverse and longitudinal positions host many forms of interference effects [37, 14]. We explore these dynamics by isolating the constituent multiphoton subsystems of the propagating beam. Each multiphoton subsystem, characterized by different polarization properties, exhibits different degrees of second-order coherence [20]. In the experiment, we perform projective measurements on polarization. These measurements unveil the possibility of extracting multiphoton subsystems with attenuated quantum fluctuations below the shot-noise limit [38, 32]. In this case, we use the four detectors depicted in the experimental setup in Figure 1b to perform full characterization of polarization [31]. These measurements enable us to characterize correlations of multiphoton subsystems with different polarization properties, which are reported in Figure 3. We plot the degree of second-order mutual coherence

$$g^{(2)}_{ijkl}(\bm{r}_{1},\bm{r}_{2},z)=\frac{G^{(2)}_{ijkl}(\bm{r}_{1},\bm{r}_{2},z)}{G^{(1)}_{i,j}(\bm{r}_{1},z)G^{(1)}_{k,l}(\bm{r}_{2},z)}.$$

(7)

The propagation of the multiphoton subsystem described by $g^{(2)}_{\text{HHHH}}$ shows a modification of the quantum statistics from super-Poissonian to nearly Poissonian [40, 34]. A similar situation prevails for the multiphoton subsystem described by $g^{(2)}_{\text{HVHV}}$. It is worth noticing that the multiphoton subsystems described by $g^{(2)}_{\text{VHHV}}$ and $g^{(2)}_{\text{HHVV}}$ unveil the possibility of extracting multiphoton subsystems with sub-shot-noise properties [41]. This implies photon-number distributions narrower than the characteristic Poissonian distribution of coherent light [22, 42]. This peculiar feature might unlock novel paths towards the implementation of sensitive measurements with sub-shot-noise fluctuations [43].

We now turn our attention to describe the quantum coherence evolution of propagating multiphoton wavepackets. This is explored by projecting the polarized components of the initial thermal beam into its constituent multiphoton wavepackets [22]. In this case, we analyze wavepackets with $n_{1}+n_{2}$ number of photons. The number of photons detected in arm 1 of our experiment is described by $n_{1}$, whereas $n_{2}$ indicates the number of photons detected in arm 2. Our findings unveil that despite the fact that the degree of second-order coherence $g^{(2)}_{\text{HHHH}}$ in Figure 3 is produced by its constituent wavepackets, these show a completely different evolution of their properties of coherence. Our experimental measurements of these wavepackets can be found from Figure 4a to c. The results in Figure 4a indicate that multiphoton wavepackets, in which $n_{1}$ and $n_{2}$ are the same, show a particular evolution. In contrast, propagating wavepackets with asymmetric values of $n_{1}$ and $n_{2}$ show different trends in the modification of quantum coherence, these are shown in Figure 4b and c. The propagation of these wavepackets can be described using Eq. (5). Specifically, we can calculate the multiphoton degree of second-order mutual coherence [44]

$$\tilde{g}^{(2)}_{ijkl}(n_{1},n_{2},z)=\frac{p_{ijkl}(n_{1},n_{2},z)}{\sum_{n}p_{ijkl}(n,n_{2},z)\sum_{m}p_{ijkl}(n_{1},m,z)}.$$

(8)

Furthermore, the multiphoton wavepackets in the projected beam, characterized by $g^{(2)}_{\text{HVHV}}$, exhibit the multiphoton dynamics reported from Figure 4d to f. These results suggest that the multiphoton dynamics in Figure 4 depend on the number of photons in each of the measured wavepackets.

This quantum field theoretic approach to studying the quantum van Cittert-Zernike theorem provides us with the ability to describe the propagation dynamics of the multiphoton systems that constitute classical light beams. We used this formalism to extract propagating multiphoton subsystems, with quantum statistical properties, from unpolarized thermal light fields. While nonlinear light-matter interactions offer the possibility of engineering complex quantum systems [26, 27, 28], our scheme exploits linear propagation of multiphoton systems [29, 45]. This feature enabled us to exploit multiphoton scattering in free space to produce wavepackets with different quantum statistical properties [15]. As such, our work combines the benefits of post-selective measurements with those of multiphoton scattering in propagating light beams, and it allows us to study the modification of the quantum statistical properties of multiphoton wavepackets in free space. Although, the incoherent combination of light beams with different polarization properties can lead to the modification of the degree of second-order coherence [17, 46], we performed direct measurements of polarized multiphoton systems with propagating quantum coherence properties (see Fig. 4a to c). Interestingly, these processes are defined by the number of particles in the measured multiphoton system. Consequently, these findings have important implications for all-optical engineering of multiphoton quantum systems.

We demonstrated the possibility of modifying the excitation mode of thermal multiphoton fields through free space propagation. This modification stems from the scattering of multiphoton wavepackets in the absence of light-matter interactions [25, 26, 27, 28, 35, 17, 36, 47]. The modification of the excitation mode of a photonic system and its associated quantum fluctuations result in the formation of different light fields with distinct quantum coherence properties [22, 12, 13]. The evolution of multiphoton quantum coherence is described through the nonclassical formulation of the van Cittert-Zernike theorem, unveiling conditions for the formation of multiphoton systems with attenuated quantum fluctuations below the sub-shot-noise limit [38, 48, 43]. Notably, these quantum multiphoton systems emerge in the absence of optical nonlinearities, suggesting an all-optical approach for extracting multiphoton wavepackets with nonclassical statistics. We believe that the identification of this surprising multiphoton dynamics has important implications for multiphoton protocols quantum information [18, 21].

J.F. acknowledges funding from the National Science Foundation through Grant No. OMA 2231387. M.H., R.B.D., C.Y. and O.S.M.L. acknowledge support from the Army Research Office (ARO), through the Early Career Program (ECP) under the grant no. W911NF-22-1-0088. R.J.L.-M. thankfully acknowledges financial support by DGAPA-UNAM under the project UNAM-PAPIIT IN101623.

In our experiment, we utilized unpolarized multimode thermal light. The light is thermal such that the electric field $E^{(+)}(x)$ obeys complex-Gaussian statistics at each point $x$, and that the mean $\langle E^{(+)}(x)\rangle$ is zero. It is unpolarized such that it is an equal mixture between two orthogonal polarizations (here we choose horizontal ($H$) and vertical ($V$)). Finally, any two spatial projections on this source will be statistically independent. It is our goal to study the quantum properties of such a source as it propagates through our experimental setup. We will now present a sufficient quantum description for such a light source before propagation.

The generation of unpolarized multimode thermal light is accomplished by mixing coherent states with different amplitudes [37]. We then pixelize the source and assume that each pixel obeys independent polarization statistics. Such a source can then be written as

$$\hat{\rho}=\int d\Sigma\bigotimes_{\bm{s}}\Big{(}|\alpha\rangle\langle\alpha|_{\Sigma,H,\bm{s}}+|\alpha\rangle\langle\alpha|_{\Sigma,V,\bm{s}}\Big{)}.$$

(9)

Here, each $\bm{s}$ represents the position of a pixel, $\alpha$ is a coherent amplitude, and the coherent states $|\alpha\rangle_{\Sigma,B,\bm{s}}$ are defined by the modes

$$\hat{a}_{\Sigma,B,\bm{s}}=\int d\bm{\rho}\text{ Rect}\left[(\bm{s}-\bm{\rho})/d\right]\Sigma(\bm{\rho})\hat{a}_{B}(\bm{\rho}),$$

(10)

where $d$ is the side-length of each pixel. In the integral, $\Sigma(\bm{\rho})$ represents one instance of a random complex electric field profile. Writing $\Sigma(\bm{\rho}_{i})\equiv\Sigma_{i}$, the action of this functional integral is characterized by the formula

$$\int d\Sigma\text{ }f(\Sigma_{1},...,\Sigma_{n})=\int d^{2}\Sigma_{1}...d^{2}\Sigma_{n}\text{ }\frac{1}{(2\pi)^{n}\sqrt{|\bm{\Gamma}|}}e^{-\frac{1}{2}(\bm{r}-\bm{\mu})^{T}\bm{\Gamma}^{-1}(\bm{r}-\bm{\mu})}f(\Sigma_{1},...,\Sigma_{n}),$$

(11)

where $\bm{r}\equiv\Big{(}\text{Re}[\Sigma_{1}],\text{Im}[\Sigma_{1}],...,\text{Re}[\Sigma_{n}],\text{Im}[\Sigma_{n}]\Big{)}$, $\bm{\mu}=\langle\bm{r}\rangle$, $\bm{\Gamma}$ is the covariance matrix of $\bm{r}$, and $|\cdot|$ represents the determinant operation. Note that, in the case of thermal statistics, $\bm{\mu}=\bm{0}$. To complete our description of unpolarized multimode thermal light, we now determine the covariance matrix $\bm{\Gamma}$. It is easy to see that $\bm{\Gamma}$ is completely determined by $\langle\Sigma(\bm{\rho}_{1})\Sigma(\bm{\rho}_{2})\rangle$ and $\langle\Sigma^{*}(\bm{\rho}_{1})\Sigma(\bm{\rho}_{2})\rangle$. The term $\langle\Sigma(\bm{\rho}_{1})\Sigma(\bm{\rho}_{2})\rangle$ will always be $0$ in the case of thermal light, and $\langle\Sigma^{*}(\bm{\rho}_{1})\Sigma(\bm{\rho}_{2})\rangle$ will have the form

$$\langle\Sigma^{*}(\bm{\rho}_{1})\Sigma(\bm{\rho}_{2})\rangle=\sqrt{\bar{n}(\bm{\rho}_{1})\bar{n}(\bm{\rho}_{2})}\left[\frac{1}{\pi\sigma}e^{-\frac{|\bm{\rho}_{1}-\bm{\rho}_{2}|^{2}}{\sigma}}\right],$$

(12)

where $\sigma$ is assumed to be small so that $\frac{1}{\pi\sigma}e^{-\frac{|\bm{\rho}_{1}-\bm{\rho}_{2}|^{2}}{\sigma}}\approx\delta(\bm{\rho}_{1}-\bm{\rho}_{2})$. Normalization gives us the mean photon number at position $\bm{\rho}$ as $\bar{n}(\bm{\rho})=\pi\sigma/d^{2}$.

The temporal evolution of a photon’s spatial probability distribution obeys classical physics. This behavior is a direct consequence of the free-space Hamiltonian for the electromagnetic field being quadratic in its quadrature variables [49, 41]. Assuming a paraxial light source, we can utilize the Fresnel diffraction formula to determine the propagated mode-structure of the unpolarized multimode thermal light [50]. The Fresnel kernel at propagation distance $z$ is given by [39]

$$K(\bm{r},\bm{\rho},z)=\frac{e^{ikz}}{i\lambda z}e^{\frac{ik}{2z}(\bm{r}-\bm{\rho})},$$

(13)

where $\lambda$ is the wavelength of the light source, $k=\frac{2\pi}{\lambda}$, and $\bm{r},\bm{\rho}$ represent positions in the measurement plane and source plane respectively. We can calculate the propagated mode structure as

$$E^{(+)}(\bm{r},z)=\int d\bm{\rho}K(\bm{r},\bm{\rho},z)E^{(+)}(\bm{r},0).$$

(14)

Given a mode $\hat{a}_{0}=\int d\bm{r}\text{ }f(\bm{r})\hat{a}_{i}(\bm{r})$ with arbitrary polarization $i$, the resulting mode in the far-field is given by

$$\hat{a}_{z}=\int d\bm{r}\left[\int d\bm{\rho}\text{ }K^{*}(\bm{r},\bm{\rho},z)f(\bm{\rho})\right]\hat{a}_{i}(\bm{r}).$$

(15)

Importantly, we note that the operator-valued distribution $\hat{a}_{i}(\bm{r})$ is not integrated against the Fresnel kernel.

In our experiment, a linear polarization grating with a position-dependent polarization angle filters the multimode unpolarized light, and the beam is propagated in free-space [9]. In the far-field, we make measurements with two point-detectors which are able to post-select on a particular configuration of polarizations. To be explicit, the operator representing a measurement of the first-order correlation at position $\bm{\rho}$ is given by $\hat{a}^{\dagger}_{i}(\bm{\rho})\hat{a}_{j}(\bm{\rho})$ where $i,j\in\{H,V\}$. We represent the transformation of the polarization grating with the matrix

$$\bm{P}(x)=\begin{pmatrix}\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)\\ \sin\left(\frac{\pi x}{L}\right)&\cos\left(\frac{\pi x}{L}\right)\end{pmatrix},$$

(16)

where $L$ is the width of the polarization grating and the represented transformation is given by

$$\hat{a}_{B}(\bm{\rho})=P_{HB}(\rho_{x})\hat{a}_{H}(\bm{\rho})+P_{VB}(\rho_{x})\hat{a}_{V}(\bm{\rho})+P_{\emptyset B}(\rho_{x})\hat{a}_{\emptyset}(\bm{\rho}),$$

(17)

for $B\in\{H,V\}$ and $\bm{\rho}=\rho_{x}\hat{\bm{x}}+\rho_{y}\hat{\bm{y}}$. The $\hat{a}_{\emptyset}(\bm{\rho})$ mode represents photon-loss at the polarization grating. Therefore, immediately after the linear polarizer, the mode structure is given by

$$\hat{a}_{\Sigma,B,\bm{s}}=\int d\bm{\rho}\text{ Rect}\left[(\bm{s}-\bm{\rho})/d\right]\Sigma(\bm{\rho})\Big{[}P_{HB}(\rho_{x})\hat{a}_{H}(\bm{\rho})+P_{VB}(\rho_{x})\hat{a}_{V}(\bm{\rho})+P_{\emptyset B}(\rho_{x})\hat{a}_{\emptyset}(\bm{\rho})\Big{]}.$$

(18)

Consequently, in the far-field, it is given by

$$\hat{a}_{\Sigma,B,\bm{s},z}=\int d\bm{r}\text{ }\left(\int d\bm{\rho}K^{*}(\bm{r},\bm{\rho},z)\text{ Rect}\left[(\bm{s}-\bm{\rho})/d\right]\Sigma(\bm{\rho})\Big{[}P_{HB}(\rho_{x})\hat{a}_{H}(\bm{r})+P_{VB}(\rho_{x})\hat{a}_{V}(\bm{r})+P_{\emptyset B}(\rho_{x})\hat{a}_{\emptyset}(\bm{r})\Big{]}\right).$$

(19)

We now make a couple of approximations to simplify our calculations. First, we suppose that $\text{ Rect}\left[(\bm{s}-\bm{\rho})/d\right]\approx\text{ Rect}\left[\bm{\rho}/L\right]$ for all pixel positions $\bm{s}$. Then, we assume that light at each position in the far-field had originated primarily from a single pixel. With these approximations, we are able to write the propagated state in the form

$$\hat{\rho}_{z}=\int d\Sigma\bigotimes_{\bm{r}}\Big{(}|\alpha\rangle\langle\alpha|_{\Sigma,H,\bm{r},z}+|\alpha\rangle\langle\alpha|_{\Sigma,V,\bm{r},z}\Big{)},$$

(20)

where the mode structure for each $\bm{r}$ is now given by

	$\displaystyle\hat{a}_{\Sigma,B,\bm{r},z}\approx\int d\bm{r}^{\prime}$	$\displaystyle\text{ Rect}\left[\frac{(\bm{r}-\bm{r}^{\prime})}{d^{\prime}}\right]\Bigg{(}\int d\bm{\rho}K^{*}(\bm{r}^{\prime},\bm{\rho},z)\text{ Rect}\left[\frac{\bm{\rho}}{d}\right]\Sigma(\bm{\rho})$		(21)
		$\displaystyle\times\Big{[}P_{HB}(\rho_{x})\hat{a}_{H}(\bm{r})+P_{VB}(\rho_{x})\hat{a}_{V}(\bm{r})+P_{\emptyset B}(\rho_{x})\hat{a}_{\emptyset}(\bm{r})\Big{]}\Bigg{)},$

where $d^{\prime}$ is the width of each pixel in the measurement plane. From here, we can compute the second-order correlation functions for various polarization projections as

	$\displaystyle G^{(2)}_{ijkl}(\bm{r}_{1},\bm{r}_{2},z)$	$\displaystyle=\int d\Sigma\text{ }\sum_{A,B}\Big{(}\langle\alpha|_{\Sigma,A,\bm{r}_{1},z}\langle\alpha|_{\Sigma,B,\bm{r}_{2},z}\Big{)}\hat{a}^{\dagger}_{i}(\bm{r}_{1})\hat{a}^{\dagger}_{j}(\bm{r}_{2})\hat{a}_{k}(\bm{r}_{1})\hat{a}_{l}(\bm{r}_{2})\Big{(}|\alpha\rangle_{\Sigma,A,\bm{r}_{1},z}|\alpha\rangle_{\Sigma,B,\bm{r}_{2},z}\Big{)}$		(22)
		$\displaystyle=\int d\Sigma\text{ }|\alpha|^{4}\int d\bm{\rho}_{1}d\bm{\rho}_{2}d\bm{\rho}_{3}d\bm{\rho}_{4}K^{*}(\bm{r}_{1},\bm{\rho}_{1},z)K^{*}(\bm{r}_{2},\bm{\rho}_{2},z)$
		$\displaystyle\text{ }\text{ }\text{ }\times K(\bm{r}_{1},\bm{\rho}_{3},z)K(\bm{r}_{2},\bm{\rho}_{4},z)\Sigma(\bm{\rho}_{1})\Sigma(\bm{\rho}_{2})\Sigma^{*}(\bm{\rho}_{3})\Sigma^{*}(\bm{\rho}_{4})$
		$\displaystyle\text{ }\text{ }\text{ }\times\sum_{A,B}P_{iA}(\rho_{1x})P_{jB}(\rho_{2x})P_{kA}(\rho_{3x})P_{lB}(\rho_{4x})$
		$\displaystyle\approx|\alpha|^{4}\int d\bm{\rho}_{1}d\bm{\rho}_{2}d\bm{\rho}_{3}d\bm{\rho}_{4}K^{*}(\bm{r}_{1},\bm{\rho}_{1},z)K^{*}(\bm{r}_{2},\bm{\rho}_{2},z)K(\bm{r}_{1},\bm{\rho}_{3},z)K(\bm{r}_{2},\bm{\rho}_{4},z)$
		$\displaystyle\text{ }\text{ }\text{ }\times\frac{\pi^{2}\sigma^{2}}{L^{4}}\text{Rect}(\frac{\bm{\rho}_{1}}{L})\text{Rect}(\frac{\bm{\rho}_{2}}{L})\text{Rect}(\frac{\bm{\rho}_{3}}{L})\text{Rect}(\frac{\bm{\rho}_{4}}{L})$
		$\displaystyle\text{ }\text{ }\text{ }\times\Big{[}\delta(\bm{\rho}_{1}-\bm{\rho}_{3})\delta(\bm{\rho}_{2}-\bm{\rho}_{4})+\delta(\bm{\rho}_{1}-\bm{\rho}_{4})\delta(\bm{\rho}_{2}-\bm{\rho}_{3})\Big{]}$
		$\displaystyle\text{ }\text{ }\text{ }\times\frac{1}{4}\Big{[}P_{iH}(\rho_{1x})P_{jH}(\rho_{2x})P_{kH}(\rho_{3x})P_{lH}(\rho_{4x})+P_{iH}(\rho_{1x})P_{jH}(\rho_{2x})P_{kV}(\rho_{3x})P_{lV}(\rho_{4x})$
		$\displaystyle\text{ }\text{ }\text{ }+P_{iV}(\rho_{1x})P_{jV}(\rho_{2x})P_{kH}(\rho_{3x})P_{lH}(\rho_{4x})+P_{iV}(\rho_{1x})P_{jV}(\rho_{2x})P_{kV}(\rho_{3x})P_{lV}(\rho_{4x})\Big{]}$
		$\displaystyle\equiv I_{0}\int d\bm{\rho}_{1}d\bm{\rho}_{2}d\bm{\rho}_{3}d\bm{\rho}_{4}\text{ }F(\bm{r}_{1},\bm{r}_{2},\bm{\rho}_{1},\bm{\rho}_{2},\bm{\rho}_{3},\bm{\rho}_{4},z)$
		$\displaystyle\text{ }\text{ }\text{ }\times\Big{[}\delta(\bm{\rho}_{1}-\bm{\rho}_{3})\delta(\bm{\rho}_{2}-\bm{\rho}_{4})+\delta(\bm{\rho}_{1}-\bm{\rho}_{4})\delta(\bm{\rho}_{2}-\bm{\rho}_{3})\Big{]},$

where we have defined $I_{0}=\pi^{2}\sigma^{2}|\alpha|^{4}/L^{4}$ and

	$\displaystyle F(\bm{r}_{1},\bm{r}_{2},\bm{\rho}_{1},\bm{\rho}_{2},\bm{\rho}_{3},\bm{\rho}_{4},z)=$	$\displaystyle\text{ }K^{*}(\bm{r}_{1},\bm{\rho}_{1},z)K^{*}(\bm{r}_{2},\bm{\rho}_{2},z)K(\bm{r}_{1},\bm{\rho}_{3},z)K(\bm{r}_{2},\bm{\rho}_{4},z)$		(23)
		$\displaystyle\times\text{Rect}(\frac{\bm{\rho}_{1}}{L})\text{Rect}(\frac{\bm{\rho}_{2}}{L})\text{Rect}(\frac{\bm{\rho}_{3}}{L})\text{Rect}(\frac{\bm{\rho}_{4}}{L})$
		$\displaystyle\times\frac{1}{4}\Big{[}P_{iH}(\rho_{1x})P_{jH}(\rho_{2x})P_{kH}(\rho_{3x})P_{lH}(\rho_{4x})+P_{iH}(\rho_{1x})P_{jH}(\rho_{2x})P_{kV}(\rho_{3x})P_{lV}(\rho_{4x})$
		$\displaystyle\text{ }\text{ }\text{ }\text{ }\text{ }+P_{iV}(\rho_{1x})P_{jV}(\rho_{2x})P_{kH}(\rho_{3x})P_{lH}(\rho_{4x})+P_{iV}(\rho_{1x})P_{jV}(\rho_{2x})P_{kV}(\rho_{3x})P_{lV}(\rho_{4x})\Big{]}.$

These definitions allow for a drastically simplified $G^{(2)}_{ijkl}(\bm{r}_{1},\bm{r}_{2},z)$, and they are used in the main body of the paper. From here, each $G^{(2)}_{ijkl}(\bm{r}_{1},\bm{r}_{2},z)$ can be calculated explicitly. Furthermore, we can use a similar approach to show that $G^{(1)}_{i,j}(\bm{r},z)=\langle\hat{a}^{\dagger}_{i}(\bm{r})\hat{a}_{j}(\bm{r})\rangle=\sqrt{I_{0}}L/(2z^{2}\lambda^{2})$. Using these, the normalized second-order correlation functions $g^{(2)}_{ijkl}(\bm{r}_{1},\bm{r}_{2},z)$ can be calculated, and this list is presented in the next section. Notably, these results are in agreement with our previous theoretical approach and our experimental data [9].

Here we explicitly write the relevant second-order coherence functions studied in our experiment. In this section, we are using the shorthands $\operatorname{sinc}(\nu)\equiv\sin(\pi\nu)/(\pi\nu)$ and $\nu=L(r_{1x}-r_{2x})/(\lambda z)$.

	$\displaystyle g^{(2)}_{\text{HHHH}}(\nu)=$	$\displaystyle\frac{1}{16}(0\operatorname{sinc}(\nu)^{2}+2(6\operatorname{sinc}(\nu+1)+\operatorname{sinc}(\nu+2)+6\operatorname{sinc}(1-\nu)+\operatorname{sinc}(2-\nu))\operatorname{sinc}(\nu)$		(24)
		$\displaystyle+6\operatorname{sinc}(\nu+1)^{2}+\operatorname{sinc}(\nu+2)^{2}+6\operatorname{sinc}(1-\nu)^{2}+\operatorname{sinc}(2-\nu)^{2}$
		$\displaystyle+4\operatorname{sinc}(\nu+1)\operatorname{sinc}(\nu+2)+4(\operatorname{sinc}(\nu_{1})+\operatorname{sinc}(2-\nu))\operatorname{sinc}(1-\nu)+6),$
	$\displaystyle g^{(2)}_{\text{HVHV}}(\nu)=$	$\displaystyle\frac{1}{16}(2\operatorname{sinc}(\nu)^{2}-2(\operatorname{sinc}(\nu+2)+\operatorname{sinc}(2-\nu))\operatorname{sinc}(\nu)+2(\operatorname{sinc}(1-\nu)$
		$\displaystyle-\operatorname{sinc}(\nu+1))^{2}+\operatorname{sinc}(\nu+2)^{2}+\operatorname{sinc}(2-\nu)^{2}+6),$
	$\displaystyle g^{(2)}_{\text{VHHV}}(\nu)=$	$\displaystyle\frac{1}{16}(6\operatorname{sinc}(\nu)^{2}-2(\operatorname{sinc}(\nu+2)+\operatorname{sinc}(2-\nu))\operatorname{sinc}(\nu)+2(\operatorname{sinc}(1-\nu)$
		$\displaystyle-\operatorname{sinc}(\nu+1))^{2}-\operatorname{sinc}(\nu+2)^{2}-\operatorname{sinc}(2-\nu)^{2}),$
	$\displaystyle g^{(2)}_{\text{HHVV}}(\nu)=$	$\displaystyle\frac{1}{16}(2\operatorname{sinc}(\nu)^{2}-2(\operatorname{sinc}(\nu+2)+\operatorname{sinc}(2-\nu))\operatorname{sinc}(\nu)$
		$\displaystyle+2(\operatorname{sinc}(1-\nu)-\operatorname{sinc}(\nu+1))^{2}+\operatorname{sinc}(\nu+2)^{2}+\operatorname{sinc}(2-\nu)^{2}).$

In this section, we present a method for determining the photon number distribution in different detection planes. In doing so, we can study the dynamics of multiphoton wavepackets. It will be challenging to compute the photon number distribution directly from $\hat{\rho}_{z}$, but we can avoid this difficulty by recognizing that

$$\eta_{AB}(\bm{r})=\int d\bm{\rho}K^{*}(\bm{r},\bm{\rho},z)\text{ Rect}\left[\bm{\rho}/L\right]\Sigma(\bm{\rho})P_{AB}(\rho_{x})$$

(25)

follows Gaussian statistics as a result of $\Sigma(\bm{\rho})$ obeying Gaussian statistics. Each $\eta_{AB}(\bm{r})$ represents one instance of a coherent state, and so by determining the probability distribution of the $\eta_{AB}(\bm{r})$ we can determine the effective quantum state as measured by our detectors. For post-selected polarization $ijkl$ at positions $\bm{r}_{1},\bm{r_{2}}$, we will need the probability distribution for $\eta_{iA}(\bm{r_{1}}),\eta_{jB}(\bm{r_{2}}),\eta_{kC}^{*}(\bm{r_{1}}),\eta_{lD}^{*}(\bm{r_{2}})\equiv\alpha_{iA},\alpha_{jB},\alpha_{kC},\alpha_{lD}$ where $A,B,C,D\in\{H,V\}$. Denoting $\bm{t}=(\text{Re}[\alpha_{iA}],\text{Im}[\alpha_{iA}],...,\text{Re}[\alpha_{lD}],\text{Im}[\alpha_{lD}])$, the desired probability distribution is given by

$$P_{iAjBkClD}(\bm{t})=\frac{1}{(2\pi)^{4}\sqrt{|\bm{\Gamma}|}}e^{-\frac{1}{2}(\bm{t}-\bm{\mu})^{T}\bm{\Gamma}^{-1}(\bm{t}-\bm{\mu})},$$

(26)

where $\bm{\mu}=\langle\bm{t}\rangle$ and $\Gamma_{nm}=\langle t_{n}t_{m}\rangle-\langle t_{n}\rangle\langle t_{m}\rangle$. With this, the resulting state describing these statistics is now

$$\hat{\rho}_{iAjBkClD}(z)=\int d^{2}\alpha_{iA}d^{2}\alpha_{jB}d^{2}\alpha_{kC}d^{2}\alpha_{lD}P_{iAjBkClD}(\alpha_{iA},\alpha_{jB},\alpha_{kC},\alpha_{lD})|\alpha_{kC},\alpha_{lD}\rangle\langle\alpha_{iA},\alpha_{jB}|.$$

(27)

It then follows that the total state is given by

$$\hat{\rho}_{ijkl}(z)=\frac{1}{4}\left[\hat{\rho}_{iHjHkHlH}(z)+\hat{\rho}_{iHjHkVlV}(z)+\hat{\rho}_{iVjVkHlH}(z)+\hat{\rho}_{iVjVkVlV}(z)\right],$$

(28)

and thus that the photon-number distribution $p(n_{1},n_{2})$ can be calculated via

$$p(n_{1},n_{2},z)=\text{Tr}\left[\hat{\rho}_{ijkl}(z)|n_{1},n_{2}\rangle\langle n_{1},n_{2}|\right].$$

(29)

In this section, we describe the realization of the polarization grating using polarization optics and a spatial light modulator (SLM) [30]. As shown in Fig. 5, the polarization rotation of the input beam is performed using a SLM in combination with two quarter-wave plates (QWPs). Specifically, the input beam is prepared by passing it through a polarizer aligned to the H polarization. The beam first passes through a QWP at an angle of $45^{\circ}$. Then, this beam is imprinted on the SLM, where a gray-value image is displayed. Finally, the reflected beam passes through another QWP at an angle of $-45^{\circ}$. This configuration provides the ability to rotate the polarization of the incident beam in a controlled fashion. We then characterize the relationship between the polarization rotation and the gray-value displayed on the SLM. This allows us to design a gray-scale image to implement the polarization grating. By adjusting the gray-scale values across different pixels along the $x$-axis of the SLM screen, we can control the polarization at each pixel. We can thus simulate the effect of a polarization grating on an unpolarized light source.