Quantum state features of the FEL radiation from the occupation number statistics

The statistical features of the radiation emitted by Free-Electron Lasers (FELs), either by Self-Amplified Spontaneous Emission (SASE-FELs) or by seeded emission (seeded-FELs), are attracting increasing attention because of the use of such light in probing high energy states of matter and their dynamics. While the experimental studies conducted so far have mainly concentrated on correlation functions, here we shift the focus towards reconstructing the distribution of the occupation numbers of the radiation energy states. In order to avoid the various drawbacks related to photon counting techniques when large numbers of photons are involved, we propose a Maximum Likelihood reconstruction of the diagonal elements of the FEL radiation states in the energy eigenbasis based on the statistics of no-click events. The ultimate purpose of such a novel approach to FEL radiation statistics is the experimental confirmation that SASE-FEL radiation exhibits thermal occupation number statistics, while seeded-FEL light Poissonian statistics typical of coherent states and thus of laser light. In this framework, it is interesting to note that the outcome of this work can be extended to any process of harmonic generation from a coherent light pulse, unlocking the gate to the study of the degree to which the original distinctive quantum features deduced from the statistical photon number fluctuations are preserved in non-linear optical processes.

The statistical features of the radiation emitted by Free-Electron Lasers (FELs), either by Self-Amplified Spontaneous Emission (SASE-FELs) or by seeded emission (seeded-FELs), are attracting increasing attention because of the use of such light in probing high energy states of matter and their dynamics. While the experimental studies conducted so far have mainly concentrated on correlation functions, here we shift the focus towards reconstructing the distribution of the occupation numbers of the radiation energy states. In order to avoid the various drawbacks related to photon counting techniques when large numbers of photons are involved, we propose a Maximum Likelihood reconstruction of the diagonal elements of the FEL radiation states in the energy eigenbasis based on the statistics of no-click events. The ultimate purpose of such a novel approach to FEL radiation statistics is the experimental confirmation that SASE-FEL radiation exhibits thermal occupation number statistics, while seeded-FEL light Poissonian statistics typical of coherent states and thus of laser light. In this framework, it is interesting to note that the outcome of this work can be extended to any process of harmonic generation from a coherent light pulse, unlocking the gate to the study of the degree to which the original distinctive quantum features deduced from the statistical photon number fluctuations are preserved in non-linear optical processes.

I. INTRODUCTION
Recent years have seen an increasing interest in the generation of ultra-bright and ultrashort electromagnetic pulses by means of radiation sources, widely known as Free-Electron Lasers (FELs), through self-amplification of spontaneous emission (SASE) from a relativistic electron bunch wiggling between a series of magnetic poles known as undulators [1]. Despite the growing experimental and theoretical activity in the field, the issue concerning the nature of statistical fluctuations in the emitted radiation remains still largely unmet. In the extreme, such a complex issue boils down to the question: "To what degree is FEL radiation laser radiation?" The experimental intensity, angular and spectral distributions of the radiation emitted by SASE-FELs are well predicted by classical electrodynamics. Not surprisingly, also a quantum field theory approach, such as that firstly proposed in [2], can recover the well-known exponential growth of the electromagnetic field intensity along with the exponential regime yielding a satisfactory description of the behavior of a single passage SASE-FEL. Instead, more surprising theoretical predictions appear in [3] and [4] that the radiation resulting from a SASE-FEL could show properties typical of squeezed light. Indeed, any quantum description of an electromagnetic phenomenon also contains a classical description as an approximation. On the contrary, a phenomenon which can be described by a classical theory, and yet exhibits purely quantum features, points to a weakness in the classical expla- * corresponding authors: Fabio Benatti (benatti@ts.infn.it) theory and Fulvio Parmigiani (fulvio.parmigiani@elettra.eu) experiment nation and requires a quantum theory of SASE-FEL radiation that is able to reproduce not only its classical features, but also the quantum characteristics. Although a possible explanation might be found in the fact that the spontaneous emission of radiation originates from the fluctuations of the electromagnetic vacuum, a solid theoretical explanation of such an effect is, to the best of our knowledge, still missing.
In this respect, the study provided in [5] made perfect sense; there, in a photon counting experiment performed on the spontaneous harmonic radiation generated by an infrared, SASE-FEL using a well-defined ensemble of electron pulses, the authors announced the detection of genuinely quantum sub-Poissonian intensity fluctuations. However, a more recent work [6] challenges such findings together with their consequences in that the reported quantum signatures of the SASE-FEL radiation can be explained, within standard, classical SASE-FEL theory, by combining detector dead-time effects with photon clustering arising from the FEL gain. Thus, to bridge the still existing gap between the theoretical predictions of quantum effects and the experimental evidence, one must develop further a fully consistent quantum explanation of SASE-FEL radiation production that predicts possible distinctive quantum features in the statistical fluctuations of SASE-FEL radiation.
A completely different scenario arises if the FEL process is triggered not by the background radiation of an oscillating relativistic electron bunch, but by a pulse of laser light with intensity several orders of magnitude larger than the background incoherent shot noise. At the moment, we can only conjecture that the coherent features of the "seeding" laser pulse are transferred to the electrons in the bunch that in turn convey them, at least in part, to the FEL radiation, first in the mod-FIG. 1: Left: cartoon of lasing in SASE and HGHG mode. In SASE, energy and density modulation of electrons at the scale of the resonant FEL wavelength happen simultaneously along the undulator. The spectrum of the amplified spontaneous radiation is an ensemble of spikes whose intensity varies shot to shot, each spike being individually coherent. In HGHG, the energy modulation generated by the electron-laser interaction in the modulator translates into density modulation via the magnetic chicane. The FEL spectrum is a narrow band single spike. Right: intensity profile of consecutive shots of the spectrum of the FERMI FEL, run in SASE (top) and HGHG mode (bottom), for the same setting of the accelerator. ulation phase and then in the radiation phase, thus producing amplified high harmonic radiation extending up to the X-ray domain and showing a coherent, that is a lasing state, of the electromagnetic field.
However, to be coherent and thus to describe laser light, the quantum states of seeded-FEL radiation must satisfy with specific requisites: they must be characterized by all 2npoint normalized Glauber correlation functions g n equal to 1 in modulus [1] and the occupation numbers of the available energy states must exhibit a Poisson distribution. Experimental investigations focussing upon the Glauber correlation functions g 1 and g 2 both for SASE-FEL [9,10] and seeded-FEL light [11] revealed typical chaotic light behaviour in the first case, |g 2 | = 1 + |g 1 | 2 , and values of |g 2 | being still greater than, but rather close to 1 in the second one. However, |g 1 | = |g 2 | = 1 would only confirm the possibility of lasing features, but does not guarantee them. For such a purpose, higher order g n should be investigated. In the following, we instead propose to shift the focus from Glauber correlation functions to reconstructing the statistical properties of the FEL radiation related to the photon occupation numbers.

II. CONCEPTS AND METHODS
The quantum nature of the FEL radiation can be proved either with specific experiments aimed at highlighting indirect clues of its granularity [13] or by retrieving the actual photon-number distribution of the field [14]. In the latter case, photon-counting techniques have many drawbacks, as evidenced by [6] in the case of [5]. In general, they are limited to the detection of up to few tens of photons [15][16][17][18][19]. If the average number of photons increases, as in the case of high intensity sources, photon counting becomes challenging, even for attenuated FEL radiation. For the reconstruction purposes outlined above, one has then to resort to other approaches, such as Maximum Likelihood (ML) [20][21][22] methods based on no-click probabilities instead of photon counting. Finally, within the present context, it is interesting to note that the question concerning the degree to which the statistical occupation of the modes of a laser pulse is preserved after a nonlinear conversion of the photon energy remains an unknown aspect that could find an adequate answer via the theory and detection technique reported here.
In the following, we pursue the study of the photonic structure of the light emitted by FELs using the statistics of noclick events in Geiger-like photo detectors [23][24][25]. Such an approach has already been successfully investigated and experimentally tested in the optical regime, both for continuouswave and pulsed sources [4,26,27], but still limited to a mean number of only tens of photons. Therefore, to apply the technique successfully to (attenuated) FEL light, in the following section we show by simulated experiments that it can be extended to reconstruct the statistics of hundreds of photons. Notice that such a procedure provides information about the photon number statistics of the state of the light generated by the FEL, without any assumption about the nature of the FEL process itself. From this point of view, one expects a nearly Poissonian distribution of the occupation number in the case of seeded FEL light and a substantially thermal one in the SASE case. Shifting the focus from photon correlations to the statistics of occupation numbers would then provide a first step towards a reliable investigation of the quantum signatures of FEL light, ultimately able to discriminate between Poissonian and non-Poissonian features.
Notice that the experimental evidence of classical, non-Poissonian statistics, such as a thermal distribution, would exclude lasing properties. On the other hand, the experimental confirmation of the presence of Poissonian statistics, though not sufficient, would in any case provide solid ground to attributing "lasing" to the mechanism behind the generation of radiation by relativistic electrons wiggling in a series of magnetic devices. It would also give support to the conjecture that the lasing properties are imprinted upon the seeded-FEL radiation by the coherence features of the seeding pulse. Furthermore, an efficient applicability of the reconstruction mechanism to the diagonal occupation number statistics of the FEL radiation would boost further steps towards the full reconstruction [25] of the quantum states of the FEL radiation and prompt the investigation of which quantum features the emitted radiation would exhibit, should one use squeezed, namely genuinely quantum seeding pulses.
From a more general perspective, the consequences of the confirmation of the lasing properties of the seeded-FEL light will be of two types: it will open the way to experiments in quantum optics with X-ray radiation and it will foster the experimental investigation of light-matter interactions in the high energy quantum regime.

FEL radiation: quantum states
We start with an overview of the quantum states of light that are most suited to describe SASE and seeded-FEL radiation.
All light is in principle describable by non-commuting creation and annihilation operators a and a † of electromagnetic energy modes of frequency ω such that [a , a † ] = 1. The number operator is N = a † a and the so-called Fock number states with n photons, |n = (n!) − 1 2 (a † ) n |0 , satisfy N |n = n|n and constitute an orthonormal basis in the Fock Hilbert space associated with the chosen mode.
Let X = Tr ρ X denote the expectation value of an electromagnetic field operator X with respect to the field state (density matrix) ρ, which, in the number state representation, is in general a non-diagonal positive matrix of trace 1, where ρ k = k|ρ| . Let us consider a monochromatic field polarized along the z axis and propagating along the x direction, with energy ω ( = 1), wave-vector k = ω/c, and right, respectively, left moving field components At a fixed spatial position x, the Glauber correlation functions of order n [1], depend on n varying instants of times and are given by with where we have disregarded the spatial dependence. By suitably choosing the times t i , one also eliminates the timedependence and obtains [14] Monochromatic coherent states |α with complex amplitudes α ∈ C such that a|α = α|α , exhibit Poissonian occupation number distributions. They are generated by displacing the vacuum state |α = D(α)|0 by means of the displacement operator D(α) = exp αa † − α * a (see the Supplementary Material for some of their properties). Coherent states yield mean-values α|(a † ) n a n |α = |α| 2n whence |g n (0)| = 1 for all n ≥ 1. The same do the uniform averages of projectors onto coherent states over ϕ in α = |α|e iϕ , which are commonly used to describe the quantum states of laser light for which the phases are not accessible.
Instead, thermal states at temperature T setting the Boltzamnn constant κ B = 1), give N 2 − N = 2 N 2 and thus g 2 (0) = 2 corresponding to the incoherent character of such states of radiation. The randomness of thermal states can be lessened by displacing them: whence N = Tr ρ T,α a † a = n T + |α| 2 . Thus, reaches the minimum value 1 when n T = 0, namely for purely coherent states and the maximum value 2 when α = 0, that is for purely thermal states, or in the limit of infinite n T .
Displaced thermal states are natural candidates to describe the radiation emitted by FELs, which is due to aggregation in micro-bunches of relativistic electrons wiggling transversally while passing through a magnetic undulator. The microbunching is induced by the electrons interacting with the electric field they generate. SASE-FELs undulator radiation is the result of interaction between the electrons and the virtual photons of the static magnetic field. The spiky behavior is due to the lack of the causal connection between segments separated by cooperation length due to the different forward velocity between the electron bunch and the e.m. radiation field [32]. To lessen the initial randomness of the electron bunch, seeded-FELs instead use an external laser radiation source to "order" the electrons before the undulator, with the effect that the various micro-bunches emit in phase within a very narrow spectral band-width, thus introducing coherence with respect to the thermal case.

No-click quantum state reconstruction
As anticipated in the Introduction, this work proposes to shift focus from the Glauber correlation functions to the occupation number statistics. Such statistics are given by the diagonal entries ρ k := ρ kk of the quantum state of light ρ in the Fock number state representation (see (A5)). The statistics will be reconstructed by means of the empirically measured probability of no-click events, thus avoiding the drawbacks of photon counting. The principal tool is a Maximal Likelihood algorithm, as proposed first by E. Fermi [33], based on the no-click probabilities where η is the detector efficiency, namely the probability that it clicks when impinged by a single photon. The no-click probability for k impinging photons is given by (1 − η) k , whence (9) gives the total no-click probability.
Since higher photon-number states have reduced populations, one truncates the sum at N such that that can be used to infer the unknown parameters ρ n via a Maximum Estimation method (see the Supplementary Material). The number of efficiencies L being in general much smaller than N , the ensuing underestimation and the existence of more than one solution can be overcome by applying an external energy constraint β (see the Supplementary Material for a short discussion of this point).
In the non-monochromatic case with M > 1 different independent modes, the quantum state of the beam is of the form ρ = M j=1 ρ (j) , where ρ (j) denotes the state of the j-th mode. Then, the joint probability forn impinging photons and no clicks reads where n j is the occupation number of mode j and ρ (j) nj is the probability of having n j photons given the state ρ (j) . Then, the total multi-mode no-click probabilities become The reconstruction algorithm described in the Supplementary Material is such that, in the monochromatic case, once fed with the no-click probabilities P ({ρ n } n ) measured for sufficiently many detector efficiencies η , 1 ≤ ≤ L, it returns the populations ρ k of the k-photon Fock states of light, from k = 0 up to k = N (see the right plot in Figure 2). Instead, in the case of multi-mode radiation, from the no-click probabilities P M (η) the algorithm is able to reconstruct the probabilities in (B17), denoted by Π n on the y axis of the right plots in  To retrieve the single-mode probabilities ρ (j) nj and thus the diagonal elements of the multi-mode state ρ, further considerations are needed based on suitable physical assumptions about the structure of the radiation quantum state.

III. NUMERICAL SIMULATIONS
To check the reconstruction algorithm, one starts from known population distributions {ρ in n } n , them simulates the no-click probabilities for chosen efficiencies, and feeds them to the algorithm thus obtaining the reconstructed {ρ out n } n , from which the fidelity can be evaluated. In the following we focus on the simulation of two types of radiation quantum states: a single-mode displaced-thermal state and a multi-mode displaced-thermal state of the form with a common temperature T and different amplitudes α j . As previously discussed, the first kind of quantum state should represent well the light produced by a seeded FEL. For SASE-FELs instead, the common physical context at the origin of the micro-bunches suggests describing the corresponding quantum state of light by the multi-mode state (14)    same degree of randomness and hence the same temperature T , but with different coherent biases and thus with different displacement amplitudes α j .
Due to the one-dimensional feature of the problem, the generation of a vector of no-click probabilities can be obtained by generating events according to the probabilities P ({ρ n } n ) by means of the rejection sampling technique [30]. Then, the occurrences of the generated points, normalized to the total number of data, coincide with the desired no-click probabilities that, for sake of simplicity, are denoted by P on the y axis of the left plots in To mimic a realistic experiment, random perturbations have been added to the no-click probabilities. The stability of the algorithm has also been tested by corrupting the simulated noclick probabilities with systematic perturbations. The iteration procedure described in the Supplementary Material has then been fed with the no-click probabilities and the initial occupation numbers ρ n (0) taken all equal: ρ n (0) = 1/N . Finally, the consistency of the photon distribution reconstructed by the algorithm with the theoretical one is evaluated by means of the fidelity (13). As mentioned above, the algorithm has been applied to two distinct dataset scenarios. The dataset parameters have been  chosen in adherence to the performance of a short wavelength high gain FEL operating in the SASE [34,35] and in an external seeding [36] regime, respectively. The FEL pulse photon distribution has been assumed to be taken close to power saturation in both regimes.
Extreme ultraviolet seeded FELs are ideally capable of producing Fourier-transform limited light pulses, hence a single transverse and longitudinal mode. The natural relative spectral bandwidth can be as low as 10 −4 [37,38]. Thus, micronscale modulations of the electron energy distribution, accumulated during the acceleration process, can add a statistically significant pedestal to the FEL spectrum [39][40][41][42].
The ideal and the perturbed seeded FEL performance is represented in our algorithm through a single mode-displaced thermal state (M = 1, α = 16), and a multi-mode-displaced thermal state (M = 3, α = 7, 8,12). In both cases, the mean occupation number has been chosen n T = 30. The reconstructed no-click frequencies and photon distributions are compared with the theoretical and simulated data in Figures 2  and 3.
SASE-FELs are driven by shot noise in the electron distribution. Spontaneous undulator radiation is amplified and, eventually, a large number of longitudinal modes -each individually coherent -is emitted. A single transverse mode is instead selected by the amplification process close to saturation [43]. The number of longitudinal modes is basically given [32] by the portion of bunch length contribution to FEL emission (∼ 100 µm) divided by the FEL coherence length (∼ 1 − 3 µ in the extreme ultraviolet regime).
The SASE-FEL radiation has therefore been modelled choosing M = 30 and randomly distributed displacements. Moreover, two thermal configurations were considered to test the algorithm: in the first example the amplitudes α j are chosen in the interval 1 − 8 with n T = 10 and in the interval 1 − 5 with n T = 100 in the second one. Results are shown in Figures 4 and 5.
The overall procedure requires less than 200 iterations, to reach a fidelity level, 1−F , smaller than 10 −4 . Similarly high fidelities have been obtained for a larger number of datasets, thus demonstrating the capability of the algorithm to retrieve the input diagonal statistics in a wide range of configurations. In Figures 2-5, the left plots display a comparison between simulated, reconstructed and theoretical no-click frequencies, while the right ones compare reconstructed and theoretical photon distributions. More specifically, the left plots in each figure show the dependence of no-click probabilities on the efficiency η where the blue dots denote the simulated probabilities obtained with the sampling rejection technique. The red dots are outputs of the reconstruction algorithm; the solid line corresponds to the expression (B18) with chosen values for the parameters M and α. In the right plots of each figure, the red line shows the reconstructed photon distribution while the black line shows the photon distribution given by (B17). The fidelity F of the reconstruction is also reported.
The figures show the ability of the no-click reconstruction algorithm to efficiently reproduce the occupation number state statistics not only for low average photon numbers, as already demonstrated in the range of visible light, but also for higher ones as required by the FEL scenario. This numerical evidence points to the feasibility of an experiment where, without incurring the difficulties inherent in photon counting with many photons, one may access the diagonal entries of the quantum state of FEL radiation in the energy representation.

IV. DISCUSSION AND CONCLUSION
We have illustrated a novel approach to the study of the laser properties of the radiation emitted by an FEL by accessing not the electromagnetic correlation functions, but rather the diagonal photon statistics of the quantum state of the radiation in the energy representation. The occupation numbers are expected to follow thermal statistics in the case of SASE-FELs and that of a displaced thermal distribution in the case of seeded-FELs. The latter statistics are closer to that of laser light, with its the higher coherence with respect to its thermal features. Due to the large average number of photons expected in concrete experimental set-ups, the occupation number statistics can be hardly accessed via photon counting. Therefore, we suggest the use of a Maximum Likelihood reconstruction algorithm based on no-click probabilities at multiple, variable detector efficiencies. The reconstruction power of the algo-rithm have been tested in numerical experiments and proved to extend well to the high photon numbers typical of FEL light.
The numerical evidence points to the feasibility of an experimental implementation of the strategy described in the previous sections. As regards this future possibility, here we only observe (some more details are provided in the Supplementary Material) that monolithically integrated, complementary metal-oxide-semiconductor (CMOS) single photon avalanche photodiodes (SPADs) appear particularly well-suited to measure no-click probabilities due to their fast response and short dead-times, especially in the case of FELs with high repetition rates. Furthermore, arrays of photon counting pixels can be easily integrated on chip allowing pixels to be selectively turned off thus directly controlling the detection efficiency.
The experimental confirmation of a thermal distribution of the occupation numbers of the energy states of the SASE-FEL radiation would confirm its classical features in agreement with the fact that its purported squeezing is an artefact due to photon counting drawbacks. On the other hand, the experimental evidence of a Poissonian distribution of seeded-FEL radiation would have two main consequences. In the first place it would represent a first step towards the development of an X-ray quantum optics in connection with high-energy light-matter interactions. In the second place, it would spur the investigation of the very mechanism behind the seeded-FEL radiation, whether and how the seeding laser imprints its coherence on the emitted light and whether a sub-Poissonian seeding might induce the appearance of genuinely quantum features in it.  In this section we present quantum states of light that are relevant for describing the SASE and seeded-FEL radiation. A monocromatic, linearly polarized plane-wave of frequency ω is described by single-mode creation and annihilation operators a and a † satisfying the commutation relations [a , a † ] = 1. The absence of photons of frequency ω is associated with the vacuum state |0 which is annihilated by the annihilation operator: a|0 = 0, whereas the so-called number state, or Fock state, describing n photons with frequency ω is obtained by acting n times on the vacuum with the creation operator a † : The number states are orthogonal n|m = δ nm and such that a |n = √ n |n , a † |n = √ n + 1 |n + 1 , whence they possess definite photon number N = a † a and thus definite energy E = ωN ( = 1): To the other extreme with respect to Fock states, one finds the so-called coherent states |α , with complex amplitudes α ∈ C; they are eigenstates of the annihilation operator and their intensities are given by the their mean photon number: Any quantum state of light described by mode-operators a and a † is representable as a density matrix ρ acting on the Fock space space spanned by the number states, which in that representation reads: In turn, any such density matrix can always be written as [1]: where P ρ (α) is the so-called Glauber-Sudarshan P-function (see e.g. [2] for details) associated with the quantum state ρ and the integration is performed with respect to the real and imaginary parts of α = α x + iα y , d 2 α = dα x dα y . Therefore any single-mode state of light reads as a linear combination of projectors |α α| onto coherent states with a function P ρ (α) that, though being normalized, need not necessarily amount to a smooth probability distribution over the complex plane. Indeed, P ρ (α) may not only degenerate into a highly singular distribution, but even become non-positive. As a consequence, a quantum state of light is said to be • classical if the Glauber-Sudarshan P-function is a bona fide distribution with P ρ (α) a smooth function or no more singular than a Dirac delta.
For instance, as we shall see in Section A, for coherent states, that is when ρ = |β β|, P ρ (α) reduces to a Dirac delta δ 2 (α) = δ(α x )δ(α y ) on the complex plane, while, for thermal states P ρ (α) reduces to a Gaussian distribution. Instead, a quantum state of light is called • non-classical or, in other words, genuinely quantum, if the function P ρ (α) is either non-positive or more singular than a Dirac delta distribution.
Indeed, the non-classical behaviour of Fock number states is accompanied by a highly singular distributional Glauber-Sudarshan P ρ function, as shown in Appendix D, The P-function being an increasingly (with n) singular distribution, its behaviour has to be tested by integration against suitably differentiable trial functions g(α, α * ) vanishing at infinity, that is by computing Evidently, these integrations can yield negative values even for positive trial functions, thus the P-function P n (α) cannot be a positive distribution, whence the Fock number states are highly non-classical.
We shall now examine some quantum states associated with different degrees of coherence and discuss swhich ones among them are good candidates for the quantum description of the radiation emitted by SASE and seeded-FELs. Technical details on how to derive their mathematical properties are given in the Appendices.
Coherent states and laser radiation: As we have seen in (A4), coherent states are eigenstates of the annihilation operator and have thus the following expansion over the orthonormal basis of number states (A1): whence the associated number or energy distribution is Possonian: Though their quantum granularity is embodied by the Poissonian distribution (A12), coherent states ρ γ = |γ γ| are nevertheless classical saccording to the definition comprising equation (A7). Indeed, (see Appendix D), their Glauber-Sudarshan P-function P γ (α) (A6) is a Dirac delta at the point γ in the complex plane: Given the complex amplitude α = α x + iα y , its phase-angle ϕ = arctan α y /α x cannot be determined; therefore, the states of monocromatic laser light are obtained by averaging uniformly over all possible ϕ. Such a procedure yields a uniform mixture of coherent states, but also a Poissonian mixture of number state projectors [3]. Writing α ϕ = α e iϕ , α ≥ 0, Single mode thermal light: A single mode thermal state of light at temperature T is described by the Gibbs density matrix (setting the Boltzmann constant κ B = 1) Then, one computes As regards the thermal Glauber-Sudarshan P-function, one gets a P-function which is a Gaussian distribution (see (D12) in Appendix D) so that the corresponding quantum state can thus be interpreted as a classical one: Indeed, thermal light can be seen as a mixture of coherent states with respect to a Gaussian distribution in the amplitude modulus |α| = α 2 x + α 2 y of the amplitude α = |α| e iϕ = α x + iα y and a uniform distribution with respect to the phase-angle ϕ = arctan α y /α x : Displaced thermal states: The randomness of thermal states expressed by their Gaussian Glauber-Sudarshan P-function can be diminished by displacing them via the displacement operators D(α) in (D4) of Appendix D: Indeed, coherent states can also be obtained by displacing the vacuum state, |α = D(α)|0 ; therefore, by means of D(α), thermal states acquire a coherent contribution; this is particularly evident for vanishing temperatures, when n T vanishes as well and thermal states behave as the vacuum state which is turned into a fully coherent state by D(α) as in (D3). A more formal characterizaton of this interpretation follows by using (A18) and (D5); indeed, one gets where one sees that the displacement introduces a coherent bias into the Gaussian distribution that originates chaotic thermal light.
A concrete consequence of this coherent bias emerges when one looks at the mean energy which also gets a coherent contribution; indeed, using (D4), one computes a † a = Tr ρ T,α a † a = n T + |α| 2 . (A21)

Appendix B: No-click quantum state reconstruction
We now illustrate how the structure of a photon state diagonal with respect to the number representation can be reconstructed by means of the empirically measured probability of no-click events from photodetectors in on/off experiments with varying detector efficiencies η.

Monochromatic light
We first consider the case of a single mode radiation, a monochromatic beam of photons of frequency ω described by a density matrix as in (A5). If photons are detected by L detectors with different efficiencies η 1 , η 2 , . . . , η L , one collects L linear equations that, as we shall presently show, can be used to infer the N unknown parameters ρ n via a Maximum Likelihood method.
Remark 1. The number L of efficiencies to be considered has to be evaluated in reference with the specific reconstruction problem: in the case of visible light, as well as in the case of the simulations discussed in the main text, L = 20 appears to be sufficient. The fact that the number N of unknowns is larger than the number of equations L makes the problem underestimated with a non-unique solution, in general. This issue can be overcome with an external energy constraint as discussed below.
For each efficiency η one counts the no-click events N out of a total number N of impinging photons coming from a radiation beam described by a state ρ, thus collecting L empirical frequencies Usually, once the empirical frequencies f of a set of independent events x are measured, the event probabilities p := p(x ) are derived from the log-likelihood function where the Lagrange multiplier λ ensures normalization. Then, maximizing with respect to the unknown p one gets whence, from f = 1 = p , one retrieves λ = −1 and p = f , as expected. In the case of (B1), the unknown parameters to be found by optimization are the populations ρ n of the n-photon states, given the overall radiation state ρ. Then, one firstly introduces the normalized no-click probabilities and then the log-likelihood Setting ∂ ρn L ({ρ n } n , {f n } n , λ) = 0, from N k=0 ρ k = 1, one obtains λ = 0 and the following N equalities Setting where the quantities ρ k are positive but not normalized: Notice that the equality (B7) can be turned into a fixed point relation whence the probabilities can be achieved by the iteration procedure where are the probabilities at the h-th iteration step.
Due to the fact that in general the number of efficiencies L is much smaller than the number N of probabilities to be determined, the linear problem (B1) is underestimated. Especially in the multi-mode case to be treated in the next section, it could happen that more than one diagonal distribution {ρ n } n may fit the same experimental no-click probabilities. In order to sort the true distribution out, it proves convenient to constrain the energy through one more Lagrange multiplier β, changing the functional L ({ρ n } n , {f n } n , λ) in (B6) into The optimization procedure leads to the following constrained version of (B7) Then, the reconstruction algorithm proceeds as above, iteratively tuning β starting from β = 0 as suggested in [4], in order to eliminate configurations not compatible with the mean energy E.

Multi-mode case
In the non-monocromatic case, let us suppose the beam contains M > 1 different independent modes corresponding to photons with M different frequencies ω 1 , ω 2 , . . . , ω M . The quantum state of the beam is then factorized: where ρ (j) denotes the state associated with the j-th mode. Given a same detector efficiency η for all modes, the joint probability of havingn photons distributed over the M modes together with no clicks reads where n j , j = 1, 2, . . . , M , is the occupation number of the mode of frequency ω j , ρ (j) nj is the probability of having n j photons with that frequency in the state ρ (j) and the summation is over all possible occupation numbers of those modes conditioned on the total number of photons beingn. Then, the total probability of no-clicks is given by Notice that normalization asks for whence as in the one-mode case, one chooses a suitable accuracy parameter and truncates the summation atn max such that Then, varying over a number of quantum efficiencies η 1 , η 2 . . . , η L , L ≥n max , one obtains the following multimode system of equations similar to (B1): Πn(M, η j ) , j = 1, 2, . . . , L .
The recursive algorithm can then reproduce the probabilities Πn(M, η j ) by means of the no-click probabilities P M (η j ) after renormalization as in (B5). Notice that, by recursion, the probabilities Πn(η) can be turned into discrete nested convolutions Multi-mode coherent states As a concrete case, let us imagine that the photons are in coherent states of amplitudes α j relative to the modes in the beam; namely, from (A12), the occupation numbers read In such a case, from (B22) and (B18) one obtains namely a damped Poisson distribution, whence the no-click probabilities are the multimode Gaussians It is important to emphasize that the above result also holds when the coherent state projections ρ (j) are replaced by the uniform phase-averages ρ L in (A14); indeed, both these states have the same diagonal elements in the Fock number state representation. Therefore, P M (η) in (B26) represents the no-click probability for a multi-mode laser beam.

Multimode-displaced thermal states
Let us suppose that the M -mode state consist of thermal states ρ T,αj at a same temperature T but displaced by different amplitudes α j . The single-mode occupation probability of the n-th Fock number state is (see Appendix E) where L n (x) = L (0) n (x) are the Laguerre polynomials. Insert the displaced thermal states ρ T,αj in the place of all ρ (j) in (B21) and consider the last sum over n M −1 : Using (B27), one obtains .
Using the sum rule [5] L (α+β+1) one finally gets Iterating this procedure yields Thence the no-click probability at efficiency η results The last equality follows from using the associated Laguerre polynomials generating equation (B32)

Appendix C: Detectors
When considering the detectors that should measure no-click probabilities, several requirements need be satisfied. The detector should have fast response (< 1 ns) and a relatively small dead time (ideally < 100 ns) to accommodate for the laser dynamics, especially in the case of FELs with high repetition rates. The quantum efficiency should be high (> 10). One of the most employed single photon detectors are photomultiplier tubes [6]. The figures of merit for photomultipliers are close to those specified above, even if response times are often larger than 1 ns and dead times in the range between 100 ns and 1 µs. One disadvantage of photomultiplier tubes is the need to bias the devices at voltages of the order of 1 kV , making gated operation extremely difficult.
Recent years have seen the development of monolithically integrated CMOS detectors. Monolithic integration of the photosensitive area with the readout electronics offers unprecedented accuracy in the processing of the detected signal and engineering freedom on the sensor characteristics in term of size and electronic response. Fully depleted p-n junctions have been used for X-ray detection and have sensitivities that can reach the single photon level [7,8]. Different materials can be employed for the sensor area via wafer bonding [9].
Monolithic CMOS detectors can operate in Geiger mode in the case of CMOS single photon avalanche photodiodes (SPADs) [10]. CMOS SPADS have very fast response times (generally < 100 ps) with short dead times (of the order of some tens of ns) and can be operated at low voltage of the order of 10 V . Furthermore their integration can ensure that arrays of photon counting pixels can be easily integrated on chip to the side of the electronic circuits used for counting [11]. Such pixels can be as small as 10 µm in size, giving the possibility of engineering arrays comprising tens or hundreds of single photon counting pixels within the FEL spot diameter. The counts from such an array can then be summed by the on-chip electronics (in a configuration similar to silicon photomultipliers [12]) and pixels can be selectively turned off thus directly controlling the detection efficiency.
Appendix D: P-function: coherent, thermal and number states Coherent states form a continuous non-orthogonal family, which is however over-complete in the sense that the corresponding projectors integrate to the identity operator Coherent states are obtained by displacing the vacuum |0 , by means of the displacement operators D(α) = e α a † − α * a = e −|α| 2 /2 e α a † e −α * a , which are such that Let ρ be the density matrix of a monocromatic light beam; its Glauber-Sudarshan representation reads as a continuous linear combination of coherent-state projectors ρ = C d 2 α P ρ (α) |α α| .
Appendix E: Displaced thermal states representation Sometimes, instead of their Glauber-Sudarshan representation in terms of projectors onto coherent states, it proves convenient to represent quantum states by means of their characteristic function Tr(ρ D(−α)) and of the displacement operators. In the case of thermal states, this representation reads whence, using the algebraic relation (D5), it yields for displaced thermal states. These representations are obtained by firstly noticing that any single-mode density matrix ρ can be written as This follows from using D † (α) D(β) = e (α * β−αβ * )/2 D(β − α) , and the overcompleteness relation (D2) to prove that Tr D † (β) D(α) = C d 2 γ π γ|D † (β) D(α)|γ = π δ 2 (α − β) .