Super-resolution imaging based on active optical intensity interferometry
Lu-Chuan Liu Cheng Wu Wei Li Yu-Ao Chen Frank Wilczek Xiao-Peng Shao Feihu Xu Qiang Zhang Jian-Wei Pan




0.0cm
0.2cm
16cm 
21cm
1.0cm
    




adja-page
#1 #2 (previous page)














 
















 
 







0.0cm 0.2cm 16cm 21cm 1.0cm



















sciabstract
	 
	












24pt



Super-resolution imaging based on active optical intensity interferometry
    Lu-Chuan Liu^1,2,3*, Cheng Wu^1,2,3*, Wei Li^1,4,5*, Yu-Ao Chen^1,2,3,

 Frank Wilczek^6,7,8,9,10, Xiao-Peng Shao^4,5, Feihu Xu^1,2,3,

Qiang Zhang^1,2,3,11 and Jian-Wei Pan^1,2,3

	

	^1Hefei National Research Center for Physical Sciences at the Microscale and School of Physical Sciences,

	University of Science and Technology of China, Hefei, 230026, China

	^2Shanghai Research Center for Quantum Science and CAS Center for 
	

	Excellence in Quantum Information and Quantum Physics,
	

	University of Science and Technology of China, Shanghai, 201315, China

	^3Hefei National Laboratory, University of Science and Technology of China, Hefei, 230026, China

	^4Hangzhou Institute of Technology, Xidian University, Hangzhou, 311200, China

	^5School of Optoelectronic Engineering, Xidian University, Xi’an, 710071, China

	^6Center for Theoretical Physics, MIT, Cambridge, MA 02139, USA

	^7T. D. Lee Institute, Shanghai Jiao Tong University, Shanghai, 200240, China

	^8Wilczek Quantum Center, School of Physics and Astronomy,
	

	Shanghai Jiao Tong University, Shanghai, 200240, China

	^9Department of Physics, Stockholm University, Stockholm, SE-106 91, Sweden

	^10Department of Physics, Arizona State University, Tempe, AZ 85287, USA
	
	^11Key Laboratory of Space Active Opto-Electronic Technology, Shanghai Institute of Technical Physics,
	

	Chinese Academy of Sciences, Shanghai, 200083, China

	^∗ These authors contributed equally to this work.

		

    
=============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





Long baseline diffraction-limited optical aperture synthesis technology by interferometry plays an important role in scientific study and practical application. In contrast to amplitude (phase) interferometry, intensity interferometry — which exploits the quantum nature of light to measure the photon bunching effect in thermal light — is robust against atmospheric turbulence and optical defects. However, a thermal light source typically has a significant divergence angle and a low average photon number per mode, forestalling the applicability over long ranges. Here, we propose and demonstrate active intensity interferometry for super-resolution imaging over the kilometer range. Our scheme exploits phase-independent multiple laser emitters to produce the thermal illumination and uses an elaborate computational algorithm to reconstruct the image. In outdoor environments, we image two-dimension millimeter-level targets over 1.36 kilometers at a resolution of 14 times the diffraction limit of a single telescope. High-resolution optical imaging and sensing are anticipated by applying long-baseline active intensity interferometry in general branches of physics and metrology.





§ INTRODUCTION


Hanbury Brown and Twiss (HBT) first proposed the seminal idea of high-angle resolution intensity interferometry for the measurement of stellar diameter <cit.>, which is a technique based on the measurement of temporal correlations of arrival time between photons recorded in different detectors. The intensity interferometry does not require a high-precision phase stabilization system, and it is robust against optical imperfections and atmospheric turbulence. Hence, it is much easier to implement than amplitude (phase) interferometry. In astrophysics, the forthcoming air Cherenkov telescope arrays, which consist of almost 100 telescopes, are the current largest optical intensity interference project <cit.>. Following its start in astronomy, intensity interferometry has been developed as a versatile lab tool in general research fields, including the probe of interactions in high-energy particle physics <cit.>, studying the photon propagation in nonlinear media <cit.> and curved space <cit.>, measurements of quantum correlations in ultracold bosonic and fermionic systems <cit.>, and identification of single photon source <cit.>. 

Although intensity interferometry has been widely implemented in passive imaging scenes <cit.>, it remains challenging to apply to the field of active imaging, i.e., light detection and ranging (LiDAR). Benefiting from active illumination, LiDAR has emerged as a powerful tool to image non-self-luminous targets <cit.>. To probe high resolution and resist atmospheric turbulence, active intensity interferometry has emerges as a great candidate. However, the lack of collimated a narrow-spectrum thermal light source and a robust imaging retrieval algorithm imposes practical difficulties in its applications in active super-resolution imaging. 

Intensity interference originates from the quantum nature of thermal light sources, which can be understood as the photon bunching effect <cit.>. The shape of the bunching peak measured by the intensity interferometer is squared proportional to the Fourier transform of the shape of the source <cit.>. In experiments, intensity interference is obtained by cross-correlating intensity fluctuations measured in different pairs of telescopes to achieve aperture synthesis and yield resolution enhancement. In general, unlike lasers, thermal light sources cannot simultaneously meet the requirements of intense light power, narrow spectra and small divergence angles. A better alternative solution involves utilizing a composite approach integrating a laser source with a spatial-temporal modulation device, such as a rotating ground glass <cit.>, a spatial light modulator <cit.>, and a projector <cit.> to synthesize a pseudo-thermal light source. However, these pseudothermal light sources are still incompetent for long-range illumination.

In this letter, we show that a phase-independent multiple laser emitter array can act as a pseudothermal light source, satisfying both the thermal nature of the source and the requirement of long-range illumination for intensity-interferometry imaging. In the experiment, we illuminate the targets with eight independent laser emitters modulated by the atmosphere and receive the back-reflected photons using two telescopes with configurable transverse range. For image retrieval, we develop an elaborate data preprocessing technique equipped with a modified phase retrieval algorithm, each tailored to adapt the intensity interferometric data. Our algorithm copes with sparse sampling in the interferometric plane, the nonlinearity of the imaging model, and the loss of Fourier phase information. Together, we experimentally realize super-resolution imaging over 1.36 km in an outdoor urban environment. The imaging resolution is 3 mm, which is 14 times higher than the diffraction limit of a single telescope. 




§ THEORY OF ACTIVE OPTICAL INTENSITY INTERFEROMETRY



An optical intensity interferometry in a classical imaging system typically necessitates that the imaging target to be a thermal light source to measure the non-shot-noise-based intensity variations of the two light detectors and compute their correlation. For ideal monochromatic polarized thermal light, Hanbury Brown and Twiss  showed that the second-order intensity correlation function is proportional to the square of the modulus of the Fourier transform of the intensity distribution on the target surface. We now express it in mathematical form. Suppose we have two light detectors D_a and D_b, whose intensities measured at time t are I_a(t) and I_b(t), respectively, and let k_a and k_b be the wave vectors from the target to the detectors D_a and D_b, respectively. Then, according to the intensity interferometry theory <cit.>, the normalized intensity correlation function c_ab^(2) can be expressed as

    c_ab^(2)=⟨Δ I_a(t)Δ I_b(t)⟩/⟨ I_a(t) ⟩⟨ I_b(t) ⟩=|f(k_a-k_b)|^2,

where ⟨...⟩ is the time averaged notation, Δ I(t)≡ I(t)-⟨ I(t)⟩ and f is the normalized Fourier function of the target surface intensity distribution ρ(r), that is

    f(Δk)=∫ρ(r)e^-iΔk·rdr/∫ρ(r)dr.

Let L be the distance from the target to the interferometer and λ be the wavelength of light. 
Since Δk and the projected interferometer baseline B have a proportional relationship that is  Δk=2πB/(λ L), super-resolution can be achieved by changing B to cover a sufficiently large  area of the u - v plane to obtain enough information about the target in the Fourier domain.  The angular resolution of the intensity interferometer is about λ/ B_max at the wavelength λ, which is larger than that of a single telescope <cit.>. 



However, if intensity interferometry is used for active imaging, it is hard to illuminate distant targets with poorly collimated thermal light sources. On the other hand, if we illuminate the target with well-collimated coherent light, such as a laser, intensity interference cannot be observed because the light detector cannot measure intensity signal fluctuations other than the shot noise, which is not correlated with the shot noise in another detector. 

To overcome the above difficulties, a basic idea is to superimpose a collection of phase-independent coherent lights to achieve pseudothermal illumination. The working principle of active intensity interferometry is schematized in Fig. <ref>. To simply verify the feasibility of applying super-resolution optical intensity interferometry to active imaging, as indicated in Fig. <ref>, we bypass the need for phase modulators by instead using atmospheric turbulence to randomly modulate the phases of multiple laser emitters. In the case of such non-ideal pseudothermal illumination, speckle-like noise occurs when measuring the intensity correlation function, so we propose a statistical optics theory (see supplementary materials), in which Eq. <ref> needs to be modified to a form of the ensemble average, namely


    ⟨ c_ab^(2)⟩_e=c_0+c_1|f(k_a-k_b)|^2,


where c_0 and c_1 are the coefficients jointly determined by the light intensity, the autocorrelation coefficient of different emitters, and the optical memory effect <cit.>. If we ignore the optical memory effect and consider a case where all laser emitters are symmetrically equivalent, then c_0=c/n and c_1=(n-1+c)/n are simple functions of the number of laser emitters n and the autocorrelation coefficient c.

The above statistical linear relationship allows us to extract information about the target in the Fourier domain from the measurement of the intensity correlation function. Specifically, considering that the ensemble average expression predicts the average result of multiple parallel measurements, for a single shot measurement of c_ab^(2), we can write c_ab^(2)=⟨ c_ab^(2)⟩_e+ϵ, where ϵ is the speckle-like noise term that we wish to suppress as much as possible. We calculate an expression for the noise intensity (i.e., the standard deviation of c_ab^(2)) and find that it decreases to 0 at a rate of O(1/√(n)) as the number of laser emitters n increases (see supplementary materials). This means that one way to suppress the experimental noise is to increase the number of laser emitters so that they are physically closer to thermal illumination. We also discuss the boundary cases of single (n=1) and infinite (n=∞) laser emitters (fig. <ref>), and show that the former will lead to a very low signal-to-noise ratio (SNR), while the latter will give predictions that are consistent with classical HBT theory based on ideal thermal light sources.



§ INTENSITY INTERFEROMETRY IMAGING SETUP





We implement an intensity interferometry imaging experiment in an urban atmospheric environment 1.36 kilometers away. Figure  <ref> shows an aerial view of the experiment configuration, with the imaging system placed at location S, facing the imaging target located at T over a distance of 1.36 km (determined by laser ranging measurement). The target is put on the 3-axis motorized rotation stages, as shown in Fig. <ref>B. 

The active optical intensity interferometer consists of phase-independent multiple laser emitters and two identical receiving telescope systems, as shown in Fig. <ref>A. A λ = 1064 nm CW laser (∼100 mW) is depolarized by a single-mode optical fiber depolarizer to reduce the intensity fluctuation caused by different polarization reflectivities of the target. The light is coupled into free space by a collimator with a focal length of 40 mm and is divided into 8 beams by a beam splitter group with equal spatial separation and power. Eight beams are adjusted to point at the target. The size of each beam after 1.36 km is approximately 0.2 m, more than the size of the target. Since light is subject to phase fluctuations through atmospheric turbulence, the phase shifts of each beam are randomly modulated <cit.>. The distance between each beam at the laser emitter array is set to 0.15 m, more than the outer scale of atmospheric turbulence (the atmospheric coherence length typically measured 0.02∼0.05 m), to ensure that each beam has an independent random phase shift. In the theoretical analysis, a collection of phase-independent multiple coherent light can be equivalent to a pseudothermal light source. 

Two identical receiver systems D_a and D_b are separately placed on two 0.4-m-long linear translation stages to change the baseline in the range of [0.07, 0.87] m. The receiving system is designed for rotational symmetry, collecting circularly polarized light through a quarter-wave plate and a polarizer. A tracking mirror mounted on a piezo platform is used to adjust the receiving field of view (FOV), ensuring that the target remains within the FOV while the telescope system moves. A dichroic mirror is introduced spectrally to separate the collected light. The signal light with a wavelength of 1064 nm transmitting to the dichroic mirror is coupled to a 62.5-μm multi-mode optical fiber with a telescope that has a diameter of 42.5 mm and a focal length of 80 mm, guided to the silicon single-photon avalanche diode (SPAD). Other background light is received by the camera with a 300 mm focal length lens for feedback on the tracking mirror pointing. 

The arrival time information of photons at two SPADs is recorded using a time-to-digital converter (TDC) with a time resolution of 8 ps, from which the normalized intensity correlation function c_ab^(2) can be calculated. The optical delay difference between the two receiving systems is controlled below 1 ns. Considering the optical coherence time of the atmosphere (typical value about 27 ms at 1064 nm, fig. <ref>), the time bin and the integration time are set to 2 ms and 2 s, respectively. The count rate of each SPAD is about 10^4 photons per second and the SNR of coincidence counting measurement is about 10 (see supplementary materials for more details).



§ IMAGE RECONSTRUCTION




 §.§ Data acquisition.


In the measurement, we change the baseline distance between two receiving telescope systems from 0.07 m to 0.87 m per 0.04 m steps. When the baseline is set to a specific length, the target is rotated every 6^∘ in the range of 0 - 360^∘, and the intensity correlation function c_ab^(2) is obtained from the coincidence counts detected by D_a and D_b. A total of 60 × 21 = 1260 coincidence count sets are gathered, and the corresponding Fourier (u,v) coverage is depicted in Fig. <ref>A. 



 §.§ Reconstruction.


Intensity interferometers do not directly provide images but intensity correlation function that is proportional to the squared-magnitude of the Fourier transform of the source brightness <cit.>. To recover an image, the phase of the Fourier transform must be determined in addition to its magnitude <cit.>. In real observations, however, there remain three main obstacles for existing image reconstruction techniques. First, the limited number of telescopes and observations yields very sparse frequency coverage. Second, since the square root operation needs to be applied to the intensity correlation function to obtain the Fourier magnitude, the imaging model is non-linear and especially sensitive to noise near zero values. Third, deploying the active intensity interferometer outdoors has to face strong urban environment noise, and the imperfection comes from the mimic pseudothermal light source. Hence robust phase retrieval techniques are needed. The reconstruction pipeline is shown in Fig. <ref>.

To address the nonconvex and ill-conditioned inverse problem, elaborate data preprocessing is designed to extract the target Fourier magnitude from the raw measurements. Before entering the phase retrieval procedure, the raw intensity correlation function c_ab^(2) should pass through five successive processes. The data in polar coordinates were first averaged with adjacent multiple angles, and then central low-frequency estimation, target Fourier magnitude estimation, and high-frequency noise suppression were performed. Finally, the data were converted into Cartesian coordinates and produced by two-dimensional Fourier interpolation to unveil the target Fourier magnitude. All data preprocessing is detailed in the supplementary materials.

After data preprocessing, we finally obtained the retrieved target Fourier magnitude y =|F(ρ(𝐫))| (Fig. <ref>B). As another key ingredient in the reconstruction process, phase retrieval is designed to reconstruct an image from the Fourier-plane magnitude and object-plane constraint. Let x̂ = ρ(𝐫)∈ℂ^n be the target surface intensity distribution and y ∈ℝ_+^N be the acquired data of |F x̂|. The phase retrieval problem is to find an approximate solution to the equation,

    y =|F x|+w,    x ∈ℂ^n

where w ∈ℝ^N represents the unknown discrepancy between the data predicted according to the Fourier transform model F and the actually acquired data. Mathematically, Eq. <ref> can be modeled as

    x̂=min _x Φ(x)+λ R(x),

where Φ(x) = 1 / 2y-|F x|_2^2 is a data-fidelity regularizer that ensures consistency between the reconstructed result and measurements, and R(x) is a regularizer that imposes a certain prior constraint. 

To solve the above problem (Eq. <ref>), we altered the iterative alternating projection approach <cit.> which incorporates a cascaded phase retrieval procedure composed of Relaxed Averaged Alternating Reflections (RAAR) <cit.> and Hybrid Input and Output (HIO) algorithm <cit.> with multiple constraints (sparsity, positivity, shrink support, and oversampling) (Fig. <ref>D). Putting all together, the elaborate data preprocessing followed by the modified phase retrieval strategy explicitly improves image recovery performance and stability (Fig. <ref>, fig. <ref>, and movie S1). More detailed mathematical treatments of the imaging process appear in the supplementary materials.





§ IMAGING RESULTS






We demonstrate the effectiveness of the proposed framework with both synthetic and experimental data. The resolution of our system is first confirmed by imaging a double slit. As shown in Fig. <ref>, the result clearly shows that an image with a 3 mm cross-range resolution at approximately 1.36 km can be achieved, which is approximately 42.5 mm / 3 mm ≈ 14 times resolution enhancement over the single receiving aperture (≈ 42.5 mm). The theoretical resolution predicted by the Rayleigh diffraction limit is 1.22λ L/B ≈ 2.03 mm, where λ, L, and B represent the illumination wavelength, imaging distance and the size of the synthetic aperture, respectively. 

We produce super-resolution imaging results in Fig. <ref>. The targets (Fig. <ref>A) were crafted from surface-blackened aluminum sheets, hollowed out to create the letters "USTC", which are covered with retroreflective sheeting. The size of each individual letter is 8 mm by 9 mm, and the width of the character is 1.5 mm. The measured intensity correlation points at sampled frequencies are rebinned to match the grid of frequels with a size of 512 × 512 points. After applying the special data preprocessing pipeline (see fig. <ref> and supplementary materials), the targets' Fourier magnitude |F(ρ(𝐫))| can be retrieved (Fig. <ref>C) and then used as input for the modified phase retrieval to obtain the two-dimensional target images (Fig. <ref>D). Shape reconstructions with PSNRs of 14 ∼ 18 dB demonstrate the feasibility of the proposed method. A superposition of 1000 individual trials of phase retrieval results is shown in Fig. <ref>E, which clearly demonstrates high image restoration fidelity and stability. Potential image quality improvement may be achieved by increasing the number of laser emitters n, the capture frames (fig. <ref>) as well as the sampling rate in the interferometric plane (figs. <ref> and  <ref>).



§ DISCUSSIONS


In summary, we have shown active intensity interferometry for super-resolution imaging over the kilometer range. 
In the experiment, phase-independent multiple laser emitters based on the phase fluctuations of atmosphere turbulence can act as an elegant pseudothermal light source to illuminate long-range objects, but the exposure time of imaging is limited to atmospheric coherence time. With the introduction of an active phase modulation scheme and a confocal receiving-emitting telescope array, high-speed active intensity interferometry imaging is an interesting avenue for future work.
On the other hand, higher-fidelity reconstruction is anticipated by applying more advanced methods, such as the optimization strategy and deep learning modules. Furthermore, by using different speckle interferometry analysis techniques, such as phase closure analysis and amplitude interference, additional phase information may be obtained from the data in some cases, reducing the ambiguity of recovery and further improving the accuracy and fidelity of phase retrieval. In the future, efforts are being devoted to joint measuring the first and second-order coherence of light to uncover unrivaled imaging resolution. Further combining the active intensity interferometry with deep space exploration and microscopy will open a new path towards super-resolution imaging from macroscopic scale to microscopic scale. 

Science




§ ACKNOWLEDGMENTS

We thank Zhen-Tao Liu and Shen-Sheng Han from Shanghai Institute of Optics and Fine Mechanics for their suggestions on this work. 



  
*Funding:This research was supported by National Natural Science Foundation of China (grants T2125010 and 12204466); Natural Science Foundation of Shanghai (grant 22ZR1468200); Anhui Initiative in Quantum Information Technologies (grant AHY010100); Shanghai Municipal Science and Technology Major Project (grant 2019SHZDZX01); Innovation Programme for Quantum Science and Technology (grants 2021ZD0300100 and 2021ZD0300300). Y.-A. C., F.-H. X. and Q. Z. were supported by the XPLORER prize from New Corner Stone Science Foundation. F. W. was supported by the U.S. Department of Energy under grant Contract Number DE-SC0012567, by the European Research Council under grant 742104, and by the Swedish Research Council under Contract No. 335-2014-7424. 



  
*Author contributions:Q. Z. and J.-W. P. conceived the research. L.-C. L. carried out the work of the theoretical derivation. L.-C. L. and C. W. built the experimental setup and collected the experimental data. W. L. worked on experimental simulations and developed the image reconstruction algorithms. All authors contributed to data processing, data analysis, and writing of the manuscript. 



  
*Competing interests:The authors declare no competing interests. 



  
*Data and materials availability:The measured intensity correlation data and the processing code supporting the findings of this study will be available online upon publication. Additional data and code are available from the corresponding authors upon request.



§ SUPPLEMENTARY MATERIALS

Materials and Methods

Figs. S1 to S6

References (36-42)

Movie S1








Supplementary Materials for Super-resolution imaging based on active optical intensity interferometry
    Lu-Chuan Liu^1,2,3*, Cheng Wu^1,2,3*, Wei Li^1,4,5*, Yu-Ao Chen^1,2,3,

 Frank Wilczek^6,7,8,9,10, Xiao-Peng Shao^4,5, Feihu Xu^1,2,3,

Qiang Zhang^1,2,3,11 and Jian-Wei Pan^1,2,3
    
=======================================================================================================================================================================================



The PDF file includes:

Materials and Methods

Figs. S1 to S6

References (36-42)

Other Supplementary material for this manuscript includes the following:

Movie S1



§ MATERIALS & METHODS




§ SIMULATIONS.
 


To gain further insight into the mechanism of active intensity interferometry, we performed numerical simulations of the whole process. The entire imaging processing and simulations are presented in Fig. <ref>, with the simulation process, data preprocessing and retrieval process indicated as red, orange, and blue backgrounds, respectively. To simulate the HBT sampling measurement, a target (Fig. <ref>A) is first transformed into the Fourier domain to achieve its Fourier magnitude (Fig. <ref>B). Then, the Fourier magnitude is sampled corresponding to the experimental measurement setup (see HBT imaging setup, data acquisition and Fig. <ref>C). The sparse Fourier sampling is then remapped using the linear transform to obtain the HBT sampling data i.e., the intensity correlation function c_ab^(2) (Fig. <ref>D). The sampled c_ab^(2) data were first converted into Polar coordinates (Fig. <ref>E) to simulate the real experiment acquisition. 

We describe our HBT simulation procedure in Algorithm <ref>. Please see Section <ref> for a detailed explanation of the simulation process.



Sfigure







§ COMPUTATIONAL RETRIEVAL METHOD.




 §.§ Image processing.



As described in the main text, to extract the target information from the measured intensity correlation function c_ab^(2), we have to face three main challenges. First, the interferometric sampling remains very sparse (for the 1260 experiment data points fall in 512 × 512 grids, leading to a 0.48% sampling rate) (see Fig. <ref>B), and the central crucial low-frequency region (frequency components lie under the 0.07 m circle) is not measured due to the physical limitations of the distance between two telescopes. Second, to extract the target Fourier magnitude |F(ρ(𝐫))| = √(c_1)|f(k_a-k_b)| with the Fourier transform operator defined as F(·) from c_ab^(2), we need to calculate Eq. <ref>, which makes the model nonlinearity and noise sensitive to the close to zero values caused by the square root operation. Third, recovering the missing Fourier phase part is difficult due to the high background noise originating from the outdoor environment and source imperfections. 

To address this problem, the raw captured intensity correlation function c_ab^(2) was first subtracted from the background and then polarization elimination was performed. The central symmetric average and adjacent five angle average were employed to make the data much approach to ⟨ c_ab^(2)⟩_e, which obeys the Gaussian process (Section <ref>). The central null part is preserved for low-frequency component estimation (Fig. <ref>E). Virtual frequency components at the positions of [-4 -2 0 2 4] cm were predicted by summing and averaging after 1-D cubic spline interpolation for each angle (Fig. <ref>F). The estimated target Fourier magnitude |F(ρ(𝐫))| (Fig. <ref>G) was calculated using the Eq. <ref> (c_0 = 0.06 in Fig. <ref> and 0.15 in Fig. <ref>) and then smoothed by a Ramp-Lank filter to remove high-frequency oscillations (Fig. <ref>H). The filtered data were transformed into Cartesian coordinates and interpolated in two dimensions to fill the voids in the sparse sampled spatial frequencies to obtain the target Fourier magnitude. The resulting Fourier magnitude was oversampled with one-fold zero values, and a two-dimensional Hanning window was applied before entering the modified phase retrieval algorithm.



 §.§ Retrieval algorithm.
 


As described in the Reconstruction section, to solve the highly ill-posed phase retrieval problem (Eq. <ref>), we leverage nonconvex optimization methods inspired by the alternating projections and reflections schemes <cit.>. Apart from normalization factors, we denote the target with reflectivity ρ(𝐫), 𝐫 is the coordinates in the object (or real) space, and the corresponding Fourier Transform ρ̃(𝐤) = F(ρ(𝐫)), with 𝐤 representing the coordinate in the Fourier (or Reciprocal) space. Four operators linking two sets 𝐒 (support) and 𝐌 (modulus) are introduced. Given a support set 𝐒, the support projection operator 𝐏_𝐒 involves setting to 0 the components outside the support while leaving the rest of the values unchanged:




    𝐏_𝐒(𝐫)= ρ(𝐫),     if 𝐫∈𝐒, 
     0,     otherwise


Similarly, one can extend this definition to enforce nonnegativity, 𝐏_𝐒+ρ(𝐫)=max(𝐏_𝐒ρ(𝐫), 0). The modulus projector 𝐏_𝐌 sets the modulus to the retrieved one m(𝐤), and leaving the phase unchanged. 


    𝐏_𝐌ρ̃(𝐤)=𝐏_𝐌|ρ̃(𝐤)| e^i φ(𝐤)=m(𝐤) e^i φ(𝐤)


The reflection operators 𝐑_𝐒 =2 𝐏_𝐒-𝐈 and 𝐑_𝐌=2 𝐏_𝐌-𝐈 apply the same step as the projector but move twice as far, where I means the unit matrix. The phase retrieval algorithm requires a starting point ρ^0(𝐫), generated by assigning a random phase to the retrieved target Fourier magnitude iteratively produce the iterate ρ^(n)(𝐫) after n steps. With the notation above, the HIO algorithm can be formulated as


    ρ^(n+1)(𝐫)=[𝐏_𝐒𝐏_𝐌+𝐏_𝐒(𝐈-β𝐏_𝐌)] ρ^(n)(𝐫)


Where 𝐏_𝐒=(𝐈-𝐏_𝐒) denotes the complement of the projector 𝐒 and the β(0<β<1) here and below means a negative feedback parameter. Empirically, small β produces better stability while β is better for escaping local minima. A value of β somewhere between 0.5 and 1.0, say 0.8, usually works well. In all experiments and simulations, the β is set as 0.8. Similarly, the RAAR algorithm can be written as


    ρ^(n+1)(𝐫)=[1/2β(𝐑_𝐒𝐑_𝐌+𝐈)+(1-β) 𝐏_𝐌] ρ^(n)(𝐫).


The retrieved target Fourier magnitude (Fig. <ref>B) was one-fold zero-padded before entering into the phase retrieval algorithm. For better image performance, a combination of RAAR (60 iterations) and HIO (one iteration) algorithms is employed, with the RAAR escaping local minima and HIO refining the solution. Along with the iterations, two sets of constraints are equipped. The first object constraint was that it had to be real and non-negative. The second object constraint was a shrink-wrap support <cit.> to refine the support region 𝐒 during the iteration. Specifically, the support region 𝐒 periodically convolves the iteration with a Gaussian of width σ and then applies a threshold ϵ to smear the support; i.e., the updated support for the two-dimensional case is taken to be 


    S={𝐫:1/2 πσ^2∫_ℝ^2ρ^(n)(𝐬) e^-|𝐫-𝐬|^2 / 2 σ^2 d 𝐬≥ϵ}


The third object constraint used was a sparsity constraint by setting 20 % of the maximum entries of the original structure. The Fourier space constraint is the magnitude constraint that sets the modulus to the retrieved one. Accompanied by the iteration, the adaptive support region is controlled by an alternative Gaussian convolution and a threshold operation, which act as dilation and erosion. With the iterative number increases, the exported support region will change from a loose one to a tight one, which can significantly reduce the solution space and exclude spurious solutions. It is worth emphasizing that using a non-centrosymmetric initial support estimation, for instance, an equilateral triangle as initial support considerably speeds up convergence and effectively avoids the twin-image problem <cit.>.

All experiments were performed with MATLAB 2018b in Windows 10 running on a dual-core chip Intel Core i9-10900K and 128G memory. A single run of the algorithm on a 1024 × 1024 pixel image retrieval (composed of 60 (RAAR) +1 (HIO) iterations) on this CPU took approximately 1 s. Notably, all results displayed in experiments and simulations are unthresholded and without color bar alternation.



§ IMAGING PERFORMANCE VS. NUMBER OF LASER EMITTERS.



To evaluate the impact of different laser emitters on image performance, we performed a double slit imaging under different laser emitters with a slit-to-slit separation of approximately 2 mm. The experimental results under n = 1, 2, 4, and 8 emitters as well as the quality performance are shown in Fig. <ref>A and B. Simulations with the same target under n = 1, 2, 4, 8, 16, 32, 64 and ‘∞’ (the ground truth Fourier domain sampling) are carried out and the corresponding results are depicted in Fig. <ref>A and C. In the simulation, we generate 50 sets of data with random initialization for each emission group and perform 50 trials of phase retrieval for each set of data, yielding a total of 50 × 50 = 2500 reconstructions. For each emission group, the average peak SNRs of the 2500 reconstructions are counted, with the error bars indicating the standard deviation of the mean peak SNR. The thermal features of the light source can be improved by increasing the number of emitters, transforming the double slit from intractable to clearly separated. The experiments and simulations embrace with each other and both indicate that at least four emitters can correctly achieve the remote target. The averaged error bars are fairly modest, especially when the emitters are greater than four. It is worth emphasizing that the experimental results were obtained by single-shot measurements.





§ DIFFERENT EMISSIONS VS ACQUISITION FRAMES.


In the current experimental setup, the incoherent source is modeled by a collection of emitters. Measuring the intensity correlation function c_ab^(2) with pseudothermal light inevitably leads to speckle-like noise due to interference of the signal with scattered light. By increasing the number of emitters or merging numerous collection frames, speckle-like noise can be reduced while also improving the quality of the returned Fourier magnitude. Two shapes of targets are tested and the reconstruction results with different emissions and acquisition frames from numerical simulation are presented in Fig. <ref>. Both results indicated that increasing the number of laser emissions brings a major enhancement in image quality when compared to accumulating the acquisition frames, demonstrating that the thermal properties of the light source are more influential to the imaging performance. To reduce the speckle-like noise, increasing the number of laser emitters n has a convergence rate of O(1 / √(n)) and this can be further accelerated when the value of the Fourier spectrum approaches to 0 (O(1 / n)). However, by accumulating the acquisition frames k has a convergence rate of 1 / √(k) (Section <ref>). The variation in convergence rate causes disparities in the effectiveness of the two techniques for image quality enhancement. However, if emissions are limited, increasing the number of acquisitions can improve imaging quality.





§ COMPARISONS AMONG DIFFERENT PHASE RETRIEVAL ALGORITHMS.


To demonstrate the superiority of the proposed phase retrieval scheme, we conducted an experimental comparison by imaging the letter ‘S’ (Fig. <ref>F and Fig. <ref>D). The letter ‘S' was chosen because curve shapes contain more complex features in Fourier space. However, the retrieved target Fourier magnitude (Fig. <ref>C) is still noise-stained and has fainter high-frequency information (Fig. <ref>C). The requirements for applying the Gerchberg-Saxton (GS) <cit.> iterative phase retrieval algorithm and its derivations such as Error Reduction (ER) <cit.> and Hybrid Input Output (HIO), are not satisfied and the resulting images (Fig. <ref>A to C) are of poor quality. More advanced alternating projection algorithms such as Hybrid Projection Reflection (HPR) <cit.> (Fig. <ref>D) and Relaxed Averaged Alternating Reflectors (RAAR) (Fig. <ref>E), are also unable to recover the target. In stark contrast, the image reconstructed from the proposed method (Fig. <ref>G) significantly outperforms the others, with a diffraction-limited resolution of 1.5 mm (limited by the virtual aperture of 87 cm). High-fidelity reconstructions are attributed to the majority voting and reduced support constraint as well as the cascade scheme. For the algorithms compared, the number of iterations was set to 60, the initial support was used as a non-centrosymmetric triangle, and the negative feedback parameter β for the single-cycle HIO, HPR, and RAAR algorithms was fixed at 0.8.





§ COMPLEX TARGETS IMAGING VIA ACTIVE INTENSITY INTERFEROMETRY.


In this proof of concept, we have demonstrated two-dimensional superresolution imaging over kilometers range with millimeters precision via the proposed active intensity interferometry. Higher-fidelity complex target reconstruction is anticipated by improving the sampling rate on the interferometric plane and multiplying the number of baselines/telescopes. To demonstrate this, complex targets (the targets ‘satellite’ and ‘UFO’) with varying reflectivity levels imaged under a 247 cm baseline with a 1 cm sampling interval were tested, and the corresponding results are depicted in Fig. <ref>. Since the intensity interferometry connects telescopes only electronically, it is practically insensitive to atmospheric turbulence and optical imperfections, permitting observations over very long baselines and through large airmasses, as well as at short optical wavelengths. With kilometer baselines, sufficient sampling frequencies, and advanced algorithms, tantalizing results are expected to open the window toward new scientific results.






§ THEORY PART I. THEORY OF MULTIPLE LASER EMITTERS WITH THERMAL NATURE.



We model our imaging system based on statistical optics. As shown in Fig. <ref> in the main text, suppose our optical system consists of n laser emitters with angular frequency ω numbered 1∼ n in the source plane, a target that is a scattering sample of a particular structure in the object plane, and two movable detectors D_a and D_b in the detection plane. Detectors D_a and D_b always use the same polarization measurement in the experiment, then we can denote the electric fields they measure on this polarization component at time t as E_a(t) and E_b(t), respectively. For simplicity, we only consider the case of far-field illumination and detection, that is, the distance between the target and the emitters or detectors is much greater than the Fraunhofer distance 2D^2/λ, where λ is the laser wavelength and D is the maximum value among laser beam waists, detector apertures and target size. 

If D is much smaller than the atmospheric coherence length, we can approximate each emitter's illumination light and its scattered light scattered by the target as an intensity- and phase-modulated plane wave. To describe these plane waves, we denote k_i as the wave vector from the i-th emitter to the target, and denote k_a and k_b as the wave vectors from the target to the detector D_a and D_b, respectively. Combining the above approximations and further ignoring all CW laser delays due to the time-of-flight which is much smaller than the atmospheric coherence time, we can approximately express E_a(t) and E_b(t) as

    E_a,b(t)=∑_i=1^n√(I_i(t))e^i(ϕ_i(t)-ω t)T(k_i,k_a,b)T_a,b(t).

Here, I_i(t) and ϕ_i(t) represent the intensity and phase of the plane wavefront formed at the target at time t by the illumination light from the ith emitter, T(k_in,k_out) is the target's transmission matrix between the k_in input field and k_out output field and T_a(t), T_b(t) are the field propagation coefficients from the target to the detectors D_a and D_b at time t, respectively. All the I_i(t),ϕ_i(t) and T_a(t), T_b(t) are treated as time-varying random variables which incorporate the whole process of atmospheric disturbance. 

Considering that it is difficult to make a complete statistical process modeling of atmospheric disturbance, we approximate the impact of atmospheric disturbance as the following two simple statistical assumptions. The first assumption is that all the I_i(t),ϕ_i(t) and T_a(t), T_b(t) are independent of each other. Note that this assumption clearly does not hold when the distances between different emitters or detectors are smaller than the atmospheric coherence length, so we need to avoid this in the experiment. The second assumption is that all the ϕ_i(t) are Gaussian random processes. Let the notation ⟨ ...⟩_t denote the time average over a sufficiently long time. Define ϕ_i≡⟨ϕ_i(t)⟩_t and Δϕ_i(t)≡ϕ_i(t)-ϕ_i, using the properties of Gaussian distribution it's easy to prove that

    ⟨e^i(ϕ_i(t)-ϕ_j(t))⟩_t=e^i(ϕ_i-ϕ_j)-1/2⟨(Δϕ_i(t)-Δϕ_j(t))^2⟩_t.

In this way, if ⟨(Δϕ_i(t)-Δϕ_j(t))^2⟩_t≫ 1 when i≠ j, we can further have an approximation of

    ⟨e^i(ϕ_i(t)-ϕ_j(t))⟩_t=δ_ij≡{[    1  i=j;    0 i≠ j ].
    .

This approximation is valid when the distance between emitters is much greater than the atmospheric coherence length. In fact, the ⟨(Δϕ_i(t)-Δϕ_j(t))^2⟩_t≈ 1 condition can be used as a definition of atmospheric coherence length. For example, the Kolmogorov turbulence model<cit.> expected that

    ⟨(Δϕ_i(t)-Δϕ_j(t))^2⟩_t= 6.88(r_ij/r_0)^5/3,

where r_ij is the distance between the i-th emitter and the j-th emitter, r_0 is the Fried parameter, which can be used as a representation of the atmospheric coherence length.

Define I_a(t)≡|E_a(t)|^2 and I_b(t)≡|E_b(t)|^2 as the intensities measured by detectors D_a and D_b at time t, respectively, after using the above assumptions, the time-averaged intensity can be expressed as

    ⟨ I_a,b(t)⟩_t=⟨|E_a,b(t)|^2⟩_t=η_a,b∑_i=1^nI_i|T(k_i,k_a,b)|^2,

where η_a,b≡⟨|T_a,b(t)|^2⟩_t and I_i≡⟨ I_i(t)⟩_t. Define c_i^(2)≡⟨ I_i^2(t)⟩_t/I_i^2-1, we further have the intensity correlation function

    ⟨Δ I_a(t)Δ I_b(t)⟩_t   ≡⟨ I_a(t)I_b(t)⟩_t-⟨ I_a(t)⟩_t⟨ I_b(t)⟩_t
       =η_aη_b[∑_i≠ jI_i I_jT(k_i,k_a)T(k_j,k_b)T^*(k_i,k_b)T^*(k_j,k_a)
       +∑_i=1^nc_i^(2)I_i^2|T(k_i,k_a)|^2|T(k_i,k_b)|^2],


It can be seen that for a particular target, both the time-averaged intensity and the intensity correlation function have a particular value determined by the transmission matrix T(k_in,k_out) of the scattering sample. Now we consider an ensemble of scattering samples with the same geometric surface but different internal microstructure. The fully developed speckle model proposed by Goodman et al. predicts that if the sample surface is sufficiently rough compared to the wavelength of light, for each (k_in,k_out) pair, the transmission coefficient T(k_in,k_out) behaves as a circular complex Gaussian random variable under different ensembles <cit.>. 

We next demonstrate that this statistical property of the transmission matrix is the key factor causing speckle-like noise in the experimental data. Let the notation ⟨...⟩_e denote the ensemble average, and the notation ⟨...⟩_t,e≡⟨⟨...⟩_t⟩_e denote the ensemble average of a time average quantity. Using the moment theorem <cit.> for complex Gaussian random variables, we have the ensemble-averaged intensity correlation function

    ⟨Δ I_a(t)Δ I_b(t)⟩_t,e
    =η_aη_b[∑_i≠ jI_i I_j⟨ T(k_i,k_a)T^*(k_i,k_b)⟩_e⟨ T(k_j,k_b)T^*(k_j,k_a)⟩_e
       +∑_i≠ jI_i I_j⟨ T(k_i,k_a)T^*(k_j,k_a)⟩_e ⟨ T(k_j,k_b)T^*(k_i,k_b)⟩_e
       +∑_i=1^nc_i^(2)I_i^2(⟨|T(k_i,k_a)|^2⟩_e ⟨|T(k_i,k_b)|^2⟩_e+|⟨ T(k_i,k_a)T^*(k_i,k_b)⟩_e|^2)].


To further obtain an explicit expression of the above formula, we need to study the coherence function of the form ⟨ T(k_in1,k_out1)T^*(k_in2,k_out2)⟩_e, where k_in1=k_in2 or k_out1=k_out2. For the first case, the classical coherence theory based on the van Cittert-Zernike theorem gives <cit.>

    ⟨ T(k_in,k_out1)T^*(k_in,k_out2)⟩_e
    ∝∫ρ(r)e^-i(k_out1-k_out2).rdr,

which is the Fourier transform of the intensity distribution on the output surface of the target. So let

    ⟨ T(k_in,k_out1)T^*(k_in,k_out2)⟩_e=η f(k_out1-k_out2),

where η is a real constant and f is the normalized Fourier function, that is,

    f(Δk)=∫ρ(r)e^-iΔk.rdr/∫ρ(r)dr.

For the second case, we assume that

    ⟨ T(k_in1,k_out)T^*(k_in2,k_out)⟩_e=η h(k_in1-k_in2),

where h is a function that satisfies h(0)=f(0)=1. For surface scattering, classical coherence theory simply predicts h=f, but for volume scattering, the physical mechanism here is much more complicated. Some modern studies on speckle have found that the function h is related to the optical memory effect. Both theory <cit.> and experiment <cit.> show that h will decay to 0 when |k_i-k_j| is large, and we will not discuss the relevant details here. Substituting Eq. (<ref>) and Eq. (<ref>) into Eq. (<ref>) and Eq. (<ref>), we obtain

    ⟨ I_a,b(t)⟩_t,e=η_a,b η∑_i=1^nI_i

and

    ⟨Δ I_a(t)Δ I_b(t)⟩_t,e   =η_aη_b η^2[∑_i=1^nc_i^(2)I_i^2+∑_i≠ jI_i I_j|h(k_i-k_j)|^2
       +(∑_i=1^nc_i^(2)I_i^2+∑_i≠ jI_i I_j)|f(k_a-k_b)|^2],

respectively. In our experiment, the quantity we measure is the normalized intensity correlation function c_ab^(2), which is defined as

    c_ab^(2)≡⟨Δ I_a(t)Δ I_b(t)⟩_t/⟨ I_a(t)⟩_t,e⟨ I_b(t)⟩_t,e.

According to Eq. (<ref>) and Eq. (<ref>), the ensemble average of c_ab^(2) is

    ⟨ c_ab^(2)⟩_e=c_0+c_1|f(k_a-k_b)|^2,

where

    c_0=∑_i=1^nc_i^(2)I_i^2+∑_i≠ jI_i I_j|h(k_i-k_j)|^2/(∑_i=1^nI_i)^2,

and

    c_1=∑_i=1^nc_i^(2)I_i^2+∑_i≠ jI_i I_j/(∑_i=1^nI_i)^2.


To quantify the level of speckle-like noise, we further calculate σ(c_ab^(2)), which is the standard deviation of c_ab^(2), defined as


    σ(c_ab^(2))≡√(⟨(c_ab^(2)-⟨ c_ab^(2)⟩_e)^2⟩_e) .

If the optical memory effect of the scattering sample is neglected, we can show that

    σ(c_ab^(2))=√(d_0+d_1|f(k_a-k_b)|^2+d_2|f(k_a-k_b)|^4),

where

    d_0=∑_i=1^n3(c_i^(2)I_i^2)^2+∑_i≠ j(I_i I_j)^2/(∑_i=1^nI_i)^4,


    d_1=∑_i=1^n14(c_i^(2)I_i^2)^2+∑_i≠ j(12c_i^(2)I_i^3I_j+2I_i^2I_j^2)+∑_i≠ j,j≠ k,k≠ i2I_i^2I_j I_k/(∑_i=1^nI_i)^4,

and 

    d_2=∑_i=1^n3(c_i^(2)I_i^2)^2+∑_i≠ j(4c_i^(2)I_i^3I_j+3I_i^2I_j^2)+∑_i≠ j,j≠ k,k≠ i2I_i^2I_j I_k/(∑_i=1^nI_i)^4.

For the derivation details of the above formula, see Section <ref>.

To better understand these results from a physical perspective, we consider two simplified cases. The first case is to use only one laser emitter for illumination, that is, n=1, then the ensemble average and the standard deviation of c_ab^(2) can be simplified as

    ⟨ c_ab^(2)⟩_e=c_1^(2)(1+|f(k_a-k_b)|^2)

and 

    σ(c_ab^(2))=c_1^(2)(3+14|f(k_a-k_b)|^2+3|f(k_a-k_b)|^4)^1/2,

respectively. In the ideal situation without atmospheric disturbance, c_1^(2)=0, then c_ab^(2)=⟨ c_ab^(2)⟩_e=σ(c_ab^(2))=0 can be further obtained, which is in line with the expectation that no intensity interference signal can be detected under ideal coherent light illumination. When atmospheric disturbance can cause an autocorrelation coefficient c_1^(2) that cannot be ignored, an intensity interference signal proportional to c_1^(2) on average can actually be measured. However, since σ(c_ab^(2)) is also proportional to c_1^(2), the level of speckle-like noise is almost the same as the signal (easy to prove that √(3)<σ(c_ab^(2))/⟨ c_ab^(2)⟩_e<√(5)), resulting in a very low SNR.

The second case is when there is more than one emitter but all emitters have identical properties and the optical memory effect of the scattering sample is neglected. We assume I_1=I_2=...=I_n and c_1^(2)=c_2^(2)=...=c_n^(2)=c is an autocorrelation coefficient, then we have the expressions

    ⟨ c_ab^(2)⟩_e=c/n+n-1+c/n|f(k_a-k_b)|^2

and 

    σ(c_ab^(2))   =[n-1+3c^2/n^3+2n^2+(12c-4)n+14c^2-12c+2/n^3|f(k_a-k_b)|^2
       +2n^2+(4c-3)n+3c^2-4c+1/n^3|f(k_a-k_b)|^4]^1/2,

respectively. We can see that as the number of laser emitters n increases, ⟨ c_ab^(2)⟩_e will approach |f(k_a-k_b)|^2, and σ(c_ab^(2)) will decrease to 0 at a rate of O(1/√(n)), that is, the speckle-like noise will gradually disappear. When n approaches infinity, our theory and the classical HBT theory based on ideal thermal light sources will give consistent predictions.



§ THEORY PART II. DETAILED DERIVATION OF THE Σ(C_AB^(2)).



To simplify the calculation of σ(c_ab^(2)), we first introduce the notation X_1∼ X_n, which are defined as

    X_i≡I_i T(k_i,k_a)T^*(k_i,k_b)  (∀ i ∈{1,2,...,n}).

Then c_ab^(2) can be rewritten as

    c_ab^(2)=∑_i≠ jX_i X_j^*+∑_i=1^nc_i^(2)|X_i|^2/η^2(∑_i=1^nI_i)^2.


We also introduce notations D and Cov to represent variance and covariance under ensemble statistics, that is, if X and Y are random variables, we define

    D(X)   ≡⟨ X^2⟩_e-⟨ X⟩_e^2,
    
    Cov(X,Y)   ≡⟨ XY⟩_e-⟨ X⟩_e⟨ Y⟩_e.

In this way, σ(c_ab^(2)) can be expressed as

    σ(c_ab^(2))   =√(D(∑_i≠ jX_i X_j^*+∑_i=1^nc_i^(2)|X_i|^2))/η^2(∑_i=1^nI_i)^2
       =√(D(∑_i≠ jX_i X_j^*)+2Cov(∑_i≠ jX_i X_j^*,∑_i=1^nc_i^(2)|X_i|^2)+D(∑_i=1^nc_i^(2)|X_i|^2))/η^2(∑_i=1^nI_i)^2


Using the momentum theorem <cit.>, we can simply prove that for any i=1∼ n, X_i has the following properties:

    ⟨ X_i⟩_e   =ηI_if(k_a-k_b),
    ⟨ X_i^2⟩_e   =2η^2I_i^2f^2(k_a-k_b),
    ⟨ |X_i|^2⟩_e   =η^2I_i^2(1+|f(k_a-k_b)|^2),
    ⟨ X_i|X_i|^2⟩_e   =2η^3I_i^3f(k_a-k_b)(2+|f(k_a-k_b)|^2),
    ⟨ |X_i|^4⟩_e   =4η^4I_i^4(1+4|f(k_a-k_b)|^2+|f(k_a-k_b)|^4).

On the other hand, if the optical memory effect of the scattering sample is neglected, it means

    ⟨ T(k_in1,k_out1)T^*(k_in2,k_out2)⟩_e=0  (∀k_in1≠k_in2).

It can be immediately proved that X_1∼ X_n are independent of each other, that is, any multivariate moment can always be decomposed into the form of the product of univariate moments, and its mathematical expression is


    ⟨∏_i=1^nX_i^α_iX_i^*β_i⟩_e=∏_i=1^n⟨ X_i^α_iX_i^*β_i⟩_e,

where α_1∼α_n and β_1∼β_n are any non-negative integers.

Using the above conclusions, we can do the following calculations:

    D(∑_i≠ jX_i X_j^*)   =∑_i≠ j(D(X_i X_j^*)+Cov(X_i X_j^*,X_j X_i^*))
    
    +   ∑_i≠ j,j≠ k,k≠ iCov(X_i X_j^*,(X_i+X_j)X_k^*+X_k(X_i^*+X_j^*))
       =∑_i≠ jη^4I_i^2I_j^2(1+2|f(k_a-k_b)|^2+3|f(k_a-k_b)|^4)
    
    +   ∑_i≠ j,j≠ k,k≠ iη^4I_i I_j(I_i+I_j)I_k(|f(k_a-k_b)|^2+|f(k_a-k_b)|^4),
    
    Cov(∑_i≠ jX_i X_j^*,∑_i=1^nc_i^(2)|X_i|^2)   =∑_i≠ jCov(X_i X_j^*,c_i^(2)|X_i|^2+c_j^(2)|X_j|^2)
       =∑_i≠ jη^4I_i I_j(c_i^(2)I_i^2+c_j^(2)I_j^2)(3|f(k_a-k_b)|^2+|f(k_a-k_b)|^4),
    
    D(∑_i=1^nc_i^(2)|X_i|^2)   =∑_i=1^nD(c_i^(2)|X_i|^2)
       =∑_i=1^nη^4(c_i^(2)I_i^2)^2(3+14|f(k_a-k_b)|^2+3|f(k_a-k_b)|^4)

Substituting them into Eq. (<ref>) yields the expression of σ(c_ab^(2)).



§ THEORY PART III. DATA ACQUISITION SIMULATION OF THE ACTIVE INTENSITY INTERFEROMETRY.



This chapter will discuss how to construct an algorithm for simulating c_ab^(2) sampling based on the aformentioned statistical optical model. For simplicity, we take all physical quantities to be dimensionless quantities, and only simulate the case where all laser emitters are symmetrically equivalent, that is, η=1 (dimensionless), I_1=I_2=...=I_n=1 (dimensionless) and c_1^(2)=c_2^(2)=...=c_n^(2)=c (autocorrelation coefficent). Substituting them into Eq. (<ref>), we have


    c_ab^(2)=|∑_i=1^nX_i|^2+(c-1)∑_i=1^n|X_i|^2/n^2,


where X_i=T(k_i,k_a)T^*(k_i,k_b) for any i=1∼ n and all the T(k_i,k_a,b)∼𝒞𝒩(0,1) are complex standard normal random variables. We futher ignore the optical memory effect of the scattering sample, then according to the discussion in the previous chapter, all the X_i are independent random variables. In this case, the algorithm only needs to independently sample the random variable X a total of n times to obtain X_1∼ X_n, and then generate a simulating sample of the c_ab^(2).

Now we discuss how to construct an algorithm to sample X. We rewrite X as X=c_1c_2^*, where c_1∼𝒞𝒩(0,1) and c_2∼𝒞𝒩(0,1) satisfy ⟨ c_1 c_2^*⟩_e=f(k_a-k_b). Using the linear combination of circular complex Gaussian random variables still satisfying the properties of circular complex Gaussian distribution, we can sample c_1 and c_2 in the following way:


    {[ c_1   =√(1-|f(k_a-k_b)|^2)z_1+f(k_a-k_b)z_2 
    
    c_2   =z_2 ].
    .


Here, z_1∼𝒞𝒩(0,1) and z_2∼𝒞𝒩(0,1) are independent of each other. On a computer that can only executes the programs with real numbers, the program can independently sample x∼𝒩(0,1/2) and y∼𝒩(0,1/2) and let z=x+iy be used as a sampling of a complex standard normal random variable. Combining all of the above, we finally construct an algorithm for sampling  c_ab^(2).



§ ESTIMATION AND OPTIMIZATION OF SNR IN COINCIDENCE COUNTING MEASUREMENT


In addition to the speckle-like noise discussed in the previous sections, for a practical intensity interferometer, the randomness of the coincidence counting measurement will also introduce noise into the intensity interference signal. For classical stellar intensity interferometers, the SNR has been shown to be estimated as<cit.>

    SNR∼ rΔτ√(T/Δ t),

where r is the count rate of the photon detectors (number of photons per unit time), Δτ is the coherence time of the light source, Δ t is the time bin for coincidence counting (or the time resolution of the photon detectors) and T is the integration time. Eq. (<ref>) uses two important assumptions in the derivation process, one is that the time resolution of the photon detectors are much larger than the coherence time of the light source, that is, Δ t≫Δτ, and the other is the Poisson noise approximation when the photon count rate is low, that is, rΔτ≪1. These two assumptions are generally valid for stellar intensity interferometers whose observation targets are weak and broad-spectrum light sources.

However, for active intensity interferometry, both assumptions may not hold because of the use of pseudothermal illumination. Let us first consider only the case where the first assumption does not hold. When Δ t≫Δτ, the visibility of the intensity interference signal is of Δτ/Δ t magnitude, but when Δ t≪Δτ, the maximum visility will be exactly 1. So reconsidering the derivation process of Eq. (<ref>) in the origin article<cit.>, when Δ t≪Δτ, the SNR need to be estimated as


    SNR∼ r√(TΔ t).


Combining the above two different boundary conditions, it can be found that the SNR is maximized when Δ t and Δτ are of the same order of magnitude.

In addition, when the count rate increases, the SNR will also increase proportionally, but it cannot be infinitely improved. When r^2ΔτΔ t≫ 1, on average, r^2ΔτΔ t coincidence events will occur every Δτ time, and they are almost all correlated. That is to say, every r^2ΔτΔ t consecutive coincidence events can only be approximately combined into one “effective” coincidence event. Therefore, the SNR given by Eq. (<ref>) or Eq. (<ref>) is overestimated by a factor of r√(ΔτΔ t) considering that the SNR is proportional to the square root of the total coincidence rate. Ultimately, the theoretical maximum value of SNR can be expressed as


    max(SNR)∼{[                   r√(TΔτ)         Δ t≈Δτ when rΔτ≪1;                  √(T/Δτ). 1/r^2Δτ≪Δ t≪Δτ when rΔτ≫1 ].
    .


In the experiment, we have time bin Δ t=2 ms, integration time T=2 s, count rate r≈ 10^4 Hz and Δτ is the atmospheric coherence time, which is measured to be about 27 ms as shown in Fig. <ref>. The value of Δ t satisfies the condition of the maximizing the SNR given by Eq. (<ref>), and it is calculated that SNR∼10.




