Accurate detection of arbitrary photon statistics
Abstract
We report a measurement workflow free of systematic errors consisting of a reconfigurable photonnumberresolving detector, custom electronic circuitry, and faithful dataprocessing algorithm. We achieve unprecedentedly accurate measurement of various photonnumber distributions going beyond the number of detection channels with average fidelity 0.998, where the error is contributed primarily by the sources themselves. Mean numbers of photons cover values up to 20 and faithful autocorrelation measurements range from to 2. We successfully detect chaotic, classical, nonclassical, nonGaussian, and negativeWignerfunction light. Our results open new paths for optical technologies by providing full access to the photonnumber information.
pacs:
The probability distribution of the number of photons in an optical mode carries a great deal of information about physical processes that generate or transform the optical signal. Along with modal structure and coherence, the statistics provides full description of light. Precise characterization of photon statistics is a crucial requirement for many applications in the field of photonic quantum technology O’Brien et al. (2009) such as quantum metrology Matthews et al. (2016); Slussarenko et al. (2017), nonclassical light preparation Harder et al. (2016); Straka et al. (2018a), quantum secure communication Dynes et al. (2018), and photonic quantum simulations Vidrighin et al. (2016); Cottet et al. (2017). Measurement of statistical properties and nonclassical feautures of light also represents enabling technology for many emerging biomedical imaging and particletracking techniques Ta et al. (2015); Israel et al. (2017); Sahl et al. (2017). Intensity correlation is routinely applied to quantify the nonclassicality of light Kimble et al. (1977); Kondakci et al. (2016). Obtaining photon statistics requires repeated measurements using a photonnumberresolving detector (PNRD). The important parameters of PNRDs are dynamic range, speed and accuracy.
The main result of our work is the PNRD design with dataprocessing workflow based on expectationmaximizationentropy, which yields unprecedented accuracy across dozens of tested optical signals ranging from a highly subPoissonian singlephoton state to superPoissonian thermal light, with nonnegligible multiphoton content up to . The accuracy is achieved despite leaving all systematic errors uncorrected and operating with raw data. The proposed method also provides faithful values for states, where the commonly used Hanbury Brown–Twiss measurement would fail due to high multiphoton content Stensson and Björk (2018).
We demonstrate the accuracy of the reported PNRD by performing photon statistics measurement for many different states of light, from which 25 states are shown in Fig. 1, covering various mean photon numbers and values. Furthermore, the reconfigurability of the presented PNRD also allows for direct measurement of correlation functions and nonclassicality witnesses Sperling et al. (2017); Straka et al. (2018a).
Contemporary PNRD technologies all rely on multiplexing with the exception of transitionedge detectors Gerrits et al. (2016); Harder et al. (2016) that require temperatures below 100 mK, offer rates of 10100 kHz and suffer from rangeversuscrosstalk compromise. Photonnumber resolution using a superconducting nanowire singlephoton detector (SNSPD) represents an emerging technology Cahall et al. (2017).
The multiplexing approach is based on dividing an input optical signal into multiple onoff detectors Paul et al. (1996); Sperling et al. (2012). Many schemes of temporal and spatial multiplexing have been reported using a few onoff detectors Banaszek and Walmsley (2003); Achilles et al. (2003); Fitch et al. (2003); Řeháček et al. (2003); Tiedau et al. (2018), hundreds of integrated onoff pixels Allevi et al. (2010a); Kalashnikov et al. (2011); KrÃ¶ger et al. (2017), or even a few photonnumber resolving detectors Calkins et al. (2013); Harder et al. (2016). Though being economical in respect to the number of onoff detectors employed, the temporal multiplexed scheme trades off a decrease of the detector speed for an increase in a number of the detection channels. Decreasing the losses and the balancing of temporal multiplexers require a great deal of optimization Mičuda et al. (2008) or even active signal switching Tiedau et al. (2018). On the other hand, multiplepixel PNRDs typically suffer from strong crosstalk effects, which demands for an experimental characterization of the detector Lundeen et al. (2008) and advanced numerical postprocessing to correct for the imperfections Kalashnikov et al. (2011); KrÃ¶ger et al. (2017). Also, the multipixel detectors offer very limited reconfigurability and complicate channel balancing. Recent onchip integration of independent onoff cryogenic detectors represents a promising direction Rath et al. (2015); Mattioli et al. (2016); Zhu et al. (2018), which has yet to be tested for various classical and, particularly, nonclassical sources.
The reported photonnumberresolving detector is based on spatial multiplexing of the input photonic signal in a reconfigurable optical network as depicted in Fig. 2(a). The multiport network consists of cascaded tunable beam splitters, which allow for accurate balancing of the output ports or, if needed, changing their number so there is no need to physically add or remove detectors. The whole network works as 1to splitter balanced to a few per mill. To measure the multiplexed signal we use singlephoton avalanche photodiodes (SPADs) with efficiency close to 70%, 250 ps timing jitter, and 25 ns recovery time. The electronic outputs of the SPADs are processed by a custom coincidence logic consisting of fast discriminators, delay lines, a summation circuit, and analogtodigital sampling. Each of the resulting distinct voltage levels corresponds to the particular number of fold coincidences. Repeated measurements give rise to click statistics.
It is important to stress here that the PNRD operates in real time and yields a result for every single input pulse with a latency lower than 30 ns (including the response of the SPADs), which allows its application also as a communication receiver, quantum discrimination device, or for a feedback operation. The use of independent detectors and wellbalanced coincidence circuitry completely removes any crosstalk between the histogram channels, see Fig. 2(a). The effects of dark counts and afterpulses are virtually eliminated by operating the detector in pulse regime with the repetition rate below 5 MHz. The period between individual measurement runs can be ultimately decreased to be only slightly longer than the recovery time of the constituent singlephoton detectors, provided that afterpulsing is low enough. Furthermore, differences in SPAD efficiencies and other optical imperfections or imbalances of the PNRD can be arbitrarily compensated by adjusting the splitting ratios of the optical network. This means that all systematic errors are eliminated either by design or by a sufficiently precise adjustment.
The theoretical click statistics of the balanced PNRD with the overall efficiency and output ports can be computed from the photon statistics of the input optical signal using
(1) 
where is a conditional probability matrix, describing the probability of detected simultaneous clicks conditioned by incident photons, assuming a fixed efficiency that was established by a separate measurement Paul et al. (1996); Fitch et al. (2003); Sperling et al. (2012).
Finding the photon statistics , , that satisfies the system of equations (1) for a measured click statistics , , represents the core problem of photon statistics retrieval. This generally illposed problem suffers from underdetermination and sampling error. Fortunately, we have additional constraints facilitating the retrieval, i.e. the photonnumber probabilities are nonnegative, normalized, and typically nonnegligible only within a finite range.
Here we present a novel method, termed expectationmaximizationentropy (EME) method, based on expectationmaximization iterative algorithm weakly regularized by maximumentropy principle. The th iteration of EME is
(2) 
(3) 
Parameter scales the entropy regularization relative to the likelihood maximization, with the optimum value of independent of the photon statistics. The initial iteration is chosen to contain no prior information about the statistics, , with sufficiently large . The process is stopped when two subsequent iterations are practically identical. The retrieved statistics does not change for different initial iterations.
To show the robustness of the novel EME method, a numerical analysis of the retrieval accuracy and convergence was performed for several different photonic states. We compared the EME method with other frequently used algorithms—direct inversion and maximum likelihood approach. EME was found to be a unique estimator that guarantees nonnegativity and the absence of numerical artifacts in the retrieved photon statistics. Total variation distance between the retrieved distribution and the true one is , one order of magnitude smaller than in the case of direct inversion and maximumlikelihood approaches. Numerical simulations yield average fidelity values using the EME algorithm and using the maximumlikelihood approach. The fidelity, defined as , cannot be evaluated for direct inversion due to negative values of estimated photon statistics. Furthermore, the EME algorithm converges 10–1000 faster compared to the maximum likelihood approach, while yielding significantly better results. More details on the EME method are included in the supplemental material.
In our experimental demonstration, we used a balanced 11channel configuration of the detector. We analyzed coherent states, thermal states, multimode thermal states, singlephoton and multiplephoton subtracted thermal states, and nonclassical multiphoton states. Furthermore, we have varied the mean number of photons, the number of modes, and the number of subtracted or superimposed photons. For each retrieved photon statistics we computed , , and other quantities presented in detail in the supplemental material.
The measurements were performed using 1nslong optical pulses with the repetition rate of few MHz. We prepared the initial coherent signal by using a gainswitched laser diode (810 nm) driven by a custom circuitry with fewpermill power stability. The resulting coherent pulses measured by the PNRD show almost perfect Poissonian statistics with up to . The thermal state is generated by temporal intensity modulation of the initial coherent light by a rotating ground glass. The scattered light is collected using a singlemode optical fiber. We measured almost ideal BoseEinstein photon statistics depicted in Fig. 3(a) with up to , . We varied the number of the collected thermal modes, which yielded a signal governed by MandelRice statistics, going from BoseEinstein to Poisson distribution as the number of modes increased. Multiplephoton subtraction from the thermal state was implemented using a beam splitter with a reflectance, see Fig. 2(b). When a (multi)coincidence was detected by a multichannel singlephoton detector in the reflected port, the heralded optical signal in the transmitted port was analyzed by the reported PNRD. A typical result of 2fold subtraction is shown in Fig. 3(b). Increasing the number of subtracted photons results in a transition from superPoissonian chaotic light to Poissonian coherent state Allevi et al. (2010b); Zhai et al. (2013).
Furthermore, we generated multiphoton states by mixing incoherently several singlephoton states from spontaneous parametric downconversion using time multiplexing, see Fig. 2(c). successive time windows, where a single photon was heralded, were merged into a single temporal detection mode. This source emulates the collective emission from identical independent single emitters Ta et al. (2015); Sahl et al. (2017); Straka et al. (2018a). The resulting photon statistics measured for these highlynonclassical multiphoton states corresponds extremely well to the ideal attenuated photon states, see Fig. 3(c,d) for and . Also the parameter computed from the measured photon statistics perfectly agrees with the theoretical model , see Fig. 1.
We utilize fidelity and total variation distance to compare the measured distribution with the ideal one. The worst and the best fidelities and are reached across all the tested sources with average fidelity being . The average distances are for all the tested sources, for nonclassical multiphoton states, and for coherent states up to . For detailed data, see supplemental material. The errors are caused by slight imbalances of splitting ratios in the PNRD, variations in PNRD efficiency , and imperfections of the tested sources, which renders the actual accuracy of the PNRD even higher. Particularly, accurate preparation and characterization of thermal and superchaotic states are highly nontrivial tasks subject to ongoing research Allevi and Bondani (2015); Kondakci et al. (2016); Mika et al. (2018); Manceau et al. (2018).
To conclude, we have reported a fully reconfigurable nearideal photonnumberresolving detection scheme with custom electronic processing and a novel EME photon statistics retrieval method. The PNRD design is free of systematic errors, which are either negligible or can be arbitrarily decreased by the user. We have demonstrated exceptional accuracy of detected photon statistics that goes beyond the conventional limit of the number of PNRD channels. We measured dozens of various photonic sources ranging from highly nonclassical quantum states of light to chaotic optical signals. The results were obtained from raw data with no other processing than EME, presuming a fixed efficiency . Despite uncorrected systematic errors and significant variability of the input signal, our approach shows superior fidelity across the board with typical values exceeding for mean photon numbers up to 20 and the parameter reaching down to a fraction of a percent. Though having been demonstrated with common singlephoton avalanche diodes, the reported measurement workflow can also accommodate superconducting singlephoton detectors ready for onchip integration to achieve a further increase in speed and efficiency Mattioli et al. (2016); Zhu et al. (2018); Münzberg et al. (2018).
Acknowledgments
We acknowledge the support from the Czech Science Foundation under the project 1726143S. This work has received national funding from the MEYS and the funding from European Union’s Horizon 2020 (20142020) research and innovation framework programme under grant agreement No 731473. Project HYPERUPS has received funding from the QuantERA ERANET Cofund in Quantum Technologies implemented within the European Union’s Horizon 2020 Programme. We also acknowledge the support from Palacký University (project IGAPrF2018010).
Supplementary Material
Appendix A Characterization of quantum states of light with focus on their photon statistics
The photon statistics retrieved from a click statistics recorded by the photonnumberresolving detector (PNRD) can be analyzed and compared with a theoretical model but very often only a limited number of characteristics are evaluated and utilized to witness a particular feature of the source under the test. The features crucial for fundamental research as well as many photonic applications are nonclassicality and a deviation from Poisson statistics. The normalized secondorder intensity correlation is routinely applied to quantify the nonclassicality of light Kimble et al. (1977); Kondakci et al. (2016). The departure of the photon statistics from Poisson distribution is frequently characterized by the Mandel parameter Mandel (1979); O’Reilly and OlayaCastro (2014). Evaluation of these parameters requires the full knowledge of the photon statistics, though an approximate value of can be measured directly for Grangier et al. (1986); Stensson and Björk (2018). Both characteristics depend on the first two moments of photon statistics, particularly the mean photon number and the second moment , but behave differently in the presence of losses. The normalized secondorder correlation is completely loss independent while the Mandel parameter changes with the applied losses. For example, the multiphoton states emitted by clusters of singlephoton emitters would ideally display , however, the real sources with limited collection efficiency show the value closer to zero. In our case of the emission emulated by merging of several heralded photons from parametric downconversion process, the measured value is given by heralding efficiency of the source.
Quantum nonGaussian (QNG) states represent another class of strongly nonclassical states, which cannot be expressed as a mixture of arbitrary Gaussian states. We were able to prove QNG of the measured multiphoton states up to nine photons Straka et al. (2018a) using reconfigurability of our PNRD detector and applying recently devised clickbased QNG criteria Lachman and Filip (2016). Having the full photon statistics retrieved, we can also apply QNG tests based on state properties in phase space Genoni et al. (2013); Park et al. (2015, 2017); Happ et al. (2018). Unfortunately, they rely on a singlemode assumption, which does not hold for our multiphoton states emulating a superimposed emission from clusters of singlephoton emitters. Anyway, we have successfully proved QNG property for the emission from the cluster of 1, 2, and 3 emitters using the criterion by Genoni et al. Genoni et al. (2013). QNG criteria feature strong sensitivity to losses compared to nonclassicality measures Straka et al. (2014).
Even higher loss sensitivity can be observed using the mean value of parity operator and the Wigner function , where and is a wave function of the th Fock state. The negativity of follows the negativity of Wigner function at the origin of phase space as , which we have detected for the singlephoton state (emission from a single emitter) and the cluster of three single emitters.
In the main text, we focus mainly on the discrepancy between the measured and the corresponding ideal photon statistics for dozens of measured photonic sources, characterized by the fidelity and the total variational distance , both ranging from 0 to 1. Higher fidelity does not necessarily correspond to a stronger nonclassical feature so other characteristics should also be evaluated Mandarino et al. (2016). We demonstrate the accuracy and wide applicability of the PNRD by plotting parameter computed from the measured photon statistics and from the corresponding ideal statistics for 25 various optical signals, see Fig. 1 of the main text. Here in the Fig. 4 of Supplementary material we demonstrate a similar characterization utilizing the Mandel parameter. We can see the exceptionally accurate experimental characterization of the deviation from Poissonian statistics with the Mandel parameter ranging from 0.55 to 5 and the mean photon number exceeding 20. Furthermore, we show all the mentioned moments, nonclassicality parameters, and Wignerfunctionnegativity characteristics for five selected photonic signals, namely the four signals, the photon statistics of which is plotted in Fig. 3 of the main text, and the initial coherent state, see Table 1.
coherent  thermal 

single photon  9photon cluster  
0.9984(9)  0.9978(5)  0.9990(4)  0.9991160(1)  0.999296(2)  
0.021(9)  0.033(3)  0.019(6)  0.0008839(1)  0.00407(1)  
4.95(2)  4.93(4)  5.77(2)  0.554675(2)  5.00786(2)  
4.99(5)  29.4(5)  17.0(2)  0.2487836(3)  2.27821(7)  
1.002(3)  2.01(1)  1.336(3)  0.0057627(9)  0.891156(3)  
0.01(1)  4.97(7)  1.94(2)  0.551478(2)  0.54507(1)  
0.003(8)  0.089(7)  0.011(5)  0.105814(4)  0.000048(2)  
0.001(3)  0.028(2)  0.004(2)  0.033682(1)  0.0000152(5) 
Appendix B Experimental setup of the detector, data acquisition and processing
The reported photonnumberresolving detector (PNRD) is based on spatial multiplexing in a reconfigurable optical network as depicted in Fig. 2(a) of the main text. The multiport optical network consists of tunable beam splitters composed of a halfwave plate and a polarizing beam splitter for accurate adjustments of the splitting ratio. The halfwave plates can be adjusted to split the light equally to any number of output channels, ten in our case, so there is no need to physically add or remove detectors. To measure the multiplexed optical signal we use 10 singlephoton avalanche photodiodes (SPADs, Excelitas) with system efficiency ranging from 60 to 70% at 810 nm, 200300 ps timing jitter, and 2030 ns dead time. Different efficiency of the SPADs are taken into account during balancing, when the overall efficiency of the channel (the product of the transmittance of the particular port and the efficiency of the SPAD sitting in that port) is balanced across all the channels.
The electronic outputs of the SPADs are processed by a custom coincidence logic consisting of 300 MHz discriminators, delay lines, and a summation circuit. We have tested two implementations of the summation circuit with the virtually same results. The first solution employs a sequence of active 250 MHz linear fanin/out units (Phillips Scientific) to produce the total output. While the response of the discriminators can be easily scaled up above 1 GHz, the bandwidth of the active fanin circuits is more challenging to improve. The second tested approach utilizes passive RF summation circuitry (Mini Circuits) with ability to sum dozens of electronic signals within a GHz bandwidth. Each of the resulting distinct voltage levels (eleven in our case) corresponds to the particular number of multicoincidences. Finally, the electronic signal is digitized with 20 GSa/s by a 1.5 GHz oscilloscope operating in a memorysegmentation regime (Teledyne LeCroy). Repeating the measurement many times for the same state of light, one can build a histogram of the coincidence events and the corresponding normalized clickstatistics. The number of channels of the PNRD can be increased to hundreds while keeping the same analog electronic signal processing technique. Alternatively, a fully digital coincidence device can be employed sporting few ns resolution and hundreds of channels using TTL logic in field programmable gate arrays or even sub100 ps resolution using discrete ECL logic with limited number of channels (verified up to 16 channels).
The use of independent detectors and wellbalanced coincidence circuitry removes completely any crosstalk between the histogram channels yielding the perfect energy quanta resolution up to number of the channels used, see Fig. 5. Furthermore, the effects of dark counts and afterpulses are virtually eliminated by operating the detector in pulse regime with the repetition rate below approx. 5 MHz Wayne et al. (2017); Ziarkash et al. (2018). The period between individual measurement runs can be ultimately decreased to be only slightly longer than the recovery time of the constituent singlephoton detectors, as far as afterpulses are negligible or fast decaying like in the case of superconducting nanowire singlephoton detectors Fujiwara et al. (2011); Burenkov et al. (2013); Münzberg et al. (2018). Consequently, the presented PNRD measurements are free of systematic errors such as the channel crosstalk and temporal correlations.
Appendix C Photonic sources characterized using the PNRD detector
An initial coherent nanosecond pulsed light is generated by gainswitched semiconductor laser diode (810 nm) driven by subnanosecond electronic pulse generator with a repetition rate of 2 MHz. We have reached few per mille pulsetopulse stability and similar longterm stability of the mean power by employing custom lownoise lowjitter pulse generator, laser diode selection, its thermal stabilization, and optimization of driving pulse duration and shape. The resulting coherent pulses measured by the PNRD detector display virtually perfect Poissonian statistics.
Thermal states are generated by temporal intensity modulation of the initial coherent light by rotating ground glass (RGG) with a random spatial distribution of speckles Martienssen and Spiller (1964). A singlemode optical fiber is used to collect the scattered light and select a single spatial mode. We have measured the corresponding photon statistics to be nearly ideal BoseEinstein distribution. Let us note that it is quite challenging to reach an ideal thermal state, see Mika et al. (2018) and reference therein. In the initial stage of our experiment, we observed a discrepancy between measured photon statistics of RGG modulated light and the ideal BoseEinstein statistics. It appeared later that the error was caused by a small inhomogeneity of the RGG disk. More consistent results can be attained using a direct programmable modulation of light intensity, which allows for preparation of nearideal thermal state and also an arbitrary photon statistics Straka et al. (2018b).
We have altered the effective number M of thermal modes by changing the size of the speckles collected via singlemode fiber. This was achieved by changing the diameter of the laser spot on the RGG disk and the distance between the disk and the fiber tip. The resulting MandelRice statistics changes from BoseEinstein to Poisson distribution with increasing number of modes. The multimode thermal state shows the largest discrepancy between the measured photon statistics and the corresponding ideal one, which is caused by its relative complicate preparation. The intensity of the initial coherent state, its focusing on the RGG, and the fiber coupling are changed to simultaneously reach the required mean photon number and the variance compatible with the MandelRice statistics, basically verifying , where is number of modes. Also, the RGG used has to be checked first for its roughness homogeneity by producing the ideal chaotic light with the BoseEinstein photon statistics. If the RGG allows generating a superchaotic statistics for any combination of the mentioned hardware parameters, it cannot be straightforwardly used to produce a proper transition from BoseEinstein to Poisson statistics via the multimode MandelRice statistics. Indeed, several modes of superchaotic light can yield the total statistics different from the BoseEinstein distribution but displaying .
Multiphoton subtraction from the thermal state is implemented by splitting the thermal state at a beam splitter with a reflectance . When a (multi)coincidence is detected by a reconfigurable multichannel detector (multiple SPADs) in the reflected port, the heralded optical signal in the transmitted port is analyzed by the reported PNRD. With increasing number of subtractions, the mean photon number of the conditioned output state is times larger than of the original thermal state and the resulting photon statistics converges to a Poisson distribution. The transition follows a similar path in – and – diagrams as multimode thermal light, however, the multiphoton subtracted thermal states are singlemode states with many applications in quantum metrology and quantum thermodynamics Allevi et al. (2010b); Zhai et al. (2013); Vidrighin et al. (2016); Hloušek et al. (2017); Bogdanov et al. (2017); Barnett et al. (2018); Fedorov et al. (2015); Katamadze et al. (2018).
Furthermore, we generate multiphoton states by mixing singlephoton states incoherently, using the process of spontaneous parametric downconversion in a PPKTP crystal and time multiplexing. We took successive time windows, when a single photon was heralded, and joined them into a single temporal detection mode. The resulting photon statistics measured for these highlynonclassical multiphoton states corresponds extremely well with the ideal attenuated photon states up to . The measured mean photon number is lower than the number of superimposed photons due to nonunity efficiency of the source, , which is contributed the efficiency of the heralding detector (65%), singlemodefiber collection efficiency (90%), and spectral filter transmission and other inefficiencies (94%). One might conclude from the – diagram shown in Fig. 4 that nonclassicality of the multiphoton state is reduced with the increasing number of photons superimposed. Photon statistics of these states are very close to a binomial distribution for all and so they are strongly nonclassical and nonGausian, as displayed by the Mandel Q parameter or other advanced criteria Straka et al. (2018a); Qi et al. (2018); Lachman et al. (2018).
We have performed approximately measurement runs to build a click statistics for classical photon states. The measurement uncertainty has been evaluated by repeating the full acquisition ten times. For nonclassical multiphoton states, we have performed the single acquisition with measurement runs and use Monte Carlo simulation for uncertainty evaluation. Monte Carlo method has also been used to quantify the statistical errors of retrieved photon statistics, fidelities, and other parameters of interest.
Appendix D Theoretical description of the PNRD and photon statistics retrieval
A perfectly balanced port PNRD with efficency and free of dark counts, crosstalk, and other imperfections is descibed by a conditional matrix,
(4) 
transforming the photon statistics , , to theoretical click statistics , ,
(5) 
The element describes the probability of detected simultaneous clicks conditioned by incident photons, for Paul et al. (1996); Fitch et al. (2003); Sperling et al. (2012). The structure of the matrix can be easily understood by representing the actual port PNRD with nonunity detection efficiency as port device with a virtual “sink” output with the probability of and equiprobable outputs with the overall probability of . The formulation (4)(5) is valid for any input state of light. The efficiency is inferred from an independent measurements, and is assumed to be constant during the PNRD operation. The simple matrix (4) with the same efficiency is used for photon statistics retrieval for all the sources characterized, without the need for a tomographic characterization of the detector Luis and SánchezSoto (1999); Fiurášek (2001); Lundeen et al. (2008); Řeháček et al. (2010); Natarajan et al. (2013); Altorio et al. (2016); Cooper et al. (2014).
Finding the photon statistics , , that satisfies the system of equations (5) for a particular measured click statistics , , represents a core problem of photon statistics measurement. This inverse problem is illposed because 1. it is obviously underdetermined, 2. the theoretical click probabilities are not available in real measurement as we acquire relative frequencies instead, which sample the true probabilities (sampling noise), and 3. for PNRDs that are not free of systematic errors, other imperfection can be present like imbalance, crosstalk, and temporal correlations. The reported PNRD implementation is almost free of these technical imperfections. The first two issues, however, remain for any PNRD detector, and the photon statistics retrieval has to take them into account. Fortunately, we have additional information facilitating the retrieval, i.e. the photon number probabilities are nonnegative and normalized. The elements of photon statistics are also typically non negligible only within a finite range of photon numbers. Indeed, the classical states of light possess a quickly decaying tail, and nonclassical states such as singlephoton states are actually defined on a finite support. There are many techniques for photon statistics retrieval, direct inversion and maximumlikelihood approach being probably the most frequently employed. In what follows we present the basic ideas of these techniques and present a novel method based on an iteration technique known as expectationmaximization algorithm weakly regularized by maximumentropy principle.
Direct inverse. The retrieval technique based on the direct inversion of the system of linear equations (5) requires setting a cutoff – the maximum photon number , typically equal to the number of PNRD ports. The truncated problem possesses a single solution
(6) 
the nonnegativity of which is not guaranteed, hence not representing a physically sound photon statistics. Here represents the inverse matrix to the conditional matrix . The solution (6) often reaches negative values and artifical oscillations, see Fig. 7. These adverse effects are particularly noticeable in the practical case of nonunity efficiency with the limited number of output ports and the mean photon number of the incident light comparable or higher than . The nonnegativity constraint can be incorporated using linear programming, for instance, which reduces the volume of the domain by the factor of . Also, the cutoff can be increased rendering the problem underdetermined, a pseudoinverse Rao and Mitra (1971); BenIsrael and Greville (2003) of which often diverges or, at least, amplifies a sampling noise. Various regularization methods are used to make these issues less pronounced Tikhonov (1943); Hoerl (1962); Tikhonov et al. (1995); Starkov et al. (2009). Despite all the mentioned issues the direct inverse methods are frequently used due to their speed and widespread implementation in many numerical libraries and computing systems.
Maximumlikelihood method and expectationmaximization algorithm. Another technique to achieve the inversion of the conditional probability matrix is well known maximumlikelihood (ML) principle and the expectationmaximization (EM) algorithm, which provides a robust method for finding a solution (ML estimate) Dempster et al. (1977); Vardi and Lee (1993). The likelihood of measuring the particular data distribution given the input photon statistics and measurement device is is given by the multinomial distribution
(7) 
which is a convex functional defined on a convex set of distributions. The maximization of the likelihood functional yields a single global maximum in the case of or a single plateau of maxima in the case of underdetermined problems. A logarithm of the likelihood is often used instead, which does not change the convexity feature. Also, the normalization condition is incorporated with the help of a Lagrange multiplier ,
(8) 
The zero variation is a necessary condition for an extreme of the likelihood functional,
(9) 
for each , which is equivalent to
(10) 
except at the boundary of the domain where . To include this boundary condition, the extremal equation is formulated as
(11) 
A summation over yields
(12) 
where the constraint and the normalization of the click data have been applied. The functional (8) can be maximized over variables using downhill simplex method or some other standard numerical method Banaszek (1998a); Marsili et al. (2009). To keep the nonnegativity constraint the variables can be parametrized as , the downside of which is even more complicated structure of the loglikelihood function. This approach is straightforward but numerically demanding as the dimension of the problem increases. Alternatively, an iterative solution of the extremal equation (11), which is a form of the EM algorithm, can be carried out as was suggested by Banaszek for the first time Banaszek (1998b); Zambra et al. (2005),
(13) 
The iteration process is started with an initial positive statistics, typically chosen as the uniform distribution, . Then the kernel of the map (13) is evaluated for the initial iteration step, and the first iteration is obtained by the application of the kernel on the initial statistics. The normalization has to be performed if the data are not properly normalized. The iteration process is repeated until the distance between th and th iteration is smaller than some given value,
(14) 
Throughout this work the value is used for all the performed photon statistics retrievals. When sufficient mathematical conditions are fulfilled, the procedure converges to the fixed point of the map (11), i.e. to the maximumlikelihood estimate of photon statistics Dempster et al. (1977); Vardi and Lee (1993).
Expectationmaximizationentropy algorithm. For underdetermined problems, when , the EM algorithm will converge to a particular solution depending on the initial distribution . All the possible solutions reaches the same value of the likelihood (given the data ) and cannot be distinguished by ML principle itself. In such cases, the common strategy is to allow for some kind of regularization or dumping to select the most “simple” solution from the plateau of all ML solutions or, in other words, to prevent overfitting of the data. Entropy characterizes the solution complexity and its maximization reflects minimum prior information. Entropy maximization is frequently used for regularization of inverse problems in various applications like image reconstruction, seismology, and electromagnetic theory Smith and Grandy (1985); Parker (1994), and also in machine learning and quantum state estimation Grandvalet and Bengio (2005); Teo et al. (2011). Adopting this idea, we have applied entropy maximization to EM algorithm to obtain the mostlikely estimate of photon statistics with the largest entropy. The resulting strategy not only offers improved fidelity of the retrieved statistics but also makes the iteration process faster. The derivation of the expectationmaximizationentropy (EME) algorithm is analogous to derivation (8)(13) but the regularized functional is used instead of simple loglikelihood,
(15) 
Parameter scales the entropy regularization relative to the likelihood maximization. Performing variation of the loglikelihoodentropy functional , eliminating the Lagrange multiplier , and rewriting the extremal equation in the iterative form leed us to the EME algorithm
(16) 
The initial iteration is chosen to contain no prior information about the statistics, , and the process is terminated based on the distance (14). We have performed hundreds of photon statistics retrievals using measured data and thousands retrievals based on Monte Carlo simulated data with not a single failure of the EME algorithm convergence. We have also verified that the retrieved photon statistics does not depend on the initial iteration.
Furthermore, we have performed a detailed analysis of accuracy and convergence speed as a functions of the regularization parameter for various photon statistics including strongly nonclassical subPoissonian states. In case of small values of the EME approaches the common EM algorithm and the accuracy and repeatability of the solution decrease. For large the entropy regularization prevails and the solution is less likely to reproduce the data  the accuracy drops. The optimum value of the regularization parameter is found to be independent of the retrieved statistics. The same optimum value is used throughout this work for all the performed photon statistics retrievals (for all tested sources).
Appendix E Comparison of the retrieval algorithms
We performed a numerical analysis comparing EME to other photon statistics retrieval methods (Fig. 7). We numerically simulated click statistics of several known initial states and then applied direct inversion, EM algorithm and EME algorithm. To quantify the match between the real and estimated photon statistics, we used total variation distance and fidelity. The direct inversion method proved to be unsatisfactory, because nonnegativity of the result is not guaranteed and therefore, some results do not represent a valid photon statistics. Those that do, exhibit the distance , which is close to the results of the EM method. The EM algorithm guarantees positivesemidefinite results with average fidelity . The total variation distances are similar to those obtained by direct inversion. An average distance is reached for all tested sources. Finally, the presented EME method gives the best match while always maintaining nonnegativity. The average fidelity and average distance is . Particularly, the total variation distance of this method is smaller by an order of magnitude across all states. The EME therefore significantly improves the results for all kinds of simulated statistics.
The convergence speed analysis of retrieval methods with respect to the measured data is shown in Fig. 8. We compare the convergence of EM and EME in the case of a coherent state, a thermal state, a two photonsubtracted thermal state and a singlephoton emitter. Fig. 8 shows that EME converges faster by orders of magnitude than EM. Only for a singlephoton emitter, both methods are on par (the green lines overlap). While EM usually requires at least iterations, EME can do with less than .
Appendix F Loss budget, multiple beam generation, and onchip integration
The presented measurements have been performed with PNRD detector not optimized for detection efficiency. The overall system efficiency is estimated to be . The efficiency parameter incorporated in photon statistics retrieval is chosen to be to assure that the efficiency of the PNRD model is the same or higher than of the actual PNRD detector used. The same value of the efficiency is used throughout this work for all the performed photon statistics retrievals (for all tested sources).
The nonunity system efficiency is caused by a sequence of five half wave plates and polarizing beam splitters with the total trasmittance of , two lenses and two fiber couplings with the transmittance of , and the efficiency of SPAD detectors ranging from to with average value of . Hence, . The efficiency can be improved by employing lowloss optics (especially polarizing beam splitters), antireflection coated fibers (transmittance 99%), and superconducting nanowire singlephoton detectors (system efficiency 90%). The improved efficiency can reach . Based on the performed numerical simulations we expect that the resulting retrieved photon statistics will be nearly identical to the ones retrieved using the current version of the PNRD detector. The highefficiency detector would find its application mainly in the case of low number of measurement runs and as heralding detector for a preparation of highlynonclassical quantum states.
The efficency can be improved significantly by a coherent detection such as homodyning with closetounity quantum efficiency photodiodes Lvovsky and Raymer (2009). Also, homodyne detection provides full information about quantum state of light, including phase information. Unfortunately, the homodyning requires a proper (frequency adjusted) local oscillator, which is not accessible in many applications and for many sources like solidstate emitters, biomedical samples, and generally all multimode source.
The discrete optical network employed in the reported PNRD feautures full reconfigurability and continuous tunability of splitting ratios, but extents over dozens square decimeters and limits the overal efficiency of the PNRD. There are other ways of producing multiple beams of uniform intensity: diffraction gratings (diffractive beam splitter) Golub (2004); Hermerschmidt et al. (2007), multiplebeam plate splitters Rajadhyaksha and Webb (1995), and fiber splitters and fanouts. These solutions possess limited efficiency and no tunability and reconfigurability. Onchip integration offers a significant reduction in size Heilmann et al. (2016), however, the limited transmittance of a waveguide network and input/output coupling losses represent an issue. The tunability can be reached using interferometer networks with adjustable phase shifters Matthews et al. (2009); Smith et al. (2009); Flamini et al. (2015); Carolan et al. (2015).
References
 O’Brien et al. (2009) J. L. O’Brien, A. Furusawa, and J. Vučković, Nat. Photonics 3, 687 (2009).
 Matthews et al. (2016) J. C. F. Matthews, X.Q. Zhou, H. Cable, P. J. Shadbolt, D. J. Saunders, G. A. Durkin, G. J. Pryde, and J. L. O’Brien, npj Quantum Information 2, 16023 (2016).
 Slussarenko et al. (2017) S. Slussarenko, M. M. Weston, H. M. Chrzanowski, L. K. Shalm, V. B. Verma, S. W. Nam, and G. J. Pryde, Nat. Photonics 11, 700 (2017).
 Harder et al. (2016) G. Harder, T. J. Bartley, A. E. Lita, S. W. Nam, T. Gerrits, and C. Silberhorn, Phys. Rev. Lett. 116 (2016), 10.1103/physrevlett.116.143601.
 Straka et al. (2018a) I. Straka, L. Lachman, J. Hloušek, M. M. M. Miková, M. Ježek, and R. Filip, npj Quantum Inf. 4, 4 (2018a).
 Dynes et al. (2018) J. F. Dynes, M. Lucamarini, K. A. Patel, A. W. Sharpe, M. B. Ward, Z. L. Yuan, and A. J. Shields, Optics Express 26, 22733 (2018).
 Vidrighin et al. (2016) M. D. Vidrighin, O. Dahlsten, M. Barbieri, M. S. Kim, V. Vedral, and I. A. Walmsley, Phys. Rev. Lett. 116 (2016), 10.1103/physrevlett.116.050401.
 Cottet et al. (2017) N. Cottet, S. Jezouin, L. Bretheau, P. CampagneIbarcq, Q. Ficheux, J. Anders, A. Auffèves, R. Azouit, P. Rouchon, and B. Huard, Proc. Natl. Acad. Sci. U.S.A. 114, 7561 (2017).
 Ta et al. (2015) H. Ta, J. Keller, M. Haltmeier, S. K. Saka, J. Schmied, F. Opazo, P. Tinnefeld, A. Munk, and S. W. Hell, Nat. Commun. 6, 7977 (2015).
 Israel et al. (2017) Y. Israel, R. Tenne, D. Oron, and Y. Silberberg, Nat. Commun. 8, 14786 (2017).
 Sahl et al. (2017) S. J. Sahl, S. W. Hell, and S. Jakobs, Nat. Rev. Mol. Cell Biol. 18, 685 (2017).
 Kimble et al. (1977) H. J. Kimble, M. Dagenais, and L. Mandel, Phys. Rev. Lett. 39, 691 (1977).
 Kondakci et al. (2016) H. E. Kondakci, A. Szameit, A. F. Abouraddy, D. N. Christodoulides, and B. E. A. Saleh, Optica 3, 477 (2016).
 Stensson and Björk (2018) K. Stensson and G. Björk, Phys. Rev. A 98, 033812 (2018).
 Sperling et al. (2017) J. Sperling, W. R. Clements, A. Eckstein, M. Moore, J. J. Renema, W. S. Kolthammer, S. W. Nam, A. Lita, T. Gerrits, W. Vogel, G. S. Agarwal, and I. A. Walmsley, Phys. Rev. Lett. 118, 163602 (2017).
 Gerrits et al. (2016) T. Gerrits, A. Lita, B. Calkins, and S. W. Nam, in Quantum Science and Technology (Springer International Publishing, 2016) pp. 31–60.
 Cahall et al. (2017) C. Cahall, K. L. Nicolich, N. T. Islam, G. P. Lafyatis, A. J. Miller, D. J. Gauthier, and J. Kim, Optica 4, 1534 (2017).
 Paul et al. (1996) H. Paul, P. Törmä, T. Kiss, and I. Jex, Phys. Rev. Lett. 76, 2464 (1996).
 Sperling et al. (2012) J. Sperling, W. Vogel, and G. S. Agarwal, Phys. Rev. A 85 (2012), 10.1103/physreva.85.023820.
 Banaszek and Walmsley (2003) K. Banaszek and I. A. Walmsley, Opt. Lett. 28, 52 (2003).
 Achilles et al. (2003) D. Achilles, C. Silberhorn, C. Śliwa, K. Banaszek, and I. A. Walmsley, Opt. Lett. 28, 2387 (2003).
 Fitch et al. (2003) M. J. Fitch, B. C. Jacobs, T. B. Pittman, and J. D. Franson, Phys. Rev. A 68 (2003), 10.1103/physreva.68.043814.
 Řeháček et al. (2003) J. Řeháček, Z. Hradil, O. Haderka, J. J. Peřina, and M. Hamar, Phys. Rev. A 67, 061801(R) (2003).
 Tiedau et al. (2018) J. Tiedau, E. MeyerScott, T. Nitsche, S. Barkhofen, T. J. Bartley, and C. Silberhorn, arXiv:1805.05881 [quantph] (2018).
 Allevi et al. (2010a) A. Allevi, M. Bondani, and A. Andreoni, Optics Letters 35, 1707 (2010a).
 Kalashnikov et al. (2011) D. A. Kalashnikov, S. H. Tan, M. Chekhova, and L. A. Krivitsky, Opt. Express 19, 9352 (2011).
 KrÃ¶ger et al. (2017) J. KrÃ¶ger, T. Ahrens, J. Sperling, W. Vogel, H. Stolz, and B. Hage, J. Phys. B 50, 214003 (2017).
 Calkins et al. (2013) B. Calkins, P. L. Mennea, A. E. Lita, B. J. Metcalf, W. S. Kolthammer, A. LamasLinares, J. B. Spring, P. C. Humphreys, R. P. Mirin, J. C. Gates, P. G. R. Smith, I. A. Walmsley, T. Gerrits, and S. W. Nam, Opt. Express 21, 22657 (2013).
 Mičuda et al. (2008) M. Mičuda, O. Haderka, and M. Ježek, Physical Review A 78 (2008), 10.1103/physreva.78.025804.
 Lundeen et al. (2008) J. S. Lundeen, A. Feito, H. ColdenstrodtRonge, K. L. Pregnell, Ch. Silberhorn, T. C. Ralph, J. Eisert, M. B. Plenio, and I. A. Walmsley, Nat. Phys. 5, 27 (2008).
 Rath et al. (2015) P. Rath, O. Kahl, S. Ferrari, F. Sproll, G. LewesMalandrakis, D. Brink, K. Ilin, M. Siegel, C. Nebel, and W. Pernice, Light Sci. Appl. 4, e338 (2015).
 Mattioli et al. (2016) F. Mattioli, Z. Zhou, A. Gaggero, R. Gaudio, R. Leoni, and A. Fiore, Opt. Express 24, 9067 (2016).
 Zhu et al. (2018) D. Zhu, Q.Y. Zhao, H. Choi, T.J. Lu, A. E. Dane, D. Englund, and K. K. Berggren, Nature Nanotechnology 13, 596 (2018).
 Allevi et al. (2010b) A. Allevi, A. Andreoni, M. Bondani, M. G. Genoni, and S. Olivares, Phys. Rev. A 82, 013816 (2010b).
 Zhai et al. (2013) Y. Zhai, F. E. Becerra, B. L. Glebov, J. Wen, A. E. Lita, B. Calkins, T. Gerrits, J. Fan, S. W. Nam, and A. Migdall, Opt. Lett. 38, 2171 (2013).
 Allevi and Bondani (2015) A. Allevi and M. Bondani, Opt. Lett. 40, 3089 (2015).
 Mika et al. (2018) J. Mika, L. Podhora, L. Lachman, P. Obšil, J. Hloušek, M. Ježek, R. Filip, and L. Slodička, New J. Phys. 20, 093002 (2018).
 Manceau et al. (2018) M. Manceau, K. Y. Spasibko, G. Leuchs, R. Filip, and M. Chekhova, arXiv:1807.11784 [quantph] (2018).
 Münzberg et al. (2018) J. Münzberg, A. Vetter, F. Beutel, W. Hartmann, S. Ferrari, W. H. P. Pernice, and C. Rockstuhl, Optica 5, 658 (2018).
 Mandel (1979) L. Mandel, Opt. Lett. 4, 205 (1979).
 O’Reilly and OlayaCastro (2014) E. J. O’Reilly and A. OlayaCastro, Nat. Commun. 5 (2014), 10.1038/ncomms4012.
 Grangier et al. (1986) P. Grangier, G. Roger, and A. Aspect, EPL 1, 173 (1986).
 Lachman and Filip (2016) L. Lachman and R. Filip, Opt. Express 24, 27352 (2016).
 Genoni et al. (2013) M. G. Genoni, M. L. Palma, T. Tufarelli, S. Olivares, M. S. Kim, and M. G. A. Paris, Phys. Rev. A 87, 062104 (2013).
 Park et al. (2015) J. Park, J. Zhang, J. Lee, S.W. Ji, M. Um, D. Lv, K. Kim, and H. Nha, Phys. Rev. Lett. 114, 190402 (2015).
 Park et al. (2017) J. Park, Y. Lu, J. Lee, Y. Shen, K. Zhang, S. Zhang, M. S. Zubairy, K. Kim, and H. Nha, Proc. Natl. Acad. Sci. U.S.A. 114, 891 (2017).
 Happ et al. (2018) L. Happ, M. A. Efremov, H. Nha, and W. P. Schleich, New J. Phys. 20, 023046 (2018).
 Straka et al. (2014) I. Straka, A. Predojević, T. Huber, L. Lachman, L. Butschek, M. Miková, M. Mičuda, G. S. Solomon, G. Weihs, M. Ježek, and R. Filip, Phys. Rev. Lett. 113, 223603 (2014).
 Mandarino et al. (2016) A. Mandarino, M. Bina, C. Porto, S. Cialdi, S. Olivares, and M. G. A. Paris, Phys. Rev. A 93, 062118 (2016).
 Wayne et al. (2017) M. A. Wayne, J. C. Bienfang, and S. V. Polyakov, Opt. Express 25, 20352 (2017).
 Ziarkash et al. (2018) A. W. Ziarkash, S. K. Joshi, M. Stipčević, and R. Ursin, Scientific Reports 8 (2018), 10.1038/s4159801823398z.
 Fujiwara et al. (2011) M. Fujiwara, A. Tanaka, S. Takahashi, K. Yoshino, Y. Nambu, A. Tajima, S. Miki, T. Yamashita, Z. Wang, A. Tomita, and M. Sasaki, Optics Express 19, 19562 (2011).
 Burenkov et al. (2013) V. Burenkov, H. Xu, B. Qi, R. H. Hadfield, and H.K. Lo, Journal of Applied Physics 113, 213102 (2013).
 Martienssen and Spiller (1964) W. Martienssen and E. Spiller, American Journal of Physics 32, 919 (1964).
 Straka et al. (2018b) I. Straka, J. Mika, and M. Ježek, Opt. Express 26, 8998 (2018b).
 Hloušek et al. (2017) J. Hloušek, M. Ježek, and R. Filip, Sci. Rep. 7, 13046 (2017).
 Bogdanov et al. (2017) Yu. I. Bogdanov, K. G. Katamadze, G. V. Avosopiants, L. V. Belinsky, N. A. Bogdanova, A. A. Kalinkin, and S. P. Kulik, Phys. Rev. A 96, 063803 (2017).
 Barnett et al. (2018) S. M. Barnett, G. Ferenczi, C. R. Gilson, and F. C. Speirits, Phys. Rev. A 98, 013809 (2018).
 Fedorov et al. (2015) I. A. Fedorov, A. E. Ulanov, Y. V. Kurochkin, and A. I. Lvovsky, Optica 2, 112 (2015).
 Katamadze et al. (2018) K. G. Katamadze, G. V. Avosopiants, Y. Bogdanov, and S. P. Kulik, Optica 5, 723 (2018).
 Qi et al. (2018) L. Qi, M. Manceau, A. Cavanna, F. Gumpert, L. Carbone, M. de Vittorio, A. Bramati, E. Giacobino, L. Lachman, R. Filip, and M. Chekhova, New J. Phys. 20, 073013 (2018).
 Lachman et al. (2018) L. Lachman, I. Straka, J. Hloušek, M. Ježek, and R. Filip, arXiv (2018), 1810.02546 .
 Luis and SánchezSoto (1999) A. Luis and L. L. SánchezSoto, Phys. Rev. Lett. 83, 3573 (1999).
 Fiurášek (2001) J. Fiurášek, Phys. Rev. A 64, 024102 (2001).
 Řeháček et al. (2010) J. Řeháček, D. Mogilevtsev, and Z. Hradil, Phys. Rev. Lett. 105, 010402 (2010).
 Natarajan et al. (2013) C. M. Natarajan, L. Zhang, H. ColdenstrodtRonge, G. Donati, S. N. Dorenbos, V. Zwiller, I. A. Walmsley, and R. H. Hadfield, Opt. Express 21, 893 (2013).
 Altorio et al. (2016) M. Altorio, M. G. Genoni, F. Somma, and M. Barbieri, Phys. Rev. Lett. 116, 100802 (2016).
 Cooper et al. (2014) M. Cooper, M. Karpiński, and B. J. Smith, Nat. Commun. 5, 4332 (2014).
 Rao and Mitra (1971) C. R. Rao and S. K. Mitra, Generalized Inverse of Matrices and Its Applications (Wiley, 1971).
 BenIsrael and Greville (2003) A. BenIsrael and T. N. E. Greville, Generalized Inverses – Theory and Applications (SpringerVerlag, 2003).
 Tikhonov (1943) A. N. Tikhonov, Doklady Akademii Nauk SSSR 39, 195 (1943).
 Hoerl (1962) A. E. Hoerl, Chemical Engineering Progress 58, 54 (1962).
 Tikhonov et al. (1995) A. N. Tikhonov, A. Goncharsky, V. V. Stepanov, and A. G. Yagola, Numerical Methods for the Solution of IllPosed Problems (SpringerVerlag, 1995).
 Starkov et al. (2009) V. N. Starkov, A. A. Semenov, and H. V. Gomonay, Phys. Rev. A 80, 013813 (2009).
 Dempster et al. (1977) A. P. Dempster, N. M. Laird, and D. B. Rubin, J. Royal Stat. Soc. B 39, 1 (1977).
 Vardi and Lee (1993) Y. Vardi and D. Lee, J. Royal Stat. Soc. B 55, 569 (1993).
 Banaszek (1998a) K. Banaszek, Phys. Rev. A 57, 5013 (1998a).
 Marsili et al. (2009) F. Marsili, D. Bitauld, A. Gaggero, S. Jahanmirinejad, R. Leoni, F. Mattioli, and A. Fiore, New J. Phys. 11, 045022 (2009).
 Banaszek (1998b) K. Banaszek, Acta Phys. Slov. 48, 185 (1998b).
 Zambra et al. (2005) G. Zambra, A. Andreoni, M. Bondani, M. Gramegna, M. Genovese, G. Brida, A. Rossi, and M. G. A. Paris, Phys. Rev. Lett. 95, 063602 (2005).
 Smith and Grandy (1985) C. R. Smith and W. T. Grandy, eds., MaximumEntropy and Bayesian Methods in Inverse Problems (Springer Netherlands, 1985).
 Parker (1994) R. L. Parker, Geophysical Inverse Theory (Princeton University Press, 1994).
 Grandvalet and Bengio (2005) Y. Grandvalet and Y. Bengio, in Advances in Neural Information Processing Systems 17, edited by L. K. Saul, Y. Weiss, and L. Bottou (MIT Press, 2005) pp. 529–536.
 Teo et al. (2011) Y. S. Teo, H. Zhu, B.G. Englert, J. Řeháček, and Z. Hradil, Phys. Rev. Lett. 107, 020404 (2011).
 Lvovsky and Raymer (2009) A. I. Lvovsky and M. G. Raymer, Rev. Mod. Phys. 81, 299 (2009).
 Golub (2004) M. A. Golub, Opt. & Photon. News , 37 (2004).
 Hermerschmidt et al. (2007) A. Hermerschmidt, S. Krüger, and G. Wernicke, Opt. Lett. 32, 448 (2007).
 Rajadhyaksha and Webb (1995) M. M. Rajadhyaksha and R. H. Webb, Appl. Opt. 34, 8066 (1995).
 Heilmann et al. (2016) R. Heilmann, J. Sperling, A. PerezLeija, M. Gräfe, M. Heinrich, S. Nolte, W. Vogel, and A. Szameit, Sci. Rep. 6, 19489 (2016).
 Matthews et al. (2009) J. C. F. Matthews, A. Politi, A. Stefanov, and J. L. O’Brien, Nature Photon. 3, 346 (2009).
 Smith et al. (2009) B. J. Smith, D. Kundys, N. ThomasPeter, P. G. R. Smith, and I. A. Walmsley, Opt. Express 17, 13516 (2009).
 Flamini et al. (2015) F. Flamini, L. Magrini, A. S. Rab, N. Spagnolo, V. D’Ambrosio, P. Mataloni, F. Sciarrino, T. Zandrini, A. Crespi, R. Ramponi, and R. Osellame, Light Sci. Appl. 4, e354 (2015).
 Carolan et al. (2015) J. Carolan, C. Harrold, C. Sparrow, E. MartinLopez, N. J. Russell, J. W. Silverstone, P. J. Shadbolt, N. Matsuda, M. Oguma, M. Itoh, G. D. Marshall, M. G. Thompson, J. C. F. Matthews, T. Hashimoto, J. L. O’Brien, and A. Laing, Science 349, 711 (2015).