The optical coherence structure, as one intrinsic property of a light beam, plays a pivotal role in not only optical interference but also a wide range of applications. The optical coherence theory provides a viable platform for studying partially coherent beams (PCBs)^{1, 2}. Such beams are characterized by a cross-spectral density (CSD) function in the spatial-frequency domain or a mutual coherence function (MCF) in the spatial-time domain. PCBs, with the prescribed optical coherence structures, have various applications in optical coherence tomography^3-5, optical communication⁶, phase imaging^7-9, incoherent source reconstruction¹⁰, and super-resolution optical imaging^{11, 12}; notably, they can be imaged through scatter media¹³. We recently proposed optical coherence imaging, wherein optical coherence structure, as the information carrier, could realize robust far-field optical imaging^14-17. It can be immune to atmospheric turbulence and opaque obstacles. On the other hand, there are multifarious theoretical and experimental methods proposed to generate predetermined PCBs, with two techniques being the most popular. The first operates based on the van Cittert–Zernike theorem^{18, 19}, employing a fully incoherent beam with the desired intensity profile via propagation. The other exploits the mode-superposition principle^20-24, involving the coherent-mode, pseudo-mode (PM), and random-mode (RM) representations. Measuring the complex optical coherence structure is of great interest from both fundamental and application perspectives.

The optical coherence structure has traditionally been measured using Young’s interference experiment, wherein its magnitude and phase can be predicted based on the visibility and position of the fringes, respectively. However, this experiment only considers two light beam points. Full characterization of the optical coherence structure requires that each point be scanned independently across the beam plane, which requires significant time and effort. Improvements to this experiment have been proposed, such as non-parallel slits²⁵, non-redundant arrays of pinholes²⁶, and programmable slits/apertures based on spatial light modulators or digital micromirror devices²⁷. Furthermore, signal interferometers²⁸, wavefront-folding interferometers²⁹, and other interference methods³⁰ have been developed to measure optical coherence structures. These interferometric methods involve a complicated and misalignment- and vibration-sensitive setup and are difficult or infeasible to implement in some spectral regions. These problems can be solved using diffraction methods such as the phase-space approach³¹, diffraction from an opaque obstacle or a transparent mask with a phase discontinuity^{32, 33}, or self-referencing holography^{8, 34, 35}. However, the zero-order approximation of the Taylor series of illuminated PCBs is involved in these methods^31-33. The authors^{34, 35} used a Dirac phase point as the perturbation point. Additionally, other methods such as the Hanbury Brown and Twiss (HBT) effect and generalized HBT effect have been proposed^{36, 37}; they are only valid for beams with Gaussian statistics. Despite all the efforts, measuring both the real and imaginary parts of the optical coherence structure still remains an open and great challenge. Thus, in this letter, we have provided rigorous mathematical solutions and employed the classical optical diffraction experiment–the generalized Arago spot experiment³⁸, wherein the transmittance function of the obstacle is complex, to fully recover the optical coherence structure of PCBs.

Let us consider a scalar, quasi-monochromatic, and statistically stationary random optical field V(r,ω) . The statistical properties of this random optical field are characterized by the CSD function in the spatial-frequency domain, defined as:

1W(r1,r2,ω)=〈V(r1,ω)V*(r2,ω)〉,

where r1≡(x1,y1) and r2≡(x2,y2) are arbitrary spatial position vectors in the source plane; the asterisk and angular brackets denote the complex conjugate and ensemble average of the random electric fields, respectively; and ω is the frequency of the light beam, which is omitted hereafter for brevity. The above equation can be rewritten as^{18, 19}:

2W(r1,r2)=τ(r1)τ*(r2)μ(r1,r2),

where τ(r) and μ(r1,r2) denote the amplitude function and the optical coherence structure, respectively. Our aim is to measure this complex structure; we have considered τ(r)=1 hereafter for brevity.

In our protocol, we adopted the classical optical diffraction experiment— the Arago spot experiment, wherein the beam is obstructed by an obstacle in the source plane and captured in the Fraunhofer zone. We used the Huygens-Fresnel principle to study the beam intensity in the observation plane:

3I(u)=∬μ(r1,r2)O(r1)O*(r2)⋅h(r1,u)h*(r2,u)d2r1d2r2,

where u denotes the position vector in the observation plane; and O(r) and h(r,u) denote the transmittance function of the obstacle and response function of the optical system, respectively. To precisely recover the optical coherence structure μ(r1,r2) , we needed to consider obstacles with three different transmittance functions, namely O1(r)=1−q(r−r0) , O2(r)=exp(iπ)−q(r−r0) , and O3(r)=exp(iπ/2)−q(r−r0) , where q(r) describes the obstacle and is equal to 1 inside the obstacle and zero otherwise, and r0 indicates the center position of the obstacle.

By substituting Oi(r),i=1,2,3 in Eq. (3), we obtained the intensity Ii(u) . We subtracted the intensity with i=1 from that with i=2 and i=3 to obtain:

4ΔI21(u)=I2(u)−I1(u)=4Re[Λ(u)],

and

5ΔI31(u)=I3(u)−I1(u)=2Re[Λ(u)]−2Im[Λ(u)],

respectively, with

6Λ(u)=∬μ(r1,r2)q(r1−r0)h(r1,u)h*(r2,u)d2r1d2r2.

From Eqs. (4–6), we derived:

7Λ(u)=ΔI21/4+i(ΔI21/4−ΔI31/2).

The observation plane was in the Fraunhofer zone, so the response function of the optical system had a Fourier structure; h(r,u)=exp(−i2πu⋅r) , when coordinate scaling was neglected. We substituted this into Eq. (6) and performed a Fourier transform on both sides; Eq. (7) was then reduced to:

8Λ˜(r)=∫μ(r',r+r')q(r'−r0)d2r',

where the tilde denotes the Fourier transform; Eq. (8) is the main result in our protocol. Using the function q(r') and capturing the far-zone intensity thrice by refreshing the obstacles, we can recover the optical coherence structure of any genuine PCB.

To simplify Eq. (8), we classified all the genuine PCBs into two types: Schell-model and non-Schell-model types. Regarding the former, the optical coherence structure only depended on the difference between the two position vectors, namely μ(r',r+r')=μ(r) . Further simplifying this equation, the optical coherence structure could be reconstructed as follows:

9μ(r)=Λ˜(r)/S0,

where S0 denotes the obstacle area. For the non-Schell model type, the optical coherence structure depended on the position of the reference point. Therefore, to accurately measure the optical coherence structure, we adopted a Dirac pinhole as the obstacle, that is, q(r−r0)=δ(r−r0) . By further simplifying Eq. (8), the optical coherence structure could be reconstructed as:

10μ(r0,r+r0)=Λ˜(r).

We concluded from Eqs. (9) and (10) that for all genuine PCBs, the far-zone intensity only needed to be captured thrice and after simple addition and subtraction calculations, the function Λ(u) [see Eqs. (4), (5), and (7)] could be achieved. To recover the complex optical coherence structure, only the Fourier transform needed to be applied to the function Λ(u) , showing the simplicity and high efficiency of our proposed method. Further, for Schell-model PCBs, the protocol operated quite independently of the obstacle, regardless of its shape, size, and location. For the non-Schell-model PCBs, the Dirac pinhole was required as the obstacle, and the location of this obstacle determined the reference point position of the optical coherence structure. It is worth mentioning that the Schell-model PCBs were much more common than the non-Schell-model ones, especially in experiments conducted in laboratories. As mentioned in the Introduction, there are two very popular ways to produce PCBs. The first one, based on the van Cittert-Zernike theorem, can only be used for generating Schell-model PCBs whereas the second one, which is based on the mode-superposition principle, can be used for customizing Schell-model PCBs as well.

To generate the PCBs with the specific optical coherence structures in the experiment and simulation, we employed the mode- superposition principle to realize it in this study. First, we briefly introduced the RM method³⁹. According to the coherence theory, the optical coherence structure can be rewritten as¹⁸:

11μ(r1,r2)=∬P(κ)K(r1,κ)K*(r2,κ)d2κ,

where P(κ) is the power spectrum density function that satisfies the non-negativity condition, and K(r,κ) is an arbitrary kernel. Under the RM method³⁹, the above equation can be rearranged as:

12μ(r1,r2)=〈T(r1)T*(r2)〉,

with

13T(r)=∬P(κ)K(r,κ)Cn(κ)d2κ,

where Cn(κ) characterizes the electric field of the white noise and obeys Gaussian statistics. Hence, the RM described by T(r) also obeys Gaussian statistics. More details could be found in ref.³⁹.

To demonstrate the feasibility of our protocol for PCBs with different optical statistics, we followed^{23, 40} and adopted the PM method. The optical coherence structure is expressed below, which is the discrete form of Eq. (11):

14μ(r1,r2)=∑m=1M∑n=1NP(κmn)K(r1,κmn)K*(r2,κmn),

where κmn characterizes the sampling point across the area κ ; and M and N are the numbers of the sampling points of the function P(κ) in the horizontal and vertical directions, respectively—to achieve high precision, M and N should be large enough; P(κmn) is the mode weight; and K(r,κmn) is treated as the PM. This method does not follow Gaussian statistics²³.

In both methods, K(r,κ) determines the beam type of the PCB. The general form is given by K(r,κ)=exp[i2πH(r)−i2πϵ0a−b(xbκ⊥a+ybκ∥a)+i2πG(κ)] , where r≡(x,y) and κ≡(κ⊥,κ∥) , respectively; H(r) and G(κ) are arbitrary real functions of r and κ , respectively; ϵ0 is a constant in meters, and the parameters a and b are integer constants. For the Schell-model type, a=1 and b=1.

In the simulation, we could first create a set of modes in the source plane by the above equations. Then, the mode obstructed by an obstacle is captured in the Fraunhofer zone. As illustrated by Eq. (3), the electric field of the mode in the observation plane is given by

15Ef(u)=∫Ei(r)O(r)h(r,u)d2r,

where the Ei(r) and Ef(u) denote the electric fields of the modes in the source plane and the Fraunhofer plane, respectively. With the help of optical wave propagation principle and the commercial software Matlab⁴¹, we can achieve the modes in the observation plane. Established by the mode-superposition principle described above, we are able to simulate the recovery of optical coherent structures.

Figure 1 shows the experimental setup for the generation (in Part 1) and measurement (in Part 2) of the complex optical coherence structures of random optical fields. In part 1, a linearly polarized beam with a wavelength λ=632.8 nm is emitted from a He-Ne laser and passes through a half-wave plate (HP). Then, it is expanded by a beam expander (BE) and illuminates a reflective phase-only spatial light modulator (SLM1). We rotated the HP to guarantee that the incident beam on the SLM had horizontal polarization because the SLM1 only responds to this direction. In our experiment, we used the RM or PM superposition principle to generate PCBs with a customized optical coherence structure. To realize the desired mode in the laboratory, the complex amplitude of the RM/PM needed to be encoded [described by T(r) in Eq. (13) or K(r,κmn) in Eq. (14)] into the phase-only SLM. We followed the protocol described by Rosales-Guzmãn and Frobes^{42, 43} to design the computer-generated holograms (CGHs):

16ΦSLM(r)=FEsin[Arg[E(r)]+2πfxx],

where E(r) denotes a pre-encoded complex function; and FE is obtained by the numerical inversion of J1(FE)=Abs[E(r)] , where J₁ is the first-order Bessel function of the first kind; we applied “Abs” and “Arg” to the function to achieve its amplitude and phase, respectively; and f_x represents the grating frequency in the x-direction. The beam reflected by the SLM1 passed through a 4f optical imaging system consisting of the thin lenses L1 and L2. We used an iris in the Fourier plane to select the positive (or negative) 1^st diffraction beam, which was the desired mode; we refreshed the CGHs to achieve multiple modes. Based on Eqs. (12–14), we could produce PCBs belonging to either the Schell-model or non-Schell-model type and those with different optical statistics. The beam transmitted through the second beam splitter (BS2) arrived at the charge-coupled device (CCD1). When a mass of the modes were captured, the desired PCBs could be recovered via the mode-superposition principle. In Part 2, the beam reflected by BS2 arrived at SLM2. We followed the same protocol [Eq. (16)] to encode an irregular obstacle, as described by Oi(r),i=1,2,3 , into SLM2. The modulated beam was focused by the thin lens L3 and captured by CCD2. The distance between SLM2 and L3 and that between L3 and CCD2 were the focal lengths of L3. Hence, the CCD2 plane was treated as a Fourier plane. We used an iris to select the positive (or negative) 1^st diffraction beams. We captured three intensity profiles using CCD2 via refreshing the CGHs (i.e., the obstacles) in SLM2 and employed Eqs. (4–10) to recover the optical coherence structure. In the following section, we present the obtained simulation and experimental results that demonstrate the validity of our protocol.

We first considered Schell-model PCBs subjected to Gaussian statistics. To illustrate the robustness of this protocol, we followed the RM method [see Eqs. (11–13)] and selected a complex PCB without any analytical expression, where we set the letter “A” as the power spectrum density function P(κ) and K(r,κ)=exp(−i2πr⋅κ) . We chose an irregular obstacle [see Fig. 1(b)]. Figure 2 shows the real and imaginary parts of μ(r) and its square of the modulus; the simulation and experimental results are presented in the upper and lower rows, respectively. The perfect agreement between both results demonstrate the validity of our protocol. More importantly, the robust applicability of our method to complex beams and arbitrary obstacles well demonstrates the robustness of our method. Further, to quantify the precision of our proposed protocol, we adopt the widely used structural similarity (SSIM)⁴⁴ to evaluate the similarity between experimental/simulation results and theoretical results (not shown here to save space). The SSIM falls into the interval [0, 1], and the larger the value, the better the experimental/simulation results. The SSIMs of the real and imaginary parts of simulation and experimental results are very close to 1, as shown in the figures. It demonstrates that our protocol has high precision.

Next, we shall describe how this protocol is also applicable to Schell-model PCBs with non-Gaussian statistics. As suggested by the PM method in Eq. (14), we employed a complex amplitude-modulation encoding algorithm [see Eq. (16)] to produce each PM and then synthesize the PCB. In such a PCB, we adopted an off-axis ellipse band as the power spectrum density function, P(κ) , where it is equal to 1 inside the band and 0 otherwise. The major-to minor-axis ratio was set to 3∶2. The width of the ellipse band was 0.16 mm, and the displacement was 0.8 mm along the horizontal axis. Furthermore, K(r,κ) could be obtained as: K(r,κ)=exp(−i2πr⋅κ) . In order to show the differences and similarities between the same PCB but with different optical statistics, we display the experimental results from the above PCBs with Gaussian (produced by the RM method) and non-Gaussian (produced by the PM method) statistics in the upper and lower rows of Fig. 3, respectively. We expected that their individual mode patterns would be different. The patterns of the RM and PM are shown in Fig. 3(a) and 3(e), respectively. Under the RM method, the generated RM obeyed Gaussian statistics, and hence, the intensity probability density function (PDF) obeyed a negative exponent. The experimental results for both space and time are shown in Fig. 3(b). The experimentally obtained intensity PDFs of the PM are shown in Fig. 3(f); they do not yield Gaussian statistics in either space or time. We employed our protocol to recover the optical coherence of the structure; the real and imaginary parts of the optical coherence structures are shown in Fig. 3(c), 3(d), 3(g), and 3(h). Although they had different optical statistics, the results were identical for both beams. This demonstrated that the proposed method would be applicable to PCBs with different optical statistics.

For a non-Schell-model PCB, that is, a PCB with a non-uniformly correlated function⁴⁵, the optical coherence structure is highly sensitive to the reference point position. We applied Eq. (10) and tried to recover the optical coherence structure of the non-Schell-model PCB. Such a beam subjected to non-Gaussian statistics, is chosen as one example. Here, the adopted power spectrum density function is the same as that in Fig. 3, and the kernel function K(r,κ) is given by:

17K(r,κ)=exp[−i2π(x2κ⊥+y2κ∥)/4a0],

where a0 is constant and assumed to be a0=1m . We followed PM method to produce such a beam. The experimental results for the non-Schell model PCB are shown in Fig. 4. We adopted different reference point positions to show the features of the optical coherence structure with the non-Schell-model type. The simulation and experimental results for the real and imaginary parts of the optical coherence structures are displayed in the upper and lower rows, respectively, in Fig. 4. We adopted the reference point positions r₀=(0 mm, 0 mm) and r₀ =(1 mm, 0 mm) in Fig. 4(a–d) and 4(e–h), respectively. The distributions of both parts of the optical coherence structures are indeed sensitive to the reference point position r₀, showing non-uniform correlation features. Overall, the simulation and experimental results agreed extremely well.

Finally, we shall show our protocol can be highly used in the optical coherence structure related applications, and can be conducive to rapid improvement of application efficiency. We took one example - optical coherence encryption protocol proposed by Peng et al¹⁷. In that work, they used a generalized HBT method for this application, which must take a lot of time for doing speckles statistics to recover a single image. Here we displayed that the proposed protocol can decrypt the dynamic images (or a video). Regarding the encryption technique, they followed Eq. (11); the P(κ) function was plaintext, the K(r,κ) function was the encryption key, and the optical coherence structure was ciphertext. Here, we have reconducted the encryption experiment, setting the dynamic letter “OES” as the P(κ) function and K(r,κ)=exp[i10π(x2+y2)−i2π(xκ⊥+yκ∥)] . We used our protocol to recover the optical coherence structures in the experiment; the obtained three instantaneous real and imaginary parts are shown in the first and second rows in Fig. 5, respectively. The images can be recovered via an inverse fractional Fourier transform of the optical coherence structure¹⁷, and the corresponding instantaneous reconstructed images are shown in the bottom row of Fig. 5. The images rotated and shifted, proving that we could decrypt the dynamic images (or videos) encrypted using the optical coherence encryption technique. The processing time of mathematical calculation by software Matlab can be ignored, and the speed of image recovery depends on the refresh rate of the SLM and frame rate of the CCD. In our experiment, the refresh rate of the SLM was 60 Hz, and the frame rate of the CCD was 60 fps. Therefore, we can decrypt the video at ideal 20 fps. We expected to improve the speed of image recovery by adopting a SLM and a CCD with a higher refresh rate and frame rate. Regarding current commercial optical devices, our protocol can recover videos at frame rates of up to 100 fps.

In this study, we propose a robust, convenient, and fast technique to measure the optical coherence structures of random optical fields using a classical optical diffraction experiment—the generalized Arago spot experiment. We have provided rigorous mathematical solutions and presented them in extremely simplified analytical expressions based on different beam types; these solutions are applicable to any genuine PCB subjected to different optical statistics. Regarding the Schell-model PCBs common in laboratories, our protocol is valid for any obstacle, regardless of its shape, size, and location. Regarding the non-Schell-model type (i.e., non-uniformly correlated PCBs), the Dirac pinhole was adopted as the obstacle, and the location of the obstacle determined the reference point position of the degree of coherence function. In the experiment, it was necessary to capture the far-zone intensity only three times, and after simple calculations, we could precisely recover the real and imaginary parts of the complex optical coherence structures. The experimental results perfectly agreed with the simulation ones. They well indicate our protocol has simplicity, high efficiency, and precision. We believe that our proposal will shed new light on phase measurement, optical imaging, and other applications.

This work was supported by the National Key Research and Development Program of China (2022YFA1404800, 2019YFA0705000), National Natural Science Foundation of China (11974218, 12004220, 12104264, 12192254, 92250304), Regional Science and Technology Development Project of the Central Government (YDZX20203700001766), China Postdoctoral Science Foundation (2022T150392), and the Natural Science Foundation of Shandong Province (ZR2021QA014, ZR2021ZD02).

X. Liu and Q. Chen completed the theoretical model and performed the experiments. J. Zeng conducted the analysis and prepared the manuscript. C. H. Liang proposed the original idea. Y. J. Cai supervised the work. All the authors discussed and analyzed the results and commented on the manuscript at all stages.

The authors declare no competing financial interests.