Efficiently reconstructing compound objects by quantum imaging with higher-order correlation functions

Quantum imaging has a potential of enhancing the precision of objects reconstruction by exploiting quantum correlations of the imaging field, in particular for imaging with low-intensity fields up to the level of a few photons. However, it generally leads to nonlinear estimation problems. The complexity of these problems rapidly increases with the number of parameters describing the object and the correlation order. Here we propose a way to drastically reduce the complexity for a wide class of problems. The key point of our approach is to connect the features of the Fisher information with the parametric locality of the problem, and to reconstruct the whole set of parameters stepwise by an efficient iterative inference scheme that is linear on the total number of parameters. This general inference procedure is experimentally applied to quantum near-field imaging with higher-order correlated light sources, resulting in super-resolving reconstruction of grey compound transmission objects. Quantum information provides a powerful method for improving imaging beyond classical limits, but the complexity typically scales non-linearly with the number of parameters. Here, the authors demonstrate a method to reduce the number of parameters required and identify an optimal degree of photon correlation for imaging.

Q uantum metrology exploits the quantum features of measurement schemes to infer the parameters of interest with enhanced precision, compared with classical schemes 1 . In some cases, such as interferometry 2,3 , the number of parameters to be determined is small, and inferring them from measurement results is straightforward. However, for a large number of parameters, the problem of inferring them from measurement results can be quite demanding, even when it is linear. It might require a prohibitive amount of measurement and computational effort to be solved. For example, to reconstruct the state of a moderate number (say, a few dozen) of the simplest quantum objects, qubits, one already needs some simplifying assumptions, such as a low rank of the state [4][5][6][7] , the possibility to approximate the state by a matrix product 8 , or by a permutationally invariant state [9][10][11] . For nonlinear problems, the task is even more difficult 12,13 . A typical case of multiparameter determination is quantum imaging 14 . Images with improved optical resolution can be obtained from higher-order correlation functions of the light field, thanks to correlation in the illumination sources [15][16][17] . For this problem, the number of parameters can be related to the number of pixels in the image to be reconstructed.
Here, we present an efficient method for nonlinear estimation problems that applies for the important class of parametrically localized measurements. For this class, the result of a particular measurement is dependent on a limited subset of parameters only. Such measurements are common for objects consisting of components well separated in physical or phase space. For example, measurements on individual systems in ion traps 18 or optical lattices 19 , direct 20 and near-field imaging 14 , or datapattern tomography 21 fall into this category. For such measurements, we develop an iterative sliding-window method (SWM) by reconstructing on each step only a subset of parameters that can be much smaller than the total number of parameters. The complexity for such an approach depends linearly on the number of times one needs to shift the window to cover the whole parameter set. To establish the use of the SWM, we develop an informational approach for the analysis of the measurement scheme. We apply here the Fisher information matrix (FIM) for the analysis of the problem and for designing the SWM. Nowadays, Fisher information analysis is firmly establishing itself as an operational tool in quantum tomography and imaging schemes [22][23][24][25][26] . We show how the structure of the FIM can be exploited for estimating the size and structure of the parameter subset of the SWM iterations. We demonstrate the efficiency of our method with the practically important problem of imaging with correlated photons by measuring a second-or higher-order correlation function for position-momentum-entangled photons generated by spontaneous parametric downconversion (SPDC) and pseudo-thermal light. We predict the existence of an optimal degree of photon correlations in the imaging field to achieve the best resolution for a given object. The FIM analysis allows us also to uncover the possibility to increase the resolution by using biased estimation with a bias stemming from physical limitations on the set of the problem parameters.

Results
Theoretical background. To elucidate our approach, let us start with the simplest linear measurement model described with the probabilities p k ¼ P M j¼1 C kj θ j , where θ j are the parameters in question and C kj is the square Hermitian measurement matrix. We call the measurement strictly parametrically l-local, if l < M and C kj = 0 for |k − j| > l, i.e., the matrix C is l-banded, and the (l + 1)th and other side diagonals are equal to zero. It means that each probability p j depends on no more than 2l neighboring parameters. The key observation here is the possibility to approximate the inverses of banded matrices with approximately banded matrices 27,28 . It would mean that the estimator of the parameter θ j depends only on probabilities in the vicinity of p j . Such a locality provides the possibility of getting an accurate estimate for some θ h , for example, by minimization of the distance between a set of experimentally measured frequencies, f k , k ∈ [h − J, h + J], where J ≥ l is an interval around h, and the probabilities estimated as p k % P h þ J j¼h À J C kj θ j 8 . Estimation can be performed for a sequence of h, thus, shifting the estimation window along the whole set of parameters. The complexity of the SWM is linear on the number of shifts required to cover the whole set of parameters. Below we elaborate on this possibility. Notice that the consideration given above holds also for nonstrictly parametrically local measurements (Supplementary Note 1). Now let us consider the general nonlinear parametric measurement model p k = C k (θ 1 , …θ M ). Strict parametric l-locality for the case would mean ∂p j ∂θ m ¼ 0 for |m − j| > l. Our suggestion is to estimate the influence of a given change of a particular parameter on the other parameters with help of the FIM and the Cramer-Rao bound (CRB). By assuming the completeness of the measurement set, P k p k ¼ 1, the FIM for this case reads For the unbiased estimate, the CRB connects the elements of the inverse FIM with the variance of the estimators, where N gives the total number of events. A banded structure of the FIM would mean that an error estimate for a particular parameter, θ h , can be influenced by variations of the parameters only in some vicinity [h − J, h + J]. Our suggestion is to use this clue for designing the SWM as it was described above for the linear case. Also, FIM and CRB can be used for optimization of the measurement scheme aiming at lowering the error bounds per given number of measured events, N. One can minimize the bound for the total measurement error described by the trace of the inverse FIM. A banded structure of the FIM gives a clue to the connection between the width of the FIM and the total error: generally, for given diagonal elements of the FIM, increasing the width (i.e., a number and value of bands) leads to an increase in the inverse trace and the total error (Supplementary Note 2). One can also define an empirical Rayleigh criterion for the parameter resolution from the banded structure: when the FIM is strongly diagonally dominant, F jj ) P k≠j jF jk j, then statistical errors for estimation of the parameter θ j are defined mainly by the measured f j , and the parameters can be well estimated by individual measurements. Notice that the diagonal dominance is a quite strong property imposing locality (Supplementary Note 1). For example, for a strictly one-banded diagonally dominant FIM, a lower bound on the variance of θ j is defined by elements F jk with |j − k| ≤ 2 29 .
Imaging with higher-order correlation functions. We illustrate the previous discussion by applying the SWM to practically relevant examples of quantum near-field imaging by means of higher-order correlation functions. Measuring higher-order correlation functions 14,30-32 is one of the ways to increase the resolution of imaging [33][34][35][36][37] and to go beyond the empirical Rayleigh limit 38 . The scheme of the measurement setup is depicted in Fig. 1. The source produces the linearly polarized field with the density matrix ρ. This field passes through the object described by the transmission function Að s ! Þ and then through the imaging system characterized by the point-spread function hð s ! ; r ! Þ.
Afterward, the field goes to the detectors. The field amplitude operator at the object plane, E o ð s ! Þ, relates to the field amplitude operator at the image plane, Eð r ! Þ, in the following way: where hð s ! ; r ! Þ is a PSF describing the field propagation between the object and the image plane 14 . Integration in Eq.
(2) is over the object plane O. We represent the object as a superposition of M pixels Að s where the function d j ð s ! Þ describes the unit transmission through the jth pixel, and x j is the value of the transmission assigned to the jth pixel. We measure the nth-order intensity correlation function, The coefficients D (k) (l 1 …l n ; m 1 …m n ) are defined by the imaging system and the state of the source. In Supplementary Note 3, D (k) are derived for SPDC-entangled photons, and pseudo-thermal states used for experimental implementations.
Sliding-window method for imaging. The problem of the object inference is to find a set of transmission values {x j } fitting the measured data described by the set of frequencies f k in the best way. For the realization of the SWM, we implement the following iterative scheme: in the first step, we define the pixels (the functions d j ð s ! Þ) in such a way that the FIM, i.e., Eq. (1), is strongly diagonally dominant and infers the initial approximation. Then, we divide each initial pixel in a subgroup of smaller pixels, assign to each of them the transmittance of the parent pixel, and calculate the FIM. Next, we define the window to be shifted as some set of adjacent "core" pixels and some "border" pixels around the "core". Thereby, we use the number and relative value of major bands of the FIM for defining the size of the "border" and perform the fitting. Notice that for building the procedure, one does not need to know the object beforehand or to perform some preliminary estimation. The pixel size and the window structure can be defined for the model object and the used imaging setup. The details of the SWM method and the pseudocode are provided in the "Methods" section.
Object inference. Let us illustrate the mechanics of the SWM with transmitting 1D and 2D objects. We take for our examples common "workhorses" of the quantum imaging field: a pseudothermal state 39 and a position-momentum-entangled state produced by SPDC 40 . In Fig. 2, one can see an illustration of the SWM for a 1D object for the simulated G (2) (Supplementary Note 3). Figure 2a, b shows an example of the typical strongly diagonally dominant FIM for a pixel larger than the Rayleigh limit (Fig. 2a), and the FIM with a pixel smaller than the Rayleigh limit (Fig. 2b). However, the FIM of Fig. 2b is still narrowly banded, and thus the inference problem is treatable by the SWM. The rule of thumb here is to choose the size of the border region larger than the number of the major bands of the FIM. Figure 2c shows the object. Figure 2d shows the simulated G (2) for the object of Fig. 2c for the thermal source. The image of the object (Fig. 2c) is shown in Fig. 2d, e. The process of reconstruction by moving the "window" is depicted in Fig. 2f (see the Methods section). The result of the reconstruction is shown in Fig. 2g, h. Notice that for this case, we have come beyond the Rayleigh limit Δl, shown with the red bar in Fig. 2h: the reconstruction result 2g is close to the original object shown in Fig. 2c, while the diagonal part of the image 2e looks differently. The object inference for the higher-order correlation functions can be realized similarly to the procedure described above.
Experiment. The experimental verification of the SWM was done with the particular realizations of a generic measurement scheme depicted in Fig. 1 for both, a pseudo-thermal and a spontaneous downconversion (SPDC) source (see the Methods section). To produce pseudo-thermal light, a rotating ground glass disk was illuminated by a monochromatic laser 39,41 operating at 405 nm. Type-0 position-momentum-entangled two-photon states were generated by a SPDC source 40 . Thereby, we use a 12-mm-long periodically poled potassium titanyl phosphate (PPKTP) nonlinear crystal pumped by a continuous-wave laser centered at 405 nm. The entangled photons are then emitted at 810 nm. Detection at the image plane was done for the pseudo-thermal light by using SuperEllen, a single-photon-sensitive 32 × 32-pixel single-photon avalanche diode (SPAD) array detector manufactured in complementary metal-oxide-semiconductor (CMOS) technology 17,42 ; and by scanning single photon counters for the entangled photons 40 . Figure 3 shows the results of the SWM for experimental data. Figure 3b shows the reconstructed 2D object inferred from the measurement of G Optimization of the imaging state. The informational approach allows us to predict the optimal correlation width of the used illumination source for the object resolution in the superresolution regime. Intuitively, it seems that the smaller the correlation width is, the better the resolution should be. However, the analysis of the collected information shows that for the object inference, perfectly correlated photons might not be the best choice. It follows from an optimization of the lower bound on the total reconstruction error. This prediction is valid for an arbitrary reconstruction method for the measurement of the second-order correlation function with both twin-photon and quasi-thermal imaging source. In Fig. 4a, an example of the optimization is shown for the image reconstruction from a G (2) function for a pseudo-thermal state. The trace of the inverse FIM and the infidelity of the reconstruction are shown for different correlation widths w c . There is an optimal w c allowing to increase the reconstruction quality for the same number of detector counts in the super-resolution regime. One can describe the most optimal state with the following rule of thumb: the correlation width should be close to the smallest object details to be resolved, i.e., to the pixel size. In the super-resolution regime, the measurement of A similar relation between the optimal photon correlation width and the size of the object features also holds for a SPDC source (Supplementary Note 4). For pixel size exceeding the Rayleigh limit, this effect disappears; decreasing the correlation width brings about enhancement of the resolution (Fig. 4c).
Inference bias. The informational approach and the SWM can capture the possibility of a considerable improvement of resolution stemming from constraints imposed on the parameters. For the special case of parameters being on the borders of the allowed regions, the estimation is generally biased. The bias can significantly modify the error bounds 43,44 (see Supplementary Note 5 and Fig. 4d). For the imaging of binary black-and-white objects (i.e., for x j being either 0 or 1; see bottom inset and thick lines in Fig. 4d), one can go far beyond the resolution limit found for gray images (top inset and thin lines in Fig. 4d) even without any prior assumption of the binary object structure. The reason for it is the dependence of the error bounds on the bias derivative with respect to the parameters 43 . Generally, the SWM shifts the estimators near borders. The closer the estimated value is to the border, the larger is the respective shift and the error-bound deviation. Notice that for the object inference demonstrated in Fig. 3, this bias effect was actually seen.

Discussion
We developed an inference method for nonlinear parametrically local problems and showed how the analysis of the information allows one to develop an estimation scheme making the complexity of the problem linear on the total number of parameters. Then, the scheme was applied to the experimental data for super-  Thick lines correspond to black-and-white objects and thin lines to gray objects resolution imaging based on higher-order correlation measurements with nonclassical two-photon and pseudo-thermal states. It was shown how the FIM can be applied to optimize the imaging state for better resolution, in particular, for the correlation width of the twin-photon or pseudo-thermal imaging fields. Generally, the correlation width should be close to the smallest details to be resolved. This prediction is experimentally confirmed for measurements with pseudo-thermal light. It was also demonstrated that bias due to marginal values of estimated parameters can improve the resolution. We believe that the suggested SWM and an information approach for nonlinear inference problems will find applications for the design and optimization of inference schemes in imaging, quantum diagnostics, and tomography.

Methods
The sliding-window method. The practical application of the SWM to the quantum imaging problem consists of the following steps. First, an initial rough estimate of the object transmission amplitude is found. The pixel size is chosen in such a way that the inverse of the FIM for the reconstruction of the object, expressed in terms of pixels of this size, d 0 ð Þ 0 , is diagonally dominant. The problem is strongly local, and a single run of the SWM is sufficient for getting the initial estimate.
Then, the estimate is refined by representing the object in terms of smaller pixels (size d) and applying the reconstruction algorithm again. The pixel size d limits the size of the object features that can be successfully reconstructed, and therefore, determines the achievable resolution.
Here, the pseudocode for the iterative reconstruction is presented. Both algorithms (the one for the first approximation inference and the one for refinement) are mathematically very similar, the only difference is the role of the window border: the first algorithm implies that pixels inside the border are unknown and not reliable (at each step): they need to be reconstructed but then discarded; the second algorithm implies that pixels inside the border are known (at each step) and thus can be included in the optimization problem as known constants (Fig. 5).
The pseudocode for the first approximation inference algorithm reads as follows: For given core window dimensions, compute all core window positions that lead to the object being fully covered by non overlapping core windows. For each position of the core window: 1. For a given border size, build a complete window: core + border. 2. Map the resulting full window to the image plane. 3. Find all detectors inside the obtained window in the image plane. 4. Combine the obtained detectors according to the order of correlations and other constraints if present. 5. For the obtained detector combinations, load the corresponding experimental joint detection frequencies. 6. Build a function that maps pixels inside the full window to residuals between theoretical joint detection probabilities and experimental frequencies. The theoretical probabilities are computed by assuming all pixels but those inside the window to be zero. 7. Perform numerical minimization of the function thus obtained. Pixels inside the full window are subjected to physical constraints. 8. Update object pixels inside the core window with the corresponding values obtained from the minimization procedure. Discard other pixel values. a c b d Fig. 5 An example of two adjacent reconstruction steps for a general two-dimensional case. Panels a, b correspond to obtaining the first approximation, and c, d correspond to iterative refinement. Although pixels are divided after the first approximation is obtained, here the pixel size is chosen to be the same for both cases in order to show mathematical similarity of the first approximation inference and refinement. Letter "X" denotes pixels for which the optimization problem will be stated ("unknown" pixels), digit "0" denotes pixels that will be completely ignored at the iteration ("irrelevant" pixels), and letter "V" denotes pixels whose values will be included in the optimization problem as known constants ("known" pixels). Dashed frame shows the core of the window, solid frame shows the whole window, the core, and the border Because the algorithm first computes all core window positions and then applies one iteration per core window position, it is easily paralleled.
The pseudocode for the iterative refinement algorithm reads as follows: 1. For a given core window position and a given border size, build a complete window: core + border. 2. Map the resulting window to the image plane. 3. Find all detectors inside the obtained window in the image plane. 4. Combine the obtained detectors according to the order of correlations and other constraints if present. 5. For the obtained detector combinations, load the corresponding experimental joint detection frequencies. 6. Build a function that maps pixels inside the core window to the residuals between theoretical joint detection probabilities and experimental frequencies. The theoretical probabilities are computed by assuming pixels outside the core window but inside the full window to be known, they are set to constant values. Pixels outside the full window are set to zero. 7. Perform numerical minimization of the function thus obtained. Pixels inside the core window are subject to physical constraints. 8. Update pixels inside the core window with the values obtained. Move the core window one step further.
Experiments. In the first experiment, we illuminate a rotating ground glass disk (RGGD) with an attenuated, monochromatic laser operating at λ = 405 nm (Fig. 6). An additional lens (L1) in front of the disk allows to vary the beam waist radius at the position of the RGGD. Subsequently, we insert a far-field lens (L2) in a 2f setting (f = 75 mm) in order to collimate the light and remove the spherical wave front given by the point-like source. (The latter has shown to induce distortions in the subsequent imaging setup.) The object plane (OP) is then located in the far field of the source. An object is then imaged onto the image plane (IP) by means of L3 (f = 150 mm), which is additionally endowed with a variable-size pinhole (PH) to control the resolution, i.e., the Rayleigh limit of the setup. The diameter of the PH was fixed to 1.7 mm. The magnification factor m = s i /s o of the imaging system is m = 1.94, whereas the object distance is s o = 234 mm, and the imaging distance is given by s i = 454 mm. At the image plane, photons are detected by SuperEllen, a single-photon-sensitive 32 × 32-pixel SPAD array detector manufactured in CMOS technology with a pixel pitch of 44.64 μm and a fill factor of 19.7% 17,42 . SuperEllen is able to provide frames with a data acquisition window of 30 ns, and a readout time of 10 μs at a frame rate of 800 kHz. The spatial correlations between pixels were evaluated between consecutive frames with a resolution of 10 μs given by the frame separation. This procedure allows for the resolution of the coherence time of the speckles of the order of microseconds. Second-and thirdorder correlation functions are measured with SuperEllen.
The here-presented pseudo-thermal light setup was used to obtain the following two results: first, the digit "5" (Group 2) from a negative U.S. Air Force (USAF) test chart was imaged and then reconstructed from the data of a G (3) function measurement. This is shown in Fig. 3a, b of the main text. Second, a negative USAF chart three-slit pattern (Group 3, Element 2) was imaged from the object plane to the image plane for various correlation widths w c . The latter was modified by changing the distance between the focusing lens L1 and the RGGD, and therefore the beam waist radius. Based on a G (2) measurement, this allowed to demonstrate the dependence of the image reconstruction quality on the correlation width of the source shown in Fig. 4b of the main text.
The setup for imaging with entangled photons is shown in Fig. 7. Our source generates type-0 position-momentum-entangled photon states by pumping a 12-mm-long PPKTP nonlinear crystal (NLC) with a continuous-wave (CW) laser centered at 405 nm 40 . The entangled photons are then emitted at 810 nm. The residual pump beam is subsequently blocked by a long-pass filter (LF), and the subsequent band-pass filter (BF) transmits photons at 810 nm with a spectral full width at half maximum (FWHM) of 10 nm to the detectors. The experimental setup contains two imaging systems: the first system consists of a 4-f image by using lenses L1 and L2 both with focal length f = 50 mm. This configuration maps the entangled photon states transverse momentum distribution from the OP1 at the center of the NLC to the OP with a magnification factor of m = 1. The OP is then imaged with a single-lens system onto the fiber tips of two multimode fibers (MMFs). Thereby, we have a magnification factor of m = 12 for s o = 65 mm and s i = 780 mm. Both detection stages can be scanned in horizontal direction. This setup was used to record the image of a three-slit pattern of a positive USAF resolution chart (Group 4, Element 1) by measuring a second-order correlation function of the photons. The experimental correlation map and the reconstructed object can be seen in Fig. 3c, d of the main text.

Data availability
The datasets generated and analyzed during the current study are available from the depository of the Center for Quantum Optics and Quantum Information of B. I. Stepanov Institute of Physics, National Academy of Sciences of Belarus http://master. basnet.by/Informational-approach-for.data.rar.

Code availability
The source codes used for the sliding-window method are available from the corresponding author upon reasonable request. IP SuperEllen Fig. 6 Pseudo-thermal light imaging setup. A monochromatic laser is focused onto a rotating ground glass disk (RGGD) by means of lens L1. The subsequent lens L2 provides the far-field speckle pattern at the object plane (OP). The resolution of the single-lens (L3) imaging system can be modified by a variable-size pinhole (PH). Single photons are detected at the image plane (IP) by SuperEllen