Ultrabroadband direct detection of nonclassical photon statistics at telecom wavelength

Broadband light sources play essential roles in diverse fields, such as high-capacity optical communications, optical coherence tomography, optical spectroscopy, and spectrograph calibration. Although a nonclassical state from spontaneous parametric down-conversion may serve as a quantum counterpart, its detection and characterization have been a challenging task. Here we demonstrate the direct detection of photon numbers of an ultrabroadband (110 nm FWHM) squeezed state in the telecom band centred at 1535 nm wavelength, using a superconducting transition-edge sensor. The observed photon-number distributions violate Klyshko's criterion for the nonclassicality. From the observed photon-number distribution, we evaluate the second- and third-order correlation functions, and characterize a multimode structure, which implies that several tens of orthonormal modes of squeezing exist in the single optical pulse. Our results and techniques open up a new possibility to generate and characterize frequency-multiplexed nonclassical light sources for quantum info-communications technology.

Multimode nonclassical light sources and their characterization have been studied extensively these years. In the frequency domain, quantum continuous variable (CV) correlations in the frequency comb were generated using a single optical parametric oscillator (OPO) and characterized for many sideband modes by homodyne detectors with designated local oscillator (LO) modes 21,22 . In the time domain, CV entanglement over massive temporal modes, was created by using two OPOs and a delayed interferometer 23 . In the spatial domain, CV entanglement over multiple transverse modes was generated from an OPO and characterized by homodyning 24,25 .
Besides these studies with CVs, discrete variable (DV) nature, i.e., photon-number statistics, of multimode quantum states have been investigated and characterized in this decade, by using PNRDs. The direct observation of nonclassical photon-number oscillation from a type-I SPDC source was demonstrated using a visible light photon counter (VLPC) 26,27 . The higher order photon-number correlations between the signal and idler beams of a twin beam from a type-II SPDC source were observed, and the influence of multimode structure to photonnumber statistics was analysed, using the APD-based time-multiplexing PNRD 28,29 . Extending this idea, reconstruction of the mode weight structure from the photon-number correlations was also demonstrated for the states containing three modes 30 . All of these works with CVs and DVs have been performed at the visible or nearinfrared wavelengths.
In this paper, we demonstrate the direct detection of nonclassical photon statistics of ultrabroadband squeezed states at the telecom wavelength. The squeezed states are generated from a SPDC source and their broadband multimode photon-number statistics is measured by a highly efficient PNRD based on a titanium TES 11,12 . The centre wavelength of the generated state is 1535 nm and the measured bandwidth is as broad as 110 nm (13.4 THz), covering the S-, C-, and L-bands which are major bands for optical fibre communication. We characterize the nonclassical property by the following two observations: first, we directly measure the photon-number distribution of a whole multimode state and confirm its nonclassicality by Klyshko's criterion 26,31 . Second, we apply the method proposed by Christ et al. 18 and reconstruct the mode weight distribution, i.e. the degrees of squeezing in each decomposed mode from the higher order photon-number correlations, which implies that several tens of independent squeezers are contained in our source. The results and techniques pave a new simple way of characterizing ultrabroadband nonclassical states and we anticipate it will accelerate further research and development on quantum information technologies for the telecom fibre infrastructures.

Results
Broadband squeezed light source. Our light source is a single-beam squeezer in a collinear configuration where the pump, signal, and idler beams are in the same spatiotemporal and polarization mode. The generation process is mathematically described by the following unitary transformation 18,32 S~exp Here A denotes the overall efficiency of the squeezing, f (v s , v i ) is the joint-spectral distribution (JSD) in terms of the signal and idler angular frequencies, v s and v i , the creation operator of the signal (idler) field in the continuous spectrum of v s (v i ). These continuous modes are, however, not necessarily a convenient representation to deal with the multimode characteristics, because they cannot be decoupled from each other for most of practical JSDs. For a JSD which distributes in an effectively finite range, one can find a more convenient discrete set of decoupled modes [32][33][34][35][36] .
In fact, when the JSD is engineered to be symmetric, which is the case here, it can then be decomposed into the following form in terms of a complete orthonormal set {w k (v)} 18 . Note that each w k (v) represents the shape of the frequency distribution for signal and idler. Here r k is a modal amplitude corresponding to a squeezing parameter for mode k. Introducing the field operators for the broadband basis modes 36 the above transformation can be expressed aŝ The generated state is a multimode squeezed vacuum state as Thus this state only includes even number of photons. In practice, however, inevitable losses in the generation, propagation, and detection processes cause the vacuum invasion, making the pure state jrae a mixed one, whose statistics has odd number components as well. If losses are below a certain level, one can still observe even-rich photon-number statistics. Direct observation of this even-odd number oscillation and multimode analysis on the distribution of r k are our tasks here.
Experimental setup. Figure 1 shows the experimental setup. The fundamental light source is a single-longitudinal-mode, passively Q-switched, diode-pumped solid-state laser (Cobolt, Tango), operating at 1535 nm with a pulse width of 4 ns and a repetition rate of 3.1 kHz. This is used as the fundamental light, and is injected into a 10 mm long, periodically poled KTiOPO 4 (PPKTP) crystal for second harmonic generation (SHG) at 767.5 nm. This second  harmonic light is then used as a pump for squeezing. An unconverted fundamental light after SHG is used as a guiding beam for fibrecoupling to the TES, but effectively eliminated by dichroic mirrors when photon counting experiment is implemented. The pump is guided into the second type-0 PPKTP crystal employed as the squeezer. The PPKTP crystals for SHG and the squeezer are mounted on copper stages and temperature-stabilized. The multimode squeezed vacuum state is fibre-coupled and guided to the TES or to a spectrometer. Our TES is based on a titanium superconductor whose size is 10 3 10 mm 2 as shown in the reference 11,12 , and is set inside a dilution refrigerator (TS-3H100-PT, Taiyo Nippon Sanso). Output waveform from the TES is amplified by a SQUID device and room-temperature electronics (Magnicon). Each waveform was sampled using a digital oscilloscope (DPO7104, Tektronix), and sent to a host computer for noise filtering and data analysis. Figure 2a shows raw waveforms from the TES. Total 1 3 10 7 waveforms are stored and overlaid here. The gating time for each detection event was 10 microseconds. Within this gating period, an average photon number of TES's dark counts was measured to be 5 3 10 24 per pulse, which corresponds to 50 counts per second. The wavelength dependence of our TES's DE is shown in Fig. 2b (also see Methods). It is almost flat from 1512 nm to 1580 nm, and then gradually drops from 1580 nm to 1620 nm. The absolute DE is 5 90.9 6 2.4% at 1535 nm. Unfortunately, we were not able to measure DE at wavelengths shorter than 1510 nm nor longer than 1620 nm, because of a limited wavelength range of the probe light.
The red line in Fig. 2c is a spectrum of the fundamental light, whose linewidth is ,0.05 nm. The green dotted line is a spectrum of the SPDC measured with a spectrometer (Acton) for a 30 sec accumulation time. Measured full-width at half maximum (FWHM) is roughly 130 nm. The blue line is a theoretically derived spectrum using the present experimental parameters (see Methods), which is symmetric about the fundamental light spectrum. The asymmetry seen in the measured spectrum (green dotted line) is an artifact caused by the DE fall-off of an InGaAs photo detector in the spectrometer. Actually the spectrometer DE gradually drops when the wavelength exceeds 1500 nm. Then the DE starts steeply falling at around 1580 nm and reaches almost 0% at around 1650 nm. The blue theoretical spectrum captures the actual SPDC spectrum whose FWHM extends over almost 150 nm in the telecom window, covering the S-, C-, and L-bands.
Nonclassical photon statistics of multimode squeezed states. Figure 3 shows the experimental results for the multimode squeezed states. Through 3a, b, and c, the histograms in red and blue correspond to a pump pulse energy of ,1 pJ and ,10 pJ, respectively. Figure 3a shows the pulse height distribution of our Ti-TES output waveforms. A voltage range of 220 to 60 mV is divided into 800 bins. The total number of events was 1 3 10 7 for each data set. The vertical axes in the main graphs are numbers of counts shown in log scale, while those in the insets are in linear scale. Every peak in each distribution was clearly separated, and was fitted with a Gaussian function in order to determine thresholds for calculating photonnumber probabilities. The FWHM of the each Gaussian peak corresponds to an energy resolution and was calculated to be about 0.2 eV, which is smaller enough compared with a photon energy (,0.8 eV) of the monochromatic fundamental light at 1535 nm. Figure 3b insets show the photon-number distributions obtained from the pulse height distributions without any correction of losses. In each condition, the results clearly show super-Poissonian distribution, whose feature could be confirmed by the Fano factors, which are greater than the unity, yielding 1.462 6 0.002 (red) and 1.494 6 0.001 (blue), as well as by g (2) values, 42.746 6 0.195 (red) and 4.978 6 0.008 (blue). Then, the overall DE including the fibre-coupling efficiency can be evaluated by comparing the single-and two-photon probabilities, according to the formula, g DE~2 P 2 =P 1 1z2P 2 =P 1 27 . The overall DE was thus estimated to be g DE 5 48.4 6 0.2%. Note the dark counts slightly changed P 1 , and we therefore subtracted its contribution from the raw statistics for calculating DE.
The main graphs in Fig. 3b are reconstructed photon-number statistics using the iterative maximum-likelihood (ML) estimation algorithm 37 . Through the ML estimation process, we used the positive operator-valued measure (POVM) of our TES described as follows 38 , The POVM elements are normalized in order to fulfill S m P m 5 I. TES is thus modelled as a linear photon counter which detects m photons from n input ones with the finite detection efficiency g DE , and with the dark counts per pulse d. M is a cutoff photon-number chosen to be large enough in order to prevent estimation errors by truncating higher photonnumbers. We chose M 5 10 here. The reconstructed statistics were fully converged after 1 3 10 6 iterations of the ML estimation algorithm. Note, in our case, the linear inversion of the binomial distribution B (m2i)k resulted in unphysical, negative photon probabilities because of alternating signs in B {1 m{i ð Þ k components 39 . The reconstructed two-photon probability P 2 in red (upper panel) is about 1%. On the other hand, the reconstructed two-and four-photon probabilities for the stronger pumping condition (lower panel in blue) resulted P 2 , 10% and P 4 , 0.6%, respectively. Thus the reconstructed probabilities exhibit even-rich photon-number statistics due to the pairwise photon generation of squeezing. We also evaluated the nonclassicality using Klyshko's criterion 26,31 .
Klyshko defined the figure of merit K n~n z1 ð Þ P n{1 P nz1 nP 2 n n~1, 2, . . . ð Þ and showed that for any {P n } originating from classical states, such as coherent state or thermal state, K n $ 1 should be satisfied for all n 31 . In other words, if this inequality is violated for any of K n , the distribution {P n } is necessarily from a nonclassical state. In Fig. 3c, K n were calculated from the raw photon distribution (Fig. 3b insets), and the classical limit (green horizontal line) is clearly violated (K n , 1) for even-photon numbers (n 5 2, 4, and 6), while all the odd photon-numbers do not violate the limit. Ideally, squeezed source consists of even-photon numbers and therefore all the K n for odd numbers must go infinite. There were in practice some contributions of odd-photon numbers caused by the finite DE, resulting in the finite K n for odd photon-numbers.
Mode analysis of multimode squeezed states. Using the above PNRD data, we evaluate experimental g (2) and g (3) , and analyse multimode structures of our squeezer by applying the method in the ref. 18. Figure 4a shows experimental g (2) and g (3) calculated from the photon-number statistics using mth factorial moment n (m) given by g m . Because g (m) is a losstolerant measure, we use AEnae 5 AEn exp ae/g DE in order to compare the experimental g (2) and g (3) to the theory. As expected from the fact that our squeezed states are highly multimode, our experimental g (2) and g (3) show clear deviations from the single-mode theory. This can be contrasted to an experimental result on single-mode squeezed states from a PPKTP, which showed that the experimental g (2) agreed well with the single-mode theory 41 . It is generally a nontrivial task to carry out the Schmidt decomposition analytically. For certain classes of states, however, it can be made. A typical class is a two-dimensional (2D) JSD f(v s , v i ) in a Gaussian form, which is actually the case here. In this case, the eigen functions and eigen-value distribution of Eq. (2) are the Hermite polynomials and the thermal distribution, respectively, as shown in Appendix in the ref. 17. Then g (2) and g (3) can be expressed as follows 18 (also see Methods), with L 4 m ð Þ~1 The red crosses in Fig. 4b show the g (3) values as a function of g (2) . By fitting the slope of g (3) -vs-g (2) to S m ð Þ using the linear least squares fitting technique, we derived S m exp ~3 :241+0:151. This corresponds to m exp 5 0.961 6 0.024. Once m exp is determined, each B exp is obtained by eqs. (6) and (7) using corresponding g (2) and g (3) , yielding 0.504 (for , 10 pJ) and 0.157 (for ,1 pJ). Note that B exp can be  independently determined by the measured average photon-number via the relation B exp~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n exp g DE q which implies 0.506 (for , 10 pJ) and 0.148 (for , 1 pJ) agreeing well with the above results. Thus we can finally obtain the mode distribution of the squeezing parameter r k (5B exp l k ). The amount of squeezing for each r k can be calculated via 210 log 10 (exp(22r k )) [dB]. Reconstructed squeezinglevel distributions (k 5 0, 1, …, 40) are shown in Fig. 4c for a pump energy of ,10 pJ. The effective mode-number K of the squeezed vacua contained in the reconstructed distribution can be calculated . It is worth to compare m exp with its theoretical predictions. From the JSD characterized by the phase-matching condition of our nonlinear crystal and an envelope of the pump pulse, we can calculate the theoretical m as m JSD , 0.99844 (see Methods for detailed calcuation) which is slightly deviated from m exp . We also estimated the origin of the uncertainty of m exp by the Monte Carlo simulations with 1 3 10 7 trials (same as that of the experiment. See Methods for details of the simulation). Figure 4d compares the experimental and simulated g (3) -vs-g (2) plots, showing that those agree very well. Then from the simulated result, we obtained S m sim ð Þ by linear least squares fitting which yields S m sim ð Þ~3:0164+0:1288 and m sim 5 0.9973 6 0.0212. We conclude that the uncertainty of our experimental data is mainly due to the statistical errors. For further investigation of the discrepancy between m exp and m JSD , the setup should be further stabilized to allow a longer data acquisition for reducing the statistical errors.

Discussion
In this work the nonclassical photon-number distributions as the even-odd photon-number oscillations were observed in an extremely broadband of 110 nm FWHM. This bandwidth is, to our knowledge, the broadest one whose nonclassicality was directly detected. The bands are technologically important ones, namely the telecom S-, C-, and L-bands. The mode reconstruction technique theoretically developed in the ref. 18 was applied, which inferred that the several tens of orthonormal modes were contained in every single optical pulse.
There are interesting open issues along with this work. First, the method used here is just for evaluating the mode weight distributions under certain assumptions on the SPDC spectral correlation. More sophisticated methods on mode analysis using direct detection with PNRDs are desired to be developed, as can be implemented by using homodyne detection 21,22 . Second, low-loss mode-separating devices should also be developed, which may include passive low-loss bandpass filters or recent proposals of quantum frequency conversion 42,43 , and combined with highly efficient PNRDs to detect and analyse mode characteristics of nonclassical states. Third, when a shaped broadband LO pulse 44 is introduced in front of a PNRD, possibly  (2) and g (3) . They show clear deviations from the single-mode cases (solid lines).
(b) Experimental g (3) are plotted as a function of experimental g (2) (red crosses). The blue solid line is a result of fitting the data. From a slope of the fitted curve, we obtained m exp , that determines a shape of the squeezer distribution. The green solid line is the single-mode case, g (3) 5 9g (2) 2 12, for reference. (c) is a reconstructed squeezing-level distribution using m exp and B exp for a pump energy of ,10 pJ. (d) Result of Monte Carlo simulations. The blue crosses are g (2) and g (3) obtained by the simulations. The red crosses are the experimental data replotted for comparison. The magenta line is multimode limit shown as a reference.
www.nature.com/scientificreports SCIENTIFIC REPORTS | 4 : 4535 | DOI: 10.1038/srep04535 combined with mode-separating devices, one can implement a displacement operation before a PNRD which could realize a phasesensitive PNRD. A displaced squeezed state is known to show unique phase-dependent oscillation in photon-number distributions 45,46 , which has never been demonstrated yet. It also provides a means to implement quantum receivers to overcome the conventional limits of optical communications [4][5][6]47 .
Finally broadband light sources play essential roles in diverse fields in the classical domain, such as high-capacity optical communications 48 , optical coherence tomography 49 , optical spectroscopy 50 , and spectrograph calibration 51 . Supercontinuum sources in the vicinity of the 1550 nm telecommunication wavelength 52 are of particular importance for wavelength division multiplexing in fibre-optic communications 53 . Our ultrabroadband source and detection scheme shown here could also provide a way to utilize the multimode nonclassical states for realizing such capabilities in the quantum domain.

Methods
Calibration of TES's detection efficiency depending on the wavelength. We used a continuous-wave (CW), wavelength-tunable ECDL (TSL-510, Santec) in order to measure the DE's wavelength dependence. The CW light beam from the ECDL was connected to the fibre-coupler. One of the outputs from the coupler was guided to our TES after passing through the variable attenuator. This power to the TES was heavily attenuated (,90 dB) such that the average power was set to be below 1 fW, and was monitored via the other output port of the coupler using a power meter, whose wavelength sensitivity was precisely calibrated. The coupler's splitting ratio and the degree of the attenuation were also precisely measured using the same power meter. Because photon arrival was random for the CW source, total counts from the TES was collected by the threshold detector (SR400, Stanford Research Systems). Photon counts from the TES were a few kHz, and the raw waveforms were not overlapped each other. We scanned the wavelength at each 3 nm step from 1512 to 1620 nm, and obtained the wavelength dependence of our TES's DE as shown in Fig. 2b.
The second-and third-order correlation functions for the multimode squeezer. In this multimode ''single-beam squeezer'' case, the second-and third-order correlation functions are described in the ref. 18 as follows, When a JSD takes a 2D Gaussian shape, l k forms a thermal distribution defined by a single parameter m, as l k~ffi ffiffiffiffiffiffiffiffiffiffiffi 1{m 2 p m k . When the r k is small enough (r k =1), that is valid in our experimental conditions, the following approximations are valid as Then, eqs. (9) and (10)  Specification of PPKTP crystal and SPDC spectrum. The JSD can be a product of two quantities as where a(v s , v i ) describes pump-frequency envelop, and w(v s , v i ) describes phase-matching function for the SPDC in a nonlinear crystal 54 . Explicitly a(v s , v i ) 5 sinc(LDk/2) with the crystal length L and phase mismatch Dk 5 k p 2 k s 2 k i 2 2p/L. L denotes a grating period of the PPKTP crystals and 24.2 mm in our case.
Assuming that the transversal phase-matching is perfect, the longitudinal one can be described as Dk 5 k p,z 2 k s,z 2 k i,z 2 2p/L, because we used collinear configuration where pump, signal, and idler are all in the same polarization. Note that the beam propagating direction is taken to be along z-axis and k m (v m ) 5 v m n z (v m )/c is the wave vectors. The frequencies are in a relation of v i 5 v p 2 v s 5 2v 0 2 v s . The corresponding Sellmeier equation for n z can be found in the ref. 55. The blue line in Fig. 2c was theoretically derived using the above formulas and parameters. The theoretical K JSD can be calculated from the JSD using the singular value decomposition (SVD) technique 29 . The JSD was obtained from the above parameters and the pump linewidth of , 0.05 nm, and then discretized into a 16000 3 16000 square matrix with a frequency resolution of 2 GHz and with a frequency range of 179.3 THz (, 1672 nm) to 211.3 THz (, 1419 nm). We carried out SVD for the JSD matrix using MATLAB, and obtained K JSD , 640. The corresponding m is m JSD , 0.99844.
Theoretical photon probability of multimode squeezed states for Monte Carlo simulations. In order to analyse the discrepancy between the experimentally obtained m exp and the theoretically obtained m JSD , we carried out Monte Carlo simulations to illustrate how finite data size affects S(m)-fitting errors. Before the simulations, we first calculate theoretical photon probabilities for the multimode squeezed vacua. The density matrix of the multimode vacuum state should be expressed as a tensor product of the decomposed squeezed vacuum states, whose squeezing parameters are characterized by r k (5B exp l k ). From this density matrix, one can derive the n-photon generation probabilities of the source. Though n is distributed to infinity, we truncated it at six photons which is a reasonable approximation to simulate our experimental situation in Fig 3b. The n-photon generation probabilities up to n 5 6 are given by where P k ð Þ 0~2 2zr 2 k À Á , P k ð Þ 2~r 2 k 2zr 2 k À Á , and N~X n P tot n n~0,2,4,6 ð Þ . Here we assume that the r k is small enough and each mode only emits vacuum or two-photons. This assumption is valid when the mode number K9 is large enough. For example, in the case of four-photon emission, the number of cases for two-photon generations from the two different modes ( K9 C 2 , O(K9 2 )) are much greater than those for fourphoton emission from only one mode ( K9 C 1 , O(K9)). Then the Monte Carlo simulation was performed along with eq. (16).