Precisely determining photon-number in real time

Superconducting transition-edge sensors (TES) are sensitive microcalorimeters used as photon detectors with unparalleled energy resolution. They have found applications from measuring astronomical spectra through to determining the quantum property of photon-number, ˆ n =ˆ a † ˆ a , for energies from 0.6–2.33 eV. However, achieving optimal energy resolution requires considerable data acquisition—on the order of 1 GB/min—followed by post-processing, preventing real-time access to energy information. We report a custom hardware processor to process TES pulses while new detection events are still being registered, allowing the photon-number to be measured in real time. We resolve photon-number up to n =16 —achieving up to parts-per-billion discrimination for low photon-numbers on the fly—providing transformational capacity for TES applications from astronomy to quantum technology.


Introduction
Superconducting transition-edge sensors (TES) are currently used in a wide range of disciplines, with applications ranging from astronomy, where they are used to count photons at the output of a spectrometer [1][2][3], through to quantum photonics, where output photons are counted in applications that have a known source spectrum [4][5][6][7][8][9].TES are single-photon detectors with a set of exquisite properties.These detectors achieve energy resolution of 0.15 eV from the infrared to Leonardo Assis Morais: leoassisfisica@gmail.com ultraviolet spectral range [10,11]; possess an exceedingly low dark count rate-even at low photon flux [12]; and operate with high efficiencyup to 98% [13,14].These characteristics combined make TES a photon detector with unique photon-number-resolving (PNR) attributes compared to other platforms like superconducting nanowire single-photon detectors (SNSPDs).
SNSPDs working principle differs fundamentally from that of TES, allowing them to achieve high detection efficiencies with fast detection rates, but they cannot intrinsically discriminate larger photon numbers [15].Different routes are being explored towards number resolution such as pulse analysis [16] or parallelisation [15,17,18].Nevertheless, both approaches have significant limitations regarding PNR and efficiency and pose significant challenges to electronic readout (see Table A1 Appendix A for a comparison of PNR implementations).
A central challenge for TES is their slower recovery (1/e decay) time, typically on the order of µs.Recent developments have improved this to less than 250 ns [19], allowing photon detection rates in the MHz range.For TES there exists a fundamental trade-off between achievable rates and PNR, as faster detectors either require higher operation temperatures and reduced detector area [20] or additional gold heat sinks [19,21], with both methods degrading the energy resolution.However parallelisation similar to SNSPDs [22] has been demonstrated as a promising path forward with less trade-offs.
Typically, energy information-be it spectrum or photon-number depending on the implementation-is extracted from a TES by measuring the time varying voltage signal pro-duced by each detection event.Usually, the entire data stream from the TES detector output is recorded and post-processed in software [23,24].Here, we introduce a different approach: a reconfigurable hardware processor, implemented with a Field-Programmable Gate Array (FPGA), that measures and stores a handful of physically significant pulse properties for every pulse, before discarding the voltage signal.
We have implemented a flexible and powerful platform to determine photon-number.Its capabilities include: digitisation of the TES analogue output; identification of pulses related to photon detection events in the digitised data stream; time stamping of detection events; measurement of relevant pulse characteristics to perform photon-number-assignment tasks; and, optionally, recording of the entire pulse, for calibration or cases where optimal information about detection events is required.

TES Signal and Analysis
The top of Figure 1 shows 200 detection events for a weak coherent pulse at 820 nm, registered using our hardware processor.After a detection is performed, a sharp rise is observed in the TES output, followed by an exponential-like decay whose constant will depend on the thin film properties [25,26].The sharp rises are of the order of hundreds of nanoseconds, and the decay time varies from 2 to 8 µs, depending on the energy of the absorbed photon wavepacket.The higher the number of photons in a detection event, the longer it takes for the detector to cool down to its operating temperature T 0 .
The bottom of Figure 1 highlights the input parameters used to register a valid detection event and the key characteristics subsequently recorded (see Appendix B for additional characteristics).The user sets: a detection baseline (solid black line) as the reference level with respect to the height measurements; the height threshold (bluedashed line); and the slope threshold (see Figure A1 in Appendix B for further details).To be considered a valid detection event, a TES pulse must have both height and slope higher than their respective thresholds.Only when this condition is met are the event packetsstructures that gather all the information about the measurements-created and recorded.FPGA showing some of the measurements our system can deliver.The solid black line is the detection baseline, and the dash-dotted blue line is the height detection threshold.
To assist the configuration of the baseline and detection thresholds, our processor contains a built-in multi-channel analyser (MCA).The MCA broadcasts histograms of a single quantity of the digitised output from the TES detectors in real time.We use the MCA to analyse the noise levels appropriately and set the baseline and detection thresholds.To determine the position of the detection baseline, we monitor the noise distribution, i.e. the detector output when no light is sent to the detectors.The detection baseline is chosen to be on the mode of the noise distribution.To determine the detection thresholds, we use the same light source which will be used in the

Counts
Figure 2: Histograms for photon detection events using a strongly attenuated pulsed diode laser at 820 nm for different aspects of the TES pulse.(Left) Height.Up to three peaks are clearly resolved.For a higher number of photons, the detector transitions to the normal phase, making it impossible to distinguish the photon number using height.(Middle) Maximum slope.Only 2 peaks are clearly resolved.We believe that the bandwidth filters used to remove high frequency noise compromise our ability to use pulse slope to distinguish photon number.(Right) Length.Due to signal noise, there is poor differentiation between events with one or more photons.
experiment.These thresholds are positioned between the end of the noise distribution and the first peak related to actual photon detections.Since anything lower than these thresholds will not be considered a valid photon detection event, these parameters will set the effective efficiency of the performed detections.Their appropriate values will depend on the experiment conditions, and a trade-off between detection efficiency and background counts must be considered when finding their optimal values.After setting the baseline and detection thresholds, our system is ready to perform measurements over different characteristics of the TES pulse.By storing only the pulse characteristics instead of the entire pulse, we achieve a 40 times reduction in transferred and stored data.The measured values are available in the FPGA less than 500 ns after the pulse falls below the height threshold.An Ethernet connection transfers pulse data to the computer, where photon number information is obtained within milliseconds after registering a pulse event.As TES are limited by their response time to the kHz counting regime, our photon number identification is fast enough to detect events in real time.

Detector Calibration
To precisely determine the best way to assign a photon number to a detection event, we measured and analysed four TES pulse characteristics: area, height, length, and maximum slope, as shown in the bottom of Figure 1.Height is measured as the distance between the baseline and the maximum point in the TES pulse.Length is the distance between the rising pulse crossing the detection threshold and the next point where the decaying pulse crosses this threshold again.Area is calculated by summing the heights between the detection threshold and the pulse every 4 ns, as shown in the shaded region.Slope is determined by differentiating the pulse, with the maximum slope being the highest value achieved by the slope signal (see Figure A1 in Appendix B).We show typical histograms for these quantities in Figures 2 and 3. We found that the area provides, by far, the best discrimination (as seen in the figures), as well as the highest precision.This reflects the fact that for TES, the pulse area is proportional to the energy absorbed during the detection event-even after the TES undergo the transition to the normal phase occurs [27]-and TES takes longer to cool down after absorbing a photon wavepacket with higher photon number (see top in Figure 1).
The important question is: given that a detection event is recorded, what is the chance of assigning a photon number m to a detection event with a true number of n photons?This probability is called discrimination.To evaluate it, we tested our system using a strongly attenuated pulsed diode laser with emission centred at 820 nm, described by the weak coherent state: where |α| 2 is proportional to the intensity of the beam.
A series of neutral density filters was employed to attenuate the optical power of our light source.16 Gaussian curves (see Appendix E for details).There are no peaks associated with zero photons due to the triggering requirement, as explained in the text.Note that the last two Gaussian curves are not visible due to their low amplitude.
Fine light intensity control was achieved by employing a pair of linear polarisers.To guarantee that our source was quasi-monochromatic, we used an 820 nm narrow-band filter with 2 nm bandwidth.The light was coupled to a single mode fibre which was directly connected to the TES detectors (see Appendix C for more details).Our optical setup was optically isolated to prevent noise contamination.To avoid multiabsorption events (as shown e.g. in Figure A3, Appendix D), we used a pulsed light source.The laser was driven with electrical pulses of 2.6 V amplitude and 50 ns width, using a repetition rate of 10 kHz.This guarantees a time separation of at least 100 µs between subsequent detection events.The TES pulse was amplified at room temperature using a pre-amplifier and differential amplifier, which also applied a 1 MHz bandwidth filter to remove high frequency noise.The filtered signal was digitised by a 14-bit analogue-to-digital converter at a rate of 250 MHz.The digitised signal was the input for our FPGA hardware processor.Our current implementation uses up to two channels simultaneously (see Appendix B for additional details on the FPGA).For further information about the hardware processor implementation, see Ref. [28].We measured a background count rate of 0.26 ± 0.04 counts/s.
A typical histogram of TES pulse areas is shown in Figure 3 (see Figure A5 Appendix E for histograms with different intensities).We de-termined the number of bins in the histogram by analysing how the reduced chi-square χ 2 R varies with the bin number (see Appendix E).Each peak on the TES pulse area histogram is associated with a different photon number.To describe the observed peaks, we used a model composed of a sum of 16 Gaussian curves.
To assign a detection event with pulse area a to a photon number m, we need to determine a set of counting thresholds t m such that if we assign a photon-number m to this event.
These counting thresholds t m were obtained after normalising each distribution obtained through the area histogram analysis.The counting thresholds t m are positioned in the intersection point between two adjacent normalised distributions, as shown in the left panel of Figure 4. Once these t m are established, assigning a photon-number to a detection event becomes a trivial task, which allows us to perform photon-number-resolving measurements in real time.Note that the distributions are wider for larger areas, thus the uncertainty increases when assigning a photon-number to larger areas.
To evaluate how precisely we can perform these photon-number assignments based on area measurements, we use the normalised distributions as the underlying distributions of photon numbers n registered by the TES in a given detection event.We are interested in determining p n m , the probability of correctly assigning the photon-number m to a detection event where n photons were detected.We calculate this probability using where g n (a) is the normalised distribution associated with the detection of photon-number n.
In the insets in Figure 4, we select two numbers (n=2 and n=15) to highlight the overlap between adjacent distributions.Note that for n=2 we obtain a precision of parts-per-billion, p 2 2 =0.99999999527, with the probability of incorrectly assigning a two event to one photon being p 2  1 =8.5×10 −10 ; and for assigning it to 3 photons being p 2  3 =3.88×10−9 .For higher photon numbers the discrimination probabilities, p i i , gradually reduce, and for 16 photons, we have p 16  16 =90% (all values summarised in Table A4 in Appendix F).
(2) (number in parenthesis corresponds to one standard deviation), and p 2  1 =8.5×10 −10 , p 2 3 =3.88×10−9 .The overlaps between these distributions are too small to be seen in the inset.(Bottom Right) Centred at n=15, we see p 15  15 =0.936 (4) (7) , and The decrease in discrimination reflects the decreasing counts available: we used a weak coherent pulse such that for n≥13 there are a relatively small number of counts.

Detector Tomography Routine
Our last step to characterise the performance of our detection system was to perform a quantum detector tomography routine [29].Quantum detector tomography aims to reconstruct the positive operator-valued measure (POVM) that completely describes the measurement apparatus [29,30].Each POVM element E m describes a possible detection outcome m.In our case, each E m describes a measurement result that we associate with a photon-number detection event of m photons.The probability of observing outcome m is given by where ρ is the probe state.To perform detector tomography, we need to use a tomographically complete set of input states.For photon-numberresolving measurements, a set of coherent states is tomographically complete [29,30].Our analysis used 4 strongly attenuated coherent states |α⟩ with α varying from 1.38 to 2.16 as input states.
The values for α were obtained from fitting coherent states to the counts obtained using the calibration method described in Section 3.
The experimental setup employed is the same described in Section 3, with a 820 nm pulsed diode laser being used as the light source.For all intensities, subsequent pulses are separated by 100 µs.For each intensity, we acquired approximately 1.7 × 10 7 counts.To measure the number of vacuum detections, we used the signal from the electric pulse generator as a reference.We chose a coincidence interval with 600 ns width (from 100 ns to 700 ns, see Figure A6 Appendix G) between the generator and TES pulses to incorporate the majority of photons from the laser while minimising accidental background photons.A generator pulse without a corresponding TES pulse was considered a vacuum count.
We note that, since TES measures the energy of the photon wavepacket, it is a phase insensitive detector [24,31].Therefore, only the diagonal of a POVM element will have non-zero terms, and only the intensity |α| 2 of the coherent states is sufficient to characterise the input states.For a coherent state |α⟩, the probability of observing n photons in a detection event is given by The tomographic routine consists in inverting Equation (4) using the recorded data for Pr(m) and the known set of input states {ρ i } to determine the set of POVM elements {E m }.To obtain the POVM, we performed a two-norm minimisation using the convex optimisation package for Python (CVXPY) [32,33].The estimated POVM is the solution to: where F is a matrix where the rows represent different α, and the columns represent their projection for different values of m used.E = E m for m = [0, N ] is the detector POVM, where N is the maximum number of photons the detector can distinguish.The POVM element E N is defined as where 1 is the identity matrix.E N accounts for all detection with N or more photons.The constraints over the POVM elements E m are to ensure the POVM is physically realisable.The estimated POVM elements E est,m are defined by where Pr(m|n) is the probability of reporting m photons given that there were n photons at the input.For an ideal POVM element, E n = |n⟩⟨n|, we would have Pr(m|n) = δ m,n in Equation (9).We show our results using a confusion matrix, see Figure 5, which takes advantage of the fact that all POVM elements are diagonal.Confusion matrices summarise the probabilities of reporting a photon-number m given that you had a detection event with photon-number n, i.e., Pr(reported m|input n).In Figure 5, each column corresponds to an estimated POVM element E est,m , and each row corresponds to an input photon number n.In the ideal case, the confusion matrix would be the identity matrix, with all terms in the diagonal equal to 1, and all nondiagonal terms 0. We note that the reconstruction quality falls considerably for n > 11.That is the case even when attempting to reconstruct an ideal POVM with the same input states used in the experimental reconstruction.The small range of intensities used for the input states and the low number of counts for high photon number limit the reconstruction quality in this range.
The readout fidelity f read provides the average probability of correctly reporting the number of photons in the range analysed, and it is calculated using: where N is the maximum number of photons in the interval analysed.When analysing the range from 0 to 11 or more photons f read = 0.99.For the full range shown in Figure 5, the readout fidelity drops to f read = 0.80, a consequence of the low number of counts in the region where n ≥ 11 1 .Additionally, we investigated the necessary intensity |α| 2 of a coherent state to perform quantum detector tomography with readout fidelities >0.99 over the entire photon number range studied here.We applied our routine in a POVM that describes a detector with 0.01 probability of incorrectly labelling a Fock state |n⟩ over the range between 1 to 16 photons: for example, the detector had 0.005 probability of labelling |10⟩ as |9⟩ or |11⟩, respectively.In this scenario, we would need coherent states with |α| 2 ranging from 0 up to at least 9 to observe the desired fidelity.In a realistic scenario, the probability of correctly assigning a photon number to a Fock state will depend on the Fock state itself and decrease for larger Fock states, further increasing the required |α| 2 .

Conclusion
Our FPGA-based hardware processor introduces features desirable across a wide variety of applications.Its parts-per-billion photon-number resolution allows Gottesman-Kitaev-Preskill (GKP) type optical quantum computing, where the detected photon-number output pattern is critical [34].In addition, number resolution shrinks the logical qubit overheads by minimising error rates, pushing towards thresholds for faulttolerant quantum computing [7][8][9], and is highly desirable in Gaussian BosonSampling [35][36][37].The reduction in collected and analysed data from the TES will enable new avenues, such as quantum-enhanced imaging where the ability to have even modest arrays of number-resolving detectors allows imaging maximal information extraction at the quantum limit [38].Possible applications range from astronomy for imaging large, remote bodies [38] to microbiology, where biomolecules can be probed at the single-photon level [39].A single FPGA-hardware processor could handle the readout of a small array, significantly reducing data transfer and storage, as 1 Due to technical problems with our detectors, we could not take data for coherent states with larger α under the same conditions as the data presented here.Nevertheless, we do not expect saturation effects to be relevant within the shown photon-number range.A single TES can distinguish photon numbers up to 37 photons, as seen e.g. in Ref. [22].
well as the number of data channels.
Finally, we dedicate this paper to our friend and co-author, the late Sae Woo Nam.We hope we did you proud.

Appendix A Comparison with superconducting nanowire single-photon detectors (SNSPDs)
Besides TES, SNSPDs are a promising technology for delivering photon-number resolving measurements.One of the main techniques to obtain photon number information with SNSPD is to use a series of detectors in a grid in such a way that the probability of two photons arriving in the same detector is negligible.In Table A1, we present a comparison of TES with these detectors, taking into account important characteristics, such as: detection efficiency η, dynamic range, wavelength λ range, operating temperature T O , and count rate.
We recently became aware of a very detailed performance comparison in Ref. [40], however, detectors are often tailored to excel in one metric and hence the best performance for different figures of merit rarely come from the same citation, as pointed out by the authors.

B FPGA Measurements
The FPGA performing the measurements on the digitised TES signal is a Xilinx 6 series implemented using the ISE toolchain.From the digitised signal, our FPGA derives two signals, as shown in Fig- ure A1: the TES signal (solid blue curve) and slope signal (solid red curve).First, the TES signal is derived directly from the ADC output after it passes a digital low-pass filter, which is used to remove high frequency noise present in the TES pulses.The slope signal is derived from the TES signal using a digital filter configured as a differentiator.
Rises in the data stream are identified using the slope signal: they start when the slope signal becomes positive and finish when it becomes negative.Once a rise has been identified, it needs to satisfy two conditions to be considered a valid detection.First, the maximum slope achieved during a rise must be larger than the slope threshold (dashed red line).Second, the height of the TES signal must be larger than the pulse threshold (dashed blue line).These detection thresholds were introduced to allow the hardware processor to distinguish between rises associated with an actual photon detection from smaller rises due to electronically associated signal fluctuations.Both detection thresholds are chosen by the user based on the information gathered and displayed by the multi-channel analyser, as described in the main text.
We also measure the pulse rise time from the TES pulse, which is the time taken from the moment the timestamp is given until it reaches the maximum height (Figure A1).Rise time must be considered when considering appropriate time windows, e.g. for counting coincidences.Additionally, we measure the number of rises in each TES pulse, which can be used to indicate the presence of background counts, for example.Additional rises indicate that two or more absorption events happened before the TES could cool down to its operating temperature, as shown in Figure A3.Based on this information, the user can perform a more careful analysis to check if it is possible to assign photon number to different events or to discard the pulse, if the information available is insufficient to assign a photon number with reasonable certainty.We note that our current implementation cannot perform height and maximum slope measurements simultaneously.

C TES Detectors
Our TES detectors were fabricated at the National Institute of Standards and Technology (NIST) Boulder and are thin films of tungsten, with a surface area of 25 µm 2 , thickness of 20 nm, fabricated on silicon.NIST reported efficiency of at least 95% for these detectors [13]; in-situ efficiency measurements were not performed.The tungsten thin films are embedded in an optical structure to optimise the Figure A2: TES bias and readout system.A stable-room temperature current source and a shunt resistance R sh (≈ 20 mΩ) are used to provide the TES voltage bias.When a photon absorption event occurs, the variation in current is amplified using an array of DC superconducting quantum interference devices (SQUID), which are inductively coupled to the TES circuit.These low-noise amplifiers are maintained at 4 K and amplify the TES signal 100 times.The SQUID output voltage is again amplified 100 times at room temperature and then sent to the differential amplifier.Finally, it passes through an analogue-to-digital converter and is fed to the FPGA, where all time and pulse shape measurements are performed.The voltage feedback, V F B , is used to bias the SQUIDS.
absorption at 820 nm wavelength.We use an adiabatic demagnetisation refrigerator (ADR) to maintain the detectors at an operating temperature of 80 mK, below the critical temperature of the tungsten film.The detectors are voltage biased, and after a detection event occurs their initial temperature is restored due to negative electro-thermal feedback.The detection circuitry is displayed in Figure A2.

D Multi-absorption events
TES have no dead time, and events with different photons being absorbed at different times can be recorded, as shown in Figure A3.A second detection event occurs at around 10 µs, before the detector had enough time to cool down after the first detection.Our FPGA can identify such events by identifying multiple rises in a single detection event and storing the height for different rises.
Since multi-absorption events add an extra layer of complexity when performing the photon-number assignment task, in the main paper we used a pulsed light source much slower (10 kHz) than the longest cool down times (≈ 10 µs) to avoid multi-absorption events.Our FPGA can identify the multiple rises in a detection event and measure the height of each rise.
The user can choose to develop a method to analyse these cases or to remove the events where multiabsorption events occur.

E Area Histogram
The model M we used to fit the area histogram-the dashed red line in Figure 3-is composed of a sum of Gaussian curves G n (a), defined by where with A n being the amplitude, µ n the mean, and σ n is the standard deviation of the Gaussian used to describe peak n.To evaluate the goodness of the fit, we used the reduced chi-squared, defined as where C i is the number of counts in the i th bin, ϵ i is the error in the i th bin, and k is the number of degrees of freedom.We define k as where n bins is the number of bins in the histogram and n par is the number of parameters in the fit.To determine n bin we observed how the χ 2 R varied as function of the number of bins used, see left panel of Figure A4.To calculate χ 2 R , we assumed our data follows Poissonian statistics, therefore ϵ i = √ C i .We divide the analysed area interval of [0, 10 × 10 6 ] into 45000 bins, as we obtained χ 2 R ≈ 1 in this case.The number of detection events used was ≈ 1.7 × 10 7 .The residual Res, shown in the right panel of Figure A4, was calculated using To determine the counting threshold position t m , we used the normalised distributions g n (a) given by which has an area equal to 1.
Once the counting thresholds t m were established, we calculated the χ 2 R for each g n (a) to quantify how well a single Gaussian could describe the data points inside the area interval [t m , t m+1 ], see Table A2.Note that for high photon numbers (n > 14), the χ 2 R is smaller than 1.This is an expected result, since this is a low count region due to the relatively low intensity of the light source used.In Figure A5 we present histograms for different intensities, showing that the source's intensity can limit the highest number of photons we can distinguish using this calibration process.
To show that this calibration procedure is robust within the intensity range, we estimated α for 4 different coherent states using different calibrations.Besides the calibration reported in the main paper, we also calibrated the detector using the coherent state with α = 1.32 (blue dots in Figure A5).In this case, our model was composed of 10 Gaussian curves due to the lower intensity of the coherent state used for calibration.The counting thresholds t n were obtained using the same method described before.In Table A3, we compare the estimated values for α for 4 different coherent states using both calibrations.For the left and right tables, the values of α for the coherent state used for calibration were 1.32 and 2.16, respectively.Note that the α values for the entire set agree within one standard deviation, showing that our calibration provides consistent estimation for the coherent states used in this low intensity range.We highlight that the standard deviation presented here is the value provided by the fitting algorithm and does not represent the experimental uncertainty with the coherent state used for calibration.

F Photon-Number Assignment Probabilities
The discrimination problem: once a photon detection event has been recorded, we want to determine the number of photons in that event.Table A4 lists the probabilities of correctly, m=n, or incorrectly, m̸ =n, assigning a photon number m to a detection event with n photons, p n m .Since only the cases where m=n−1 and m=n+1 have non-negligible values, we present only these three columns.Note that for 17 photons or more, we can no longer discriminate the photon number, therefore p n 17+ should be read as the probability of assigning 17 photons or more to a detection event where n photons were detected.The uncertainties in the probabilities of making a correct assignment p n m=n were calculated by considering the uncertainties in the position of the detection thresholds t m .Code available on request.

G Coincidence window for tomography
As stated in the main paper, we set a coincidence window of 600 ns for the tomography analysis.The window for the detection events to be considered valid spanned from 100 ns to 700 ns after a detection pulse was registered.In Figure A6, we show the relative time the detection events were reported after an electric pulse was registered in the FPGA.We note that 525 detection events were outside this time range, out of approximately 1.8 × 10 7 counts.From these, 518 were 1-photon detection events, most likely stray thermal background photons.

Figure 1 :
Figure1: (Top) 200 TES pulses captured by our hardware processor using 820 nm pulsed light from a strongly attenuated diode laser.The voltage reading is proportional to the current change due to a photon wavepacket absorption.By keeping the sensor in the transition region between its superconducting and normal phases, we can measure energy with enough precision to distinguish the pulse associated with different photon numbers.(Bottom) Single TES pulse digitised by our FPGA showing some of the measurements our system can deliver.The solid black line is the detection baseline, and the dash-dotted blue line is the height detection threshold.

Figure 3 :
Figure 3: TES Pulse area histogram (black circles) for a strongly attenuated pulsed diode laser at 820 nm.Dashed red line is the model comprised of a sum of16 Gaussian curves (see Appendix E for details).There are no peaks associated with zero photons due to the triggering requirement, as explained in the text.Note that the last two Gaussian curves are not visible due to their low amplitude.

Figure 4 :
Figure 4: (Left) Solid curves are normalised area distributions, g n (a).The counting thresholds t m are positioned in the intersection between adjacent distributions.The probability of assigning m photons to a n photon detection event p n m is the percentage of area of g n (a) within a coloured interval [t m , t m+1 ]. (Right) Insets highlighting the overlaps between three adjacent distributions.(Top Right) Centred at n=2, we see p 2 2 =0.99999999527

Figure 5 :
Figure 5: Confusion matrix showing the results for the detector tomography routine.POVM elements are the columns of the confusion matrix, see Equation (9).Colour bar at figure right shows the probability of reporting a photon-number m given an input photon-number n, Pr(m|n).The quality of reconstruction degrades for n > 11 due to the low intensity of the coherent states used in the tomographic routine.We see high quality reconstruction for n ≤ 11, with Pr(m = n|n) ≥ 0.99 in this range.

Figure A1 :
FigureA1: TES pulse detailed with measurements that our FPGA can perform.The solid blue line is the TES pulse and the solid red line is the TES pulse slope (amplified 1000× for clarity).

Figure A3 :
FigureA3: TES analogue output showing a multi-absorption event obtained using an oscilloscope.Our FPGA can identify the multiple rises in a detection event and measure the height of each rise.The user can choose to develop a method to analyse these cases or to remove the events where multiabsorption events occur.

χ 2 R
Figure A4: (Left) χ 2 R for our model M as a function of the number of bins used in the area histogram.(Right) Residual calculated from Equation (15).

Figure A5 :
FigureA5: Area histogram for different intensities registered with our detection system.Each colour corresponds to a different intensity of the light source.Note that the peaks-which correspond to a given photon number-occur in the same position independent of the intensity used.

6 )Figure A6 :
Figure A6: Histogram of relative times of detection events in a 600 ns window with respect to the registration of an initial electric pulse from the laser.

Table A2 :
Reduced chi-squared, χ 2 R , between model and data for all peaks.They were calculated by considering the Gaussian distribution g n (a) used to describe a peak in the interval [t n , t n+1 ].

Table A3 :
Comparison between estimated values for α using two different coherent states for calibration.For the left table, the coherent state used for calibration had α = 1.32.For the right table, α = 2.16.We note that all values agree within one standard deviation σ.

Table A4 :
Probabilities of assigning a photon number m to a detection event with n photons.The middle column is the probability of a correct assignment, m=n; the left and right columns are, respectively, probabilities of the incorrect assignments m=n−1 and m=n+1.The upper and lower numbers in parenthesis indicate the error (one standard deviation) in the reported values for correct probability assignments.