Accurate detection of arbitrary photon statistics

Accurate detection of arbitrary photon statistics

Josef Hloušek Department of Optics, Palacký University, 17. listopadu 12, 77146 Olomouc, Czech Republic    Michal Dudka Department of Optics, Palacký University, 17. listopadu 12, 77146 Olomouc, Czech Republic    Ivo Straka Department of Optics, Palacký University, 17. listopadu 12, 77146 Olomouc, Czech Republic    Miroslav Ježek Department of Optics, Palacký University, 17. listopadu 12, 77146 Olomouc, Czech Republic

We report a measurement workflow free of systematic errors consisting of a reconfigurable photon-number-resolving detector, custom electronic circuitry, and faithful data-processing algorithm. We achieve unprecedentedly accurate measurement of various photon-number 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, non-classical, non-Gaussian, and negative-Wigner-function light. Our results open new paths for optical technologies by providing full access to the photon-number information.


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), non-classical 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 non-classical feautures of light also represents enabling technology for many emerging biomedical imaging and particle-tracking techniques Ta et al. (2015); Israel et al. (2017); Sahl et al. (2017). Intensity correlation is routinely applied to quantify the non-classicality of light Kimble et al. (1977); Kondakci et al. (2016). Obtaining photon statistics requires repeated measurements using a photon-number-resolving detector (PNRD). The important parameters of PNRDs are dynamic range, speed and accuracy.

The main result of our work is the PNRD design with data-processing workflow based on expectation-maximization-entropy, which yields unprecedented accuracy across dozens of tested optical signals ranging from a highly sub-Poissonian single-photon state to super-Poissonian thermal light, with non-negligible multi-photon 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).

Figure 1: The normalized second-order intensity correlation evaluated from the measured photon statistics (solid marker) and the corresponding ideal statistics (empty marker) of various optical signals with mean photon number . Shown are: coherent states with (blue triangle up), thermal states (also termed chaotic light) with (red circle), -mode thermal states with (violet triangle down), and -photon-subtracted termal states for (black circle). The cases of and coincide with thermal state. Furthermore, the emission from a cluster of single-photon emitters is shown for with .

Contemporary PNRD technologies all rely on multiplexing with the exception of transition-edge detectors Gerrits et al. (2016); Harder et al. (2016) that require temperatures below 100 mK, offer rates of 10-100 kHz and suffer from range-versus-crosstalk compromise. Photon-number resolution using a superconducting nanowire single-photon detector (SNSPD) represents an emerging technology Cahall et al. (2017).

The multiplexing approach is based on dividing an input optical signal into multiple on-off detectors Paul et al. (1996); Sperling et al. (2012). Many schemes of temporal and spatial multiplexing have been reported using a few on-off 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 on-off pixels Allevi et al. (2010a); Kalashnikov et al. (2011); Kröger et al. (2017), or even a few photon-number resolving detectors Calkins et al. (2013); Harder et al. (2016). Though being economical in respect to the number of on-off 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, multiple-pixel PNRDs typically suffer from strong crosstalk effects, which demands for an experimental characterization of the detector Lundeen et al. (2008) and advanced numerical post-processing to correct for the imperfections Kalashnikov et al. (2011); Kröger et al. (2017). Also, the multi-pixel detectors offer very limited reconfigurability and complicate channel balancing. Recent on-chip integration of independent on-off 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, non-classical sources.

The reported photon-number-resolving 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 1-to- splitter balanced to a few per mill. To measure the multiplexed signal we use single-photon 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 analog-to-digital 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 well-balanced 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 single-photon 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.

Figure 2: Experimental schemes of (a) the PNRD and photon statistics retrieval, (b) preparation of multiple-photon subtracted thermal state, and (c) multi-photon light source.

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


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 ill-posed problem suffers from underdetermination and sampling error. Fortunately, we have additional constraints facilitating the retrieval, i.e. the photon-number probabilities are non-negative, normalized, and typically non-negligible only within a finite range.

Here we present a novel method, termed expectation-maximization-entropy (EME) method, based on expectation-maximization iterative algorithm weakly regularized by maximum-entropy principle. The th iteration of EME is


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 non-negativity 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 maximum-likelihood approaches. Numerical simulations yield average fidelity values using the EME algorithm and using the maximum-likelihood 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 11-channel configuration of the detector. We analyzed coherent states, thermal states, multi-mode thermal states, single-photon and multiple-photon subtracted thermal states, and non-classical 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 1-ns-long optical pulses with the repetition rate of few MHz. We prepared the initial coherent signal by using a gain-switched laser diode (810 nm) driven by a custom circuitry with few-per-mill 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 single-mode optical fiber. We measured almost ideal Bose-Einstein 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 Mandel-Rice statistics, going from Bose-Einstein to Poisson distribution as the number of modes increased. Multiple-photon 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 single-photon detector in the reflected port, the heralded optical signal in the transmitted port was analyzed by the reported PNRD. A typical result of 2-fold subtraction is shown in Fig. 3(b). Increasing the number of subtracted photons results in a transition from super-Poissonian chaotic light to Poissonian coherent state Allevi et al. (2010b); Zhai et al. (2013).

Furthermore, we generated multi-photon states by mixing incoherently several single-photon states from spontaneous parametric down-conversion 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 highly-nonclassical multi-photon 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 non-classical multi-photon 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 super-chaotic 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 near-ideal photon-number-resolving 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 non-classical 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 single-photon avalanche diodes, the reported measurement workflow can also accommodate superconducting single-photon detectors ready for on-chip integration to achieve a further increase in speed and efficiency Mattioli et al. (2016); Zhu et al. (2018); Münzberg et al. (2018).

Figure 3: Measured (blue bar) and the corresponding theoretical photon statistics (green dot) for (a) thermal state with , (b) 2-photon subtracted thermal state, (c) single-photon and (d) 9-photon states emulating emission from clusters of single-photon emitters. Inset: Wigner function evaluated from the measured statistics. Note that data agree with theory even beyond the number of channels of the PNRD.


We acknowledge the support from the Czech Science Foundation under the project 17-26143S. This work has received national funding from the MEYS and the funding from European Union’s Horizon 2020 (2014-2020) research and innovation framework programme under grant agreement No 731473. Project HYPER-U-P-S has received funding from the QuantERA ERA-NET Cofund in Quantum Technologies implemented within the European Union’s Horizon 2020 Programme. We also acknowledge the support from Palacký University (project IGA-PrF-2018-010).

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 photon-number-resolving 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 non-classicality and a deviation from Poisson statistics. The normalized second-order intensity correlation is routinely applied to quantify the non-classicality 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 Olaya-Castro (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 second-order correlation is completely loss independent while the Mandel parameter changes with the applied losses. For example, the multi-photon states emitted by clusters of single-photon 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 non-Gaussian (QNG) states represent another class of strongly non-classical states, which cannot be expressed as a mixture of arbitrary Gaussian states. We were able to prove QNG of the measured multi-photon states up to nine photons Straka et al. (2018a) using reconfigurability of our PNRD detector and applying recently devised click-based 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 single-mode assumption, which does not hold for our multi-photon states emulating a superimposed emission from clusters of single-photon 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 non-classicality 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 single-photon 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 non-classical 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, non-classicality parameters, and Wigner-function-negativity 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.

Figure 4: The left sub-plot is reproduced from the main text: the normalized second-order intensity correlation evaluated from the measured photon statistics (solid marker) and the corresponding ideal statistics (empty marker) of various optical signals with mean photon number . Shown are: coherent states with (blue triangle up), thermal states (also termed chaotic light) with (red circle), -mode thermal states with (violet triangle down), and -photon-subtracted termal states for (black circle). The cases of and coincides with thermal state. Furthermore, the emission from a cluster of single-photon emitters is shown for with . The right sub-plot shows the Mandel parameter evaluated from the measured photon statistics of the same optical signals with the same marker/color coding. Coherent states are compatible with Poisson statistics , thermal states show super-Poissonian statistics with , multi-mode thermal states and multi-photon-subtracted termal states converges to Poisson distribution with increased number of modes/subtractions. Highly nonclassical -photon states reach given by the limited efficiency of the single-photon source employed for their generation. The error bars are typically smaller than the symbol size.
coherent thermal
subtracted thermal
single photon 9-photon 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)
Table 1: Characteristics of measured photon statistics for selected states: coherent state, thermal state, 2-photon subtracted thermal state, 1-photon state, and 9-photon mixture. Standard deviations are evaluated by repeating ten times the whole process of the PNRD measurement, photon statistics retrieval, and characteristics evaluation. The amount of data acquired for single-emitter clusters were several orders of magnitude larger that for other states, which yields the corresponding error bars much smaller.

Appendix B Experimental setup of the detector, data acquisition and processing

The reported photon-number-resolving 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 half-wave plate and a polarizing beam splitter for accurate adjustments of the splitting ratio. The half-wave 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 single-photon avalanche photodiodes (SPADs, Excelitas) with system efficiency ranging from 60 to 70% at 810 nm, 200-300 ps timing jitter, and 20-30 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 fan-in/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 fan-in 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 multi-coincidences. Finally, the electronic signal is digitized with 20 GSa/s by a 1.5 GHz oscilloscope operating in a memory-segmentation 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 click-statistics. 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 sub-100 ps resolution using discrete ECL logic with limited number of channels (verified up to 16 channels).

The use of independent detectors and well-balanced 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 single-photon detectors, as far as afterpulses are negligible or fast decaying like in the case of superconducting nanowire single-photon 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.

Figure 5: Pulse-height spectrum of the analog output of the reported PNRD for coherent states of various mean photon number 1 (a), 5 (b), and 20 (c). The spectra are plotted in log-scale to emphasize the perfect photon-number resolution and the absence of any crosstalk effects or background noise.

Appendix C Photonic sources characterized using the PNRD detector

An initial coherent nanosecond pulsed light is generated by gain-switched semiconductor laser diode (810 nm) driven by sub-nanosecond electronic pulse generator with a repetition rate of 2 MHz. We have reached few per mille pulse-to-pulse stability and similar long-term stability of the mean power by employing custom low-noise low-jitter 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 single-mode 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 Bose-Einstein 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 Bose-Einstein 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 near-ideal 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 single-mode 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 Mandel-Rice statistics changes from Bose-Einstein to Poisson distribution with increasing number of modes. The multi-mode 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 Mandel-Rice 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 Bose-Einstein photon statistics. If the RGG allows generating a super-chaotic statistics for any combination of the mentioned hardware parameters, it cannot be straightforwardly used to produce a proper transition from Bose-Einstein to Poisson statistics via the multi-mode Mandel-Rice statistics. Indeed, several modes of super-chaotic light can yield the total statistics different from the Bose-Einstein distribution but displaying .

Multi-photon 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 multi-channel 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 multi-mode thermal light, however, the multi-photon subtracted thermal states are single-mode 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 multi-photon states by mixing single-photon states incoherently, using the process of spontaneous parametric down-conversion 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 highly-nonclassical multi-photon 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 non-unity efficiency of the source, , which is contributed the efficiency of the heralding detector (65%), single-mode-fiber collection efficiency (90%), and spectral filter transmission and other inefficiencies (94%). One might conclude from the diagram shown in Fig. 4 that non-classicality of the multi-photon 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 non-classical and non-Gausian, 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 non-classical multi-photon 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,


transforming the photon statistics , , to theoretical click statistics , ,


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 non-unity 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ánchez-Soto (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 ill-posed 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 non-negative 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 single-photon states are actually defined on a finite support. There are many techniques for photon statistics retrieval, direct inversion and maximum-likelihood 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 expectation-maximization algorithm weakly regularized by maximum-entropy principle.

Direct inverse. The retrieval technique based on the direct inversion of the system of linear equations (5) requires setting a cut-off – the maximum photon number , typically equal to the number of PNRD ports. The truncated problem possesses a single solution


the non-negativity 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 non-unity efficiency with the limited number of output ports and the mean photon number of the incident light comparable or higher than . The non-negativity constraint can be incorporated using linear programming, for instance, which reduces the volume of the domain by the factor of . Also, the cut-off can be increased rendering the problem underdetermined, a pseudoinverse Rao and Mitra (1971); Ben-Israel 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.

Maximum-likelihood method and expectation-maximization algorithm. Another technique to achieve the inversion of the conditional probability matrix is well known maximum-likelihood (ML) principle and the expectation-maximization (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


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 ,


The zero variation is a necessary condition for an extreme of the likelihood functional,


for each , which is equivalent to


except at the boundary of the domain where . To include this boundary condition, the extremal equation is formulated as


A summation over yields


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 non-negativity constraint the variables can be parametrized as , the downside of which is even more complicated structure of the log-likelihood 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),


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,


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 maximum-likelihood estimate of photon statistics Dempster et al. (1977); Vardi and Lee (1993).

Expectation-maximization-entropy 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 most-likely 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 expectation-maximization-entropy (EME) algorithm is analogous to derivation (8)-(13) but the regularized functional is used instead of simple log-likelihood,


Parameter scales the entropy regularization relative to the likelihood maximization. Performing variation of the log-likelihood-entropy functional , eliminating the Lagrange multiplier , and rewriting the extremal equation in the iterative form leed us to the EME algorithm


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.

Figure 6: Accuracy and speed of photon statistics retrieval using EME algorithm versus the strength of entropy regularization characterized by the parameter . The total variational distance characterizing the accuracy (black) and the number of EME iterations required (blue) are plotted for four different photon statistics. The data are simulated numerically (from true photon statistics) using Monte Carlo approach. For each value of the statistics retrieval is performed several times for different data sets to evaluate the repeability, which is represented by error bars of .

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 non-classical sub-Poissonian 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).

Figure 7: Photon statistics retrieval of various different states of light from numerically simulated click statistics. Each column shows a different retrieval method: direct inversion (Moore-Penrose pseudoinverse), EM algorithm, and EME algorithm, from left to right. The same amount of data (number of measurement runs) is used for each retrieval method to facilitate the comparison. The green points represent the corresponding true photon statistics.

Appendix E Comparison of the retrieval algorithms

Figure 8: Photon statistics retrieval convergence demonstration of selected states: coherent state (blue), thermal state (red), 2-photon subtracted thermal state (black), and single-photon emitter (green). Shown are results for EM iterative process (full lines) based on Eq. (13) and for EME algorithm (dashed lines) based on Eq. (16). For each state five runs of individual retrievals were done.

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 non-negativity 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 positive-semidefinite 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 non-negativity. 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 photon-subtracted thermal state and a single-photon emitter. Fig. 8 shows that EME converges faster by orders of magnitude than EM. Only for a single-photon 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 on-chip 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 non-unity 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 low-loss optics (especially polarizing beam splitters), anti-reflection coated fibers (transmittance 99%), and super-conducting nanowire single-photon 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 high-efficiency detector would find its application mainly in the case of low number of measurement runs and as heralding detector for a preparation of highly-nonclassical quantum states.

The efficency can be improved significantly by a coherent detection such as homodyning with close-to-unity 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 solid-state emitters, biomedical samples, and generally all multi-mode 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), multiple-beam plate splitters Rajadhyaksha and Webb (1995), and fiber splitters and fan-outs. These solutions possess limited efficiency and no tunability and reconfigurability. On-chip 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).


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description