Recently, there has been a growing interest in nonlinear optics [1–3] with intensely structured light [4], which allows for the mixing and manipulation of light beams via light-light interactions. In addition to frequency conversion [5], the generated spatial structure is derived from the product of input modes, going well beyond the linear framework based on the superposition principle of interference [6]. Early pioneering investigations revealed a doubling effect of orbital angular momentum (OAM) [7] of light in the second-harmonic generation [8]. In terms of the multiple spatial degrees of freedom of light, the coupling, transfer, and conservation of angular momentum among light beams have been widely studied in various nonlinear processes [9–17], including spontaneous parametric down-conversion (SPDC) [18,19] and four-wave mixing (FWM) [20].

Entangled photon-pair sources generated through wave mixing lie at the heart of many quantum optical applications [21]. In two-photon optics, the standard observation method shifts from integral exposure to coincidence counting or intensity correlation measurements [22–25]. Spatial entanglement [26] between twin photons has been experimentally observed and certified in high-dimensional momentum spaces [27–32]. In this context, quantum ghost imaging (GI) and interference [33,34] depict techniques for acquiring images (or diffraction patterns) using idler photons, where only the spatially separated twin photons interact with targets. In Klyshko’s advanced-wave picture [35], the formation mechanism of a quantum image is intuitively interpreted using geometric optics: the nonlinear crystal can be analogized as a planar or spherical mirror [36,37], depending on the spatial structure of the pump beam, whether it is a planar wave or a converging spherical wave. Subsequently, Monken et al. demonstrated that the angular spectrum information of the pump light in SPDC can be directly transferred to the spatial profiles of the generated biphoton state [18,38]. As a result, when SPDC meets structured pump light, a batch of quantum optical experiments and schemes is triggered, including OAM entanglement and conservation [27,29], biphoton states engineering [39–44], and even pattern recognition [45].

SPDC is fundamentally a three-wave mixing process from a second-order nonlinear crystal, and the rich nonlocal optical information interactions among the injected pump beam and the generated pair of correlated beams can be aptly described as triple correlations by Abouraddy et al. [46]. Conventionally, such triple correlations were deemed exclusive to nonlinear optical regimes, lacking viable linear analogs, despite the existence of classical counterparts in GI that exploit Hanbury Brown-Twiss (HBT)-type spatial correlations [22,23] in thermal light [34,47–52]. In prior work, we reported generating a class of correlated twin beams in a classical linear system based on the principle of holography [53], which fills the significant gap between biphoton sources and classical thermal sources, exhibiting unique dual properties: it retains the quantum-mimic spatial correlation characteristics while possessing the advantages of classical illumination, such as high brightness, low-cost fabrication, and convenient detection.

Here, we reveal, for the first time, to our knowledge, such triple correlations in the random hologram (RH) source produced by a commercial spatial light modulator (SLM) and low-power continuous-wave light. Within this platform, we demonstrate a variety of nonlocal tripartite optical information processing functionalities, including pattern recognition and matching, intelligent image filtering, and phase-sensitive metrology. Intriguingly, recent research trends show that some nonlinear mapping functions [54–56] can be found in some linear optical configurations equipped with structural nonlinearity, such as recurrent scattering. From a correlation optics perspective, our findings likewise unveil a profound linkage between the intensity correlation in spatially separated linear configurations and wave mixing in nonlinear media. This connection not only provides linear analogs for selecting nonlinear processes but also pioneers novel methodologies for multiparty optical information processing using linear optics.

In correlation optics, the observation is based on two-detector intensity correlation measurement rather than single-detector intensity distribution in conventional optics. In this section, we will generalize intensity correlation functions G(2) [24] in different classes of light sources (Section 2.A) and then illustrate in principle how to use the triple correlation in RH to execute various optical tasks using two examples (Section 2.B). The RH source is inherently a type of incoherent light and can thus be described by the well-known van Cittert-Zernike theorem by introducing the two-field amplitude correlation defined in Ref. [53]. This work follows the theoretical framework established in Refs. [25,46] and Fourier optics [57]. We assume throughout a planar source and a 1D geometry in the transverse plane for the sake of simplicity but without loss of generality. The dimension of interest throughout this paper is the second-order spatial coherence of light or spatial correlations [34], so the temporal dimension is not considered here. Besides, all physical models proposed here can be well interpreted by the modified Klyshko’s advanced-wave picture shown in Section 5.A, which intuitively maps the results of two-detector joint correlation detection to those of traditional single-detector intensity observations.

The two illustrations in Fig. 1(a) show simplified beam-emission models, with the left one typically depicting a spatially incoherent thermal light source made up of a large number of independent emitters stochastically emitting uncorrelated photons or secondary waves. In an HBT intensity interferometer [22,23] with thermal light, a beam splitter is usually used to split the thermal light into two copies. Intensity correlation measurements are performed using two scanning point-like detectors, yielding G(2)(x1,x2) with a width inversely proportional to the source size. Because true thermal light has an ultra-short coherence time on the order of sub-nanoseconds [58], pseudothermal light with a long coherence time is made in the lab by illuminating time-varying random phase masks, such as rotating ground glass [48,59] and SLMs [50,51], with the incident coherent beam Ep(x). Now inserting two linear optical systems, possibly including lenses and objects, into the two arms of the HBT interferometer, the correlation function can be yielded as GHBT(2)(x1,x2)∝|∫dxIp(x)h1*(x1,x)h2(x2,x)|2,(1)where Ip=EpEp* represents the intensity profile of the incident “pump” light at the entrance to the phase modulator, and hj(xj,x) represents the impulse response function for propagating from a point x on the source to a point xj (j=1,2) in the detector plane.

When it comes to the SPDC process in Fig. 1(b), an intense laser beam (pump) is impinged onto a nonlinear crystal with quadratic nonlinear susceptibility, and only a very small number of pump photons disintegrate into pairs of photons that can be spatially separated by some distinguishing features, known as signal and idler photons. Such biphoton states are very close to the EPR state conceived by Einstein et al. in 1935 [26]. Considering that the signal and idler photons pass through the same two linear optical systems as in Eq. (1), the coincidence rate of biphotons at the two detectors is proportional to [25,46] GSPDC(2)(x1,x2)∝|∫dxEp(x)h1(x1,x)h2(x2,x)|2.(2)According to Eq. (2), the spatial profile of the pump beam is directly encoded into the biphoton amplitude known as angular spectrum transfer [38]; however, in Eq. (1), the field degenerates into intensity, losing the most critical part—phase, thus being unable to flexibly customize classical correlations like in biphoton states.

The protagonist of this paper is newly found classically correlated beams produced by using a coherent beam to illuminate varying RHs [53], which exhibit EPR-like classical correlations. The RH is the interference pattern between a tilted plane wave and a chaotic wavefront and can be viewed as an array of independent sub-gratings with the same period but random initial phases; when the RH is illuminated by a coherent beam, each sub-grating therein emits mainly three diffracted beams (±1 and 0 orders) as illustrated in the right panel of Fig. 1(a), which is similar to Ref. [60]. Upon propagation, the yielded speckle pattern in Fig. 1(c) is the result of the superposition of all illuminated sub-gratings. The RH source behaves similarly to the SPDC process, at least in terms of geometrical optics [Fig. 1(b)]. Now the correlation function in Fig. 1(d) between the ±1-order diffracted beams can be expressed as GRH(2)(x1,x2)∝|∫dxEp2(x)h1(x1,x)h2(x2,x)|2.(3)Compared to Eqs. (1) and (2), the quadratic field of the incident “pump” beam can now be transferred into the spatial correlations. In other words, by subtly manipulating the complex amplitude of the pump field, it is promising to say that the correlation function of arbitrary shapes can be customized. Such a nontrivial transfer mechanism of the angular spectrum is more like twin beams produced in FWM [61–63] than in SPDC. So, of great interest here is the design of the structured pump light, and it is certainly possible to spatially shape the incident coherent beam using a linear optical system via Ep(x)=∫dx0E0(x0)h0(x,x0) in Fig. 1(d).

In the previous subsection, we have unveiled the triple correlations among the structured pump and the converted twin beams in Eq. (3), which opens up a new dimension for information encoding beyond the traditional case in Eq. (1). One may alternatively reallocate a target and imaging components in any possible combinations among the three beams, and the lens-demanded imaging equation [53] has already been proven to be identical to two-photon geometric optics [36]. Also, one can place multiple targets and filtering components in the three beams coupled in the intensity correlation via multiplication or convolution, thus showing “nonlocal” or distributed image-processing potentials.

Here, we showcase principally two specific configurations shown in Fig. 1(e) (Case 1) and Fig. 1(f) (Case 2). In Fig. 1(e), two targets, T1(x′) and T2(x′′), are placed in the image plane of the RH source via a 4f imaging system, whereupon a 2f system performs Fourier transform, and two point-scanning or array detectors are used to measure the correlation function between ±1-order diffracted beams. At this time, the impulse response functions of the two symmetric optical systems are hj(xj,x)=Tj(−x)exp(ikfxxj),(4)where k=2π/λ is the wave number, λ is the wavelength, and f is the focal length of the lens. By substituting Eq. (4) into Eq. (3), one can obtain G1(2)(x1,x2)∝|∫dxEp2(x)T1(−x)T2(−x)eikf(x1+x2)x|2.(5)For simplicity, the structured pump is yielded using a plane wave to pass through a target, Tp(xp), followed by a 4f system, so we have Ep=Tp(−x). We assume transmittance functions of the three targets are phase-only, namely, Tj(xj)=exp[iϕj(xj)], and Eq. (5) can be rewritten as $$\begin{gathered} G1(2)(x1,x2)∝F(Tp2T1T2) \\ =|∫dxei(2ϕp+ϕ1+ϕ2)e−ikf(x1+x2)x|2. \end{gathered}$$(6)According to Eq. (6), the three targets placed in three spatially separated beams are coupled into the spatial correlation of the converted twin beams in the form of the Fourier transform (F) of their product, in which the target in the pump beam is amplified to its quadratic field. This theoretical result can be well visualized by Klyshko’s advanced-wave picture in Section 5 [Fig. 7(a)]. Three unrelated phase targets will generally result in a chaotic diffraction pattern. Still, there is a special case when the sum of the phase profiles of the three targets, for as many spatial locations x as possible, is a constant value, or the three phases are well matched, namely, 2ϕp(x)+ϕ1(x)+ϕ2(x)=0,(7)and then a strong correlation peak, proportional to δ(x1+x2), will be generated wherever positions of the two point-like detectors satisfy x1+x2=0. Usually, a detector in one arm is fixed at the midpoint, for example, x1=0. In this sense, by utilizing such tripartite phase-matching conditions in Eq. (7), applications such as phase pattern recognition, matching, and encryption can be flexibly implemented, and corresponding experimental observations will be presented in Section 3.A. Besides, it should be noted that the coherent information mixing of the three targets can only be retrieved by two-detector correlation measurement, whereas any single arm, including the pump beam, cannot. Also, due to the spatial incoherence, the two detectors alone are unable to observe any valid information and usually obtain Gaussian profiles concerning the phase targets.

In Fig. 1(f), the RH source is imaged onto the detector plane with a target, T1(x′), sandwiched between two lenses, while another detector in the other arm is placed in the far field of the RH source via a 2f system (since no target is placed, this arm can be regarded as the idler or reference arm). The impulse response functions of the two sub-optical systems are respectively $$\begin{gathered} h1(x1,x)=T1˜[kf(x+x1)], \\ h2(x2,x)=exp(−ikfxx2), \end{gathered}$$(8)where T1˜=F[T1(x)]. Here we default x1 is fixed at x1=0, and by substituting Eq. (8) into Eq. (3), we have $$\begin{gathered} G2(2)(x2)∝|∫dxEp2(x)T1˜(kx/f)e−ikfxx2|2 \\ =|F[Tp2(−x2)]⊗T1(−x2)|2, \end{gathered}$$(9)where ⊗ denotes the convolution operation. Equation (9) exhibits a Fourier spatial filtering or optical convolution between a filter (placed in the pump beam) and a target (placed in one of the twin beams), and the filtered image occurs in the other beam via intensity correlation. The experimental validation is given in Section 3.B.

Following the two cases presented above, one can, of course, design arbitrary linear optical image-processing systems, featured by the “nonlocal” nature and triple correlations. This paradigm of distributed information processing in two-photon optics is first and foremost unattainable in conventional optics. Secondly, such triple correlations with the rare angular spectrum transfer that mimics the information interaction in FWM do not exist in the traditional model of thermal light. In addition, the two configurations above can be unfolded with Klyshko’s advanced-wave picture, and especially the pump light can also be equivalent to a singular mirror with a modulation function proportional to Ep2. See Section 5.A for detailed descriptions as well. Furthermore, an approach of relative phase measurement (Case 3) based on the triple correlation of the RH source is proposed.

Based on the principle in Section 2, we design several experiments here to illustrate how the manipulated triple correlations in the RH source can be harnessed to demonstrate pattern recognition, image filtering, and optical metrology. We will see that some applications with entangled photon-pair sources can now be implemented using the classical RH source rather than HBT-type thermal light. Also, some brief discussions and analyses on experimental results are provided in this section. Of course, quantum resources may have the advantages of few photons and high sensitivity, but our scheme also provides a practical and cost-effective alternative as it fills the gap under bright beams.

Section 2.A gives a brief introduction to EPR-type correlation in the RH source. In this work, a commercial liquid-crystal SLM (Holoeye-VIS-016, the pixel pitch is 8 μm) is imprinted with varying RHs, and a pair of ±1-order correlated chaotic beams is generated using the He-Ne laser light (at 632.8 nm) to illuminate onto the SLM. As for these RHs used here, each sub-grating is composed of 4×4 discrete pixels on the reflective SLM, and the period d of the sinusoidal grating pattern is set to 36 μm, so the separation angle α between ±1 orders is approximately 2° according to sinα=2λ/d. Two CMOS-based image sensors are used to measure the normalized intensity correlation function with zero time delay: g(2)(r1,r2)=⟨I1(r1)I2(r2)⟩⟨I1(r1)⟩⟨I2(r2)⟩,(10)where Ij(rj) represents the intensity profile of a spatial location rj≡(xj,yj) in one beam, and ⟨·⟩ stands for ensemble average over a large number of realizations. Usually, r1 remains fixed, meaning that only one pixel (midpoint) in the detector is activated, while another detector selects all pixels in a region of interest. Now g(2)(r1,r2) can be rewritten as g(2)(r2), for simplicity, with a single independent variable r2, and g(2)(0) denotes the normalized second-order correlation coefficient. Also see Section 5.D for details of the experimental data processing.

Figure 2(a) shows the experimental setup following Fig. 1(e). The polarized He-Ne laser beam is expanded through a telescopic system and then passes through a half-wave plate (not shown in the setup) followed by a transmissive SLM0 (Daheng Optics, GCI-770104; the pixel pitch is 12.5 μm), which spatially reshapes the incident wavefront onto the reflective SLM plane through a two-lens-based (f=30cm) imaging system. The accurate phase-only modulation in this paper is realized by imprinting computer-generated holograms (CGHs) [64] onto the SLMs along with a spatial filter (i.e., iris) that allows only the ±1-order diffracted beam to pass through. See Section 5.D for details. The reflective SLM in Fig. 2(a) is imprinted with varying computer-generated RHs at a refreshing rate of 10 Hz. The generated ±1-order diffracted beams are respectively imaged onto SLM1 and SLM2 (Daheng Optics, GCI-770104) with a two-lens-based (f=30cm) imaging system, and the two transmissive SLMs are also optionally imprinted with CGHs to realize arbitrary phase-only modulation. Then, two lenses with a focal length of 60 cm focus the two beams onto two detectors for joint correlation measurement of g(2)(r2) mentioned above. The fixed specific phase patterns loaded on SLM0, SLM1, and SLM2 are denoted as ϕp, ϕ1, and ϕ2, respectively, the same as Eq. (7).

The tripartite phase relation in Eq. (7) actually reveals the correlation generation rule, whereby we present two practical optical applications for phase pattern recognition and correlation manipulation. We consider the human face patterns for pattern recognition, as shown in Fig. 2(b), and a total of seven human faces {Φ} and their negative films {−Φ} are selected. For clarity, we utilize seven different colors and capital letters (A–G) to represent the seven human faces and the corresponding lowercase letters (a–g), as well as their complementary colors to represent their negative patterns. The notations A and a in Fig. 2(c) represent {Φ/2} and {−Φ/2}, respectively. The process of pattern recognition is to scan the pattern library and compare them with the input unclassified pattern, and optical correlations enable this process to be realized in a nonlocal fashion, implying that the two detectors can be located far apart; then, the unknown pattern can be identified by quantitatively comparing the magnitude of the correlation coefficient g(2)(0) obtained from the joint correlation measurement over 4000 realizations.

There are two manners to perform nonlocal optical recognition. The first, akin to Ref. [65], is to insert the two phase patterns in the two arms of ±1-order diffracted beams in Fig. 2(a), and the pump beam is now treated as a plane wave (ϕp=0), namely, a fixed blank pattern is imprinted onto SLM0. Each time a human-face pattern ϕ1 to be recognized is imprinted onto SLM1, 14 phase patterns ϕ2 (i.e., A–G and a–g) from the library are sequentially loaded onto SLM2 for joint correlation measurements of g(2)(r2;ϕ1,ϕ2), where a pair of conditional variables, ϕ1 and ϕ2, is added. In this way, by switching the patterns on the two SLMs, we perform a total of 16×16 correlation measurements and combine these functions into a 16×16 array in Fig. 2(d). In this case, whenever ϕ1+ϕ2=0, a correlation peak occurs (i.e., the principal diagonal of the correlation array), whereby an unknown pattern can be inferred as the negative film of the reference pattern that has the maximum correlation coefficient g(2)(0). We extract g(2)(0) from each cell in the correlation array and form a 2D confusion matrix in Fig. 2(f) (left panel), from which we can see that the diagonal elements are generally much larger than the off-diagonal elements, indicating a successful job with all targets recognized.

The second manner is to place the human-face patterns and the reference patterns in the pump beam (SLM0) and the signal beam (SLM1), respectively, and a fixed blank phase is loaded onto SLM2 (ϕ2=0). The seven human faces (a–g) to be recognized are imprinted onto SLM0 in the pump beam, and the 14 reference phases (A–G and a–g) in the library are imprinted onto SLM1; at that time, likewise, we obtain a 7×14 correlation array of g(2)(r2;ϕp,ϕ1) in Fig. 2(e) and the corresponding confusion matrix in the right panel of Fig. 2(f). As expected in Eq. (7), the correlation peak stands out whenever ϕ1+2ϕp=0, and the unknown pattern can now be inferred as half of the negative film of the reference pattern that receives the maximum correlation coefficient. It is worth noting that the proposed manner is similar to structured-pump-enabled quantum human face recognition [45] in the SPDC light, and the target pattern is placed in the pump beam, the reference phase patterns in the database are scanned sequentially in the signal beam, and the correlation signal is generated in the idler arm through joint detection.

We furthermore introduce a more general application, named correlation manipulation, to induce correlation generation between two arbitrary and uncorrelated phase modes. Actually, when the three phases are helical phases involving OAM, we have already observed the conservation of OAM in Ref. [53]. Similarly in Section 5.B, we observe the correlation matrix of radial momentum [31] and its conservation law. Back to Eq. (7), the phase-matching relation consists of three variables, namely, having two degrees of freedom; thus for two arbitrary phases, ϕ1 and ϕ2, placed in the two correlated beams of Fig. 2(a), we only need to reshape the pump beam’s wavefront to a mixed spatial profile, ϕp=−(ϕ1+ϕ2)/2. We experimentally customize the matched correlations between seven human-face patterns and seven spiral phases (OAM modes with the topological charge ℓ from −3 to 3) in Fig. 3(a) as well as seven random-grid patterns in Fig. 3(b). For instance in Fig. 3(a), the human-face set (ϕ1 in SLM1) and the helical phase set (ϕ2 in SLM2) are paired one by one to produce a mixed set consisting of another seven patterns (ϕp in SLM0) with Eq. (7). For each fixed structured pump, we can obtain a 7×7 correlation array akin to Figs. 2(d) and 2(e), so seven such phases ϕp will yield a 3D correlation cube (7×7×7) of g(2)(r2;ϕp,ϕ1,ϕ2). Figure 3(c) displays the two 3D confusion matrices corresponding to Figs. 3(a) and 3(b). We notice that the establishment of such matched correlations is so specific that each layer in the correlation cubes only has one correlation peak signal; in other words, the correlation peak signal only appears in one of the body diagonals of the cube. In this sense, we manipulate the triple correlation by engineering the spatial profile of the pump beam, so that there is only a unique pair in the 49 combinations in the library that produces a specific matching peak signal.

Of significance is that the triple correlation is completely absent in traditional HBT-type thermal sources. According to Eq. (1), the phase structure of the pump beam cannot be transferred to the correlation function, and the phase-matching relation or correlation generation rule in this case is ϕ1−ϕ2=0 regardless of ϕp, trivially showing a dual correlation. Also, there is a subtle difference between the RH source and the entangled photon-pair source in SPDC, and according to Eq. (2), the phase-matching relation is ϕp+ϕ1+ϕ2=0. However, in FWM, the phase-matching relation is identical to the RH source. Given that these applications are not demonstrated in the twin beams via FWM, our experimental results thus give insightful instructions.

Generally speaking, an image processing system transforms one image into another with enhanced features. One of the most famous optical information processing systems is the Fourier filtering [57,66], which introduces a well-designed mask (i.e., Fourier filter) in the middle of a 4f system (i.e. Fourier-transform or far-field plane). Spiral phase contrast imaging [45,67–70], developed in the last two decades, is a typical application of Fourier filtering in which a spiral phase plate with the topological charge of ±1 is added to the Fourier-transform plane; it has been found that this technology can easily enhance or extract the edge information, particularly in regions with steep phase changes, as well as map phase singularities [71]. Equation. (9) in Section 2.B [Fig. 1(f), Case 2] shows a nonlocal Fourier filtering in principle, and according to the advanced-wave picture [Fig. 7(b) in Section 5], the target T1 is imaged onto the detection plane in the idler arm through a 4f system and especially a spatial filter (proportional to Ep2) in between.

Now we demonstrate nonlocal spiral phase contrast imaging based on the RH source for edge enhancement and vortex mapping. Figure 4(a) shows the experimental setup based on Fig. 1(f), and the spiral filter is realized through SLM0 in the pump beam. It should be noted that, due to the second-harmonic effect (Ep2) of the pump beam, the spiral phase changes from 0 to π (half of 2π) every revolution (ℓ=0.5). The yielded helical wavefront is imaged onto the reflective SLM imprinted with varying RHs, and the produced ±1-order correlated beams are spatially split into two arms by a BS. In one arm, +1-order diffracted beam is imaged onto D1 via a 4f system (the focal length is 30 cm), and in between the target (SLM1) is placed in the far field of the RH source; in the other arm, D2 is directly placed in the far field of the RH source. The intensity correlation is implemented between a single fixed point’s light-intensity fluctuation of D1 and all points’ light-intensity fluctuations of D2, and it can be thought that the correlated pattern naturally occurs in the idler arm.

Figure 4(b) shows an amplitude target (a circular aperture) used for nonlocal edge-enhanced imaging. When the helical wavefront is removed (ℓ=0), Fig. 4(c) shows the target’s correlated image; by contrast, Fig. 4(d) shows the edge-enhanced correlated image when imprinting the helical phase of ℓ=−0.5 onto the SLM0. Figure 4(e) displays a phase target with two phase singularities (surrounded by black dashed boxes), and the phase function is encoded to a hologram that is imprinted onto SLM1. Figure 4(f) shows the correlated pattern when the vortex beam of ℓ=−0.5 drives the RH source, and the positions of the two yielded correlation peaks correspond to the two counterclockwise optical vortices. Through the experimental verification above, the spatial structure of the pump beam is actually manifested as a spatial filter, which performs Fourier filtering on the target in the signal beam and then “transfers” the filtered image to the idler beam through spatial correlation.

The phase-only filter in Fig. 4 is simple and exhibits a certain degree of universality for various types of targets. For specific categories of targets, using specially designed filters can more efficiently achieve some optical tasks, such as pattern recognition. The human face recognition in Fig. 2 requires strict phase matching to achieve pattern recognition. This may not be practical and intelligent, as it requires the database in the phase library to include as many poses of all users as possible, which inevitably leads to a dual reduction in recognition efficiency and performance. Modern pattern recognition technology relies on artificial intelligence driven by neural networks [72], which optimize the framework and parameters of neural networks by being fed by the database (training set); for an unknown target (testing set, not identical to any target in the training set) to be recognized, by directly inputting the target into the neural network, efficient recognition can be performed without the need to compare each target with all the images in the database one by one. It has recently been found that those electrical neural networks can be realistically executed in some optical manners, known as optical neural networks [73], among which passive cascaded diffractive layers [74–77] use the diffraction to accelerate recognition at the speed of light. The recognition here can be regarded as a special optical transformation process, and the cascaded diffraction layers are like a decoder that converts the object beam into some basic orthogonal optical modes, e.g., Gaussian spots at different spatial locations.

In this context, one can configure the target and the decoder separately in any two beams in the triple correlation among pump, signal, and idler beams. This is not uncommon, whether in the HBT-type thermal light [77] or entangled photon pairs [78], placing the target and the decoder in the signal and idler arms to perform diverse nonlocal optical information processing. When it comes to EPR-type thermal light, here we put the single-layer decoder in the pump beam and the target in the signal arm. The experimental apparatus in Fig. 5(a) is basically the same as Fig. 4(a), except for the SLM placed in the pump beam, where a higher-resolution reflective SLM is used for better modulation accuracy. Following the same method in Refs. [76,77], we train the decoder, a single-layer phase mask in Fig. 5(b), to implement mode conversion from two categories of handwritten digits (0 and 1) to two basic modes in Fig. 5(c) (refer to Section 5.C, and the recognition accuracy reaches 99%). Specifically, the beam with a “digit 1”-like shape passing through the decoder will tend to produce a focal point on the left side in Fig. 5(c), while the beam with a “digit 0”-like shape will be on the right side. In this case, we only need to compare the relative energy or intensity of the two given focus positions to complete the inference. Figures 5(d) and 5(e) showcase the experimental results of six different targets (the first column, testing set) sequentially placed in the signal arm (SLM1), and the second column displays the correlated pattern g(2)(r) over 20,000 realizations between a fixed point in D1 and all points in D2. In fact, we are only interested in the relative magnitude of the correlation coefficients corresponding to the positions of the two focus points in Fig. 5(c). For intuitive comparison, the third column in Figs. 5(d) and 5(e) thus shows the results of g(2)(r)·TM, where TM is the output function corresponding to Fig. 5(c). Consistent with the prediction, the two categories will produce a much higher correlation signal at the corresponding position, so all six targets are now successfully classified. Of course, one can, according to the method proposed here, design multi-layer diffractive elements to perform more complex optical classification tasks.

The correlation-enabled quadratic-field effect from Eq. (3) in a sense simulates the second-harmonic generation (except for frequency conversion) in a classical linear system, thus leading to some novel applications previously only observed in nonlinear optics. The so-called quadratic field naturally and primarily results in a doubling of the phase, reflected in the phase-matching condition in Eq. (7). With this motivation, we further demonstrate a linear phase amplifier akin to the nonlinear version via harmonic generation [79–81] or a multi-photon entangled state (NOON state) [82,83].

As proof of conceptualization, we measure the relative phase as an example. Considering a double slit with a phase difference Δφ (to be measured) between two single slits [Fig. 6(a); the slit width ws and slit spacing ds are respectively 162.5 μm and 475 μm], we can yield its Fourier-plane expression based on the 2f system and the plane-wave coherent illumination: Uf(xf)=A0iλf2wssinc(xfwsλf)cos(πxfdsλf−Δφ2),(11)where A0 is a complex constant value. Since ws≪πds, the intensity envelope of sinc2(xws/λf) has little effect on the position of the main maximum fringe, namely, the peak of the center stripe in the double slit’s interference pattern. As Δφ varies, the fringe pattern shifts periodically with the dependence of cos2(πxfds/λf−Δφ/2) as shown in Fig. 6(b). We define the position of the center maximum value of the interference fringes of Δφ=0 as the zero point, so for an arbitrary phase shift Δφ, the transverse offset Δx with respect to the zero point can be obtained as Δx=λf2πds·Δφ+b=kCΔφ+b,(12)where b=0 for 0≤Δφ<π and b=−λf/ds for π≤Δφ<2π. The linear relationship in Eq. (12) cleverly converts the relative phase difference that cannot be directly measured to the light-intensity offset in the transverse position, and this double-slit-based apparatus is also commonly used to measure spatial coherence by fringe visibility [84,85]. In this sense, the resolution of this measurement system depends on the smallest distinguishable displacement, and the sensitivity of this metrology system can thus be well represented by the slope kC.

Owing to the phase doubling effect through joint correlation measurement, intriguingly, the shift Δx in the correlation function caused by Δφ will become twice that directly obtained under coherent illumination, thus yielding Δx=λfπds·(Δφmodπ)+b=kE(Δφmodπ)+b,(13)where the function “Δφmodπ” returns the remainder of Δφ divided by π; b=0 for (Δφmodπ)<π/2 and b=−λf/ds for (Δφmodπ)≥π/2; and apparently the sensitivity kE is twice of kC.

The complex amplitude function of the double slit is encoded into a hologram and imprinted onto SLM0 in Fig. 6(a). The double slit is imaged onto the RH plane (SLM1) as a structured wavefront via a 4f system (the focal length is 30 cm), and two detectors are placed in the far field of the RH source via the 2f system. Also, see Fig. 7(c) for the unfolded equivalent diagram of this apparatus in Section 5. In the experimental observation, we continuously adjust Δφ between the two slits. For each Δφ to be measured, the stripe patterns of the three cases in Fig. 6(b) are obtained for a fair comparison, and the transverse shift of the center stripe is calculated with respect to Δφ=0. For the case of direct coherent detection, one fixed grating is loaded on SLM1, and the far-field diffraction patterns are directly captured in the single idler arm; for the case of indirect correlation measurement (8000 realizations), SLM1 functions as the RH source, and the HBT-type (between two identical +1-orders) or EPR-type (between ±1-orders) correlations are realized by selecting different diffraction orders of interest. Figure 6(b) displays experimental results of the three light sources (HBT-type thermal light, coherent beam, and EPR-type thermal light), including spatial profiles and their cross-sections, of three different values of Δφ=0, π/5, 2π/5. Figure 6(c) plots three curves of Δx with the change of Δφ of the three cases, and through linear fitting, the measured kC and kE are 68.2 μm/rad and 128.6 μm/rad, respectively, which generally satisfies Eqs. (12) and (13) with kE≈2kC and therefore demonstrates a doubling of sensitivity. Especially in the HBT-type correlation measurement, the diffraction pattern does not shift with Δφ (kH=0); that is, only the intensity profile, regardless of phase, can be embedded into the correlation function as previously highlighted in Eq. (1), and thus the promising turbulence-free or phase-insensitive nature can be embodied in such an HBT-type intensity interferometer [23,86]. It is not wise to declare whether HBT-type or EPR-type thermal light is generally better, because it depends on the specific purpose, such as phase super-resolution measurement or phase-insensitive intensity-profile imaging in a turbulent or scattering environment. Whatever the objective, due to our apparatus’s flexible tunability between the two extremes of optical correlations, the RH source is expected to play an important role in a variety of correlation-based applications and metrology. Since the second-order intensity correlation in EPR-type thermal light can achieve a two-fold resolution gain for relative-phase measurement, it is possible to apply the N-order intensity correlation to reach N-fold resolution gain, similar to high-order harmonic generation [79,80].

In the previous section, diverse applications based on triple correlations in the RH source have been experimentally demonstrated in three devices. The underlying physics can be concisely expressed in Eq. (3) and especially in Eq. (7), where the tripartite phase-matching relationship highly simulates the angular spectrum dependence of parametric four-wave mixing. These experimental results confirm that distributed classical linear optical systems can simulate the unique image transmission process in quantum and nonlinear optical systems by employing the two-detector intensity correlation. Compared to the twin beams produced in SPDC, the main differences lie in the background term of the thermal-fluctuation nature and the quadratic-field dependence of the incident “pump” beam. As a result, the merit of these quantum and nonlinear optical image processing schemes is the high sensitivity at the single-photon level, though constrained by the detection signal-to-noise ratio and cost. By contrast, the feature of our work is that identical image processing tasks can be achieved under high-brightness classical illumination, accompanied by a high detection signal-to-noise ratio and high cost-effectiveness. The trade-off and complementarity between classical and quantum, nonlinear and linear, will more comprehensively guide the implementation of these distributed optical information processing tasks tailored to specific needs and scenarios.

Among these applications, it is noticed that the number of realizations (or the number of acquiring frames) required to yield a visible correlation pattern varies significantly among different cases, ranging from thousands of realizations in pattern recognition (Figs. 2 and 3) to hundreds of thousands of realizations in correlated image filtering (Fig. 4). The reasons for these differences depend on the complexity of the target on the one hand and on the processing task itself on the other. For instance, the edge imaging task in Fig. 4 requires obtaining a complete image and is, therefore, often very complex and sensitive to noise, whereas, in pattern recognition (Figs. 2, 3, and 5), it is only necessary to detect the correlation peak in a specified area. It is found that only a few hundred realizations are enough to ensure good recognition accuracy. Hence, these experimental differences show that the proposed method is more suitable for image-free recognition and matching tasks. Besides, the main limitation of the proposed method lies in the prerequisite ensemble averaging, which is usually time-consuming—it usually takes several minutes to measure the intensity correlation function for the recognition task. This limitation can be completely circumvented by replacing the SLM with a digital micromirror device, which has a much higher pattern-switching speed. On the detection side, especially when only the correlation peak signal [e.g., g(2)(0) in Figs. 2(f), 3(c), and 5(d)] is of interest, the two cameras can be replaced by two point-like photodetectors with the much faster response time.

To sum up, we introduce the angular spectrum dependence in the RH source using three experimental cases and show how to manipulate the triple spatial correlations with the structured pump for nonlocal optical tasks such as pattern recognition, correlated image filtering, and phase-sensitive metrology. Except for the time-frequency component, the intensity correlation in the proposed linear system successfully simulates optical interactions with spatial degrees of freedom in multi-wave mixing. This study enables implementing nonlinear coherent optical information processing, which previously relied solely on second- or third-order nonlinearities (e.g., SPDC, four-wave mixing, and second-harmonics generation), in a simple and linear system with optical correlation. In addition to the spatial degree of freedom, we are investigating high-dimensional correlations in the wavelength and/or polarization degrees of freedom of light in the RH source, similar to the Type-II and/or non-degenerate SPDC processes, to perform all-sided nonlinear information interactions and processing in the context of correlation optics.

In this subsection, we will introduce Klyshko’s advanced-wave picture and how to use it to illustrate the results obtained from the intensity correlation in Section 3 by means of geometrical optics. The key insight is that the joint detection of g(2) using Eq. (10) between the two detectors is, in fact, perfectly equivalent to the observations made by a single detector following some suitable rules. •Replacing one detector with an equivalent light source. In this work, g(2)(r2) is established between a fixed point r1=0 in one detector and all points r2 in another. In this regard, the single-point detector can be replaced by a point-like light source that emits light towards the other detector through the three optical subsystems.•Replacing the RH source with a linear reflective optical element whose modulation function is proportional to Ep2. It should be noted that, when considering other correlated light sources, the alternatives are quite different. For example, the conventional thermal source should be replaced by a phase-conjugate mirror [49], while the nonlinear crystal in the SPDC source should be replaced by a reflective optical element proportional to Ep [36].

Following the two rules, the three setups in the left panels of Fig. 7, corresponding to Section 3, can now be unfolded to diffraction-limited optical systems shown in the right panels of Fig. 7.

In Section 3.A, we have demonstrated the pattern recognition and correlation manipulation using the triple correlations in the RH source, among diverse phase patterns including human-face patterns, OAM-carried helical phase patterns, and random-grid patterns. Here we supplement the phase patterns involving radial momentum pr [31]. It is well known that the azimuthal OAM mode exhibits a continuous phase variation from 0 to 2ℓπ along the angular direction, and similarly, the radial phase mode exhibits a phase variation from 0 to 2prπ along the radius. Figure 8(a) shows the experimental setup, identical to Fig. 2(a), to observe the radial correlation function of g(2)(ppr,p1r,p2r), where the three SLMs (SLM0, SLM1, and SLM2) are used to load the holograms corresponding to the radial phase modes of ppr, p1r, and p2r. Figure 8(b) displays a radial phase mode of pr=−3 (upper) and the corresponding hologram (bottom) imprinted onto the three SLMs. By varying the three radial phase modes, Fig. 8(c) showcases experimental observations of five radial correlation functions of g(2)(p1r,p2r) with the change of ppr in the pump beam, and satisfying experimentally 2ppr=p1r+p2r.(14)In fact, Eq. (14) can be regarded as a radial version of Eq. (7), and the OAM version is 2ℓp=ℓ1+ℓ2 [53].

In Ref. [77], we have demonstrated a type of deep diffractive neural network using the HBT-type thermal light, where the object (encoder) and the decoder are placed separately in the two arms. In the configuration shown in Figs. 1(f) and 5(a) with the RH source, the difference is that we place the phase-only decoder in the pump light, which is invalid for the traditional HBT-type thermal light. According to Klyshko’s advanced-wave picture depicted in Fig. 7(b), Fig. 9(a) shows the equivalent unfolded apparatus for the process of electrical training; this structure is similar to Fourier deep diffractive neural networks [75], except that there is only one diffractive layer to be trained. We define two ideal output modes (TM) in Fig. 5(c); when the amplitude target with the shape of handwritten digit “1” is placed and filtered by the phase mask, a bright Gaussian focus will be produced on the left side of the detection plane. Similarly, for the handwritten digit “0”, a bright focus will be produced on the right side. Hence, the binary classification can be implemented by comparing the relative energy (light intensity) at two positions on the output plane.

The optimization algorithm used here is the wavefront matching algorithm [87], which is also used in Refs. [76,77]. We take 2800 samples for each of the two classes of objects (5600 samples in total), of which 2000 samples are used as the training set and the remaining 800 samples are used as the test set. Figure 9(b) shows the optimized phase pattern of the training set after 200 iterations, and Fig. 9(c) shows the numerical classification results of the two classes of handwritten digits. The four patterns, from left to right, in each panel of Figs. 9(c) and 9(d) are the amplitude target, the complex amplitude in the Fourier-transform plane, the resultant image Iout filtered by the phase mask, and Iout·TM, respectively. Note that the complex amplitude is visualized by means of pseudo-color processing (the HSV colormap), in which the intensity and phase are mapped on the brightness (value) and hue, respectively. The well-trained Fourier filter is used to classify the samples in the test set, as shown in Fig. 9(d), and the recognition accuracy of the binary classification task reaches 99%. It should be noted that the recognition results obtained here come from the single-layer Fourier deep diffractive neural network under coherent illumination [Fig. 9(a)], rather than the nonlocal image filtering or classification in Fig. 5(a). To counteract the phase-doubling effect in the EPR-type thermal light, an extra but compulsory step is to configure the half-phase pattern of Fig. 9(b) inside the pump beam, akin to {Φ/2} in Fig. 2.

Both transmissive and reflective SLMs play a crucial role in this work, including the preparation of the RH source and the implementation of complex amplitude modulation. All patterns loaded on these SLMs are real-valued holograms normalized to [0,1]. For an arbitrary complex amplitude EH, the hologram can be produced via IH=|EH+tgEg|2, where tg is a weighting factor set to the maximum value of |EH|, Eg=exp(ikgx) represents a horizontally inclined plane wave, and the phase profile of Eg shows the structure of the blazed grating. When using a coherent beam to illuminate the SLM imprinted with the IH, the produced +1-order diffracted beam carries the information proportional to EH. By varying the value of kg, one can modulate the separation angle between ±1-order and 0-order diffracted beams. Regarding the RH source generated using the reflective SLM, EH represents a spatially random phase pattern that obeys a uniform distribution in the range of [0,2π], and the value of kg is set to ∼0.0176 to yield the separation angle of ∼2° between ±1-order beams of interest. The equivalent coherence time of the RH source is determined by the switching rate of the reflective SLM, and 10 frames of different computer-generated RHs are played in sequence per second. Compared to the rotating RH used in Ref. [53], the merit of this method is that the halo and periodic repeatability issues caused by rotation can be eliminated. On the other hand, the transmissive SLM is mainly used to achieve the desired modulation and to act as targets, including the human-face phase patterns (Fig. 2), helical phase patterns (Fig. 3), amplitude-only or phase-only targets (Figs. 4 and 5), and complex-valued targets (Fig. 6). In this case, an iris is used to filter out beams other than +1 order, and the value of kg is set to ∼0.0127 to yield the separation angle of ∼0.73° between +1-order and 0-order beams.

During data acquisition, two sets of pairwise CMOS-based cameras (Daheng Imaging; MER-133-54U3M, 960×1280pixels; the pixel pitch is 3.75 μm; MER-231-41U3M, 1200×1920pixels; the pixel pitch is 5.86 μm) are used to perform intensity correlation detection of g(2)(r). The low-resolution ones are used for these experiments in Figs. 2, 3, 6, and 8, while the high-resolution ones are used in Figs. 4 and 5 for a broader field of view. An external signal generator triggers both cameras to acquire the speckle patterns at 10 frames per second. The integration time of the detectors is set at the order of hundreds of microseconds, which is much smaller than the coherence time of the RH source. Regarding the data processing, the pixel values of a fixed midpoint in one set of speckle patterns (detector 1) are correlated, using Eq. (10), with all pixel values of the region of interest in the other set of speckle patterns (detector 2).