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

F. Benatti  Dipartimento di Fisica, Università degli Studi di Trieste, I-34151
Trieste, Italy  corresponding authors: Fabio Benatti<EMAIL_ADDRESS>theory and Fulvio Parmigiani<EMAIL_ADDRESS>experiment INFN,
Sezione di Trieste, I-34151 Trieste, Italy  S. Olivares  Dipartimento di
Fisica “Aldo Pontremoli,” Università degli Studi di Milano, I-20133 Milan,
Italy  INFN, Sezione di Milano, I-20133 Milan, Italy  G. Perosa  Dipartimento
di Fisica, Università degli Studi di Trieste, I-34151 Trieste, Italy  Elettra
Sincrotrone Trieste, I-34149 Basovizza, Trieste, Italy  D. Bajoni
Dipartimento di Ingegneria Industriale e dell’Informazione, Università degli
Studi di Pavia, I-27100 Pavia, Italy  S. Di Mitri  Dipartimento di Fisica,
Università degli Studi di Trieste, I-34151 Trieste, Italy  Elettra Sincrotrone
Trieste, I-34149 Basovizza, Trieste, Italy  R. Floreanini  INFN, Sezione di
Trieste, I-34151 Trieste, Italy  L. Ratti  Dipartimento di Ingegneria
Industriale e dell’Informazione, Università degli Studi di Pavia, I-27100
Pavia, Italy  INFN Sezione di Pavia, I-27100 Pavia, Italy  F. Parmigiani
Dipartimento di Fisica, Università degli Studi di Trieste, I-34151 Trieste,
Italy  Elettra Sincrotrone Trieste, I-34149 Basovizza, Trieste, Italy
International Faculty, University of Cologne, Albertus-Magnus-Platz, 50923
Cologne, Germany

###### Abstract

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

Figure 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.

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 Schmueser . 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 Preparata , 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 Becker
and Gjaja  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 explanation 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 Chen  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 Park
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 modulation 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 $2n$-point normalized Glauber correlation functions
$g_{n}$ equal to $1$ in modulus Glauber  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 Singer ; Roling  and seeded-FEL light Elettra1
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 Eberly  or
by retrieving the actual photon-number distribution of the field Mandel . In
the latter case, photon-counting techniques have many drawbacks, as evidenced
by Park  in the case of Chen . In general, they are limited to the detection
of up to few tens of photons Allevi:hyb ; chesi ; Perina12 ; Dovrat13 ;
Avenhaus10 . 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)
Dempster ; Vardi ; Boyles  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 no-click events in Geiger-like photo
detectors Mogilevtsev ; Brida11 ; Rossi . Such an approach has already been
successfully investigated and experimentally tested in the optical regime,
both for continuous-wave and pulsed sources Zambra ; Allevi ; Brida2 , 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 Brida11  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^{\dagger}$ of electromagnetic energy modes
of frequency $\omega$ such that $[a\,,\,a^{\dagger}]=1$. The number operator
is $N=a^{\dagger}a$ and the so-called Fock number states with $n$ photons,
$|n\rangle=(n!)^{-\frac{1}{2}}(a^{\dagger})^{n}|0\rangle$, satisfy
$N|n\rangle=n|n\rangle$ and constitute an orthonormal basis in the Fock
Hilbert space associated with the chosen mode.

Let $\langle\,X\,\rangle={\rm Tr}\Big{(}\rho\,X\Big{)}$ denote the expectation
value of an electromagnetic field operator $X$ with respect to the field state
(density matrix) $\rho$, which, in the number state representation, is in
general a non-diagonal positive matrix of trace $1$,

$\rho=\sum_{k,\ell=0}^{\infty}\rho_{k\ell}\,|k\rangle\langle\ell|\
,\quad\sum_{\ell=0}^{\infty}\rho_{\ell\ell}=1\ ,$ (1)

where $\rho_{k\ell}=\langle k|\rho|\ell\rangle$. Let us consider a
monochromatic field polarized along the $z$ axis and propagating along the $x$
direction, with energy $\omega$ ($\hbar=1$), wave-vector $k=\omega/c$, and
right, respectively, left moving field components

$E_{+}(x,t)=i\,\frac{\omega}{c}\,{\rm e}^{i(kx-\omega t)}\,a\ ,\quad
E_{-}(x,t)=E^{\dagger}_{+}(x,t)\ .$ (2)

At a fixed spatial position $x$, the Glauber correlation functions of order
$n$ Glauber , depend on $n$ varying instants of times and are given by

$g_{n}\left(\left\\{t_{i}\right\\}_{i=1}^{n}\right)\equiv\frac{G_{n}\left(\left\\{t_{i}\right\\}_{i=1}^{n}\right)}{\prod_{j_{=}1}^{2n}\sqrt{G_{1}(t_{j})}}\
,$ (3)

with

$G_{n}\left(\left\\{t_{i}\right\\}_{i=1}^{n}\right)=\left\langle\,E_{-}(t_{1})\cdots
E_{-}(t_{n})E_{+}(t_{n+1})\cdots E_{+}(t_{2n})\,\right\rangle,$ (4)

where we have disregarded the spatial dependence. By suitably choosing the
times $t_{i}$, one also eliminates the time-dependence and obtains Mandel

$g_{1}(0)=1\ ,\quad g_{2}(0)=1+\frac{\langle N^{2}\rangle-\langle
N\rangle^{2}-\langle N\rangle}{\langle N\rangle^{2}}\ .$ (5)

Monochromatic coherent states $|\alpha\rangle$ with complex amplitudes
$\alpha\in\mathbb{C}$ such that $a|\alpha\rangle=\alpha|\alpha\rangle$,
exhibit Poissonian occupation number distributions. They are generated by
displacing the vacuum state $|\alpha\rangle=D(\alpha)|0\rangle$ by means of
the displacement operator $D(\alpha)=\exp(\alpha a^{\dagger}-\alpha^{*}a)$
(see the Supplementary Material for some of their properties). Coherent states
yield mean-values
$\langle\alpha|(a^{\dagger})^{n}a^{n}|\alpha\rangle=|\alpha|^{2n}$ whence
$|g_{n}(0)|=1$ for all $n\geq 1$. The same do the uniform averages of
projectors onto coherent states over $\varphi$ in
$\displaystyle\alpha=|\alpha|{\rm e}^{i\varphi}$ , 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
$\kappa_{B}=1$),

$\hskip-8.5359pt\rho_{T}=\sum_{k=0}^{\infty}\frac{n^{k}_{T}}{(1+n_{T})^{k+1}}|k\rangle\langle
k|,\ n_{T}\equiv\langle N\rangle=\frac{1}{{\rm e}^{\omega/T}-1}\ ,$ (6)

give $\langle N^{2}\rangle-\langle N\rangle=2\,\langle N\rangle^{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:

$\rho_{T,\alpha}:=D(\alpha)\,\rho_{T}\,D^{\dagger}(\alpha)\ ,$ (7)

whence $\langle N\rangle={\rm
Tr}\left(\rho_{T,\alpha}\,a^{\dagger}a\right)=n_{T}+|\alpha|^{2}$. Thus,

$|g_{2}(0)|=1+\frac{n_{T}}{|\alpha|^{2}}\frac{1}{(1+\frac{n_{T}}{|\alpha|^{2}})^{2}}+\frac{n_{T}}{|\alpha|^{2}}\frac{1}{1+\frac{n_{T}}{|\alpha|^{2}}}\leq
2\ ,$ (8)

reaches the minimum value $1$ when $n_{T}=0$, namely for purely coherent
states and the maximum value $2$ when $\alpha=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 micro-bunching 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 Bonifacio2 .
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 $\rho_{k}:=\rho_{kk}$ of the
quantum state of light $\rho$ in the Fock number state representation (see
(19)). 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 Fermi , based on the no-click probabilities

$P(\\{\rho_{n}\\}_{n},\eta)=\sum_{k=0}^{\infty}(1-\eta)^{k}\,\rho_{k}\ ,$ (9)

where $\eta$ 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-\eta)^{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 $\sum_{k=1}^{N}\,\rho_{k}\geq 1-\epsilon$ for a
sufficiently small $\epsilon$. Photon detection at $L$ different efficiencies
$\\{\eta_{\ell}\\}_{\ell=1}^{L}$ yields a system of $L$ linear equations

$P_{\ell}(\\{\rho_{n}\\}_{n}):=P(\\{\rho_{n}\\}_{n},\eta_{\ell})=\sum_{k=0}^{N}(1-\eta_{\ell})^{k}\,\rho_{k}\
,$ (10)

that can be used to infer the unknown parameters $\rho_{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 $\beta$ (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
$\rho=\bigotimes_{j=1}^{M}\rho^{(j)}$, where $\rho^{(j)}$ denotes the state of
the $j$-th mode. Then, the joint probability for $\bar{n}$ impinging photons
and no clicks reads

$\Pi_{\bar{n}}(M,\eta)=\sum_{n_{1}+n_{2}+\cdots+n_{M}=\bar{n}}\,\left[\prod_{j=1}^{M}(1-\eta)^{n_{j}}\rho^{(j)}_{n_{j}}\right]\
,$ (11)

where $n_{j}$ is the occupation number of mode $j$ and $\rho^{(j)}_{n_{j}}$ is
the probability of having $n_{j}$ photons given the state $\rho^{(j)}$. Then,
the total multi-mode no-click probabilities become

$P_{M}(\eta)=\sum_{\bar{n}=0}^{\infty}\Pi_{\bar{n}}(M,\eta)\ .$ (12)

The reconstruction algorithm described in the Supplementary Material is such
that, in the monochromatic case, once fed with the no-click probabilities
$P_{\ell}(\\{\rho_{n}\\}_{n})$ measured for sufficiently many detector
efficiencies $\eta_{\ell}$, $1\leq\ell\leq L$, it returns the populations
$\rho_{k}$ of the $k$-photon Fock states of light, from $k=0$ up to $k=N$ (see
the right plot in Figure 3). Instead, in the case of multi-mode radiation,
from the no-click probabilities $P_{M}(\eta)$ the algorithm is able to
reconstruct the probabilities in (52), denoted by $\Pi_{n}$ on the $y$ axis of
the right plots in Figures 3-5. To retrieve the single-mode probabilities
$\rho^{(j)}_{n_{j}}$ and thus the diagonal elements of the multi-mode state
$\rho$, further considerations are needed based on suitable physical
assumptions about the structure of the radiation quantum state.

Figure 2: Single-mode displaced thermal state, $M=1$ $\alpha=16$, $n_{T}=30$.
Figure a: simulated no-click probabilities (blue dots), outputs of the
reconstruction algorithm (red dots), plot of (9) (solid line).
Figure b: reconstructed photon distribution (red line), photon distribution
from (9) (black line). Fidelity $F$ also shown.

Figure 3: Multi-mode displaced thermal state: $M=3$, $\alpha=[7,8,12]$ and
$n_{T}=30$. Figure a: simulated no-click probabilities (blue dots), outputs of
the reconstruction algorithm (red dots), plot of (53) (solid line).
Figure b: reconstructed photon distribution (red line), photon distribution
from (53) (black line). Fidelity $F$ also shown.

##  III Numerical simulations

To check the reconstruction algorithm, one starts from known population
distributions $\\{\rho^{in}_{n}\\}_{n}$, them simulates the no-click
probabilities for chosen efficiencies, and feeds them to the algorithm thus
obtaining the reconstructed $\\{\rho^{out}_{n}\\}_{n}$, from which the
fidelity

$F=\sum_{j=1}^{N}\sqrt{\rho^{in}_{j}\rho^{out}_{j}}\ ,$ (13)

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

$\rho=\bigotimes_{j=1}^{M}\rho^{(j)}_{T,\alpha_{j}}\ .$ (14)

with a common temperature $T$ and different amplitudes $\alpha_{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) with a same
degree of randomness and hence the same temperature $T$, but with different
coherent biases and thus with different displacement amplitudes $\alpha_{j}$.

Figure 4: Multimode-displaced thermal state: $M=30$, with $\alpha_{j}$ chosen
randomly in the range $[1,8]$ and $n_{T}=10$. Figure a: simulated no-click
probabilities (blue dots), outputs of the reconstruction algorithm (red dots),
plot of (53) (solid line). Figure b: reconstructed photon distribution (red
line), photon distribution from (53) (black line). Fidelity $F$ also shown.

Figure 5: Multimode-displaced thermal state: $M=30$, $\alpha_{j}$ chosen
randomly in the range $[1,5]$ and $n_{T}=100$.Figure a: simulated no-click
probabilities (blue dots), outputs of the reconstruction algorithm (red dots),
plot of (53) (solid line).
Figure b: reconstructed photon distribution (red line), photon distribution
from (53) (black line). Fidelity $F$ also shown.

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_{\ell}(\\{\rho_{n}\\}_{n})$ by means of the rejection
sampling technique SRM . 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 Figures 3-5. From datasets of thousands of values, we
extract a subset of $L$ elements with the corresponding efficiencies.
In the following, $L$ is chosen to be equal to $20$. Such choice is consistent
with previous tests of the algorithm Genovese ; Brida and with possible
experimental set-ups, as described below. The range of selected efficiencies
is chosen taking particular care to avoid regions of $\eta$ where the no-click
probabilities are flat. This condition leads to the exclusion of points with
low Fisher information.

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 no-click 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
$\rho_{n}(0)$ taken all equal: $\rho_{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
Bonifacio1 ; Kondratenko  and in an external seeding Yu  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}$
Allaria1 ; Allaria2 . Thus, micron-scale modulations of the electron energy
distribution, accumulated during the acceleration process, can add a
statistically significant pedestal to the FEL spectrum Marcus ; Hemsing ;
Zhang ; Perosa .

The ideal and the perturbed seeded FEL performance is represented in our
algorithm through a single mode-displaced thermal state $(M=1,\alpha=16)$, and
a multi-mode-displaced thermal state $(M=3,\alpha=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 3 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 Saldin . The number of longitudinal modes is basically given
Bonifacio2  by the portion of bunch length contribution to FEL emission ($\sim
100$ $\mu$m) divided by the FEL coherence length ($\sim 1-3$ $\mu$ 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
$\alpha_{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 5 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 3–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 $\eta$ 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 (53) with chosen values for the parameters $M$
and $\alpha$. 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 (52). 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 algorithm 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.

## authorship statements

Conception and design of the project: F.B., S.O., R.F., F. P.
Analysis and interpretation of research data: G.P., S.O.
Drafting and critically revising the research output: F.B.,F.P.
Contributions to the analysis and interpretation: R.F., S.D.M, D.B., L.R.
F.B. and F.P. are corresponding authors for the theoretical, respectively
experimental aspects of the manuscript.

## acknowledgments

F.P. would like to thanks William A. Barletta (MIT) for the many useful
discussions about the FEL physics relevant to the present study. F.B. and R.F.
acknowledge that their research has been conducted within the framework of the
Trieste Institute for Theoretical Quantum Technologies. S.D.M. and G.P.
acknowledge that the experimental spectra have been collected at the FERMI FEL
facility, within the framework of the FERMI Team operation.

## References

  * (1) P Schmüser, M Dohlus, J Rossbach, C Behrens, Free-Electron Lasers in the Ultraviolet and X-Ray Regime, Springer Tracts in Modern Physics 258, Springer Verlag (2014).
  * (2) G Preparata, Quantum Field Theory of the Free-electron Laser, Phys. Rev. A 38, 233 (1988).
  * (3) W Becker, MO Scully, MS Zubairy, Generation of Squeezed Coherent States via a Free-Electron Laser, Phys. Rev. Lett. 48, 475 (1982).
  * (4) I Gjaja, A Bhattacharjee, Generation of squeezed radiation from a free-electron laser, Phys. Rev. A 36, 5486 (1987).
  * (5) T Chen, J M J Madey, Observation of Sub-Poisson Fluctuations in the Intensity of the Seventh Coherent Spontaneous Harmonic Emitted by a rf Linac Free-Electron Laser, Phys. Rev. Lett. 86, 5906 (2001).
  * (6) JW Park, KJ Kim, R Lindberg, An Investigation of Possible Non-standard Photon Statistics in a Free-Electron Laser I: Experiment in 39th Free Electron Laser Conf. FEL-2019, Hamburg, Germany JACoW Publishing (2019).
  * (7) RJ Glauber, The Quantum Theory of Optical Coherence, Phys. Rev. 130, 2529 (1963).
  * (8) CC Gerry, PL Knight, Introductory quantum optics, (Cambridge University Press, 2005).
  * (9) A Singer et al., Hanbury Brown–Twiss Interferometry at a Free-Electron Laser, Phys. Rev. Lett. 111, 034802 (2013).
  * (10) S Roling et al., Temporal and spatial coherence properties of free-electron-laser pulses in the extreme ultraviolet regime, Phys. Rev. ST Accel. Beams 14, 080701 (2011).
  * (11) OY Gorobtsov et al., Seeded X-ray free-electron laser generating radiation with laser statistical properties, Nat. Comm. 9, 4498 (2018).
  * (12) OY Gorobtsov et al., Publisher correction: Seeded X-ray free-electron laser generating radiation with laser statistical properties, Nat. Comm. 10, 5217 (2019).
  * (13) JH Eberly, NB Narozhny, and JJ Sanchez-Mondragon, Periodic Spontaneous Collapse and Revival in a Simple Quantum Model, Phys. Rev. Lett. 44, 1323-1326 (1980).
  * (14) L Mandel and E Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, 1995).
  * (15) A Allevi and M Bondani, Statistics of twin-beam states by photon-number resolving detectors up to pump depletion, J. Opt. Soc. Am. B 31, B14-B19 (2014).
  * (16) G Chesi, et al, Measuring nonclassicality with silicon photomultipliers, Opt. Lett. 44, 1371-1374 (2019).
  * (17) J Peřina, Jr., M Hamar, V Michálek, and O Haderka, Photonnumber distributions of twin beams generated in spontaneous parametric down-conversion and measured by an intensified CCD camera, Phys. Rev. A 85, 023816 (2012).
  * (18) M Avenhaus, K Laiho, MV Chekhova, and C Silberhorn, Accessing higher order correlations in quantum optical states by time multiplexing, Phys. Rev. Lett. 104, 063602 (2010).
  * (19) L Dovrat et al., Direct observation of the degree of correlations using photon-number-resolving detectors, Phys. Rev. A 87, 053813 (2013).
  * (20) AP Dempster, NM Laird, and DB Rubin, Maximum Likelihood from Incomplete Data Via the EM Algorithm, J. R. Stat. Soc. Ser. B. Methodol. 39, 1-22 (1977).
  * (21) Y Vardi and D Lee, From Image Deblurring to Optimal Investments: Maximum Likelihood Solutions for Positive Linear Inverse Problems, J. R. Stat. Soc. Ser. B. Methodol. 55, 569-612 (1993).
  * (22) RA Boyles, On the convergence of the EM algorithm, J. R. Stat. Soc. Ser. B. Methodol., 45, 47–50 (1983).
  * (23) D Mogilevtsev, Diagonal element inference by direct detection, Opt. Comm. 156, 307–310 (1998).
  * (24) AR Rossi, S Olivares and MGA Paris, Photon statistics without counting photons, Phys. Rev. A 70, 05580 (2004).
  * (25) G Brida et al. Quantum State Reconstruction Using Binary Data from On/Off Photo-detection, Adv. Sci. Lett. 4, 1 (2011).
  * (26) G Zambra et al., Experimental Reconstruction of Photon Statistics without Photon Counting, Phys. Rev. Lett. 95, 063602 (2005).
  * (27) A Allevi et al., State reconstruction by on/off measurements, Phys. Rev A 80, 022114 (2009).
  * (28) G Brida et al., Constrained MaxLik reconstruction of multimode photon distributions, J. Mod. Opt. 56, 196-200 (2009).
  * (29) SJ van Enk, CA Fuchs, The Quantum State of a Propagating Laser Field, Phys. Rev. Lett. 88, 027902 (2001).
  * (30) MT Wells, G Casella, CP Robert, Generalized Accept-Reject sampling schemes, Institute of Mathematical Statistics Lecture Notes - Monograph Series 45, 342-347 (2004).
  * (31) M Genovese et al., Measuring the Photon Distribution with ON/OFF Photodetectors, Laser Physics 16, 385–392 (2006).
  * (32) R Bonifacio, BW McNeil, P Pierini, Superradiance in the high-gain free-electron laser, Phys. Rev. A 40, 4467–4475 (1989).
  * (33) J Orear, Notes on Statistics for Physicists, U. S. Atomic Energy Commission, (1958).
  * (34) R Bonifacio, C Pellegrini, LM Narducci, Collective instabilities and high-gain regime in a free electron laser, Opt. Commun. 50, 373-378 (1984).
  * (35) AM Kondratenko, EL Saldin, Generating of Coherent Radiation by a Relativistic Electron Beam in an Ondulator, Part. Accel. 10, 207–216 (1980).
  * (36) LH Yu, Generation of intense UV radiation by sub-harmonically seeded single-pass free-electron lasers, Phys. Rev. A 44, 5178 (1991).
  * (37) E Allaria et al., Highly coherent and stable pulses from the FERMI seeded free-electron laser in the extreme ultraviolet, Nat. Photonics 6, 699 (2012).
  * (38) E Allaria et al., Two-stage seeded soft-X-ray free-electron laser, Nat. Photonics 7, 913 (2013).
  * (39) G Marcus et al., Experimental observations of seed growth and accompanying pedestal contamination in a self-seeded, soft x-ray free-electron laser, Phys. Rev. Accel. Beams 22, 080702 (2019).
  * (40) E Hemsing, A Havalanau, and Z Zhang, Statistical theory of a self-seeded free electron laser with noise pedestal growth, Phys. Rev. Accel. Beams 23, 010701 (2020).
  * (41) Z Zhang, G Marcus, E Hemsing, WM Fawley, Z Huang, and A Lutman, Statistical analysis of a self-seeded x-ray free-electron laser in the presence of the microbunching instability, Phys. Rev. Accel. Beams 23, 010704 (2020).
  * (42) G Perosa et al., Linear optics control of sideband instability for improved free-electron laser spectral brightness, Phys. Rev. Accel. Beams 23, 110703 (2020).
  * (43) EL Saldin, EA Schneidmiller and MV Yurkov, Statistical and coherence properties of radiation from x-ray free-electron lasers, New J. Phys. 12, 035010 (2010).

##  Appendix A Quantum states of light

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 $\omega$ is
described by single-mode creation and annihilation operators $a$ and
$a^{\dagger}$ satisfying the commutation relations $[a\,,\,a^{\dagger}]=1$.
The absence of photons of frequency $\omega$ is associated with the vacuum
state $|0\rangle$ which is annihilated by the annihilation operator:
$a|0\rangle=0$, whereas the so-called number state, or Fock state, describing
$n$ photons with frequency $\omega$ is obtained by acting $n$ times on the
vacuum with the creation operator $a^{\dagger}$:

$|n\rangle=\frac{(a^{\dagger})^{n}}{\sqrt{n!}}|0\rangle\ .$ (15)

The number states are orthogonal $\langle n|m\rangle=\delta_{nm}$ and such
that

$a\,|n\rangle\,=\,\sqrt{n}\,|n\rangle\ ,\qquad
a^{\dagger}|n\rangle\,=\,\sqrt{n+1}\,|n+1\rangle\ ,$ (16)

whence they possess definite photon number $N=a^{\dagger}a$ and thus definite
energy $E=\omega N$ ($\hbar=1$):

$N\,|n\rangle\,=\,n\,|n\rangle\ .$ (17)

To the other extreme with respect to Fock states, one finds the so-called
coherent states $|\alpha\rangle$, with complex amplitudes
$\alpha\in\mathbb{C}$; they are eigenstates of the annihilation operator and
their intensities are given by the their mean photon number:

$a\,|\alpha\rangle=\alpha\,|\alpha\rangle\
,\quad|\alpha|^{2}=\langle\alpha|a^{\dagger}\,a|\alpha\rangle\ .$ (18)

Any quantum state of light described by mode-operators $a$ and $a^{\dagger}$
is representable as a density matrix $\rho$ acting on the Fock space space
spanned by the number states, which in that representation reads:

$\rho=\sum_{k,\ell}\rho_{k\ell}\,|k\rangle\langle\ell|\
,\quad\rho_{k\ell}=\langle k|\rho|\ell]\rangle\ .$ (19)

In turn, any such density matrix can always be written as Glauber :

$\rho=\int_{\mathbb{C}}{\rm
d}^{2}\alpha\,\mathcal{P}_{\rho}(\alpha)\,|\alpha\rangle\langle\alpha|\quad\hbox{with}\quad\int_{\mathbb{C}}{\rm
d}^{2}\alpha\,\mathcal{P}_{\rho}(\alpha)=1\,$ (20)

where $\mathcal{P}_{\rho}(\alpha)$ is the so-called Glauber-Sudarshan
P-function (see e.g. Gerry  for details) associated with the quantum state
$\rho$ and the integration is performed with respect to the real and imaginary
parts of $\alpha=\alpha_{x}+i\alpha_{y}$, ${\rm d}^{2}\alpha={\rm
d}\alpha_{x}\,{\rm d}\alpha_{y}$.

Therefore any single-mode state of light reads as a linear combination of
projectors $|\alpha\rangle\langle\alpha|$ onto coherent states with a function
$\mathcal{P}_{\rho}(\alpha)$ that, though being normalized, need not
necessarily amount to a smooth probability distribution over the complex
plane. Indeed, $\mathcal{P}_{\rho}(\alpha)$ 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

$\mathcal{P}_{\rho}(\alpha)\geq 0\ ,$ (21)

with $\mathcal{P}_{\rho}(\alpha)$ a smooth function or no more singular than a
Dirac delta.

For instance, as we shall see in Section Coherent states and laser radiation,
for coherent states, that is when $\rho=|\beta\rangle\langle\beta|$,
$\mathcal{P}_{\rho}(\alpha)$ reduces to a Dirac delta
$\delta^{2}(\alpha)=\delta(\alpha_{x})\delta(\alpha_{y})$ on the complex
plane, while, for thermal states $\mathcal{P}_{\rho}(\alpha)$ reduces to a
Gaussian distribution. Instead, a quantum state of light is called

  * •

non-classical or, in other words, genuinely quantum, if the function
$\mathcal{P}_{\rho}(\alpha)$ is either non-positive

$\mathcal{P}_{\rho}(\alpha)\ngeq 0\ ,$ (22)

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 $\mathcal{P}_{\rho}$
function, as shown in Appendix D,

$\mathcal{P}_{n}(\alpha)=\frac{\displaystyle{\rm
e}^{|\alpha|^{2}}}{n!}\frac{\partial^{2n}}{\partial_{\alpha}^{n}\partial_{\alpha^{*}}^{n}}\delta^{2}(\alpha)\
.$ (23)

The P-function being an increasingly (with $n$) singular distribution, its
behaviour has to be tested by integration against suitably differentiable
trial functions $g(\alpha,\alpha^{*})$ vanishing at infinity, that is by
computing

$\mathcal{P}_{n}[g]=\int_{\mathbb{C}^{2}}{\rm d}^{2}\alpha\,{\rm
d}^{2}\alpha^{*}\,\mathcal{P}_{n}(\alpha)\,g(\alpha,\alpha^{*})=\frac{\displaystyle{\rm
e}^{|\alpha|^{2}}}{n!}\frac{\partial^{2n}g(\alpha,\alpha^{*})}{\partial_{\alpha}^{n}\partial_{\alpha^{*}}^{n}}\
.$ (24)

Evidently, these integrations can yield negative values even for positive
trial functions, thus the P-function $\mathcal{P}_{n}(\alpha)$ 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 (18), coherent states are eigenstates of the annihilation
operator and have thus the following expansion over the orthonormal basis of
number states (15):

$\ket{\alpha}={\rm
e}^{-|\alpha|^{2}/2}\,\sum_{n=0}^{\infty}\frac{\alpha^{n}}{\sqrt{n!}}\,\ket{n}\
,\quad\alpha\in\mathbb{C}\ ,$ (25)

whence the associated number or energy distribution is Possonian:

$|\langle n|\alpha\rangle|^{2}={\rm
e}^{-|\alpha|^{2}}\frac{|\alpha|^{2n}}{n!}\ .$ (26)

Though their quantum granularity is embodied by the Poissonian distribution
(26), coherent states $\rho_{\gamma}=|\gamma\rangle\langle\gamma|$ are
nevertheless classical saccording to the definition comprising equation (21).
Indeed, (see Appendix D), their Glauber-Sudarshan P-function
$\mathcal{P}_{\gamma}(\alpha)$ (20) is a Dirac delta at the point $\gamma$ in
the complex plane:

$\mathcal{P}_{\gamma}(\alpha)=\delta^{2}(\alpha-\gamma)\ .$ (27)

Given the complex amplitude $\alpha=\alpha_{x}+i\alpha_{y}$, its phase-angle
$\varphi=\arctan\alpha_{y}/\alpha_{x}$ cannot be determined; therefore, the
states of monocromatic laser light are obtained by averaging uniformly over
all possible $\varphi$. Such a procedure yields a uniform mixture of coherent
states, but also a Poissonian mixture of number state projectors Enk . Writing
$\displaystyle\alpha_{\varphi}=\alpha\ {\rm e}^{i\varphi}$, $\alpha\geq 0$,

$\displaystyle\rho_{L}=\frac{1}{2\pi}\,\int_{0}^{2\pi}{\rm
d}\varphi\,\ket{\alpha_{\varphi}}\bra{\alpha_{\varphi}}={\rm
e}^{-\alpha^{2}}\,\sum_{m,n=0}^{\infty}\frac{\alpha^{m+n}}{\sqrt{m!n!}}\,\ket{n}\bra{m}\,\frac{1}{2\pi}\,\int_{0}^{2\pi}{\rm
d}\varphi\,{\rm e}^{i\varphi(n-m)}={\rm
e}^{-\alpha^{2}}\,\sum_{n=0}^{\infty}\frac{\alpha^{2n}}{n!}\,\ket{n}\bra{n}.$
(28)

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 $\kappa_{B}=1$)

$\rho_{T}=\left(1-{\rm e}^{-\omega/T}\right)\,{\rm
e}^{-\omega/T\,a^{\dagger}\,a}=\frac{1}{1+n_{T}}\sum_{n=0}^{\infty}\left(\frac{n_{T}}{1+n_{T}}\right)^{n}\,|n\rangle\langle
n|\ ,\ n_{T}\equiv\langle N\rangle=\frac{1}{{\rm e}^{\omega/T}-1}\ .$ (29)

Then, one computes

$\langle N^{2}\rangle=\frac{{\rm e}^{\omega/T}+1}{\Big{(}{\rm
e}^{\omega/T}-1\Big{)}^{2}}\ ,\qquad\langle N^{2}\rangle-\langle
N\rangle=2\,\langle N\rangle^{2}\ .$ (30)

As regards the thermal Glauber-Sudarshan P-function, one gets a P-function
which is a Gaussian distribution (see (79) in Appendix D) so that the
corresponding quantum state can thus be interpreted as a classical one:

$\mathcal{P}_{T}(\alpha)=\frac{1}{\pi}\frac{\displaystyle{\rm
e}^{-\frac{|\alpha|^{2}}{n_{T}}}}{n_{T}}\ .$ (31)

Indeed, thermal light can be seen as a mixture of coherent states with respect
to a Gaussian distribution in the amplitude modulus
$|\alpha|=\sqrt{\alpha_{x}^{2}+\alpha_{y}^{2}}$ of the amplitude
$\displaystyle\alpha=|\alpha|\,{\rm e}^{i\varphi}=\alpha_{x}+i\alpha_{y}$ and
a uniform distribution with respect to the phase-angle
$\varphi=\arctan\alpha_{y}/\alpha_{x}$:

$\rho_{T}=\int_{\mathbb{R}^{2}}{\rm d}\alpha_{x}{\rm
d}\alpha_{y}\,\frac{\displaystyle{\rm
e}^{-(\alpha_{x}^{2}+\alpha_{y}^{2})/n_{T}}}{\pi
n_{T}}\,|\alpha_{x}+i\alpha_{y}\rangle\langle\alpha_{x}+i\alpha_{y}|\ .$ (32)

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(\alpha)$ in (71) of Appendix D:

$\rho_{T,\alpha}:=D(\alpha)\,\rho_{T}\,D^{\dagger}(\alpha)\ .$ (33)

Indeed, coherent states can also be obtained by displacing the vacuum state,
$|\alpha\rangle=D(\alpha)|0\rangle$; therefore, by means of $D(\alpha)$,
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(\alpha)$ as in (70). A more formal characterizaton of this interpretation
follows by using (32) and (72); indeed, one gets

$\rho_{T,\alpha}=\int_{\mathbb{C}}{\rm d}^{2}\beta\frac{\displaystyle{\rm
e}^{-|\beta|^{2}/n_{T}}}{\pi
n_{T}}\,|\alpha+\beta\rangle\langle\alpha+\beta|=\int_{\mathbb{C}}{\rm
d}^{2}\beta\frac{\displaystyle{\rm e}^{-|\beta-\alpha|^{2}/n_{T}}}{\pi
n_{T}}\,|\beta\rangle\langle\beta|\ ,$ (34)

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 (71), one
computes

$\langle a^{\dagger}a\rangle={\rm
Tr}\left(\rho_{T,\alpha}\,a^{\dagger}a\right)=n_{T}+|\alpha|^{2}\ .$ (35)

##  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 $\eta$.

### Monochromatic light

We first consider the case of a single mode radiation, a monochromatic beam of
photons of frequency $\omega$ described by a density matrix as in (19). If
photons are detected by $L$ detectors with different efficiencies
$\eta_{1},\eta_{2},\ldots,\eta_{L}$, one collects $L$ linear equations

$P_{\ell}(\\{\rho_{n}\\}_{n})=\sum_{k=0}^{N}(1-\eta_{\ell})^{k}\,\rho_{k}\
,\quad\ell=1,2,\ldots,L\ ,$ (36)

that, as we shall presently show, can be used to infer the $N$ unknown
parameters $\rho_{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 $\eta_{\ell}$ one counts the no-click events
$\mathcal{N}_{\ell}$ out of a total number $\mathcal{N}$ of impinging photons
coming from a radiation beam described by a state $\rho$, thus collecting $L$
empirical frequencies

$f_{\ell}:=\frac{\mathcal{N}_{\ell}}{\mathcal{N}}\ .$ (37)

Usually, once the empirical frequencies $f_{\ell}$ of a set of independent
events $x_{\ell}$ are measured, the event probabilities
$p_{\ell}:=p(x_{\ell})$ are derived from the log-likelihood function

$L(\\{p_{n}\\}_{n},\\{f_{n}\\}_{n},\lambda)=\frac{1}{\mathcal{N}}\log\prod_{\ell}(p_{\ell})^{\mathcal{N}_{\ell}}\,+\,\lambda\Big{(}\sum_{\ell}p_{\ell}-1\Big{)}\
,$ (38)

where the Lagrange multiplier $\lambda$ ensures normalization. Then,
maximizing with respect to the unknown $p_{\ell}$ one gets

$\partial_{p_{\ell}}L(\\{p_{n}\\}_{n},\\{f_{n}\\}_{n},\lambda)=\frac{f_{\ell}}{p_{\ell}}\,+\,\lambda=0\
,$ (39)

whence, from $\sum_{\ell}f_{\ell}=1=\sum_{\ell}p_{\ell}$, one retrieves
$\lambda=-1$ and $p_{\ell}=f_{\ell}$, as expected.

In the case of (36), the unknown parameters to be found by optimization are
the populations $\rho_{n}$ of the $n$-photon states, given the overall
radiation state $\rho$. Then, one firstly introduces the normalized no-click
probabilities

$\widetilde{P}_{\ell}(\\{\rho_{n}\\}_{n}):=\frac{P_{\ell}(\\{\rho_{n}\\}_{n})}{P(\\{\rho_{n}\\}_{n})}\
,\quad\hbox{where}\quad
P(\\{\rho_{n}\\}_{n}):=\sum_{\ell=1}^{L}P_{\ell}(\\{\rho_{n}\\}_{n})\ ,$ (40)

and then the log-likelihood

$L\left(\\{\rho_{n}\\}_{n},\\{f_{n}\\}_{n},\lambda\right):=\sum_{\ell=1}^{L}\,f_{\ell}\,\log\widetilde{P}_{\ell}(\\{\rho_{n}\\}_{n})\,+\lambda\Big{(}\sum_{n=1}^{N}\rho_{n}-1\Big{)}\
.$ (41)

Setting
$\partial_{\rho_{n}}L\left(\\{\rho_{n}\\}_{n},\\{f_{n}\\}_{n},\lambda\right)=0$,
from $\sum_{k=0}^{N}\rho_{k}=1$, one obtains $\lambda=0$ and the following $N$
equalities

$\frac{P(\\{\rho_{n}\\}_{n})}{F}\,\sum_{\ell=1}^{L}\frac{(1-\eta_{\ell})^{n}}{\sum_{k=1}^{L}(1-\eta_{k})^{n}}\,\frac{f_{\ell}}{P_{\ell}(\\{\rho_{n}\\}_{n})}=1\
,\qquad F:=\sum_{\ell=1}^{L}f_{\ell}\\\ ,\quad n=1,2,\ldots,N\ .$ (42)

Setting

$\rho^{\prime}_{k}=\frac{F}{P(\\{\rho_{n}\\}_{n})}\,\rho_{k}\Longrightarrow
P_{\ell}(\\{\rho^{\prime}_{n}\\}_{n})=\frac{F}{P(\\{\rho_{n}\\}_{n})}\,P_{\ell}(\\{\rho_{n}\\}_{n})$
(43)

one recasts (42) as

$\sum_{\ell=1}^{L}\frac{(1-\eta_{\ell})^{n}}{\sum_{k=1}^{L}(1-\eta_{k})^{n}}\,\frac{f_{\ell}}{P_{\ell}(\\{\rho^{\prime}_{n}\\}_{n})}=1\
,$ (44)

where the quantities $\rho^{\prime}_{k}$ are positive but not normalized:

$\sum_{k=1}^{N}\rho^{\prime}_{k}=\frac{F}{P(\\{\rho_{n}\\}_{n})}\neq 1\ .$
(45)

Notice that the equality (42) can be turned into a fixed point relation

$\sum_{\ell=1}^{L}\frac{(1-\eta_{\ell})^{p}}{\sum_{k=1}^{L}(1-\eta_{k})^{p}}\,\frac{f_{\ell}}{P_{\ell}(\\{\rho^{\prime}_{n}\\}_{n})}\rho^{\prime}_{p}\,=\,\rho^{\prime}_{p}\
,$ (46)

whence the probabilities

$\rho_{p}=\frac{\rho^{\prime}_{p}}{\sum_{k=0}^{N}\rho^{\prime}_{k}}$ (47)

can be achieved by the iteration procedure

$\sum_{\ell=1}^{L}\frac{(1-\eta_{\ell})^{p}}{\sum_{j=1}^{L}(1-\eta_{j})^{p}}\,\frac{f_{\ell}}{P_{\ell}(\\{\rho_{n}(h)\\})_{n}}\,\frac{\rho_{p}(h)}{\sum_{k=1}^{N}\rho_{k}(h)}\,=\,\rho_{p}(h+1)\
,$ (48)

where $\displaystyle\frac{\rho_{p}(h)}{\sum_{k=0}^{N}\rho_{k}(h)}$ 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 (36)
is underestimated. Especially in the multi-mode case to be treated in the next
section, it could happen that more than one diagonal distribution
$\\{\rho_{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 $\beta$, changing the functional
$L\left(\\{\rho_{n}\\}_{n},\\{f_{n}\\}_{n},\lambda\right)$ in (41) into

$L_{\beta}\left(\\{\rho_{n}\\}_{n},\\{f_{n}\\}_{n},\lambda\right):=L\left(\\{\rho_{n}\\}_{n},\\{f_{n}\\}_{n},\lambda\right)+\beta\Big{(}\sum_{k=1}^{N}k\,\rho_{k}-E\Big{)}\
.$ (49)

The optimization procedure leads to the following constrained version of (42)

$\frac{P(\\{\rho_{n}\\}_{n})}{F}\,\sum_{\ell=1}^{L}\frac{(1-\eta_{\ell})^{p}}{\sum_{k=1}^{L}(1-\eta_{k})^{p}-\frac{P(\\{\rho_{n}\\}_{n})}{F}\beta(p-E)}\,\frac{f_{\ell}}{P_{\ell}(\\{\rho_{n}\\}_{n})}=1\
.$ (50)

Then, the reconstruction algorithm proceeds as above, iteratively tuning
$\beta$ starting from $\beta=0$ as suggested in Brida2 , 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
$\omega_{1},\omega_{2},\ldots,\omega_{M}$. The quantum state of the beam is
then factorized:

$\rho=\bigotimes_{j=1}^{M}\rho^{(j)}\ ,$ (51)

where $\rho^{(j)}$ denotes the state associated with the $j$-th mode. Given a
same detector efficiency $\eta$ for all modes, the joint probability of having
$\bar{n}$ photons distributed over the $M$ modes together with no clicks reads

$\Pi_{\bar{n}}(M,\eta)=\sum_{n_{1}+n_{2}+\cdots+n_{M}=\bar{n}}\left[\prod_{j=1}^{M}(1-\eta)^{n_{j}}\rho^{(j)}_{n_{j}}\right]\
,$ (52)

where $n_{j}$, $j=1,2,\ldots,M$, is the occupation number of the mode of
frequency $\omega_{j}$, $\rho^{(j)}_{n_{j}}$ is the probability of having
$n_{j}$ photons with that frequency in the state $\rho^{(j)}$ and the
summation is over all possible occupation numbers of those modes conditioned
on the total number of photons being $\bar{n}$. Then, the total probability of
no-clicks is given by

$P_{M}(\eta)=\sum_{\bar{n}=0}^{\infty}\Pi_{\bar{n}}(M,\eta)\ ,$ (53)

Notice that normalization asks for

$\sum_{\bar{n}=0}^{\infty}\ \sum_{n_{1}+n_{2}+\cdots+n_{M}=\bar{n}}\
\left[\prod_{j=1}^{M}\rho^{(j)}_{n_{j}}\right]=1\ ,$

whence as in the one-mode case, one chooses a suitable accuracy parameter
$\epsilon$ and truncates the summation at $\bar{n}_{max}$ such that

$\sum_{\bar{n}=0}^{\bar{n}_{max}}\sum_{n_{1}+n_{2}+\cdots+n_{M}=\bar{n}}\prod_{j=1}^{M}\rho^{(j)}_{n_{j}}\geq
1-\epsilon\ .$ (54)

Then, varying over a number of quantum efficiencies
$\eta_{1},\eta_{2}\,\ldots,\eta_{L}$, $L\geq\bar{n}_{max}$, one obtains the
following multimode system of equations similar to (36):

$P_{M}(\eta_{j})=\sum_{\bar{n}=0}^{\bar{n}_{max}}\Pi_{\bar{n}}(M,\eta_{j})\
,\quad j=1,2,\ldots,L\ .$ (55)

The recursive algorithm can then reproduce the probabilities
$\Pi_{\bar{n}}(M,\eta_{j})$ by means of the no-click probabilities
$P_{M}(\eta_{j})$ after renormalization as in (40).

Notice that, by recursion, the probabilities $\Pi_{\bar{n}}(\eta)$ can be
turned into discrete nested convolutions

$\displaystyle\Pi_{\bar{n}}(M,\eta)$ $\displaystyle=$
$\displaystyle\sum_{n_{1}=0}^{\bar{n}}(1-\eta)^{n_{1}}\rho^{(1)}_{n_{1}}\,\Pi_{\bar{n}-n_{1}}(M-1,\eta)=\sum_{n_{1}=0}^{\bar{n}}\sum_{n_{2}=0}^{\bar{n}-n_{1}}(1-\eta)^{n_{1}+n_{2}}\,\rho^{(1)}_{n_{1}}\rho^{(2)}_{n_{2}}\,\Pi_{\bar{n}-n_{1}-n_{2}}(M-2,\eta)$
(56) $\displaystyle=$
$\displaystyle(1-\eta)^{\bar{n}}\sum_{n_{1}=0}^{\bar{n}}\sum_{n_{2}=0}^{\bar{n}-n_{1}}\cdots\sum_{n_{M-1}=0}^{\bar{n}-\sum_{k=1}^{M-2}n_{k}}\,\prod_{j=1}^{M}\rho^{(j)}_{n_{j}}\
,\quad n_{M}=\bar{n}-\sum_{j=1}^{M-1}n_{j}\ .$ (57)

#### Multi-mode coherent states

As a concrete case, let us imagine that the photons are in coherent states of
amplitudes $\alpha_{j}$ relative to the modes in the beam; namely, from (26),
the occupation numbers read

$\rho^{(j)}_{n_{j}}={\rm
e}^{-|\alpha_{j}|^{2}}\,\frac{|\alpha_{j}|^{2n_{j}}}{n_{j}!}\ .$ (58)

In such a case, from (57) and (53) one obtains

$\displaystyle\Pi_{\bar{n}}(M,\eta)$ $\displaystyle=$
$\displaystyle(1-\eta)^{\bar{n}}{\rm
e}^{-\sum_{j=1}^{M}|\alpha_{j}|^{2}}\,\sum_{n_{1}=0}^{N}\sum_{n_{2}=0}^{N-n_{1}}\cdots\sum_{n_{M-1}=0}^{\bar{n}-\sum_{k=1}^{M-2}n_{k}}\,\prod_{j=1}^{M}\frac{|\alpha_{j}|^{2n_{j}}}{n_{j}!}\
,\quad n_{M}=\bar{n}-\sum_{j=1}^{M-1}\,n_{j}$ (59) $\displaystyle=$
$\displaystyle{\rm
e}^{-\sum_{j=1}^{M}|\alpha_{j}|^{2}}\,\frac{(1-\eta)^{\bar{n}}}{\bar{n}!}\,\Big{(}\sum_{k=1}^{M}|\alpha_{j}|^{2}\Big{)}^{\bar{n}}\
,$ (60)

namely a damped Poisson distribution, whence the no-click probabilities are
the multimode Gaussians

$P_{M}(\eta)=\sum_{\bar{n}=0}^{\bar{n}_{max}}\Pi_{\bar{n}}(M,\eta)\simeq{\rm
e}^{-\eta\,\sum_{j=1}^{M}|\alpha_{j}|^{2}}\ .$ (61)

It is important to emphasize that the above result also holds when the
coherent state projections $\rho^{(j)}$ are replaced by the uniform phase-
averages $\rho_{L}$ in (28); indeed, both these states have the same diagonal
elements in the Fock number state representation. Therefore, $P_{M}(\eta)$ in
(61) 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
$\rho_{T,\alpha_{j}}$ at a same temperature $T$ but displaced by different
amplitudes $\alpha_{j}$. The single-mode occupation probability of the $n$-th
Fock number state is (see Appendix E)

$\rho_{T,\alpha}(n)=\left(\frac{n_{T}}{1+n_{T}}\right)^{n}\,\frac{\displaystyle{\rm
e}^{-\frac{|\alpha|^{2}}{1+n_{T}}}}{1+n_{T}}\,L_{n}\left(-\frac{|\alpha|^{2}}{n_{T}(1+n_{T})}\right)\
,$ (62)

where $L_{n}(x)=L^{(0)}_{n}(x)$ are the Laguerre polynomials.

Insert the displaced thermal states $\rho_{T,\alpha_{j}}$ in the place of all
$\rho^{(j)}$ in (56) and consider the last sum over $n_{M-1}$:

$\Sigma_{M}:=\sum_{n_{M-1}=0}^{\bar{n}-\sum_{k=1}^{M-2}n_{k}}\,\rho^{(M-1)}_{n_{M-1}}\rho^{(M)}_{\bar{n}-\sum_{k=1}^{M-1}n_{k}}\
.$

Using (62), one obtains

$\displaystyle\Sigma_{M}$ $\displaystyle=$
$\displaystyle\frac{1}{(1+n_{T})^{2}}\,{\rm
e}^{-2\frac{|\alpha|^{2}}{1+n_{T}}}\,\left(\frac{n_{T}}{1+n_{T}}\right)^{\bar{n}-\sum_{k=1}^{M-2}n_{k}}\,\times$
$\displaystyle\times$
$\displaystyle\sum_{n_{M-1}=0}^{\bar{n}-\sum_{k=1}^{M-2}n_{k}}\,L_{n_{M-1}}\left(-\frac{|\alpha_{M-1}|^{2}}{n_{T}(1+n_{T})}\right)\,L_{\bar{n}-\sum_{k=1}^{M-2}n_{k}\,-\,n_{M-1}}\left(-\frac{|\alpha_{M}|^{2}}{n_{T}(1+n_{T})}\right)\
.$

Using the sum rule Gradshtein

$L_{n}^{(\alpha+\beta+1)}(x+y)=\sum_{k=0}^{n}L^{(\alpha)}_{k}(x)\,L^{(\beta)}_{n-k}(y)\
,$ (63)

one finally gets

$\Sigma_{M}=\frac{1}{(1+n_{T})^{2}}\,{\rm
e}^{-2\frac{|\alpha|^{2}}{1+n_{T}}}\,\left(\frac{n_{T}}{1+n_{T}}\right)^{\bar{n}-\sum_{k=1}^{M-2}n_{k}}\,L^{(1)}_{\bar{n}-\sum_{k=1}^{M-2}n_{k}\,-\,n_{M-1}}\left(-\,\frac{|\alpha_{M-1}|^{2}+|\alpha_{M}|^{2}}{n_{T}(1+n_{T})}\right)\
.$ (64)

Iterating this procedure yields

$\Pi_{\bar{n}}(M,\eta)=\frac{1}{(1+n_{T})^{M}}\,{\rm
e}^{-\frac{\sum_{j=1}^{M}\,|\alpha_{j}|^{2}}{1+n_{T}}}\,\left(\frac{(1-\eta)n_{T}}{1+n_{T}}\right)^{\bar{n}}\,L^{(M-1)}_{\bar{n}}\left(-\frac{\sum_{j=1}^{M}|\alpha_{j}|^{2}}{n_{T}(1+n_{T})}\right)\
.$ (65)

Thence the no-click probability at efficiency $\eta$ results

$P_{M}(\eta)\simeq\sum_{\bar{n}=0}^{\infty}P^{(M)}_{\bar{n}}(\eta)=\frac{1}{(1+\eta
n_{T})^{M}}\,\exp\Big(-\frac{\eta\,\sum_{j=1}^{M}\,|\alpha_{j}|^{2}}{1+\eta\,n_{T}}\Big{missing})\
.$ (66)

The last equality follows from using the associated Laguerre polynomials
generating equation

$(1-t)^{-1-\lambda}\,\exp\Big(-\frac{x\,t}{1-t}\Big{missing})=\sum_{n=0}^{\infty}L^{(\lambda)}_{n}(x)\,t^{n}\
.$ (67)

##  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 Howe2006 . 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\,\mu 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
Giampaolo2019 ; Hirono2016 . Different materials can be employed for the
sensor area via wafer bonding Ruat2012 .

Monolithic CMOS detectors can operate in Geiger mode in the case of CMOS
single photon avalanche photodiodes (SPADs) ghioni17 . 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 braga14 . Such pixels can be as small as $10\,\mu
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 Conca2020 ) 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,

$\langle\alpha|\beta\rangle={\rm
e}^{\displaystyle-\frac{1}{2}(|\alpha|^{2}+|\beta|^{2}-2\alpha^{*}\beta)}\ ,$
(68)

which is however over-complete in the sense that the corresponding projectors
integrate to the identity operator

$\int_{\mathbb{C}}\frac{{\rm
d}^{2}\alpha}{\pi}\,\ket{\alpha}\bra{\alpha}=\mathbb{I}\ .$ (69)

Coherent states are obtained by displacing the vacuum $|0\rangle$,

$D(\alpha)\ket{0}=\ket{\alpha}\ ,$ (70)

by means of the displacement operators

$D(\alpha)={\rm e}^{\alpha\,a^{\dagger}\,-\,\alpha^{*}\,a}={\rm
e}^{-|\alpha|^{2}/2}\,{\rm e}^{\alpha\,a^{\dagger}}\,{\rm e}^{-\alpha^{*}\,a}\
,$ (71)

which are such that

$D(\beta)\,D(\alpha)={\rm
e}^{(\beta^{*}\alpha-\beta\alpha^{*})/2}\,D(\alpha+\beta)\ ,\quad
D^{\dagger}(\beta)\,D(\alpha)\,D(\beta)={\rm
e}^{\beta^{*}\alpha-\beta\alpha^{*}}\,D(\alpha)\ ,\quad
D^{\dagger}(\alpha)\,a\,D(\alpha)=a+\alpha\ .$ (72)

Let $\rho$ be the density matrix of a monocromatic light beam; its Glauber-
Sudarshan representation reads as a continuous linear combination of coherent-
state projectors

$\rho=\int_{\mathbb{C}}{\rm
d}^{2}\alpha\,\mathcal{P}_{\rho}(\alpha)\,|\alpha\rangle\langle\alpha|\ .$

Using the scalar product (68), with $\alpha=x+iy$ and $\beta=u+iv$, one
obtains

$\displaystyle{\rm e}^{|\beta|^{2}}\,\langle-\beta|\rho|\beta\rangle$
$\displaystyle=$ $\displaystyle\int_{\mathbb{C}}{\rm
d}^{2}\alpha\,\mathcal{P}_{\rho}(\alpha)\,{\rm
e}^{-|\alpha|^{2}-\alpha]\beta^{*}+\alpha^{*}\beta}$ (73) $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{2}}{\rm d}x{\rm
d}y\,\mathcal{P}_{\rho}(x+iy)\,{\rm e}^{-x^{2}-y^{2}}\,{\rm e}^{2i(xv-yu)}\ ,$

whence, by anti-Fourier transformation,

$\int_{\mathbb{R}^{2}}\frac{{\rm d}(2u){\rm d}(2v)}{4\pi^{2}}\,{\rm
e}^{u^{2}+v^{2}}\,\langle-u-iv|\rho|u+iv\rangle\,{\rm e}^{2i(xv-yu)}={\rm
e}^{-x^{2}-y^{2}}\,\mathcal{P}_{\rho}(x+iy)\ .$ (74)

Then, the Glauber-Sudarshan P-function $\mathcal{P}_{\rho}(\alpha)$ can be
expressed by means of the matrix elements $\langle-\beta|\rho|\beta\rangle$:

$\mathcal{P}_{\rho}(\alpha)=\frac{\displaystyle{\rm
e}^{|\alpha|^{2}}}{\pi^{2}}\int_{\mathbb{C}}{\rm d}^{2}\beta\,{\rm
e}^{|\beta|^{2}}\,{\rm
e}^{\alpha\beta^{*}-\alpha^{*}\beta}\,\langle-\beta|\rho|\beta\rangle\ ,$ (75)

where $|\pm\beta\rangle$ are coherent states.

If $\rho$ is itself a coherent state
$\rho_{\gamma}=|\gamma\rangle\langle\gamma|$, the scalar product (68) yields

$\langle-\beta|\rho_{\gamma}|\beta\rangle={\rm
e}^{-|\beta|^{2}-|\gamma|^{2}-\gamma\beta^{*}+\gamma^{*}\beta}\ .$ (76)

Using (75) and the fact that, with $\alpha=x+iy$ and $\beta=u+iv$,

$\frac{1}{\pi^{2}}\int_{\mathbb{C}}{\rm d}^{2}\beta\,{\rm
e}^{\alpha\beta^{*}-\alpha^{*}\beta}=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\frac{{\rm
d}x}{\pi}\frac{{\rm d}y}{\pi}{\rm e}^{2i(yu-
xv)}=\delta(x)\delta(y)=\delta^{2}(\alpha)\ ,$ (77)

one then obtains

$\mathcal{P}_{\rho_{\gamma}}(\alpha)=\frac{\displaystyle{\rm
e}^{|\alpha|^{2}-|\gamma|^{2}}}{\pi^{2}}\int_{\mathbb{C}^{2}}{\rm d}^{2}\beta\
{\rm
e}^{-(\gamma-\alpha)\beta^{*}+(\gamma^{*}-\alpha^{*})\beta}=\delta^{2}(\alpha-\gamma)\
.$ (78)

For thermal states as in (29), by Gaussian integration one obtains

$\langle-\beta|\rho_{T}|\beta\rangle=\frac{\displaystyle{\rm
e}^{-|\beta|^{2}\frac{1+2n_{T}}{1+n_{T}}}}{1+n_{T}}\Longrightarrow\mathcal{P}_{T}(\alpha)=\frac{1}{\pi}\frac{\displaystyle{\rm
e}^{-\frac{|\alpha|^{2}}{n_{T}}}}{n_{T}}\ .$ (79)

On the other hand, when $\rho=|n\rangle\langle n|$, one computes

$\langle-\beta|\rho|\beta\rangle=\frac{\displaystyle{\rm
e}^{-|\beta|^{2}}}{\pi^{2}}(-)^{n}\frac{|\beta|^{2}}{n!}\ ,$ (80)

whence, using (77), the P-function becomes a highly singular distribution:

$\mathcal{P}_{n}(\alpha)=\frac{\displaystyle{\rm
e}^{|\alpha|^{2}}}{\pi^{2}}\int_{\mathbb{C}}{\rm
d}^{2}\beta\,(-)^{n}\frac{|\beta|^{2}}{n!}\,{\rm
e}^{\alpha\beta^{*}-\alpha^{*}\beta}=\frac{\displaystyle{\rm
e}^{|\alpha|^{2}}}{\pi^{2}\,n!}\frac{\partial^{2n}}{\partial_{\alpha}^{n}\partial_{\alpha^{*}}^{n}}\int_{\mathbb{C}}{\rm
d}^{2}\beta\,{\rm e}^{\alpha\beta^{*}-\alpha^{*}\beta}=\frac{\displaystyle{\rm
e}^{|\alpha|^{2}}}{n!}\frac{\partial^{2n}}{\partial_{\alpha}^{n}\partial_{\alpha^{*}}^{n}}\delta^{2}(\alpha)\
.$ (81)

##  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 ${\rm Tr}(\rho\,D(-\alpha))$
and of the displacement operators. In the case of thermal states, this
representation reads

$\rho_{T}=\int_{\mathbb{C}}{\rm d}^{2}\beta\,\frac{1}{\pi}{\rm
e}^{-|\beta|^{2}\,(\frac{1}{2}+n_{T})}\ D(\beta)\ ,$ (82)

whence, using the algebraic relation (72), it yields

$\rho_{T,\alpha}=\int_{\mathbb{C}}{\rm d}^{2}\beta\,\frac{1}{\pi}{\rm
e}^{-|\beta|^{2}\,(\frac{1}{2}+n_{T})}\ {\rm
e}^{\beta\alpha^{*}-\beta^{*}\alpha}D(\beta)$ (83)

for displaced thermal states. These representations are obtained by firstly
noticing that any single-mode density matrix $\rho$ can be written as

$\rho=\int_{\mathbb{C}}\frac{{\rm d}^{2}\beta}{\pi}\,{\rm
Tr}\left(\rho\,D^{\dagger}(\alpha)\right)\,D(\alpha)\ .$ (84)

This follows from using

$D^{\dagger}(\alpha)\ D(\beta)={\rm e}^{(\alpha^{*}\beta-\alpha\beta^{*})/2}\
D(\beta-\alpha)\ ,$

and the overcompleteness relation (69) to prove that

${\rm
Tr}\left(D^{\dagger}(\beta)\,D(\alpha)\right)=\int_{\mathbb{C}}\frac{{\rm
d}^{2}\gamma}{\pi}\,\langle\gamma|D^{\dagger}(\beta)\,D(\alpha)|\gamma\rangle=\pi\,\delta^{2}(\alpha-\beta)\
.$

Then, using the P-function (31) and the Glauber-Sudarshan representation
together with the coherent state overcompleteness and Gaussian integration,
yield

${\rm
Tr}\left(D^{\dagger}(\beta)\,\rho_{T}\right)=\frac{1}{\pi\,n_{T}}\int_{\mathbb{C}}{\rm
d}^{2}\gamma\ {\rm
e}^{-|\gamma|^{2}/n_{T}}\,\langle\gamma|D^{\dagger}(\beta)|\gamma\rangle=\exp\left(-|\beta|^{2}\Big{(}\frac{1}{2}+n_{T}\Big{)}\right)\
.$ (85)

Given a displacement operator $D(\alpha)$ and a coherent state
$|\beta\rangle=D(\beta)|0\rangle$ (see (71)), using (69) and then (25), one
computes in two ways the scalar products with Fock number states $\ket{n}$,
$a^{\dagger}a\ket{n}=n\ket{n}$:

$\displaystyle\langle n|D(\alpha)|\beta\rangle$ $\displaystyle=$
$\displaystyle{\rm e}^{(\alpha\beta^{*}-\alpha^{*}\beta)/2}\,\langle
n|D(\alpha+\beta)|0\rangle=\frac{(\alpha+\beta)^{n}}{\sqrt{n!}}\,{\rm
e}^{-(|\alpha|^{2}+|\beta|^{2}+2\alpha^{*}\beta)/2}$ $\displaystyle=$
$\displaystyle{\rm
e}^{-|\beta|^{2}/2}\,\sum_{m=0}^{\infty}\frac{\beta^{m}}{\sqrt{m!}}\,\langle
n|D(]\alpha)|m\rangle\ .$

It thus follows that

$(\alpha+\beta)^{n}\,{\rm e}^{-\alpha^{*}\beta}={\rm
e}^{|\alpha|^{2}/2}\sum_{m=0}^{\infty}\sqrt{\frac{n!}{m!}}\beta^{m}\,\langle
n|\,D(\alpha)\,|m\rangle\ .$

By setting $\beta=\alpha\,y$, the equality reads

$(1+y)^{n}\,{\rm e}^{-y|\alpha|^{2}}=\sum_{m=0}^{\infty}\left({\rm
e}^{|\alpha|^{2}/2}\,\alpha^{m-n}\,\sqrt{\frac{n!}{m!}}\,\langle
n|\,D(\alpha)\,|m\rangle\right)\,y^{m}\ ,$

whose right hand side is recognizable as one of the generating functions of
the associated Laguerre polynomials $L^{(\alpha)}_{n}(y)$ Gradshtein ,

$(1+t)^{\lambda}\,{\rm
e}^{-xt}=\sum_{m=0}^{\infty}L^{(\lambda-m)}_{m}(x)\,t^{m}\ .$

The identification

${\rm e}^{|\alpha|^{2}/2}\alpha^{m-n}\sqrt{\frac{n!}{m!}}\,\langle
n|D(\alpha)|m\rangle=L^{(n-m)}_{m}(|\alpha|^{2})$ (86)

finally yields

$\langle n|D(\alpha)|m\rangle={\rm
e}^{-|\alpha|^{2}/2}\alpha^{n-m}\sqrt{\frac{m!}{n!}}\,L^{(n-m)}_{m}(|\alpha|^{2})\
.$ (87)

From (83) and (87) one gets

$\displaystyle\rho_{T,\alpha}(n)$ $\displaystyle=$ $\displaystyle\langle
n|\rho_{T,\alpha}|n\rangle=\int_{\mathbb{C}}{\rm
d}^{2}\beta\,\frac{1}{\pi}{\rm e}^{-|\beta|^{2}\,(\frac{1}{2}+n_{T})}\ {\rm
e}^{\beta\alpha^{*}-\beta^{*}\alpha}\langle n|D(\beta)|n\rangle$ (88)
$\displaystyle=$ $\displaystyle\int_{\mathbb{C}}{\rm
d}^{2}\beta\,\frac{1}{\pi}{\rm e}^{-|\beta|^{2}\,(\frac{1}{2}+n_{T})}\ {\rm
e}^{\beta\alpha^{*}-\beta^{*}\alpha}{\rm
e}^{-|\alpha|^{2}/2}\,L_{n}(|\alpha|^{2})\ ,$

where $L_{n}(x)=L^{(0)}_{n}(x)$ are the Laguerre polynomials.

By passing to polar coordinates, ${\rm d}^{2}\beta=r{\rm d}r{\rm d}\varphi$
and setting $\displaystyle\alpha=|\alpha|{\rm e}^{i\psi}$, one rewrites

$\rho_{T,\alpha}(n)=\int_{0}^{\infty}\frac{{\rm d}}{\pi}\,r\,{\rm
e}^{-r^{2}(1+n_{T})}\,L_{n}(r^{2})\,\int_{0}^{2\pi}{\rm d}\varphi\,{\rm
e}^{2ir|\alpha|\sin(\varphi-\psi))}\ .$ (89)

The last integral equals the Bessel function $J_{0}(2r|\alpha|)$ so that one
can use the relation Gradshtein

$\int_{0}^{+\infty}{\rm d}x\,x\,{\rm e}^{-\alpha x^{2}/2}\,L_{n}(\beta
x^{2}/2)\,J_{0}(xy)=\frac{(\alpha-\beta)^{n}}{x^{n+1}}\,{\rm
e}^{y^{2}/(2\alpha)}\,L_{n}\left(\frac{\beta
y^{2}}{2\alpha(\beta-\alpha)}\right)$

and get

$\rho_{T,\alpha}(n)=\left(\frac{n_{T}}{1+n_{T}}\right)^{n}\,\frac{\displaystyle{\rm
e}^{-\frac{|\alpha|^{2}}{1+n_{T}}}}{1+n_{T}}\,L_{n}\left(-\frac{|\alpha|^{2}}{n_{T}(1+n_{T})}\right)\
.$

## References

  * (1) RJ Glauber, The Quantum Theory of Optical Coherence, Phys. Rev. 130, 2529 (1963).
  * (2) CC Gerry, PL Knight, Introductory quantum optics, (Cambridge University Press, 2005).
  * (3) SJ van Enk, CA Fuchs, The Quantum State of a Propagating Laser Field, Phys. Rev. Lett. 88, 027902 (2001).
  * (4) G Brida et al., Constrained MaxLik reconstruction of multimode photon distributions, J. Mod. Opt. 56, 196-200 (2009).
  * (5) IS Gradshteyn, IM Ryzhik, Table of Integrals, Series and Products, (Academic Press, 2015).
  * (6) J Howe et al., Comparison of film detectors, charged-coupled devices, and imaging plates in x-ray spectroscopy of hot dense plasma, Rev. Scie. Inst. 77, 036105 (2006).
  * (7) RA Giampaolo et al., Depleted MAPS on a 110 nm CMOS CIS Technology, in 26th International Conference on Electronics, Circuits and Systems (ICECS), IEEE (2019).
  * (8) T Hirono et al., Characterization of fully depleted CMOS active pixel sensors on high resistivity substrates for use in a high radiation environment, in Nuclear Science Symposium, Medical Imaging Conference and Room-Temperature Semiconductor Detector Workshop (NSS/MIC/RTSD), IEEE (2016).
  * (9) M Ruat, C Ponchut, Characterization of a Pixelated CdTe X-Ray Detector Using the Timepix Photon-Counting Readout Chip, IEEE Trans. Nucl. Scie. 59, 2392 (2012).
  * (10) M Ghioni et al., Progress in Silicon Single-Photon Avalanche Diodes, IEEE J. Sel. Top. Quant. Electr. 13, 852 (2007).
  * (11) LHC Braga et al., A Fully Digital 8 $\,\times\,$16 SiPM Array for PET Applications With Per-Pixel TDCs and Real-Time Energy Output, IEEE J. Solid-State Circuits 1, 301 (2014).
  * (12) E. Conca et al., Large-Area, Fast-Gated Digital SiPM With Integrated TDC for Portable and Wearable Time-Domain NIRS, IEEE J. Solid-State Circuits 55, 3097 (2020).