Rémi Bardenet
email: remi.bardenet@univ-lille.fr

Alexandre Feller
email: alexandre.feller@univ-lille.fr

Jérémie Bouttier

Pascal Degiovanni

Adrien Hardy

Adam Rançon

Benjamin Roussel

Grégory Schehr

Christoph I Westbrook

From point processes to quantum optics and back

this is Part I of a cross-disciplinary survey that is work in progress; All comments are welcome.] Some fifty years ago, in her seminal PhD thesis, Odile Macchi introduced permanental and determinantal point processes. Her initial motivation was to provide models for the set of detection times in fundamental bosonic or fermionic optical experiments, respectively. After two rather quiet decades, these point processes have quickly become standard examples of point processes with nontrivial, yet tractable, correlation structures. In particular, determinantal point processes have been since the 1990s a technical workhorse in random matrix theory and combinatorics, and a standard model for repulsive point patterns in machine learning and spatial statistics since the 2010s. Meanwhile, our ability to experimentally probe the correlations between detection events in bosonic and fermionic optics has progressed tremendously. In Part I of this survey, we provide a modern introduction to the concepts in Macchi's thesis and their physical motivation, under the combined eye of mathematicians, physicists, and signal processers. Our objective is to provide a shared basis of knowledge for later cross-disciplinary work on point processes in quantum optics, and reconnect with the physical roots of permanental and determinantal point processes.

Introduction

The photoelectric effect is the release of individual electrons from a metal, when light falls onto that metal. The empirical observation that no electrons are released if the light frequency is beneath a certain threshold, along with an analogy with elastic collisions, prompted Einstein to posit the existence of light quanta, also known as photons. A century later, the full explanation of the photoelectric effect is considered a success of the quantum theory of light. In parallel, experimental devices have been designed to amplify the current resulting from a small number of released electrons, eventually giving detectors of light so sensitive that they are able to detect single photons. These detectors have led to investigations on the quantum coherence properties of light. Coherence here means the ability of sources of light to generate interference patterns, such as in Young's celebrated double slit experiment [Mandel and Wolf, 1965, Sections 1-3]. One puzzling aspect of coherence was demonstrated by Hanbury [START_REF] Brown | Interferometry of the intensity fluctuations in light. II. an experimental test of the theory for partially coherent light[END_REF], and is called photon bunching or the HBT effect, after the initials of its discoverers. When placing a detector of single photons in the electric field created by a thermal light source, such as incandescent matter, Hanbury [START_REF] Brown | Interferometry of the intensity fluctuations in light. II. an experimental test of the theory for partially coherent light[END_REF] showed that the detector tends to produce clicks that are grouped in time. Everything happens as if the photons were tightly bunched together when arriving at the detector, hence the name of photon bunching. A few years later, with the invention of the laser, it was realized that this bunching effect disappeared for (so-called coherent) laser light.

At the turn of the 60s, there was intense research in finding the right mathematical framework to describe such coherence phenomena. A celebrated contribution is that of [START_REF] Glauber | Coherent and incoherent states of the radiation field[END_REF], who introduced concepts like coherent states and coherence functions, which quickly became textbook material for quantum optics. Having met Glauber at Les Houches in 1964, French signal processing pioneer Bernard Picinbono launched a team in statistical optics in Orsay, determined to explore the transition between bunching and non-bunching measurements using the formalism of stochastic processes. Picinbono assembled a small and diverse group, including PhD student Odile Macchi, who had just obtained a degree in mathematics. The starting point of Macchi's 1972 PhD thesis was to find the right stochastic object to describe the detection times in HBT and explain photon bunching. Her thesis turned out to be foundational in many respects. This justified a recent translation of Macchi's thesis (originally in French) by Hans Zessin [START_REF] Macchi | Point processes and coincidences -Contributions to the theory, with applications to statistical optics and optical communication, augmented with a scholion by[END_REF], along with a large appendix written by Hans Zessin and Suren Poghosyan.

Odile Macchi introduced what we now call correlation functions to describe point processes, i.e., random configurations of points in a generic space. She showed how common models for physical sources and detection led to point processes with closed-form correlation functions. In particular, for a given model of the source used in the HBT experiment, she showed that the resulting point process of detection times is a permanental point process, for which bunching can be fully characterized and related to properties of the electric field. Teaming up with fellow PhD student in physics Christine Bénard, they used Glauber's formalism to identify the point process describing the parallel situation of detection times of electrons. Unlike for photons, the latter point process naturally exhibits anti-bunching, with detection times being very regularly spaced. [START_REF] Bénard | Detection and "emission" processes of quantum particles in a "chaotic" state[END_REF] is one the foundational stones1 for what we call today determinantal point processes (DPPs).

Outside physics, determinantal point processes (DPPs) have since then become a cornerstone of the theory of random matrices [START_REF] Anderson | An introduction to random matrices[END_REF]Zeitouni, 2010, Johansson, 2006], with applications in combinatorics [START_REF] Borodin | Asymptotics of Plancherel measures for symmetric groups[END_REF] and number theory [START_REF] Rudnick | Zeros of principal L-functions and random matrix theory[END_REF]. DPPs are also a popular model for repulsive point pattern data in spatial statistics [START_REF] Lavancier | Determinantal point process models and statistical inference: Extended version[END_REF] and machine learning [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF]. Many DPP users in these fields 2 have little idea of the physical origin of DPPs. Given the late blooming of DPPs outside their original field, and some fifty years after Macchi's first generic formalization of DPPs, we felt it a natural endeavour to reconnect DPPs with their physical roots, as models for non-interacting fermions and tools to probe generalizations of the HBT effect. This has been the purpose of two workshops so far.

We started with a two-day event3 in 2019 in Lille, France, opened by Odile Macchi and featuring both theoretical and experimental physicists exposing their view of fermionic coherence to an audience of DPP users across mathematics, computer science, and signal processing. The workshop having been met with cross-disciplinary enthusiasm, we made two decisions. The first was to organize an ambitious, two-week follow-up to the workshop,4 which took place in 2022 in Lyon, France. The second decision was to unite forces and write a survey on the links between point processes and optics, to pin down a common ground for discussions. The current document is the first part of this survey. It is both a joint introduction to point processes and quantum optics, and organized notes from a modern reading of Odile Macchi's thesis. The novelty of our document relies in its crossdisciplinary target audience of mathematicians, physicists, and signal processers, with a solid undergraduate background in probability and functional analysis. In particular, while tackling topics in modern physics, we assume little physics knowledge from the reader beyond undergraduate exposure to wave optics. One of our leitmotivs is to sketch the thought process behind some fundamental arguments in quantum optics. Indeed, in our experience, arguments thought as basic by physicists can be hard to grasp by people trained in mathematics, mostly because the implicitly assumed lore differs across communities. In a reverse movement towards physicists, and following the spirit of [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF], we motivate most mathematical concepts by their need as models in optics, including point processes.

We have striven to maintain a balance between mathematical rigour, clarity, and brevity, giving references whenever we had to take shortcuts. We expect that every reader will find some of the material basic and some more exotic, depending on their background, and we hope that all readers will eventually learn something useful. Our objective is to make the potential barrier for crossing from one discipline to another as low as possible, so that ideas can flow more easily. A second part of the document is in preparation, presenting selected advanced topics from the Lille workshop, including experimental measurements of the HBT effect, non-interacting trapped fermions in statistical physics [START_REF] Dean | Noninteracting fermions at finite temperature in a d-dimensional trap: Universal correlations[END_REF][START_REF] Dean | Noninteracting fermions in a trap and random matrix theory[END_REF], fermions in combinatorics, and electronic quantum signals. The two parts are ultimately to be bound in a single manuscript.

The rest of the document is organized as follows. In Section 2, we introduce key examples of point processes. Poisson, Cox, and permanental point processes are motivated there by the so-called semi-classical derivation of the HBT photon bunching effect, treating only the detector as a quantum object, not the electric field. Determinantal point processes are also introduced, but their physical motivation requires to go beyond the semi-classical picture, which justifies the next three sections. Section 3 is a crash course in quantum field theory, from the basics of quantum mechanics to Wick's theorem on the average of products of ladder operators. Wick's theorem yields two very different results depending on the commutation rules of the operators it applies to, which in turn derive from modeling either bosonic particles (like photons), or fermionic particles (like electrons). Ultimately, this dichotomy is at the origin of the appearance of permanental and determinantal point processes. The section concludes with a discussion of the coherent states of [START_REF] Glauber | Coherent and incoherent states of the radiation field[END_REF] and their relation to time-frequency analysis in signal processing. Section 4 covers the modern view on photodetection, culminating in the full quantum justification of the HBT effect using permanental point processes as a model for the coincidence measurements of photons, as well as considerations on the role of the source in the bunching properties of the measurements. Following the lines of Section 4, Section 5 covers the detection of electrons, finally resulting in the appearance of determinantal point processes. The section concludes on the comparative difficulties of recovering non-quantum arguments from the quantum treatment in the case of fermions. Finally, Section 6 wraps up this first part, commenting on the generic derivation of a permanental or determinantal point process from a model of free fermions, and discussing open questions motivated by this construction.

A note on the style. Because of its cross-disciplinary objectives, the style of our document is hybrid. We mostly follow a style inspired by texts in mathematics, with definitions, theorems, assumptions, and examples, sometimes merged with the main text, sometimes fleshed out to draw the reader's attention. Assumptions, in particular, often stand out. By assumption we mean a statement for the reader to accept in order to progress in a discussion or a computation. It can be, e.g., a modeling choice or a mathematical approximation.

Examples are often borrowed from physics; to make the text self-contained, some of these examples are relatively long. To help the reader to visually isolate examples from the rest of the text, we conclude each example with a symbol. Finally, footnotes abound, and are usually meant as side remarks to one of the targeted scientific communities, e.g. to discuss notational or minor conceptual differences between different domains.

Point processes and the semiclassical picture of HBT

In a physical experiment where a detector clicks when hit by a particle, the observation consists of a set of real numbers, the times at which the detector clicks. The natural probabilistic model for such a situation, where both the number and the location of the observed points are uncertain, is a point process, i.e, a random configuration of points. In Sections 2.1 to 2.4, we introduce some of the modern vocabulary of point process theory, along with three key families of point processes: Poisson, permanental and determinantal point processes.

While optical models that feature determinantal point processes will have to wait for the fully quantum treatment of electron detection in Section 5, we already motivate Poisson and permanental point processes in this section by the so-called semi-classical (as opposed to fully quantum) derivation of the photon bunching effect. Our running examples are based on the simple setup given here as Example 0.

Example 0 (A simple photodetection setup). As introduced in Section 1, photodetection is based on the photoelectric effect: when light falls onto a metal, electrons are released, generating a current that we can measure. Light is an electromagnetic field, so typically has an electric and magnetic component, but the physics of the photoelectric effect is essentially dependent on the electric component of the field.

The electric field at a point r ∈ R 3 in space and time t is modeled by a square-integrable function E : R 3 × R → R 3 of space and time. The three components of E(r, t) correspond to what is called in physics the polarization of light. Throughout this paper, for simplicity, we assume that the field is linearly polarized. In this section,5 linear polarization amounts to assuming the existence of a unit vector u ∈ R 3 such that, for all (r, t) ∈ R 3 × R, E(r, t) = E(r, t)u, with E(r, t) ∈ R. This simplification allows us to focus here on the scalar function E : R 3 × R → R. Moreover, we consider a photodetector placed at a fixed position r ∈ R 3 . The effect of the field on the detector is assumed to only depend on the value E(r, •) of the field at the detector position, and we thus further focus on the function t → E (r, t) in our examples, which we also denote6 by E in this section.

We are interested in the times at which our detector clicks, i.e., detects a single photon. Because the measured times vary from one run of the experiment to the next, we want to model them as a random set of (distinct) real numbers.

Finally, note that we only describe in this section an idealized version of the HBT experiment, following [Macchi, 1975, Section 4.2]. For physics arguments, we refer the reader with little prior exposition to physics to [Mandel and Wolf, 1995, Chapter 9], which is a recent textbook treatment of the survey [START_REF] Mandel | Coherence properties of optical fields[END_REF] to which [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF] originally referred.

The correlation functions of a point process

Let us fix a complete metric space X, with µ a Borel measure on the Borel sets of X. A majority of the examples in this paper will deal with X = R equipped with the Lebesgue measure. As mentioned in the header of Section 2, this choice describes the ideal detection times of physical particles. In spatial statistics, we typically encounter X = R d with µ the Lebesgue measure [START_REF] Baddeley | Spatial point patterns: methodology and applications with R[END_REF]. It is also not uncommon to take X to be a compact manifold M (like a sphere; [START_REF] Beltrán | Energy of the Coulomb gas on the sphere at low temperature[END_REF]) or a discrete space (like a large dataset in machine learning [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF]), in which case µ is typically taken to be the volume form on M or the counting measure of the discrete space, respectively.

A point process γ on X is a random configuration7 of points in X. In other words, γ is a random variable taking its values in the space of locally finite subsets of X,

Conf(X) := γ ⊂ X : #(γ ∩ A) < ∞ for all compact A ⊂ X .
In terms of modeling beams of physical particles, the locally finite assumption entails that for any given time interval, there is a finite number of particles that were detected in that interval.

For any k ≥ 1, the k-point correlation function

ρ k : X k → [0, ∞] satisfies, when it exists, E     x1,...,x k ∈γ xi =xj if i =j f (x 1 , . . . , x k )     = X k f (x 1 , . . . , x k ) ρ k (x 1 , . . . , x k ) k i=1 dµ(x i ) (1)
for any bounded (or positive) and measurable test function f . Here

X k = k X × • • • × X
and the symbol E in (1) stands for the expectation8 under the law of the random variable γ. Thus, the k-point correlation function ρ k encodes the distribution of k-tuples of points from γ. Indeed, an informal rewriting of (1) reads

ρ (k) (x 1 , . . . , x k )µ(dx 1 ) . . . µ(dx k ) = P   There are at least k points in γ, one in each of the infinitesimal balls B(x i , dx i ) for i = 1, . . . , k   . ( 2
)
It is common to actually define a point process γ by a sequence (ρ k ) of compatible correlation functions; see [START_REF] Daley | An introduction to the theory of point processes[END_REF]. By compatible, we mean that not every sequence (ρ k ) actually defines a point process, and that usually a mathematical argument for existence is necessary.

Of particular practical significance are the correlation functions for k = 1 and k = 2. The first correlation function ρ 1 describes the marginal distribution of particles, and is called the intensity in probability and statistics. 9 The second correlation function ρ 2 describes pairwise correlations, and is often discussed in its normalized form

g(x, y) = ρ 2 (x, y) ρ 1 (x)ρ 1 (y) . (3)
Finally, when X = R d and the distribution of γ is translation-and rotation-invariant, g only depends on r = x -y , so that we write g(x, y) = g 0 (r). The function g 0 is called the pair correlation function of γ.

Figure 1 shows samples of three translation-invariant point processes in R, along with their pair correlation functions. The three point processes share the same (constant) firstorder correlation function. Note how the pair correlation function in Figure 1(f) is lower than 1 close to 0, indicating fewer small pairwise distances than the reference point process in Figure 1(d): this is a sign of a very regular, more grid-like distribution of the points, as seen in Figure 1(c). On the contrary, the pair correlation function of Figure 1(e) shows more small pairwise distances than the reference: points are lumped together, as confirmed by Figure 1 (b). The rest of this section describes these three examples in more detail.

Poisson and Cox point processes

Let λ : X → R + be locally integrable, that is, B λ(x)dµ(x) < ∞ for every bounded B ⊂ X. Further assume for simplicity that the measure λdµ has no atom. 10 The point process with correlation functions

ρ k (x 1 , . . . , x k ) = λ(x 1 ) . . . λ(x k ), k ≥ 1, (4)
always exists and is called the Poisson point process with parameter function11 λ. The first correlation function of the process is thus ρ 1 = λ. Moreover, the separable form of ( 4) implies the lowest level of correlation among the points of the process. In particular, the pairwise correlation function g ≡ 1 is constant, so that no pairwise distances are preferred.

E + = 2F -1 F(E) × 1 (0,∞) . ( 5
)
Taking the analytic signal [Mandel and Wolf, 1995, Section 3.1] is a partial isometry of L 2 (R) commonly used in signal processing, which has several natural properties. For instance, E + removes the symmetry in |FE|, so that each frequency is represented in an unambiguous, non-redundant manner. To a mathematician, the analytic signal is the boundary value of a particular analytic function of the upper-half plane, defined as the Cauchy transform of f ; see e.g. [Pugh, 1982, Section 2.1]. A more physical property is that when E is quasi-monochromatic, |E + | is a good approximation of the envelope of the signal E; see Figure 2. This idea of an envelope, insensitive to rapid oscillations, along with intuition from classical electromagnetism relating a physical intensity to the square of the amplitude of a wave, is at the origin of the following modeling assumption, which we flesh our for future reference. , where E + is given by (5).

The full justification for the assumption can be obtained by a more precise model of light and matter interaction using quantum theory; see [Mandel and Wolf, 1995, Chapter 9] or our Section 4. Meanwhile, we see that Poisson processes naturally result from deterministic fields. Poisson point processes are useful both as a reference point process and as a building block in statistical modeling. For instance, a Poisson point process with a random parameter function λ is known as a Cox process [Daley and Vere-Jones, 2003, Section 6.2]. A Cox process has correlation functions

ρ k (x 1 , . . . , x k ) = Eλ(x 1 ) . . . λ(x k ), k ≥ 1,
where the expectation is over λ. In particular, by Jensen's inequality,

ρ 2 (x, x) = Eλ(x) 2 ≥ (Eλ(x)) 2 = ρ 1 (x) 2 . ( 6
)
If inequality is strict and ρ 2 is continuous, one can thus expect that for x, y close to each other, the probability that there is a point near x and a point near y in the same realization is larger for a Cox process than for a Poisson process with the same first correlation function.

In other words, samples from a Cox process exhibit clusters, or bunching.

Example 2 (Random classical sources imply photon bunching). When the physical field E comes from a thermal source (i.e., incandescent matter or a gas discharge), experimentally observed detection events do not seem to form a Poisson process, but rather exhibit some form of bunching, or clustering. This is the HBT effect [START_REF] Brown | Interferometry of the intensity fluctuations in light. II. an experimental test of the theory for partially coherent light[END_REF].

Mathematically, it can be seen as a consequence of Assumption 1, as soon as one represents the field resulting from the thermal source by a stochastic process making x, y → ρ 2 (x, y) smooth near the diagonal x = y. The point process of detection times is then a Cox process that favors clusters of points, compared to a Poisson process.

The archetypal Cox processes are the so-called permanental point processes, which we now introduce.

Permanental point processes

A point process γ is said to be permanental when there exists a so-called correlation kernel K : X × X → C such that the correlation functions of γ read

ρ k (x 1 , . . . , x k ) = per K(x i , x j ) k i,j=1 , k ≥ 1 . ( 7
)
The permanent perA in (7) of a matrix A ∈ R k×k is defined, analogously to the determinant, by a sum over permutations perA =

σ∈S k k i=1 a iσ(i) ,
where S k denotes the symmetric group. Assuming existence for a moment, the point process described by ( 7) has nontrivial correlation functions. Taking k = 1 and k = 2 in (1), we obtain for instance

E x∈γ f (x) = f (x)K(x, x)µ(dx), ( 8
) E     x,y∈γ x =y f (x, y)     = f (x, y) K(x, x)K(y, y) + K(x, y)K(y, x) µ(dx)µ(dy) . (9)
When the kernel K is Hermitian, that is, K(x, y) = K(y, x), the second correlation function in (9), once normalized, becomes

g(x, y) = ρ 2 (x, y) ρ 1 (x)ρ 1 (y) = 1 + |K(x, y)| 2 K(x, x)K(y, y) ≥ 1. ( 10
)
In particular, (10) shows that permanental point processes are attractive: the larger |K(x, y)| 2 , the more likely two particles at x and y are to co-occur. Moreover, they are always more likely to co-occur than if γ were a Poisson process with the same first correlation function.

The existence of a point process satisfying (7) requires conditions on K [START_REF] Shirai | Random point fields associated with certain Fredholm determinants I: fermion, Poisson and boson point processes[END_REF]. A standard set of conditions comes from a representation of γ as a Cox process. This representation derives from the semiclassical treatment of the HBT effect for photons by [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF], which we present below as Example 3. To isolate the mathematical statement, let us simply mention that a Poisson point process with random first correlation function λ taken to be the squared modulus of a Gaussian process is a permanental point process; see e.g. [Hough et al., 2006, Proposition 35]. Making sure that the underlying Gaussian process exists in turn guarantees existence of the attached permanental point process.

Example 3 (Gaussian classical fields yield permanental point processes). Continuing Example 2 on thermal sources, quasimonochromatic thermal sources are actually represented by zero-mean stationary Gaussian processes, following the intuition that they result from the superposition of many zero-mean, independent random contributions from the source at roughly the same frequency [Macchi, 1975, Section 4.2]. Mathematically, a random function f is said to have for distribution a zero-mean Gaussian process with kernel K if, for any number n of observations and for any times t 1 , . . . , t n ∈ R, we have

(f (t 1 ), . . . , f (t n )) T ∼ N (0, . . . , 0) T , ((K(t i , t j )) 1≤i,j≤n , (11)
where N (µ, Σ) stands for the multivariate Gaussian with mean µ and covariance matrix Σ. Assuming E is such a zero-mean Gaussian process, and that its distribution is invariant under translations along the time axis, the linear transform E + in ( 5) is also a zero-mean, translation-invariant Gaussian process,12 with the specific property that EE + (t)E + (s) = 0.

As an example, the reader can think of E as having kernel

K Lorentz (t, s) = exp - |t -s| σ cos(ω(t -s)) = K Lorentz,0 (t -s), (12)
for some ω 1/σ. K Lorentz is the product of a slowly varying envelope, which we arbitrarily take to be exponential for concreteness, and a fast oscillating function. In particular, the Fourier spectrum of K Lorentz,0 is concentrated around ω, making E quasi-monochromatic in the sense of [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF]. In that case, E + has kernel

C Lorentz (t, s) = 2K + Lorentz,0 (t -s) ≈ 2 exp - |t -s| σ e iω(t-s) . ( 13
)
The first equality is a general property of analytic transforms of second-order stationary processes, and can be proved by direct computation. The approximation in ( 13) is a consequence of Bedrosian's theorem, stating that the slowly varying envelope is preserved when taking the analytic signal; see e.g. [START_REF] Picinbono | On instantaneous amplitude and phase of signals[END_REF]. Finally, we note that the covariance ( 13) is a slight modification of the Lorentz kernel example given by [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF], where we have introduced a phase factor since, stricto sensu, analytic covariance kernels cannot take only real values. A sample E with kernel K Lorentz is shown in Figure 2. The modulus of the analytic signal E + visibly plays the role of an envelope for E, averaging over local oscillations of E.

We now go back to a generic E following a zero-mean Gaussian process with kernel K. Following Assumption 1, the detection events now form a Poisson process with random parameter function, a.k.a. a Cox process; see Section 2.2. Its correlation functions read

ρ k (t 1 , . . . , t k ) = E|E + (t 1 )| 2 . . . |E + (t k )| 2 = EE + (t 1 ) . . . E + (t k )E + (t 1 ) . . . E + (t k ).
Now the expectation of a product of Gaussians can be expressed in terms of pairwise expectations using a theorem by Isserlis. 13 Further recalling that EE + (t)E + (s) = 0 for all t = s, we obtain

ρ k (t 1 , . . . , t k ) = σ∈S k k i=1 EE + (t i )E + (t σ(i) ). ( 14
)
Recognizing a permanent in ( 14), we conclude that

ρ k (t 1 , . . . , t k ) = per((C(t i , t j ))), k ≥ 1,
where C(t, s) = EE + (t)E + (s) is the covariance kernel of the Gaussian process E + . The photon detection process is thus a permanental point process, see Section 2.3, with kernel the covariance kernel of the analytic signal of the electromagnetic field. In particular, for a Hermitian and translation-invariant kernel C(t, s) = C 0 (t -s) = C 0 (s -t), we obtain, for all s, t such that t -s = r,

g(r) = ρ 2 (t, s) ρ(t)ρ(s) = 1 + |C 0 (r)| 2 . ( 15
)
The probability of coincidence of a pair of detection times is thus larger than under a Poisson process with the same first correlation function. The pair correlation function g in Figure 1(e) is actually (15) with the Lorentz kernel in (13). The corresponding samples exhibit bunching when compared to the Poisson samples of Figure 1(a).

Equation ( 15) provides a remarkably simple mathematical derivation of photon bunching in terms of correlation functions, due to [Macchi, 1975, Chapter 4]. It is noteworthy that the representation of permanental point processes as Cox processes, which has become more of a side result for probabilists [Hough et al., 2006, Proposition 35], was actually the physical motivation for the introduction of permanental point processes by [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF].

Determinantal point processes

A point process γ is said to be determinantal when there exists a so-called correlation kernel K : X × X → C such that the correlation functions of γ read, for any k ≥ 1,

ρ k (x 1 , . . . , x k ) = det K(x i , x j ) k i,j=1 . ( 16
)
We write DPP(K, µ) for the point process with correlation functions (184) with respect to the reference measure µ. Note that, by definition of ρ k , the kernel K(x, y) has to be chosen so that the right hand side of ( 184) is nonnegative for any k ≥ 1 and x 1 , . . . , x k ∈ X; like for permanental point processes, not every kernel yields a well-defined DPP.

Assuming existence for a moment, we look at the first and second correlation functions. Taking k = 1 and k = 2 in (1), we obtain for instance

E x∈γ f (x) = f (x)K(x, x)µ(dx), ( 17
) E     x,y∈γ x =y f (x, y)     = f (x, y) K(x, x)K(y, y) -K(x, y)K(y, x) µ(dx)µ(dy). ( 18
)
When the kernel K is Hermitian, that is, K(x, y) = K(y, x), the second correlation function in (18), once normalized, becomes:

g(x, y) = ρ 2 (x, y) ρ 1 (x)ρ 1 (y) = 1 - |K(x, y)| 2 K(x, x)K(y, y) ≤ 1. ( 19
)
In particular, (19) shows that DPPs with Hermitian kernels are repulsive: the larger |K(x, y)| 2 , the less likely the two particles at x and y are to co-occur. Moreover, they are always less likely to co-occur than if γ were a Poisson process with the same first correlation function. In particular, if γ were to model the clicks of a detector as particles arrive, the arrival times would exhibit antibunching. We shall see examples of antibunching arrival times in Section 5. Unlike Poisson and permanental point processes, whose appearance results from Assumption 1, to understand why DPPs appear in quantum optics, we shall need first to introduce elements of quantum field theory. Necessary and sufficient conditions for existence of a DPP are known when the operator K : f → K(•, y)f (y)dµ(y) on L 2 (µ) is Hermitian and locally trace-class. In particular, letting (λ k ) denote the eigenvalues of K, DPP(K, µ) is well-defined if and only if 0 ≤ λ k ≤ 1 for any k ∈ N, namely when K is a contraction operator. This is now known as the Macchi-Soshnikov theorem.

One constructive proof of this existence theorem relies on the following decomposition lemma [START_REF] Hough | Determinantal processes and independence[END_REF]. Informally, a DPP with kernel

K(x, y) = ∞ k=0 λ k ϕ k (x)ϕ k (y), ( 20
)
is a statistical mixture of projection DPPs, i.e., DPPs with projection kernels. More precisely, if (b k ) k∈N are independent Bernoulli random variables with

P(b k = 1) = 1 -P(b k = 0) = λ k and K b (x, y) := k b k ϕ k (x)ϕ k (y), then DPP(K b , µ) = DPP(K, µ)
in law; see again [START_REF] Hough | Determinantal processes and independence[END_REF]. In particular, if we restrict ourselves to projection kernels, i.e., if λ k ∈ {0, 1} for every k, the cardinality of γ ∼ DPP(K, µ) is thus Tr K almost surely. If we further assume Tr K = N < ∞, the corresponding projection DPP generates exactly N particles, with joint probability distribution on X N given by

1 N ! det ϕ k-1 (x j ) N j,k=1 2 N j=1 µ(dx j ). ( 21
)

Elements of quantum field theory

To deepen our understanding of the relationship between point processes and quantum physics, we now introduce a few elements of what is known as many-body quantum physics, quantum statistical physics or quantum field theory. We first introduce the basic quantum formalism to describe a single particle, its evolution in time and what information we can gather on it, mostly following [START_REF] Folland | Quantum Field Theory: A tourist guide for mathematicians[END_REF]. Then, we introduce how to deal with systems made of many or an indefinite number of particles.

The mathematical framework of quantum theory

States. The state ψ of a quantum system is represented by a non-zero element of a complex Hilbert space (H , •|• ), the space of all possible states of the quantum system of interest. Two elements of H that are equal up to multiplication by a (complex) scalar are understood to represent the same physical states,14 and we henceforth always assume that our states are normalized, i.e. ψ = ψ|ψ = 1. Finally, we will often use the so-called bra-ket notation, denoting a vector by |ψ = ψ ∈ H and the associated linear form on H by ψ| : h → ψ|h .

Observables. Any observable quantity is represented by a self-adjoint operator A : H → H . To understand in what sense, take a self-adjoint operator A : H → H and a vector |ψ ∈ H , and define an ordinary probability measure P A,ψ on the spectrum σ

(A) ⊂ R of A by P A,ψ : E → ψ|1 E (A)|ψ , E ⊂ σ(A). ( 22
)
The operator 1 E (A) in ( 22) is defined by the spectral functional calculus; see [Reed and Simon, 1980, Chapter VII]. Informally, 1 E (A) is the projection operator with the same eigenvectors as A, but replacing each eigenvalue λ by 1 if λ ∈ E, and 0 otherwise. One easy setting to understand the definition of P A,ψ is that of a compact (self-adjoint) operator

A, for in that case the support σ(A) = {λ n } of P A,ψ is discrete, with no other accumulation point than 0, and there is an orthonomal basis (|ϕ n ) of H with A |ϕ n = λ n |ϕ n ; see e.g. [Reed and Simon, 1980, Theorem VI.16]. In particular, for λ n = 0 with multiplicity one, P A,ψ ({λ n }) = | ψ|ϕ n | 2 , which indeed sums to ψ 2 = 1 in n. Without entering into details, we finally note that the definition (22) extends to unbounded operators by representing them as multiplication operators [Reed and Simon, 1980, Theorems VIII.4 and VIII.5]. Unbounded operators are necessary to represent physical observables that have an unbounded support, such as the position of a particle on the real line.

Interpreting a self-adjoint operator A as an observable means that the experimentally accessible information about the observable A is in P A,ψ . In particular, if we had a machine that could repeatedly prepare unrelated copies of the physical system in question in a state |ψ , we would model the measurements of observable A on such a sequence of states as independent draws from P A,ψ . Unless |ψ is an eigenvector of A with eigenvalue λ, in which case P A,ψ = δ λ , measuring the physical quantity encoded by A thus becomes inherently probabilistic: one can only talk about the probability that a measured observable will lie in a given E ⊂ R. Finally, we note that by the law of large numbers, the average of a large number of measurements of A is expected to be close to the expectation of a random variable with law P A,ψ , which is

λ dP A,ψ (λ) = ψ|A|ψ . ( 23
)
Again, the derivation of ( 23) is easier for a compact operator A = n λ n |ϕ n ϕ n | , say even Hilbert-Schmidt, so that λ 2 n < ∞. In that case, both sides of ( 23) are equal to i λ n | ψ|ϕ n | 2 < ∞. The general treatment of self-adjoint operators requires to carefully define the spectral measure of an operator, and we refer to [START_REF] Reed | Methods of modern mathematical physics: Functional analysis[END_REF] with p a probability measure on H . The linear combination Equation ( 24) is typically interpreted as being the output of a noisy preparation device, which outputs the state |ψ ψ| with |ψ drawn from p. When the agent runs its machine many times, they will describe the output by Equation (24).15 More abstractly, a density matrix ρ is a traceclass positive operator on H with unit trace. 16 A density matrix is often called a mixed state, as opposed to a pure state, which corresponds to a single element of H . In other words, pure states correspond to density matrices that are also projectors. Generally speaking, a density matrix ρ exhausts all the statistical content that an observer can predict about a system. For instance, for ρ defined by ( 24),

Tr [ρ1 E (A)] = ETr [|ψ s ψ s | 1 E (A)] = E ψ s | 1 E (A) |ψ s , E ⊂ R.
In particular, the map P A,ρ : E → Tr (ρ1 E (A)) is a probability distribution on R that corresponds to a statistical mixture of pure states, and simplifies to (22) for pure states. It gives the probability that a measurement of A will belong to E. By linearity, the expectation of the random variable17 describing the measurement of observable A when the system is in the state ρ, denoted by convention A ρ , is given by Born's rule

A ρ = Tr (ρA) . ( 25
)
This formula is a cornerstone of the quantum formalism. It contains Equation ( 22) as a special case and thus describes both noisy and noiseless situations.

Conditioning on measurements. To define the joint distribution of the measurement of two observables A and B when the physical system is in a (possibly mixed) state ρ, one would like to use the operator AB in a spectral formula like (22). It turns out to be very natural to define a joint distribution using the spectral calculus when A and B commute, i.e. AB = BA. In particular, one can then talk of the conditional distribution of the measurement of B on ρ given that we observed A on ρ. Without entering into details, Bayes' formula yields that evaluating the conditional amounts to evaluating P B, ρ, where

ρ = CρC † , ( 26
)
and C is related to the measurement of A; see [START_REF] Bouten | An introduction to quantum filtering[END_REF] for a precise statement that includes monitoring a state across time, i.e., filtering. The "sandwiched" updated state ρ in ( 26) is interpreted as the state of the system immediately after the measurement of A.

Incompatible observables and Heisenberg's uncertainty principle. When AB = BA, it is not even guaranteed that AB is self-adjoint, let alone that A and B have common eigenspaces. Physicists associate this mathematical difficulty to the fact that noncommuting observables are incompatible: it is not possible to obtain a joint measurement of both observables using a single prepared copy of a state ρ.

Another hint that non-commuting observables are peculiar is Heisenberg's celebrated uncertainty principle; see [Folland, 2008, Section 3.3]. In a nutshell, for two non-commuting observables A and B and a state ρ, the product of the standard deviations σ A,ρ of P A,ρ and σ B,ρ of P B,ρ is lower-bounded,

σ A,ρ σ B,ρ ≥ 1 2 |Tr [(AB -BA)ρ]| . ( 27
)
Since A and B cannot be measured simultaneously, this does not imply anything on measuring A and B on the same copy of ρ, but it is rather a property of the model of the overall experiment, i.e., of the physical system and the probes corresponding to A and B. The concrete consequences of the uncertainty relation ( 27) can be seen by repeatedly preparing a copy of the state ρ and measuring either A or B. If the measure represented by A yields measurements with a small empirical variance, then measuring B on similar repeated copies of ρ will lead to a comparatively large empirical variance.

Schrödinger's equation and Hamiltonians. To complete our description of a physical system, we need to model its evolution in time. In quantum physics, the evolution of the state of the system is given by a one-parameter group {U (t), t ∈ R} of unitary operators on H , that is,

U (t)U (s) = U (t + s) and U (t) -1 = U (t) † = U (-t).
In particular U (0) = 1 is the identity operator. After time t, the state of a system that was in state |ψ(0) at time t = 0 is considered to be

|ψ(t) = U (t) |ψ(0) . ( 28
)
Under weak assumptions on H and U (t), one can show that U (t) = expit H in the sense of the spectral functional calculus again, for some (possibly unbounded) self-adjoint operator H on H [Folland, 2008, Section 3.1]. Conversely, any choice of self-adjoint H gives a one-parameter unitary group. The operator H, thought as an observable, is called the Hamiltonian of the system. As the generator of the dynamics, H is often thought as the energy of the system. In particular, the group identity of U yields by differentiation that U (t) † HU (t) = H for all t, so that the average energy ψ(t)|H|ψ(t) of the state ψ(t) in ( 28) is constant over time. Building a model in quantum physics usually boils down to choosing a state space H , a set of self-adjoint operators as observables, and a Hamiltonian to describe the evolution of the system in time. In particular, when s 1, (28) informally yields

|ψ(t + s) = 1 - isH + O(s 2 ) |ψ(t) , so that i d dt |ψ(t) = H |ψ(t) , ( 29
)
which the reader may recognize as the celebrated Schrödinger equation. This equation of motion can be generalized to mixed states (24) by linearity.

To finish, note that instead of letting |ψ(t) depend on time through (28) and keeping observables constant (the so-called Schrödinger picture), we can obtain the same measurement probabilities ( 22) by keeping states constant and letting observables vary as A(t) = U (t) † A(0)U (t) (the so-called Heisenberg picture). The two pictures are thus equivalent as to what observations they predict. There is a third equivalent convention that is commonly used, the interaction picture. The latter is tailored to problems where the Hamiltonian H = H 0 + H I is the sum of a well-studied Hamiltonian H 0 , typically describing standard prior information on a subpart of the physical system (say, a free field), and an additional term H I , called interaction Hamiltonian, typically encoding the interaction between parts of the system. In the interaction picture, one lets both states and observables evolve in time. Observables evolve according to A(t) = U 0 (t) † A(0)U 0 (t), using the one-parameter unitary group U 0 associated to Hamiltonian H 0 . On the other hand, states evolve as |ψ(t) = V (t) |ψ(0) , where V (t) = U 0 (t) † U (t) and U (t) is the group associated to the full Hamiltonian H. While the state evolution now differs from (29), the interaction picture yields again the same probabilities ( 22). Unless otherwise specified, we use the interaction picture in this work.

Perturbation theory.

What is convenient with the interaction picture is that it singles out the role of H I as follows. To compute U (t), it is enough to compute V and then apply U (t) = U 0 (t)V (t), since U 0 (t) is assumed to be known. Now, by definition, V satisfies

d dt V (t) = i U 0 (t) † H 0 U (t) - i U 0 (t) † HU (t) = 1 i H I (t)V (t). ( 30
)
Even in seemingly simple physical situations, as when modeling the detection of photons in Section 4.1, solutions to (30) have to be approximated. Writing (30) in integral form, and noting that V (0) = 1, we obtain

V (t) = V (0) + t 0 d dt V (τ )dτ = 1 + 1 i t 0 H I (τ )V (τ )dτ . ( 31
)
Plugging (31) into itself yields

V (t) = 1 + 1 i t 0 H I (τ )dτ - 1 2 t 0 dτ τ 0 dτ H I (τ )H I (τ )V (τ ). ( 32
)
Iteratively plugging (31) into (32) yields a series representation for V , the study of which is called perturbation theory, and is at the heart of quantum field theory and its current mathematical difficulties [Folland, 2008, Chapter 6]. In practice, physicists will often assume that stopping this iterative process early on yields a good approximation to V . For instance, we shall see in Section 4.1 a physical assumption that is mathematically interpreted as keeping only the first two terms in the right-hand side of (32).

Two fundamental one-particle systems

We saw in Section 3.1 that a quantum model consists in a Hilbert space H , a Hamiltonian H specifying the dynamics, and self-adjoint operators that describe observables. To make things more concrete, we now describe two simple systems: the qubit and the harmonic oscillator.

The qubit

Consider a physical system whose state is described by one of two labels, say {e, g}, like whether a two-level atom is excited or in its ground state. One would then use H = C 2 with the usual inner product. Let now (|e , |g ) be any orthonormal basis of H , and |ψ = α |e + β |g be a state. Normalization implies |α| 2 + |β| 2 = 1.

Observables are described by self-adjoint operators of H . For instance, the observable corresponding to checking whether the system is in state |e is the projector |e e|. Following (22), the probability to obtain the result e in the state |ψ is given by

P |e e|,ψ ({e}) = e| (|ψ ψ|) |e = | e|ψ | 2 = |α| 2 = 1 -P |e e|,ψ ({g}) = 1 -|β| 2 .
More generally, every self-adjoint operator on C 2 can be expressed as a linear combination of four matrices. Expressed in the basis (|e , |g ), the four so-called Pauli matrices are

1 = 1 0 0 1 , σ x = 0 1 1 0 , σ y = 0 -i i 0 , σ z = 1 0 0 -1 . ( 33
)
All four matrices squared are the identity, so they are all diagonalizable with spectrum {±1}. For instance, an orthonormal basis of eigenvectors of σ x is given by

|± x = 2 -1/2 (|e ± |g ),
which is simply a rotation of the basis (|e , |g ). We can compute for instance the probability ( 22) of observing the outcome + x , when the state is a generic |ψ = α |e + β |g , namely

P |+x +x|,ψ ({+ x }) = ψ| |+ x + x | |ψ = | ψ|+ x | 2 = |α + β| 2 2 ≤ 1.
Similarly, the average value of the observable σ x is given by

ψ| σ x |ψ = |α + β| 2 2 - |α -β| 2 2 = 2R(αβ * ).
Finally, we give the example of a rather common and simple Hamiltonian for the qubit, called free evolution. In the basis we have chosen, it reads

H FQ = 1 2 ω eg σ z , ( 34
)
where ω eg is the transition frequency between the excited (|e ) and ground (|g ) states, and the label FQ stands for free qubit. The resulting evolution in time is easy to understand. The states |g and |e are eigenstates of H FQ . Calling U 0 the one-parameter group corresponding to H FQ , and starting e.g. at state |g at t = 0, the probability of finding the evolved state U 0 (t) |g in state |g is g|U 0 (t)|g = 1. The basis vectors |e and |g are thus stationary states. Finally, note that the form (34) of the free dynamics is mostly conventional; we could for instance add a term proportional to the identity to shift the energy reference without changing the physics. What matters is that the energy difference between the two states is ω eg . Richer physical behaviors are modeled by adding interaction terms to the free Hamiltonian. There is no general prescription as to how to write interaction Hamiltonians and they mostly depend on the physical context we wish to describe. One simple but useful model is to consider a qubit interacting with an external field B(t). Note that we mean here a classical field, in the sense that B(t) is a vector (B x (t), B y (t), B z (t)) T of functions from R to R. Concretely, the model describes a small magnet or the spin of a particle in an external magnetic field. The basic ideas behind nuclear magnetic resonance, and magnetic resonance imaging, come from this choice of interaction Hamiltonian [Cohen-Tannoudji et al., 2019a, Chapter IV, Appendix F]. Adding the free part, we obtain the total Hamiltonian

H = 1 2 ω eg σ z + gσ • B(t), (35)
where g ∈ R is a coupling constant and σ • B(t) is a shorthand notation for the linear combination of matrices

σ x B x (t) + σ y B y (t) + σ z B z (t).
When the field is time-independent, the dynamical problem can be interpreted geometrically. Indeed, the Hamiltonian is of the general form H = ωn • σ with n a unitary vector. The evolution operator is thus

U (t) = e -iωtn•σ ,
and we recognize one way of parametrizing an element of the rotation group SU(2), where ωt is the angle of the rotation around the axis n. The dynamics of a qubit in this model can be represented as a vector, representing the state |ψ , which precesses around an axis fixed by the external field; see [START_REF] Degiovanni | Physique quantique, information et calcul: des concepts aux applications[END_REF] for details.

The harmonic oscillator

To model the position of a particle living in R d , the state space is commonly taken to be H = L 2 (R d ). A natural observable is the position of the particle. We associate 18 the measurement of the jth coordinate to the multiplication operator

X j f : x → x j f (x), defined for all functions f of H = L 2 (R d ) such that x → x j f (x) ∈ L 2 (R d ).
Note that the only candidates to be common eigenvectors of all X j are delta functions, and are thus not in L 2 . When performing computations, however, it is convenient to consider the tempered

distributions x| := δ x , x ∈ R d
as generalized states, and interpret δ a (x)δ b (x)dx = δ(a -b) and f = f (x)δ x (•)dx as making |x a "generalized orthogonal basis". In particular, for a smooth state |ψ ∈ H , x|ψ = ψ(x), and the probability distribution ( 22) corresponding to X j in a smooth state |ψ becomes

P Xj ,ψ (E j ) = Rו••×R×Ej ×R×...R |ψ(x)| 2 dx, E j ⊂ R .
The reader may recognize here the common interpretation of the squared modulus |ψ(x)| 2 of the wave function as the probability density function for the position of the particle. Similar considerations allow for considering the momentum operator P j , defined for 1 ≤ j ≤ d through its Fourier transform

F(P j f )(k) = k j F(f )(k), for all f ∈ H such that k j F(f )(k) ∈ L 2 (R d ).
The "generalized basis" of tempered distributions k| : f → F(f )(k) allows talking about the momentum representation k|ψ of a smooth state |ψ . Furthermore, by the inverse Fourier transform, we can express the momentum operator P j in the position basis as a derivative operator,

P j = -i ∂ j .
A natural Hamiltonian in this setting is the so-called quantum harmonic oscillator, modelling the movement of a single particle in a quadratic potential. It writes

H = d j=1 1 2m P 2 j + 1 2 mω 2 Q 2 j , ( 36
)
18 Choosing what operator to associate to what classical quantity is a process known as quantization and is a whole research area; see [START_REF] Folland | Quantum Field Theory: A tourist guide for mathematicians[END_REF] and references therein.

where m > 0 is called the mass and ω ∈ R the angular frequency of the oscillator. For ψ a smooth (Schwartz) function, we get

x| H |ψ = - 2 2m ∆ψ(x) + 1 2 mω 2 x 2 2 ψ(x).
In particular, the eigenequation x| H |ψ = λψ(x) leads to a well-known differential equation. The solutions are tensor products of Hermite functions [Folland, 2008, Section 3.4], and form an orthonormal basis of L 2 (R d ) of (smooth) eigenfunctions of H.

The harmonic oscillator can also be diagonalized algebraically. Define the so-called ladder operators

a j := mω 2 Q j + i mω P j and a † j := mω 2 Q j - i mω P j .
They satisfy the commutation relation

a i a † j -a † j a i = δ ij .
The Hamiltonian can now be rewritten as

H = ω d j=1 a † j a j + 1 2 , ( 37
)
and its spectrum is given by that of the number operators N j = a † j a j . Using that the number operator is positive, one can easily show that its eigenvalues are n j ∈ N, with corresponding eigenstates |n j , i.e. N j |n j = n j |n j , called the Fock states [Cohen-Tannoudji et al., 2019a, Chapter V]. The corresponding wavefunctions x|n j are the Hermite functions discussed above.

The ladder operators were introduced first by Dirac to solve the harmonic oscillator of a single quantum particle. But these operators have taken a life of their own, to describe many-particles quantum states, as we discuss now.

Modelling a finite number of particles

Section 3.1 focused on describing systems made of one particle. In this section, we show how to combine such systems to describe one made of a (known) finite number of particles.

Subsystems, bosons and fermions

Let us now consider a system of N < ∞ identical particles. Concatenating several physical systems is represented by tensor products, so that one is tempted to consider the Hilbert space H ⊗N = N i=1 H, where H is the Hilbert space corresponding to a one-particle system. A generic state would thus be a linear combination of states of the form

|ψ = N i=1 |ψ i =: |ψ 1 • • • ψ N . ( 38
)
However, experiments rather suggest to use a state that encapsulates the property that the particles are indistinguishable, i.e., that permuting particles leaves the state invariant.

Arguments from (projective) representation theory of the symmetric group [START_REF] Degiovanni | Physique quantique, information et calcul: des concepts aux applications[END_REF]

∈ S N |ψ 1 , . . . , ψ N = U σ |ψ 1 , . . . , ψ N := ψ σ(1) , . . . , ψ σ(N ) ,
where we have implicitly defined the unitary representation U σ of the permutation group on the Hilbert space H ⊗N . In other words, the Hilbert space for N bosons is the symmetric subspace of H ⊗N , namely the range of the orthogonal projector

S = 1 N ! σ∈S N U σ .
The other type of particles is called fermion, and corresponds to antisymmetric states

|ψ 1 , . . . , ψ N = ε(σ)U σ |ψ 1 , . . . , ψ N , σ ∈ S N ,
where ε(σ) is the signature of σ. Systems of N fermions are thus represented by states in the antisymmetric subspace of H ⊗N , namely in the range of the orthogonal projector

A = 1 N ! σ∈S N ε (σ) U σ .
Note that observables that are permutation-invariant leave both the ranges of S and A invariant. If the Hamiltonian is permutation-invariant, bosons thus stay bosons across time, and fermions stay fermions. One particular case is that of N free (a.k.a. non-interacting) particles. Letting H be the Hamiltonian for one particle, the dynamics of N free particles is given by the Hamiltonian

H (N ) = N n=1 1 ⊗ • • • ⊗ 1 ⊗ H n-th term ⊗1 ⊗ • • • ⊗ 1.
(39)

Example N -particle states

Consider the harmonic oscillator setting of Section 3.2.2, say with d = 1 for simplicity.

Since the Hamiltonian is interpreted as measuring the energy of the system, and since the eigenstates of the Hamiltonian H in (36) form a basis of L 2 (R), the state with minimal energy is the eigenfunction of H with minimal eigenvalue. In position representation, it is the first Hermite function ϕ 0 , that is, a normalized Gaussian function. Now, if we consider N free bosons, each coming with the same 1-dimensional harmonic oscillator Hamiltonian, the ground state of H (N ) in ( 39) is easily seen to be, in position representation,

ψ Sym 0 (x 1 , . . . , x N ) = ϕ 0 (x 1 ) . . . ϕ 0 (x N ),
so that, in this ground state, the probability density of having a configuration of N free bosons at (x 1 , . . . , x N ) ∈ R N is given by

|ψ Sym 0 (x 1 , . . . , x N )| 2 = |ϕ 0 (x 1 )| 2 • • • |ϕ 0 (x N )| 2 .
This is the law of N independent random variables with identical distribution |ϕ 0 (x)| 2 dx, i.e. N i.i.d. Gaussians. Still within the harmonic oscillator setting, but this time for N free fermions, the situation is less straightforward. With similar eigendecomposition arguments, one can show that the ground state of H (N ) on AH ⊗N is given by

ψ Asym 0 (x 1 , . . . , x N ) = 1 √ N ! det    ϕ 0 (x 1 ) • • • ϕ 0 (x N ) . . . . . . ϕ N -1 (x 1 ) • • • ϕ N -1 (x N )    , ( 40
)
where

ϕ k-1 ∈ L 2 (R d
) is the (unit-norm) eigenstate associated with the k th smallest eigenvalue of H, i.e. here, the kth Hermite function. The right hand side of ( 40) is known as a Slater determinant. In particular, in this ground state, the density probability of having a configuration at (x 1 , . . . , x N ) ∈ R N d is given by

|ψ Asym 0 (x 1 , . . . , x N )| 2 = 1 N ! det    ϕ 0 (x 1 ) • • • ϕ 0 (x N ) . . . . . . ϕ N -1 (x 1 ) • • • ϕ N -1 (x N )    2 = 1 N ! det K N (x i , x j ) N i,j=1
where

K N (x, y) := N -1 k=0 ϕ k (x)ϕ k (y).
Since the ϕ k 's are orthonormal, K N is a projection kernel. In other words, the N positions of the free fermions form a projection DPP on X = R, as introduced in Section 2.4. In our setting where d = 1 and the ϕ k are Hermite functions, this DPP is known as the Gaussian unitary ensemble (GUE); see e.g. [START_REF] Dean | Noninteracting fermions in a trap and random matrix theory[END_REF]. The GUE is a fundamental example in random matrix theory, as it arises as the eigenvalues of a (A + A * )/2, where A is filled with i.i.d. complex Gaussian random variables [START_REF] Anderson | An introduction to random matrices[END_REF]. More generally, any projection DPP can be obtained as the position representation of the ground state of an (adhoc and non-necessarily physically realizable) free Hamiltonian for fermions. This is our first encounter of DPPs as arising from a fermionic construction.19

Occupation number representation

Since the Hilbert space for modeling N bosons, respectively N fermions, is SH ⊗N , respectively AH ⊗N , it is useful to have an orthonormal basis of these spaces [Folland, 2008, Section 4.5]. A convenient way of writing states of indistinguishable particles is through the occupation number basis.

The rationale behind this basis is simply to count the number of particles that are in a specific basis state of the one-particle Hilbert space, chosen to make observables simple to express. Let (e i ) be a Hilbert basis of H. Then an orthonormal basis of SH ⊗N is given by

|n 1 , n 2 , . . . := 1 Z S e 1 , • • • , e 1 n1 , e 2 , • • • , e 2 n2 , • • • , where i≥1 n i = N,
with Z the normalization constant. The label of the state |n 1 , n 2 , . . . represents occupation numbers: n i represents the number of particles in the state |e i . By construction, the inner product between two states with different occupation numbers is zero. The normalization constant Z > 0 is fixed by requiring the state to be normalized, i.e.,

1 = 1 Z 2 n 1 , n 2 , . . .|S † S|n 1 , n 2 , . . . = 1 Z 2 1 N ! ∞ i=1 n i ! .
The same construction applies for fermions. An orthonormal basis of AH ⊗N is given by

|n 1 , n 2 , • • • = √ N !A e 1 , • • • , e 1 n1 , e 2 , • • • , e 2 n2 , • • • , where i≥1 n i = N , ( 41
)
with the important difference that for fermions, n i ∈ {0, 1} for all i. Indeed, if there is some j such that n j ≥ 2, then the corresponding state of H ⊗N is in the kernel of the antisymmetrization operator A. Physicists identify this property with Pauli's exclusion principle: no two fermions of the same multi-particle state can be in the same mode, i.e., correspond to the same vector of the basis e i of the single-particle space H.

Modeling an indefinite number of particles

In order to form a state space that can accommodate any number of particles, and in particular to model uncertainty in that number, physicists introduce so-called Fock spaces.

Fock spaces for bosons and fermions

Denote by H ⊗0 ∼ = C the Hilbert space consisting of all multiples of some unit norm state labeled as |0 , called the vacuum, and representing the absence of any particle. Taking the convention Sh = Ah = h for h ∈ H ⊗0 , we define two Hilbert spaces H ∞ Bosons and H ∞ Fermions , respectively called the bosonic and fermionic Fock spaces, as orthogonal direct sums

H ∞ Bosons := k≥0 SH ⊗k , H ∞ Fermions := k≥0 AH ⊗k . ( 42
)
Several remarks are in order. Note first that in (42), by an abuse of notation, we write the symmetrization and antisymmetrization operators in any dimension as S and A, respectively. Second, the orthogonal sum in (42) makes two vectors coming from different summands have inner product zero by definition. Third, we underline that the absence of particles is modeled by the vacuum state |0 , not by the null element of either Hilbert space

H ∞ Bosons or H ∞ Fermions .
In particular, unlike any null element, the vacuum state represents a physical state, and thus has norm 1 in the Fock space.

The Fock spaces in ( 42) are defined as orthogonal sums of Hilbert spaces, and are thus Hilbert spaces themselves. For instance,

H ∞ Bosons is the set of collections h = (h k ) ≥0 , where h k ∈ SH ⊗k and h := k h k 2 < ∞.
It is sometimes convenient to work instead with the algebraic sum of the same vector spaces, meaning that all but a finite set of coefficients are constrained to be zero. This is the case, for instance, when defining the number operator through

N |n 1 , n 2 , . . . := i n i |n 1 , n 2 , . . . , (43)
which counts the particles in a state, and is defined by the same formula for bosons and fermions. N is then uniquely extended to an unbounded operator on the Fock space [Folland, 2008, Section 4.5]. In this text, we will assume that this extension step can always be done unequivocally, and work directly with the (Hilbert) Fock spaces (42). Finally, a state of the form |n 1 , n 2 , . . . with i≥1 n i < ∞ is called a Fock state, and together they form an orthonormal basis of the corresponding Fock space. Any state of a Fock space can thus be represented as a linear combination of Fock states.

Creation and annihilation operators

Besides the number operator, we now introduce two types of operators on Fock spaces that are respectively thought of as creating and destroying a particle in a given one-particle state. The choice of the symmetrization operator S or A has important consequences on the commutation relations between these operators, and we thus separate the treatment of bosons and fermions. We paraphrase here [Folland, 2008, Chapter 4.5].

Bosons.

For v ∈ H , define the operator b(v) on the algebraic sum corresponding to (42), i.e., finite linear combinations of Fock states, by

b(v)S(u 1 ⊗ • • • ⊗ u k ) = 1 k k j=1 v|u j S(u 1 ⊗ • • • ⊗ u j-1 ⊗ u j+1 ⊗ • • • ⊗ u k ),
for any k ≥ 1 and u 1 , . . . , u k ∈ H . Now, for later ease of writing, we renormalize b and define the operator a(v) on the algebraic sum corresponding to (42), i.e., finite linear combinations of Fock states, by a(v)w = √ kb(v)w whenever w ∈ SH ⊗k . Alternately, using the number operator introduced in Section 3.

4.1, a(v) = √ N + I b(v). The operator a(v) is called the annihilation operator in state v. Its adjoint a † (v) := a(v) † on the same algebraic sum is defined as a(v) † = S √ N b(v) †
, and is called the creation operator in state v. Together with a(v), they are called ladder operators. 20 21 Writing [u, v] = uv -vu for the commutator of two operators, we can check that the ladder operators satisfy the so-called canonical commutation relations (CCR),

[a(u), a(v)] = a † (u), a † (v) = 0, a(u), a † (v) = u|v I. ( 44
)
To see why the creation and annihilation operators bear such names, consider again a basis (e i ) of H , and let a i = a(e i ) for i ≥ 1. It can be checked that the ladder operators act on Fock states as

a † i |n 1 , • • • , n i , • • • = √ n i + 1|n 1 , • • • , n i + 1, • • • , ( 45a
)
a i |n 1 , • • • , n i , • • • = √ n i |n 1 , • • • , n i -1, • • • , (45b)
where, by convention, the right-hand side of (45b) is the null element of the Fock space whenever n i = 0. Physically, a † i thus models the creation of a particle in the state |e i , while a i removes a particle from the same state. Note also that a i |0 = 0 is the null element of the Fock space; we say that the annihilation operator annihilates the vacuum.

In terms of ladder operators, the number operator ( 43) can be rewritten

N = ∞ i=1 a † i a i , ( 46
)
with N |0 = 0 |0 = 0, the null element of the Fock space, as expected from an operator that counts particles. Note also that the CCRs imply [N,

a † i ] = a † i and [N, a i ] = -a i .
Finally, any Fock state can be rewritten as the action of creation operators on the vacuum as

|n 1 , • • • , n i , • • • = 1 ∞ i=1 n i ! ∞ i=1 a † i ni |0 . ( 47
)
Note that by definition of Fock states, the products in (47) consist of a finite number of terms.

Fermions. Following the same lines, we can now define the creation and annihilation operators for fermions on the Fock space H ∞ Fermions . Formally, a(v) is defined on finite linear combinations of Fock states 22 by

a(v)A(u 1 ⊗ • • • ⊗ u k ) = 1 √ k k j=1 (-1) j v|u j A(u 1 ⊗ u j-1 ⊗ u j ⊗ • • • ⊗ u k ), while one can check that its adjoint a † (v) := a(v) † satisfies a † (v)A(u 1 ⊗ • • • ⊗ u k ) = √ k + 1Au 1 ⊗ • • • ⊗ u k .
As for bosons, these definitions lead to particular commutation relations. Writing {u, v} = uv + vu for the so-called anti-commutator of two operators, the fermionic annihilation and creation operators satisfy

{a(u), a(v)} = {a † (u), a † (v)} = 0 {a(u), a † (v)} = u|v I, ( 48
)
known in the literature as the canonical anti-commutation relations (CARs). Here again, the presence of an anti-commutator is interpreted by physicists as a manifestation of Pauli's exclusion principle. For instance, the identity {a † (u), a † (u)} = 2a † (u)a † (u) = 0 translates the fact that there cannot be two fermionic particles in the same quantum state.

spaces. In the rest of the paper, we only use the ladder operators acting on the Fock space.

22 and then extended to the whole Fock space.

Again, the ladder operators a i = a(e i ) in a given basis apply straightforwardly to the corresponding Fock states. For i ≥ 1, and remembering that n i ∈ {0, 1}, we obtain

a † i |n 1 , • • • , n i , • • • = δ ni0 ε(σ)|n 1 , • • • , n i + 1, • • • (49a) a i |n 1 , • • • , n i , • • • = δ ni1 ε(σ)|n 1 , • • • , n i -1, • • • , (49b)
where n = n i < ∞, σ is the permutation

σ = 1 • • • S i + 1 S i + 2 • • • n + 1 S i + 1 • • • S i S i + 2 • • • n + 1 , ( 50
)
where

S i = i-1
j=1 n j are the partial sums of the sequence of occupation numbers. These fermionic ladder operators call for comments again. First, the presence of the delta symbols ensures that in order to create a fermion in a given mode, i.e. in a state described by one of the basis vectors, this mode must be empty; similarly, destroying a fermion in a mode requires that mode to be initially occupied by one particle. This is a natural implementation of the exclusion principle which requires that a mode can only be either empty or occupied by one and only one particle.

Finally, we note that the fermionic number operator can be written as N = ∞ i=1 a † i a i , and that the following commutators are identical to their bosonic counterparts,

N, a † i = a † i [N, a i ] = -a i . (51)

Modes

As seen in Section 3.4, a Fock space H ∞ is built starting from a single-particle Hilbert state H of arbitrary dimension, of which we single out a basis (e k ). Basis vectors of the singleparticle Hilbert space H are typically called modes in physics. Fock states |n 1 , n 2 , . . . form a basis of the Fock space, where the notation stands for the (anti-)symmetrization of the tensor product of n 1 times the basis vector (mode) e 1 , n 2 times the basis vector e 2 , etc.

There is a natural isomorphism of Hilbert spaces between the Fock space H ∞ built on H and the tensor product k (Ce k ) ∞ of the Fock spaces built on each of the one-dimensional (or "single-mode") Hilbert spaces Ce k . Indeed, taking bosons as an example, one simply needs to map the Fock state

|n 1 , n 2 , . . . ∈ H ∞ to k S(e k ⊗ • • • ⊗ e k n k times ),
where S is the symmetrization operator, and should be replaced by A for fermions.

Thinking of the Fock space as a product of Fock spaces across modes is often implicit in physics texts. 23 It is, for instance, customary to introduce sophisticated states assuming a single mode, i.e., that H is one-dimensional, and then to write tensor products of such states across modes to cover the case of an H of arbitrary dimension. A concrete example will be bosonic coherent states in Section 3.6.

Field operators

There are often several natural bases for a physical situation, for instance a natural basis to describe a source of radiation, and a natural basis to describe a measurement. The linearity of the definition of ladder operators yields easy "change of basis formulas" connecting ladder operators in two different bases. Formally, let (|u i ) i∈I and (|v i ) i∈I be two bases of the single-particle Hilbert space H. The creation and annihilation operators transform exactly as a regular change of basis, i.e.,

a † (v k ) = i a † (u i ) u i |v k and a(v k ) = i a(u i ) v k |u i . ( 52
)
Note that we implicitly assumed that the two bases were actual Hilbert bases(i.e., countable), but the formula naturally extends in a weak sense to "generalized orthogonal bases", such as the position basis for the harmonic oscillator in Section 3.2.2. The creation and annihilation position field operators are precisely the creation operator and the annihilation operator corresponding to the position basis, and are usually denoted by ψ † (x) and ψ (x), respectively. Their physical interpretation is that the operator ψ † (x) creates a particle of a certain type at the position x. They are called (position) field operators because they are the building blocks for position-dependent operators, which we shall see are the mathematical description of quantum physical fields. As ladder operators, the field operators also satisfy the commutation relations corresponding to the particles being described, say for bosons

[ψ(x), ψ(y)] = 0 , ψ † (x), ψ † (y) = 0 , ψ(x), ψ † (y) = δ(x -y) . ( 53
)
To illustrate a generalized case of the change of basis in ( 52), consider the other fundamental basis for the harmonic oscillator, the momentum basis |k . Since, in an appropriate sense k|x = e -ik•x (2π) d , Equation ( 52) then yields the Fourier transform relationship between position and momentum field operators,

ψ † (x) = a † (k)e -ik•x dk (2π) d and ψ(x) = a(k)e ik•x dk (2π) d , ( 54
)
with a † (k) the creation operator corresponding to the momentum basis, creating a particle of momentum k. Once again, the important algebraic relations are the commutation relations written in the momentum basis, say for bosons,

[a(k), a( )] = 0 , a † (k), a † ( ) = 0 , a(k), a † ( ) = (2π) 3 δ k . ( 55
)
Note that the whole section translates to fermions, with the anticommutation relations (48) replacing the commutation relations of bosons. Finally, we make a point on notation. While it is customary in physics to use ψ(x) for the field operator of bosons and fermions, to avoid confusion, we write a(x) for a position bosonic field operator, and reserve ψ(x) for a position fermionic field operator.

Observables

Ladder operators allow one to define observables that model clusters of interacting particles. Depending on the number of particles interacting, we call the corresponding operators oneparticle observables, two-particle observables, etc.

One-particle observables. An example of one-particle observable is the number operator N in (46).

A generic one-particle operator O (1) is an operator that takes as in input a one-particle state |u i and gives as an output |u j , in a given basis (u i ) of the one-particle Hilbert space H. This means that we destroy the particle in the state |u i and create one in the state |u j . The operator O (1) can actually be rewritten [Cohen-Tannoudji et al., 2019b, Section XV.B] as

O (1) = ij O (1) ij a † (u j )a(u i ) . ( 56
)
The number operator N is diagonal, in the sense that the only nonzero terms in the sum (56) are those for which i = j. A typical non-diagonal one-particle observable correspond to a hopping between states. For example, if i, j are labels for sites on a lattice, then ij J ij a † i a j with J ij a symmetric matrix implements the hopping of particles from i to j with amplitude J ij . Two-particle observables. Two-particle observables are used to describe the interaction between two particles. The intuition is that an interaction takes two particles in a state |u k ⊗ |u l and yields a new two-particle state |u i ⊗ |u j . Acting on the Fock space, a generic two-body operator O (2) can be written using ladder operators as

O (2) = ijkl O (2) ijkl a † (u i )a † (u j )a(u l )a(u k ) ;
(57) see [Cohen-Tannoudji et al., 2019b, Chapter XV, Section C] for manipulations using multiparticle observables. A simple Hamiltonian featuring one-and two-particle observables is the Bose-Hubbard model [START_REF] Fisher | Boson localization and the superfluid-insulator transition[END_REF]],

H = i ω i N i + i,j J ij a † i a j + i U i N i (N i + 1), (58)
where (ω i , J ij , U i ) are real numbers that parametrize the model, and N i = a † i a i counts the particles in mode i. The Hamiltonian ( 58) is composed of both one-particle and two-particle contributions. The two-particle term in N i (N i + 1) models the extra energetic cost to pay to add a particle on a site i when some are already present.

Wick's theorem

Wick's theorem is a cornerstone of quantum field theory and the backbone of perturbation theory. It gives rise to the famous Feynman diagrams, and importantly for us, it yields the permanents and determinants in the correlation functions of fields, later to be turned into correlation functions of point processes when we introduce physical detection in Section 4 and Section 5.

In essence, Wick's theorem is a generalization of the calculation of the moments of multivariate Gaussian distributions to the case of Gaussian density matrices. We follow the derivation from [Cohen-Tannoudji et al., 2019b, Appendix C XVI].

Gaussian density matrices

A commonly used mixed state on the Fock space H ∞ of either fermions or bosons corresponds to the so-called Gaussian density matrices, defined, when it exists, as the unit-trace, self-adjoint operator

ρ = 1 Z e - i,j a † i Mij aj , ( 59
)
with M ij = M ji . As in Section 3.1, the definition (59) relies on the spectral calculus: ρ is the operator with the same eigenvectors as the argument of the exponential, but with the exponential applied to the corresponding eigenvalues, followed by division by Z. The normalization constant

Z = Tr e - i,j a † i Mij aj
is called the partition function in statistical physics. The use of of the word Gaussian in the name is sometimes confusing to non-physicists. It is motivated by the quadratic expression in the exponential, and the fact that, as we shall see below, the moments of products of creation and annihilation operator under these density matrices behave similarly to the moments of multidimensional Gaussian distributions. Finally, we note that there are conditions on the operator appearing in the exponential in (59) for ρ to be a proper mixed state, in particular to guarantee that Z < ∞. These conditions are best discussed for any particular M .

Example 4 (Grand canonical ensemble). Consider a system described by a Hamiltonian H acting on a Fock space of either bosons or fermions, with the Hamiltonian being quadratic in the ladder operators. In other words, we require that the Hamiltonian is a one-particle observable (56), like the Bose-Hubbard Hamiltonian (58) with U i = 0. Physicists usually consider additional fluctuations in energy and number of particles, and describe the system at thermal equilibrium by the so-called grand canonical ensemble. Formally, this amounts to considering the Gaussian density matrix

ρ GC = e -β(H-ζN ) Z GC (60)
where N = i a † i a i is the number operator, β > 0 is the inverse temperature and ζ > 0 is the chemical potential, which can be adjusted to change the average number of particles in the system. The normalization Z GC = Tr e -β(H-ζN ) is the so-called grand-canonical partition function.

As with any state ρ, a natural question is to compute the expectation A ρ = Tr (ρA) of an observable A; see Section 3.1. It turns out that for Gaussian density matrices, this is a simple mechanical computation using canonical (anti)commutation relations, as soon as A is a product of linear combinations of annihilation and creation operators. The result of this computation is precisely Wick's theorem, a close parent to Isserlis' theorem in classical statistics on the computation of moments of a Gaussian distribution. Up to a change of basis, and because we only care about linear combinations of ladder operators, it is enough to treat the case of

ρ = 1 Z e - i νia † i ai . ( 61
)
In the rest of this section, we henceforth assume ρ to be in diagonal form (61).

Where permanents and determinants appear

For brevity, every average in this section is implicitly meant as under the Gaussian density ρ in ( 61), and we thus write • instead of • ρ . Additionally, to treat both bosons and fermions in a single theorem, we define the generalized commutator [u, v] η := uv -ηvu, where η = ±1. We note that both the CCR (44) and the CAR (48) can be rewritten as

[a i , a j ] η = 0, a † i , a † j η = 0, and a i , a † j η = δ ij , ( 62
)
with η = 1 for CCR and η = -1 for CAR. We note right away that the density matrix (61) conserves the number of particles, in the sense that [ρ, N ] = 0, with N = i a † i a i . This leads to the following important lemma.

Lemma 2 (Creation and annihilation numbers must match). Consider the Gaussian density matrix ρ given by (61), acting on either the bosonic or fermionic Fock space, and an operator O nc,na composed of the product, in any order, of n c creation operators a † i1 , . . . , a † in c

and n a annihilation operators a i1 , . . . , a in a . Then O nc,na = 0 unless n c = n a .

Note that this implies that if the total number of creation and annihilation operator is odd, the average necessarily vanishes.

+ p a = n a . Then [N, O nc,na+1 ] =[N, O mc,ma ]a i O pc,pa + O mc,ma [N, a i ]O pc,pa + O mc,ma a i [N, O pc,pa ], =(m c -m a )O nc,na+1 -O nc,na+1 + (p c -p a )O nc,na+1 , =(n c -(n a + 1))O nc,na+1 . (63)
The case O nc+1,na is similar. Now, for an even integer N = 2k, define a contraction of order N as a permutation σ ∈ S N such that σ(1) < σ(3) < ... < σ(2k -1), and σ(2i -1) < σ(2i) for i = 1, . . . , k. For instance, there are three contractions of order 4, namely 1 2 3 4 1 2 3 4 , 1 2 3 4 1 3 2 4 , and 1 2 3 4 1 4 2 3 .

(64)

In words, contractions are built as follows. Starting from Z = {1, . . . , N } and i 1 = 1, pair i 1 with an arbitrary i 2 ∈ Z \{i 1 }. Then select the smallest element i 3 of Z \{i 1 , i 2 } and pair it with any element of Z \ {i 1 , i 2 , i 3 }. Repeat the procedure until all integers in Z have been paired. The corresponding permutation σ is the one such that (σ(1), σ(2), . . . , σ(N )) = (i 1 , i 2 , . . . , i N ). This constructive definition shows that there are actually (N -

1)!! = 1 × 3 × 5 • • • × (N -3) × (N -
. b N = σ contraction η ε(σ) b σ(1) b σ(2) b σ(3) b σ(4) . . . b σ(N -1) b σ(N ) , ( 65
)
where η = +1 for bosons and η = -1 for fermions. The sum is over contractions of order N , with parity ε(σ).

Before proving Wick's theorem for the sake of completeness, we give two illustrative examples. First, when N = 4, remembering the three contractions (64) of order 4, we have

b 1 b 2 b 3 b 4 = b 1 b 2 b 3 b 4 + η b 1 b 3 b 2 b 4 + b 1 b 4 b 2 b 3 . ( 66
)
Second, Wick's theorem is the mathematical reason why permanents and determinants will appear when we consider point processes of detection times in Sections 4, 5, and more abstractly in Section 6. This application is so fundamental to our paper that we highlight it here.

Example 5 (Wick's theorem for coherence functions). We are interested in the expected number of k-uplets of particles simultaneously appearing in (distinct) modes 1 to k. Since the operator N i = a † i a i counts the number of particles in mode i, we aim to compute

N 1 . . . N k = a † 1 a 1 . . . a † k a k .
Note that all N i s commute, so that the order is irrelevant. We now use the (anti-)commutation relations (62) to bring all creation operators to the front,

N 1 . . . N k = η 1+2+•••+(k-1) a † 1 a † 2 . . . a † k a 1 a 2 . . . a k .
By convention, we also put the annihilation operators in decreasing order, which removes the sign,

N 1 . . . N k = a † 1 a † 2 . . . a † k a k a k-1 .
. . a 1 . Now, we are ready to apply Theorem 3. Upong noting, thanks to Lemma 2, that pairing two creation or two annihilation operators results in a zero average, the only non-zero terms in (65) result from permutations σ ∈ S 2k of the form

σ = 1 2 3 4 5 • • • 2k -1 2k 1 ν(k + 1) 2 ν(k + 2) 3 • • • k + 1 ν(2k) ,
where ν is a permutation of {k + 1, . . . , 2k}. Composing σ with 2

× (1 + 2 + • • • + (k -1)) transpositions, we obtain 1 • • • k k + 1 • • • 2k 1 • • • k ν(k + 1) • • • ν(2k) , so that ε(σ) = ε(ν).
In particular, Wick's theorem yields

N 1 . . . N k = ν∈S k η ε(ν) k i=1 a † i a ν(i) ,
where we recognize the permanent or determinant of the matrix ( a † i a j ) 1≤i,j≤k , depending on η = ±1.

We now prove Wick's theorem for the sake of completeness.

Proof. We want to compute the expectation Tr (O N ρ), where

O N = b 1 b 2 . . . b N , b i = α A i,α a † α + B i,α a α , and ρ = e - α ναa † α aα /Z,
where we have expanded the b i 's in the basis in which the density matrix is diagonal.

Let us start with two general remarks. First, by Lemma 2, we can restrict the number of b i 's to be even, since otherwise the average vanishes. We henceforth assume that N = 2k is even. We are going to repeatedly use the generalized commutation relations. In particular, we note that, for any η and i, j, [b i , b j ] η is a multiple of the identity and thus commutes with all operators.

We are now ready to compute

O N = b 1 . . . b N = Tr (b 1 . . . b N ρ).
The general idea is (i) to push b 1 through all the other operators b i using the (anti)commutation relations. Then, (ii) using the cyclicity of the trace, we bring b 1 back in front of the rest. And finally, (iii) using that ρ is Gaussian to recover the original average O N up to a constant.

Starting this programme, we have

O N = [b 1 , b 2 ] η b 3 . . . b N + η b 2 b 1 . . . b N , = [b 1 , b 2 ] η b 3 . . . b N + η[b 1 , b 3 ] η b 2 b 4 . . . b N + η 2 b 2 b 3 b 1 . . . b N .
Iterating the same steps, we write

O N = η N -1 b 2 b 3 . . . b N b 1 + N j=2 [b 1 , b j ] η η j-2 b 2 . . . b j-1 b j+1 . . . b N . ( 67
)
Now, by the cyclicity of the trace, we have

b 2 b 3 . . . b N b 1 = Tr (ρb 2 b 3 . . . b N b 1 ) , ( 68
) = Tr (b 1 ρb 2 b 3 . . . b N ) . ( 69
)
Let us for a moment assume that b 1 is either a specific creation or annihilation operator, i.e. b 1 = a ( †) α for some specific α. Now, using the explicit Gaussian form of the density matrix, we have b

1 ρ = λρb 1 , with λ = e -να /Z if b 1 = a α , and λ = e να /Z if b 1 = a †
α , as can be seen by checking the action of these operators on the basis of Fock states. We thus come to

b 2 b 3 . . . b N b 1 = λ b 1 . . . b N . ( 70
)
Since η = ±1, Equation ( 67) is simplified into:

b 1 . . . b N = N j=2 η j [b 1 , b j ] η 1 -ηλ b 2 . . . b j-1 b j+1 . . . b N . ( 71
)
In particular, note that (71) with N = 2 becomes

b 1 b j = [b 1 , b j ] η 1 -ηλ . ( 72
)
Plugging this into (71), we conclude that

b 1 . . . b N = N j=2 η j b 1 b j b 2 . . . b j-1 b j+1 . . . b N . ( 73
)
Since this last expression is explicitly linear in b 1 , it is actually valid for any linear combination of creation and annihilation operators. In particular, (73) replaces the calculation of an average of a product of N operators by the weighted sum of N -1 averages of products of N -2 operators. Repeating the procedure for each product of N -2 operators in (73), we obtain Wick's theorem. This is the simplest form of Wick's theorem. It can be further generalized to Hamiltonians with in addition the "anomalous" quadratic terms nm A nm a n a m + nm Ānm a † n a † m , which arise in mean-field theories of interacting Bose gases and superconductors, and generate a i a j = 0; or to the case of Hamiltonians with linear terms, which generate nonzero a i [Berezin, 1966, Chapter 3].

Remark 4. The proof of Wick's theorem also holds for the vacuum, and the vacuum is also usually said to be a Gaussian state. The vacuum can also be obtained as a limiting state when the temperature 1/β is taken to zero and ζ = 0; see the computations in Section 6.

Finally, note that the definition of a Gaussian state is not universal and may be community-dependent. Quantum opticians, for instance, might define a Gaussian state as one that has a Gaussian Wigner transform; see Section 4.

Bosonic coherent states model classical fields

Many quantum systems have a classical (i.e., non-quantum) description, and physicists also sometimes start from a classical description to build a quantum theory, a procedure known as quantization. But in the end, a physical system is fundamentally quantum and we should try to understand how a classical behavior emerges from a quantum description. This phenomenon is known as decoherence, and is a vast modern research programme. An important part of the answer is to build quantum states that behave as closely as possible to classical ones. [START_REF] Glauber | Coherent and incoherent states of the radiation field[END_REF] made seminal contributions in that regard, defining coherent states for systems of photons. 24 We follow the introduction of bosonic coherent states in [Mandel and Wolf, 1995, Chapter 11], before pointing out the strong connections with the subfield of signal processing called time-frequency analysis. We defer the more subtle discussion on fermionic coherent states to Section 5.4.

Definition and properties

As discussed around Equation ( 26), measuring an observable is modelled by sandwiching the state operator. After measurement, the quantum state is thus different from the initial state. This is very different from classical physics, where it is possible, in principle, to perform passive measurements, e.g. looking at a screen, or recording the intensity of the electromagnetic field. Consider henceforth the boson Fock space built using the modes of the harmonic oscillator from Section 3.2.2. If we are to build a state of this Fock space that behaves classically, we would intuitively like to take a state that remains unchanged when we measure it. Since most measurements involve absorbing bosons (think photons in optics; see Section 4.1), good candidates for classical-like states are the left eigenstates of the annihilation operator, or equivalently the right eigenstates of the creation operator.

To make things concrete, consider a single mode for simplicity, that is, a Fock space built on a one-dimensional H = Ce 1 ; see Section 3.4.3. Denote the ladder operators by a = a(e 1 ) and a † = a † (e 1 ). Our candidate classical-like states are |α ∈ H ∞ Bosons such that

a |α = α |α , α ∈ C . ( 74
)
Informally, 25 the sandwiched state (26) that appears when modeling absorption mechanisms 24 The name coherent state comes from the concept of coherence in optics. 25 For details, see the Lindblad equation in [START_REF] Bouten | An introduction to quantum filtering[END_REF].

is then a |α α| a † = |α| 2 |α α|, so that coherent states are stable under absorption of a photon.

Solving the eigenstate equation ( 74) for |α = n c n |n expressed in the Fock basis, we write

n≥1 c n √ n |n -1 = n≥0 c n |n , so that c n = αc n-1 / √ n for n ≥ 1. This leads to |α = e -|α| 2 /2 n≥0 α n √ n! |n , ( 75
)
up to a complex number of unit modulus, where we have determined the modulus of the normalizing constant by imposing α|α = 1. The states (75), labeled by nonzero complex numbers α, are called canonical coherent states. By construction, they are mixed states formed as an infinite linear combination of Fock states. If we were to measure the number of photons in |α , we would obtain n with probability | α|n | 2 ∝ |α| 2n /n!, and we recognize the Poisson distribution with mean |α| 2 . In particular, the Poisson distribution is in a sense the maximum entropy distribution for independent counts [START_REF] Harremoës | Binomial and Poisson distributions as maximum entropy distributions[END_REF]: for a fixed mean λ > 0, the maximal entropy of a sum of k independent Bernoulli variables is increasing with k, and converges to the entropy of the Poisson variable with mean λ. In that limited sense, one can think of the coherent states as states with a maximally uncertain number of particles.

An equivalent point of view on coherent states, which sometimes serves as their definition in mathematical physics [START_REF] Ali | Coherent states, wavelets and their generalizations[END_REF], is that they are naturally associated to a projective representation of the translation group in phase space,26 the complex plane indexed by α.

To see this, start from the description of Fock states |n = (a † ) n / √ n! |0 . The coherent state ( 75) is thus obtained from the vacuum by the action of an operator D(α) on the boson Fock space called a displacement operator,

|α = D(α) |0 , D(α) = e αa † -α * a . ( 76
)
This formula is a consequence of the CCRs and the so-called Baker-Campbell-Hausdorff formula expanding the exponential of a sum of (non-commuting) operators. The same formula yields properties of D like

D -1 (α) a D(α) = a + α (77) D(α) D(β) = e i Im{(α β * )} D(α + β) . ( 78
)
These relations are central to the theory of coherent states. First, they show that D is a projective representation of the translation group in phase-space. Second, they naturally generalize to other groups, yielding coherent states for different systems than the harmonic oscillator; see e.g. spin coherent states. In fact, coherent states have been defined for any locally compact Lie group [START_REF] Perelomov | Coherent states for arbitrary lie group[END_REF][START_REF] Arecchi | Atomic coherent states in quantum optics[END_REF][START_REF] Zhang | Coherent states: theory and some applications[END_REF].

Another way to justify that coherent states are almost classical is through Heisenberg's uncertainty principle (27). Remember that for the harmonic oscillator, position and momentum operators X, P are incompatible observables, and that both can be expressed as sums of ladder operators; see the discussion around Equation (37). With the notation of Section 3.1, one can check that the product of σ X,|α α| σ P,|α α| is minimal among states. One says that coherent states saturate Heisenberg's uncertainty relations.

Another quasi-classical aspect of coherent states is their dynamical evolution. Since the Fock states are eigenvectors of the Hamiltonian of a harmonic oscillator H = ωa † a, U (t) = e -itH applied to α decomposed on the Fock basis (75) yields

|α(t) = αe -iωt . ( 79
)
In phase space, i.e., the complex plane parametrized by α, the coherent state parameter α(t) simply rotates over time, with an angular frequency ω. This is precisely the phase-space behavior of a classical harmonic oscillator. Finally, we mention two more properties of coherent states. First, they are an overcomplete family of the Fock space, with a reconstruction formula

1 π C |α α| dα = 1 . ( 80
)
The fact that they are not an orthonormal basis is further seen from the overlap

β|α = exp α * β - 1 2 |α| 2 - 1 2 |β| 2 , α, β ∈ C. (81)
A good reference for such properties of coherent states of the harmonic oscillator is [START_REF] Folland | Harmonic analysis in phase space[END_REF], or, actually, books on harmonic analysis applied to time-frequency signal processing [START_REF] Gröchenig | Foundations of time-frequency analysis[END_REF], on which we say a few words in the next section. Finally, we have built here bosonic coherent states assuming a single-mode Fock state. Multi-mode coherent states are naturally obtained as tensor products of single-mode coherent states, see Section 3.4.3, and one typically writes them

|α = |α 1 ⊗ |α 2 ⊗ . . .

The Husimi distribution and time-frequency analysis

Single-mode bosonic coherent states of the harmonic oscillator are intimately linked to the subfield of signal processing called time-frequency analysis; see e.g. [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF][START_REF] Gröchenig | Foundations of time-frequency analysis[END_REF]. Consider e.g. the short-time Fourier transform

V g : L 2 (R) → L 2 (R 2 ), defined as V g (f ) : t, ω → f (τ )g(τ -t)e -2iπωτ dτ ,
where g(t) = 2 1/4 e -πt 2 is a unit-norm Gaussian window. Intuitively, if f is a signal, say representing an audio recording, then |V g (f )(t, ω)| will be large whenever frequency ω is present at time t. In other words, V g (f ) is akin to a musical score. Time-frequency transforms such as V g are thoroughly used in signal processing, for tasks such as detection or estimation of signal corrupted with noise [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF]. Mathematically, V g is a unitary linear operator, which can be inverted [START_REF] Gröchenig | Foundations of time-frequency analysis[END_REF]. Furthermore, it is tightly linked to decomposing a state into coherent states. Indeed, up to a non-vanishing factor,

V g (f )(t, -ω) ∝ α|ϕ , ( 82
)
where α = t+iω indexes a coherent state (75), and the state ϕ is described in the Fock basis by ϕ|n = f h * n dt, where h n is the nth Hermite function, a special basis of L 2 (R) that we already met in Section 3.2.2. Physicists call the squared modulus of the right-hand side of Equation ( 82) the Husimi distribution of ϕ, and think of it as a phase space representation of the state ϕ, in pretty much the same way signal processers think of the spectrogram t, ω → |V g (f )(t, ω)| 2 as a time-frequency representation of the signal f .

As a final note and to introduce another related point process, one can prove that, in a suitable sense, the STFT V g of white Gaussian noise is equal to the so-called planar Gaussian analytic function, up to an nonvanishing term [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF]Hardy, 2019]. Similarly, if one had access to a white noise-like state |ξ such that ξ|n ∼ N C (0, 1) are i.i.d. complex unit Gaussians, then its coherent-state "decomposition" α → ξ|α would be the so-called planar Gaussian analytic function, up to a nonvanishing term again. This planar Gaussian analytic function, along with the point process of its zeros, plays a role in the analysis of chaotic dynamical systems in statistical physics [START_REF] Nonnenmacher | Anatomy of quantum chaotic eigenstates[END_REF], where it is also called the chaotic analytic function [START_REF] Hannay | The chaotic analytic function[END_REF]. The zeros of the planar Gaussian analytic function are not a DPP. Yet, being the zeros of a random smooth function, they have a repulsive behaviour. They actually share many properties with the Ginibre ensemble, itself a fundamental DPP [START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF].

Photodetection and bosonic coherences

To obtain signals from a physical situation, the picture to have in mind consists in the processes of emission, propagation, and detection of radiation. The first two processes are usually described by physicists using one or several sources, classical or quantum, and a field theory, classical or quantum. Questions that arise regard, for instance, the dynamics of the source(s), how they generate radiation, and the physical properties of this radiation. The third process, detection, lies at the interface of physics and signal processing. In this section, we discuss simple standard models for the sources, the field, and the detection for bosons, having in mind photodetection. Following the footsteps of [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF], our goal is to arrive at the description of the measurement of arrival times at a simple detector as a point process. We shall see how some of the point processes introduced in Section 2 naturally appear from Glauber's coherent state decompositions. As a guiding thread, we will comment on the HBT effect in its different guises.

In Section 4.1, we decribe simple models for sources, fields, and detectors, and arrive at the marginal probability of detecting an event at a given time. In Section 4.2, we examine the correlation between more than one detection events. The central objects that encapsulate information about detection times are the coherence functions. In Section 4.3, we examine how to turn the coherence functions of a physical detection setup into the correlation functions of a point process. We discuss special particular cases, and recover the permanental point processes already announced using a semi-classical treatment in Section 2.3. Finally, in Section 4.4, we show that single-atom sources can yield antibunching detection events. This is a warning that bosons should not be identified with bunching particles in general, but that the properties of the source should be mentioned.

Modeling photodetection events

Photomultipliers are experimental devices that turn incoming radiation into a measurable electric current. The theory of photodetection aims at understanding what kind of signals can be observed as the output of a photomultiplicator. As we shall see, these signals are directly related to the so-called coherence functions27 of the field, quantum or classical.

Modeling the radiation

Before discussing how we introduce coherence functions, we need to specify how we describe the radiation we want to probe, using the objects of quantum field theory introduced in Section 3. Note that, while the framework below is enough to describe photodetection, a complete description would require to go into the details of the construction of (relativistic) quantum field theories and especially its dynamical content [START_REF] Cohen-Tannoudji | Atom-photon interactions: basic processes and applications[END_REF].

The electromagnetic field is described by a set of quantum fields satisfying some commutation relations. For photodetection, the dominant contribution comes from the electric field at the detection device. As a reminder, a field is a quantity that depends on both position and time. A bosonic quantum field is a collection of operators associated to an indefinite number of bosonic particles, called e.g. photons, on a Fock space H ∞ Bosons (see Section 3.4). Formally, the state of the field is built by acting on the vacuum |0 with creation and annihilation operators a † k,ε and a k,ε [Mandel and Wolf, 1995, Chapter 10]. The quantities k and ε index modes of the field (i.e., solutions to Maxwell's equations), and typically correspond to the momentum and the polarization of the photon.

For simplicity, and because it already contains the ingredients that relate the point processes of Section 2 to photodetection, we consider a fixed linear polarization, and henceforth drop the index ε. This corresponds to an assumption we made in the semi-classical treatment of Section 2.3. The electric field is then described by a single (instead of one per space coordinate) time-and space-dependent operator E(r, t), acting on the Fock space H ∞ Bosons .

To obtain the explicit form of the electric field operator, the stategy is to consider the field as enclosed in a (large) box and perform the so-called canonical quantization of Maxwell's equations. We refer the reader to the literature for such details, e.g. [Mandel and Wolf, 1995, Chapter 10], and consider as given the resulting form of the field, namely

E(r, t) = k N k a k e i(k•r-ω k t) + a † k e -i(k•r-ω k t) , ( 83
)
where N k is a normalization coefficient, r ∈ R d indexes space and t ∈ R indexes time. The time dependence in ( 83) is a consequence of working in the interaction (or equivalently, so far, Heisenberg) picture. The explicit plane-wave form corresponds to the evolution implied by the free electromagnetic Hamiltonian

H FF = k ω k a † k a k . (84)

Modeling the detector: the first-order coherence function

The usual image to have in mind for a detector is a two-level atom like the qubit described in Section 3.2.1, described quantum-mechanically with an eigenfrequency ω eg . The atom interacts with the electric field E(r, t). In mathematical terms, the Hilbert space of the detector is finite-dimensional, isomorphic to C 2 ; see Section 3.2.1. Since the Hilbert space of the electric field is H ∞ Bosons , the joint system of the detector and the field is C 2 ⊗ H ∞ Bosons . To describe the time evolution of the system, we use the interaction picture and decompose the Hamiltonian as a free and an interaction part

H = H 0 + H I , (85)
see Section 3.1. The free part is taken to be H 0 = H FQ ⊗ 1 + 1 ⊗ H FF , a combination of the free Hamiltonians of a qubit and the field, respectively defined in ( 34) and ( 84).

To describe the interaction Hamiltonian, let first q ∈ R model the modulus of the dipolar moment multiplied by the charge of an electron. Once gain, Maxwell's equations hint that the interaction can be taken as [START_REF] Cohen-Tannoudji | Atom-photon interactions: basic processes and applications[END_REF]]

H I (t) = -q σ x (t) ⊗ E(r, t), ( 86
)
where

σ x (t) = U FQ (t)σ x (t)U FQ (-t) = U FQ (t) 0 1 1 0 U FQ (-t),
is the observable σ x from Section 3.2.1, evolved through time using the group U FQ (t) corresponding to the free part of (85), as befits the interaction picture. 28We want to compute the probability p ωeg (r, t) of a detector at position r to be in its excited state |e after a time t, knowing that the initial state of the detector is |g and the field is any arbitrary initial state, but irrespective of the final state of the field. We assume that the initially prepared state of the detector and the field is |g ⊗ |i = |g, i , where |g is the ground state of the detector and |i is an arbitrary initial state of the field. The corresponding density matrix is thus the projector ρ(0) = |g, i g, i|. In the interaction picture described in Section 3.1, the state evolves until time t through the action of the evolution operator

V (t) = U † 0 (t)U (t) as ρ(t) = V (t)ρ(0)V † (t).
Measuring the detector in its excited state corresponds to the observable

|e e| ⊗ 1 = |e e| ⊗ f |f f | = f |e, f e, f | , ( 87
)
where f indexes Fock states, and we decomposed the identity operator onto the Fock states, 1 = f |f f |. Now the observable (87), in the interaction picture, evolves as

U † 0 (t) f |e, f e, f | U 0 (t).
We are now ready to compute the probability p ωeg (r, t) of the detector being measured in its excited state at time t. By Born's rule (see Section 3.1), this probability is the average value of the evolved observable onto the evolved state ρ(t), namely

p ωeg (r, t) = Tr   U † 0 (t)   f |e, f e, f |   U 0 (t)V (t) |g, i g, i| V † (t)   (88a) = f Tr (|e, f e, f | U (t) |g, i g, i| U (-t)) (88b) = f | e, f | U (t) |g, i | 2 (88c) = f | e, f | V (t) |g, i | 2 , ( 88d
)
where the last line comes from the definition V (t) = U † 0 (t)U (t) and the fact that |e, f is an eigenvector of H 0 (t), so that applying U 0 (t) to |e, f simply multiplies by a complex number of modulus 1. Note also that the dependence in r is hidden in V (t).

Up until now, all the discussion is exact. We now introduce a physical assumption that allows for a simple expression of V (t).

Assumption 5 (Weak Coupling). The interaction between the detector and the field is assumed to be weak, allowing us to treat the dynamics using perturbation theory. Formally, this means that we assume the existence of a series expansion of the evolution operator in terms of the interaction Hamiltonian; see the discussion around (32). To first order, this amounts to

V (t) 1 - i t 0 H I (t )dt . ( 89
p ωeg (r, t) = q 2 [0,t] 2 e iωeg(t -t ) i|E(r, t )E(r, t )|i dt dt . ( 91
)
This equation is straightforwardly generalized to the case where the initial state of the field is not a pure state but a mixed state, described by a density matrix ρ = i p i |i i| in some basis (|i ) of initial states; see Section 3.1. The probability of a detection event between 0 and t then becomes

p ωeg (r, t) = q 2 [0,t] 2 e iωeg(t -t ) E(r, t )E(r, t ) ρ dt dt . ( 92
)
Equation ( 92) underlies most of photodetection theory. It relates properties of the field, encoded in a coherence function E(r, t )E(r, t ) ρ , to the excitation probability of the detector.

The second important approximation usually made in photodetection theory is the rotating wave approximation (RWA), which translates the intuition that a proper photodetection event is an absorption of an excitation. To motivate and illustrate it, we consider first a simple example where the field is assumed to always be in a coherent state.

Example 6 (Monochromatic classical wave). We assume here that the field always remains in a coherent state as time evolves and is factorized with respect to the state of the qubit; see Section 3.6. Since the free part of the Hamiltonian preserves coherent states up to a phase, what we are assuming is that the evolution under the interaction Hamiltonian approximately preserves coherent states. Physically, this assumption corresponds to, for instance, a macroscopic source emitting a coherent state or a statistical mixture thereof, like a common lightbulb: the fact that one observes the emitted light does not change its classical character. We also assume, for simplicity, that we have an incoming monochromatic coherent state at frequency ω. Formally, we reduce the sum in (83) to a single term k such that ω k = ω, and denote by E(t) = E(r, t) the average value of this field at the fixed position r of the detector. The terms in r can be omitted by changing the time origin. More explicitly, we set E(t) = E sin ωt. We want a coherent state |α such that the average value of the electric field operator

α|E(t)|α = α|a k e -iωt + a † k e iωt |α = N k (αe -iωt + α * e iωt )
is equal to E(t). We thus choose the coherent state parameter α to be -E/2iN k .29 Furthermore, assume that we initially prepare the detector in its ground state |g . In the interaction picture, our approximation that the state of the field remains the same coherent state30 and factorized allows to write

V (t) |g ⊗ |α ≈ (V sc (t) |g ) ⊗ |α . ( 93
)
The operator V sc (t) in ( 93) is defined as the interaction picture evolution operator of the so-called semi-classical Hamiltonian H sc (t) = -qE(t) σ x (t). Using the weak coupling approximation to evaluate the time evolution, we obtain

p ωeg (r, t) = Tr U † 0 (t) |e e| ⊗ 1 U 0 (t) V (t) |g g| ⊗ |α α| V † (t) (94a) = Tr |e e| ⊗ 1 V (t) |g g| ⊗ |α α| V † (t) (94b) = Tr V sc (t) |g g| V † sc (t) |e e| (94c) = iq t 0 E(t ) e iωegt dt 2 (94d) = qE 2 i 1 -e i(ωeg-ω)t ω eg -ω - 1 -e i(ωeg+ω)t ω eg + ω 2 , ( 94e
)
where the second line is obtained because the state |e is an eigenvalue of the free evolution operator U 0 (t), the third follows from the semi-classical assumption that the state of the field always stays in a coherent state, the fourth from the weak coupling assumption and the form of H sc , and the last from the monochromatic assumption. At this point, we make an extra assumption called the rotating wave approximation (RWA). In our special case, the RWA amounts to neglecting the second term of the righthand side of (94e). Because both ω and ω eg are strictly positive, the denominator of the second term never vanishes, unlike the denominator of the first term. Around resonance ω ≈ ω eg , the first term will largely dominate the second.

We make the RWA from Example 6 a general principle that applies to arbitrary quantum fields. It echoes the use of the analytic signal in representing a real classical field; see Section 2.

Assumption 6 (Rotating wave approximation). The photodetection response is a function of the positive frequency part of the quantum field only.

With Assumption 6, the transition probability (92) now reads

p ωeg (r, t) = e q 2 t 0 e iωeg(t -t ) G (1) ρ r t , r t dt dt , ( 95
)
where 31 we introduced the first-order coherence function

G (1) ρ r t , r t = E -(r , t )E + (r, t) ρ , ( 96
)
where E + is associated to the positive part of the spectrum (annihilation operators) and E -is associated to the negative part (creation operators). The coherence function ( 96) was originally introduced by [START_REF] Glauber | Coherent and incoherent states of the radiation field[END_REF]. 32 At least informally, it is not difficult to generalize the computations leading to (92) to more sophisticated detectors. In experiments, for instance, photodetection is often performed through ionization: an electron of an atom is extracted by the incident light and this event is amplified and detected. This detector is best described by a system having, on top of its ground |g state, a countable set of excited states |e i , and a continuum of so-called diffusive states |p , p ∈ R. This leads to non-compact Hamiltonian operators, having a spectrum that is partly discrete, partly continuous. Nonetheless, the computation of the transition probability p(r, t) follows from the same lines as above, with the new observation corresponding to projecting on all states that are not the ground state, that is, ( n |e n e n | + |p p| dp) ⊗ 1 instead of |e e| ⊗ 1. Because we are working with a first-order approximation in Assumption 5, this leads to a probability of transition that is an expectation of the single-state probability of transition (92), namely

p(r, t) = p ωeg (r, t)dξ(ω eg ) = [0,t] 2 κ(t -t )G (1) ρ r t , r t dt dt , ( 97
)
where ξ has a support that is the union of a discrete set, corresponding to the states |e n , and an uncountable set, corresponding to |p , p ∈ R, and the structure function is defined as κ(τ ) := e 2 e -iωτ dξ(ω).

We can think of either ξ or κ as being a characteristic of the detector, describing its time or frequency efficiency. For instance, if ξ is supported on a single point or closely around a single point, we only collect photons at a very specific energy, and call the detector narrowband. If ξ is supported on a large interval, we collect all photons equally, and call the detector broad-band.

For completeness, we close this section by recovering the semi-classical treatment of Section 2.3 through our quantum treatment, using the coherent state formalism of Section 3.6.

Example 7 (General rotating wave). By linearity arguments, the monochromatic light case of Example 6 can be used to derive the excitation probability in the case of a more general field. Thinking of the field being enclosed in a large box, so as to be able to work with Fourier series instead of Fourier transforms for simplicity, let us decompose the electric field E α (r, t) on a set of monochromatic plane waves

E α (r, t) = k N k α k e i(k•r-ω k t) + α * k e -i(k•r-ω k t) , ( 98
)
with N k a normalization coefficient and α k ∈ C the amplitude of the mode k. Note that E is a (real-valued) function, not an operator: we are still treating the field as a classical field. The first part of the sum corresponds to the positive frequency part of the field (the "analytic" field), while the second part corresponds to the negative frequency part E -. An important observable of the electromagnetic field is its intensity. A general result of Maxwell's equations [START_REF] Mandel | Optical coherence and quantum optics[END_REF] is that the energy carried by the field is given 31 The missing commas in the arguments are a convention that will become handy when we consider higher-order coherence functions.

32 As a side remark, remember that this coherence function is defined after choosing a specific polarization of the field. In the general case, the electric field should be expanded along a basis of the three-dimensional polarization space and we should consider all combinations of coherence functions

G (1) ij (t , t) = E - j (t )E + i (t)
ρ where i, j label the basis elements.

by the squared modulus of the field. Under our rotating wave Assumption 6, the intensity I(r, t) of the field reads

I(r, t) = E - α (r, t)E + α (r, t) = |E + α (r, t)| 2 . ( 99
)
Investigating this intensity, and especially its correlation structure in the case of random sources, has led to the major discoveries of photon bunching and in the end of the general theory of optical coherence in quantum theory [START_REF] Mandel | Optical coherence and quantum optics[END_REF]. More precisely, say the field ( 98) is not perfectly known; it is then reasonable to model it as a stochastic process, i.e., a random function. A very common assumption is to assume the process to be a stationary Gaussian process. One partial justification is that the field at a given spacetime point is the sum of the contributions of many independent sources, and is well approximated by a Gaussian. Formally, the field state becomes a statistical mixture of coherent states, a mixed state in the vocabulary of Section 3.1. The probability of excitation remains of the form ( 92), with first-order coherence function

G (1) ρ r t , r t = E α α|E -(r , t )E + (r, t)|α (100) = E α E - α (r , t )E + α (r, t) (101) = p classical (α)E - α (r , t )E + α (r, t) dα . ( 102
)
The expectation is over the stochastic process E. Physicists prefer to denote it informally using a probability density p classical over the Fourier coefficients α = (α k ) of the field in ( 98), and we use that notation here for future notational compatibility with Glauber-Sudarshan decompositions. In Section 4.2, we shall use Wick's theorem to derive higherorder coherence functions, and recognize the set of correlation functions of a permanental point process.

Role of the detector structure function

In this section, we consider the influence of the detector response κ on the probability of excitation (97).

Broad-band detector. Consider first a broad-band detector, i.e. a detector with structure function that is sharply peaked in time or, equivalently, does not depend much on the frequency. Ideally, this corresponds to κ(t -t ) = κ 0 δ(t -t ) with κ 0 a real constant. The probability of excitation (97) becomes

p(r, t) = κ 0 t 0 G (1) ρ r t , r t dt . ( 103
)
Anticipating over Section 4.3, we would like to map this probability to a point process on the positive half-line R + . One way to associate the two notions is to imagine t > 0 infinitely small in ( 103), and "re-setting" the detector in its ground state at each time multiple of t, keeping track of whether the detector was found in its excited state, say in the middle of each small time interval. Loosely speaking, in the limit t → 0, these detection times would form a point process on the positive half of the real line, with first correlation function G

(1) ρ r t , r t . To obtain a more formal correspondence, one needs to further model the detection process, for instance including de-excitation of the detector, see e.g. [START_REF] Bouten | An introduction to quantum filtering[END_REF].

Narrow-band detector. A narrow-band detector is a detector that is well-resolved in frequency, like the two-level atom we used in our derivation of the photodetection signal. At the extreme, we think of the structure function as a delta distribution in Fourier space at a frequency ω 0 , so that Fκ(ω) = κ 0 δ(ω -ω 0 ), and thus κ(t) = κ 0 / √ 2πe iω0t . The excitation probability (97) now reads

p(r, t) = κ 0 √ 2π [0,t] 2 G (1) ρ r t , r t e iω0(t -t ) dt dt . , ( 104
)
The meaning of this equation is even more transparent when we switch to the frequency representation of the coherence function. Indeed, let us define the Fourier representation of the first-order coherence function through the Fourier transform(-like) equation

G (1) ρ (r t , r t ) = R 2 G (1) ρ (r ω , r ω )e -iω t e iω t dω dω 2π , ( 105
)
where we follow the physicists' convention to denote both the function and its Fourier transform by the same notation and index the function by the name of its variable. Putting everything together, we arrive at

p(r, t) = κ 0 √ 2π R 2 G (1) ρ (r ω, r ω ) 2t sin(t(ω -ω 0 )) (ω -ω 0 ) dωdω 2π . ( 106
)
In the long time limit33 , the integration window is approximately equal to tδ(ω -ω 0 ). We end up with the simple expression

p(r, t) ≈ κ 0 t √ 2π G (1) ρ (r ω 0 , r ω 0 ) . ( 107
)
The interpretation of this expression is easier now. For a narrow-band detector in frequency space centered around a frequency ω 0 , the excitation probability is proportional to the diagonal ω 0 part of the first-order coherence function which is none other than the average intensity of the mode ω 0 .

Correlation between photodetection events

We now examine the joint probability of detecting two, and then more, photons.

Second-order coherence

Consider two detectors labelled 1 and 2 and placed at two distinct positions r 1 and r 2 . They are excited by an incident radiation, and we ask for the probability that detector 1 gets excited at time t 1 and detector 2 gets excited at time t 2 . Actually, the derivation of the joint probability of detection is quite similar to the derivation for the one photodection event in Section 4.1.2, and we will only give the general principles of the method here again. The dynamics are once again modeled through a Hamiltonian, whose interaction term is given by the following generalization of ( 86),

H I (t) = q 1 σ (1) x (t) ⊗ 1 ⊗ E(r 1 , t) + q 2 1 ⊗ σ (2) x (t) ⊗ E(r 2 , t) . ( 108
)
Such a form is natural to physicists, in the sense that it models the local interaction of each of two detectors with the field at their respective position.

Computing the joint probability of excitation requires to compute first a matrix element of the evolution operator in the interaction picture, which is again done assuming weak coupling. We will see that in order to have a non-trivial result, we need to push the expansion in Assumption 5 up to second order.

Assumption 7 (Weak Coupling; second order). The interaction between the detectors and the field is assumed to be weak, allowing to treat the dynamics using perturbation theory. To second order, this amounts to

V (t) 1 - i t 0 H I (t 1 )dt 1 - 1 2 t 0 t1 0 H I (t 1 )H I (t 2 ) dt 2 dt 1 . ( 109
)
Being interested in the joint probability of excitation after a time t, we have to compute the matrix element e 1 , e 2 , f |V (t)|g 1 , g 2 , i where (|g 1 , |g 2 ) and (|e 1 , |e 2 ) are the ground and excited states of both detectors. Since they are all pairwise orthogonal, the 1 term in (109) will give a zero contribution. As for the linear term, and given the form of the interaction (108) as a sum of two contributions from each detector, it will split into two terms, each acting trivially on one of the detectors, again yielding a zero result. For instance e 1 , e 2 , f |σ (1) x (t

) ⊗ 1 ⊗ E(r 1 , t)|g 1 , g 2 , i = e 1 |σ (1) x (t)|g 1 e 2 |g 2 f |E(r 1 , t)|i = 0
by orthogonality of |e 2 and |g 2 . This is the reason why we have to go to the second order in perturbation theory to compute the joint probability distribution. In general, to compute the joint probability distribution of N detectors clicking, we will need to push the expansion of Assumption 7 to the order N , since none of the lower-order terms contributes. Going back to N = 2, the quadratic term also simplifies with only the cross terms E(r 2 , t 2 )E(r 1 , t 1 ) giving a non-zero contribution. Using again the RWA of Assumption 6, we end up with

e 1 , e 2 , f |V (t)|g 1 , g 2 , i (110) = - 1 2 t 0 t t1 e iωeg 1 t1 e iωeg 2 t2 f |E + (r 2 , t 2 )E + (r 1 , t 1 )|i dt 2 dt 1 + 1 ↔ 2 = - 1 2 t 0 t 0 e iωeg 1 t1 e iωeg 2 t2 f |E + (r 2 , t 2 )E + (r 1 , t 1 )|i dt 2 dt 1 . ( 111
)
The notation 1 ↔ 2 is a shorthand for the same term with the two indices transposed. The second line is obtained by grouping the two terms into one and changing the integration range, noting that the involved operators commute. Just like we did in Section 4.1.2, we use Born's rule and sum over the final state of the field, since we are only interested in the state of the detectors, to obtain the joint excitation probability. Once again, we express it as the Fourier transform of a second-order coherence function. Assuming a fixed polarization for the field for simplicity, we have

G (2) ρ (r 1 t 1 , r 2 t 2 , r 2 t 2 , r 1 t 1 ) := E -(r 1 , t 1 )E -(r 2 , t 2 )E + (r 2 , t 2 )E + (r 1 , t 1 ) ρ . ( 112
)
By once again introducing structure functions κ 1 and κ 2 for both detectors, we obtain our final expression for the joint excitation probability

p(r 2 t, r 1 t) = t 0 κ 1 (t 1 -t 1 )κ 2 (t 2 -t 2 )G (2) ρ (r 1 t 1 , r 2 t 2 , r 2 t 2 , r 1 t 1 )) dt 1 dt 1 dt 2 dt 2 , ( 113
)
Intensity correlations. Of particular experimental relevance is the diagonal of the second-order coherence function, namely

G (2) ρ (r 1 t 1 , r 2 t 2 , r 2 t 2 , r 1 t 1 ) = E -(r 1 , t 1 )E -(r 2 , t 2 )E + (r 2 , t 2 )E + (r 1 , t 1 ) ρ . ( 114
)
For reference, if we were treating a classical field (i.e., modeled as a real-or complex-valued stochastic process rather than a collection of noncommuting operators), all the fields would commute and from the definition of the field intensity, we would obtain

G (2) ρ (r 1 t 1 , r 2 t 2 , r 2 t 2 , r 1 t 1 ) = EI(r 2 , t 2 )I(r 1 , t 1 ) . ( 115
)
Equation ( 114) is the non-commuting counterpart to (115), which gives us information about the intensity correlations of the field. The order of the operators in (114), called normal order, is crucial for the interpretation of the photodetection signal as the absorption of excitations by the detector, but also to take properly into account the fundamental quantum fluctuations that affect the detection signal. Intensity correlations can be probed experimentally. For instance, assume a perfect temporal resolution for each detector, and that they are functioning in the time intervals [0, t i ], i = 1, 2. This is modeled by choosing for structure functions

κ(t i -t i ) = κ i Θ [0,ti] δ(t i -t i ), i = 1, 2,
with κ 1 , κ 2 > 0. The probability of joint excitation is then given by

p(r 2 t 2 , r 1 t 1 ) = κ 1 κ 2 t2 0 t1 0 G (2) ρ (r 1 t , r 2 t , r 2 t , r 1 t )) dt dt . ( 116
)
In particular, taking the derivative in t 1 and t 2 yields the so-called transition rate,

Γ(r 2 t 2 , r 1 t 1 ) = κ 1 κ 2 G (2) ρ (r 1 t 1 , r 2 t 2 , r 2 t 2 , r 1 t 1 ) , ( 117
)
which is proportional to (114), and reduces to the intensity correlation function in the classical setting.

Higher-order coherences

The previous derivation can be generalized straightforwardly to obtain the joint probability of n detection events. Denoting the spacetime coordinates by x i = (r i , t i ), we define the nth order coherence function by

G (n) ρ (x 1 , • • • , x 2n ) = E -(x 1 ) • • • E -(x n )E + (x n+1 ) • • • E + (x 2n ) ρ , n ≥ 1. ( 118
)
This definition generalizes the first-and second-order coherence functions. Upon noting that Tr (ρO † O) ≥ 0 for any operator O, we can derive general inequalities for the coherence functions (118). For instance, consider the operator O = m k=1 λ k E + (x k ), for arbitrary λ k 's. The nonnegativity of the quadratic form Tr (ρO † O) implies a contraint in the underlying determinant. For m = 2, we obtain

G (1) ρ (x, x)G (1) ρ (y, y) ≥ |G (1) ρ (x, y)| 2 ,
which we interpret as a sign of bunching. Finally, we note that for general states, higher-order coherence functions are typically hard to compute. One special case is for a Gaussian state (59), where Wick's Theorem 3 applies. The nth-order coherence function can then be expanded as

G (n) ρ (x 1 , • • • , x 2n ) = σ∈Sn n i=1 E -(x i )E + (x n+σ(i) ) ρ . ( 119
)
We now have all the tools to identify the point processes behind the detection of photons.

Turning coherence functions into correlation functions

There is a specific set of states for which the coherence functions are particularly simple, namely the coherent states of Section 3.6. Specifically, take a multi-mode coherent state |α = ⊗ k |α k . These states are eigenstates of the electric field (83) itself, with eigenvalue the classical field from (98) , i.e.,

E + (r 1 , t 1 ) |α = E + α (r 1 , t 1 ) |α = k E + α k (r 1 , t 1 ) ⊗ k |α k . ( 120
)
Physically, coherent states correspond to ideal lasers. Now, by construction, the quantum coherence functions (118) are reduced to the classical coherence functions when the field is prepared in such a state, i.e., writing

x i = (r i , t i ), G (n) |α α| (x 1 , • • • , x 2n ) = α|E -(x 1 ) • • • E -(x n )E + (x n+1 ) • • • E + (x 2n )|α = n i=1 E + α (x i )E + α (x n+i ). (121)
This factorization across coherent states is what allows us recovering the correlation functions of a Poisson point process.

Example 8 (Ideal lasers yield Poisson processes). Take x i = x n+i = (r, t i ) for some fixed position r ∈ R 3 in (121), and assume a broad-band detector placed at r, see Section 4.1.3. The joint detection probabilities yield measurement times that follow a point process with correlation functions

ρ n (t 1 , . . . , t n ) = |E + α (t 1 ) • • • E + α (t n )| 2
. Such a point process indeed exists: it is the Poisson point process of first correlation function |E + α | 2 ; see Section 2.2. Now, consider a slightly more general state ρ that can be written as a statistical mixture ρ = E α∼p classical |α α| of coherent states. Concretely, this state describes non-ideal laser light, e.g. where there is uncertainty on which light amplitude is emitted. For instance, a single-mode randomly phased laser falls into this category, where p classical is uniform over a circle in the complex plane: the mixture is over all possible phases with a fixed intensity |α| 2 . As for the coherence functions, by linearity, we obtain

G (n) ρ (x 1 , • • • , x 2n ) = E α∼p classical n i=1 E + α (x i )E + α (x n+i ).
Example 9 (Mixture states in phase space yield Cox processes). Taking again x i = x n+i = (r, t i ) for some fixed position r ∈ R 3 in ( 121) and a broad-band detector, the detection times follow a point process with correlation functions

ρ n (t 1 , . . . , t n ) = E α∼p classical |E + α (t 1 ) • • • E + α (t n )| 2 .
The detection times thus form a Cox point process; see Section 2. In particular, arrival times tend to exhibit bunching: they are more clustered together than a Poisson process. Moreover, assuming that each α k in ( 98) is Gaussian, E α is a Gaussian process, and the detection times then form a permanental point process. Permanental point processes are the archetype of a point process exhibiting bunching; see Section 2.3 and a few samples in Figure 1(b). The fact that we obtain a permanental point process validates a posteriori the semiclassical analysis carried out in Section 2.

Additionally, note that while we mention ideally resolved detectors in time for simplicity, non-broad-band detectors also yield standard point processes when the detector's characteristic functions can be normalized. Indeed, consider for simplicity the point process formed by the detection times at detectors that share the same characteristic function κ. The correlation functions then correspond to the point process formed by the marks in a marked point process with independent marks. In the vocabulary of [Daley and Vere-Jones, 2003, Chapter 6], the ground process is formed by the ideal detection times, and for each such time t i , an actual detection time is drawn from the probability density function proportional to κ(• -t i ), independently from all other ideal detection times and marks. Now we consider states that are not statistical mixtures of coherent states, the socalled non-classical or genuinely quantum states. As we saw in Section 3.6, coherent states form an overcomplete family of states. This suggests that even nonclassical states admit a decomposition as a linear combination of coherent states. The decomposition of any mixed state as a linear combination of projectors onto coherent states is known as the Glauber-Sudarshan decomposition. However, the linear decomposition is not a statistical mixture over a bona fide probability density, and needs to be interpreted in a weak sense. Informally, physicists write

ρ = P ρ (α) |α α| dα, (122)
with P ρ (α) a tempered distribution for each α; see [Mandel and Wolf, 1995, Section 11.8] and references therein. For a general state, there is no reason for P ρ to be a positive function, or even a function at all. We do not attempt to formalize the Glauber-Sudarshan decomposition here, but we study the simplest example of a nonclassical state: a Fock state.

Consider a single-mode Fock state

ρ = |k k| ⊗ I ⊗ . . . .
The associated Glauber-Sudarshan distribution, denoted for simplicity P |k , is the "derivative of a Dirac delta", namely

P |k (α) = e |α| 2 k! ∂ 2k ∂ k α∂ k α δ(α) ; (123)
see [Mandel and Wolf, 1995, Section 11.8]. This representation of a Fock state is very singular and comes with no simple statistical interpretation. To get some intuition about Equation ( 123), let us fix two arbitrary states |u and |v . We have in fact the matrix elements

u|ρ|v = C 2 ∂ 2k ∂ k α∂ k α e |α| 2 k! u|α α|v δ(α) dα (124) = ∂ 2k ∂ k α∂ k α e |α| 2 k! u|α α|v α=α=0 , ( 125
)
as long as α → u|α α|v is smooth enough. This indeed corresponds to the informal (122), when computing matrix elements and with P ρ a specific tempered distribution. In particular, if |m and |n are two single-mode Fock states, and remembering from Section 3.6

that n|α = e -|α| 2 2 α n / √ n!, we verify that (125) yields n|ρ|m = ∂ 2k ∂ k α∂ k α e |α| 2 k! e -|α| 2 α n α m √ n! √ m! α=α=0 = δ kn δ km
which is 1 if n = m = k and otherwise 0, as expected.

As for the coherence functions themselves, we can use this representation to give a form very similar to the statistical mixture we discussed above. Indeed, by linearity, we can write

G (n) ρ (x 1 , . . . , x 2n ) = P ρ (α) n i=1 E + α (x i )E + α (x n+i ) dα . ( 126
)
As an example, consider again a single-mode state in a Fock state with k photons. The coherence function can then be obtained directly from Equations ( 123) and ( 126)

G (n) |k (x 1 , . . . , x 2n ) = ∂ 2k ∂ k α∂ k α e |α| 2 k! n i=1 E + α (x i )E + α (x n+i ) α=α=0 . ( 127
)
Example 10 (A candidate point process corresponding to a Fock state.). Taking again x i = x n+i and a broad-band detector, we obtain candidate correlation functions to define a point process for a single-mode Fock state with k photons

ρ n (t 1 , . . . , t n ) = ∂ 2k ∂ k α∂ k α e |α| 2 k! |E + α (t 1 ) • • • E + α (t n )| 2 α=α=0 . ( 128
)
Explicitly, we have

E + α (t) = N αe -iωt . ( 129
)
Then (128) implies ρ n (t 1 , . . . , t n ) = 0 for k < n. For k ≥ n, we obtain

ρ n (t 1 , . . . , t n ) = N n ∂ 2k ∂ k α∂ k α e |α| 2 k! α n α n α=α=0 (130) = N n k! ∂ 2k ∂ k α∂ k α (αα) k (n -k)! α=α=0 (131) = N n k! (k -n)! . ( 132
)
We find an expression reminiscent of the correlation functions of k independent and uniformly distributed times [Johansson, 2006, Example 2.6]. However, these tentative correlation functions are not integrable, a problem linked with the fact that the plane waves in (83) or (98), while allowing simple formal derivations and a clear separation of positive and negative frequencies, are not actually in L 2 . In fact, it is possible to build a quantum field theory of radiation replacing plane waves by any orthonormal set of L 2 solutions ϕ (r, t) of the classical equation of motion [START_REF] Fabre | Modes and states in quantum optics[END_REF]. The positive part of the field can then be decomposed on these modes as r, t), where the + subscript still corresponds to the analytic signal 5 taken in the time variable. The whole framework introduced in Section 3 follows, with ladder operators and coherent states now relative to this new basis. Similarly, Equation ( 128) stills holds with α corresponding to the coherent state parameters in the chosen modes. However, in the current example of a single-mode Fock state with k photons, the electric field ( 129) becomes E + α (t) = N αϕ + (t), where we have dropped the dependence on r as in (129) since the detector remains at a fixed position. We thus obtain, in lieu of (133),

E + α (r, t) = α ϕ + (
ρ n (t 1 , . . . , t n ) = N n k! (k -n)! |ϕ + (t 1 )| 2 . . . |ϕ + (t n )| 2 , n ≥ 1. ( 133
) Since ϕ + ∈ L 2 ,

Single-photon sources can lead to anti-bunching

The fact that photons emitted by a Gaussian classical field yield bunched detection times is known as the Hanbury-Brown & Twiss (HBT) effect, and was first evidenced by Hanbury [START_REF] Brown | Interferometry of the intensity fluctuations in light. II. an experimental test of the theory for partially coherent light[END_REF]. While we shall come back to HBT-type experiments in Part II, we pause here to insist on the fact that bunching is a consequence of considering multiple independent sources. It is perfectly possible to obtain point processes that anti-bunch with photons, in the sense that the resulting second correlation function ρ 2 (x, y) of the point process is small close to the diagonal {(x, y) : x = y}.

To gain intuition without the quantum overhead, we first recast the HBT effect as a consequence of interference between classical sources. Consider S independent 34 classical sources, represented by S independent, zero-mean stochastic processes E s (r, t). The total field is E(r, t) = s E s (r, t).

As derived in Section 4.1, the first-order coherence function is

G (1) (r t , r t) = E s E - s (r , t ) s E + s (r, t) . ( 134
)
By assumption, all cross-terms in the product have expectation zero, so that

G (1) (r t , r t) = s E E - s (r , t )E + s (r, t) := s G (1) s (r t , r t) . ( 135
)
The total first-order coherence is the sum of the first-order coherence of each source. Quite generically, the sum will tend to zero when |r -r | or |t -t | are large compared to the so-called spatial coherence scale and correlation time , respectively. We can also compute physical quantities like the average field intensity, which, in this simple model, is the sum of the average field intensities of each source,

E I(r, t) = E E * (r, t)E(r, t) = G (1) (r t, r t) = s E (|E s (r, t)| 2 ) = s E I s (r, t) . ( 136
)
34 For instance, we neglect interatomic coupling through the radiation for atomic sources.

To go further and be able to discuss bunching or anti-bunching, we need to consider the second-order coherence function, which characterizes the density of pairs of detection events; see (118). It writes

G (2) (r 1 t 1 , r 2 t 2 , r 2 t 2 , r 1 t 1 ) (137) = E s E - s (r 1 , t 1 ) s E - s (r 2 , t 2 ) s E + s (r 2 , 2 t) s E + s (r 1 , t 1 ) . ( 138
)
Using our assumptions of independence and zero mean, this coherence function is a sum of one-source and two-source terms. We focus on two-source terms. Developing (138), the contribution of two-source terms to EI(r, t)I(r , t ) is

s =s E (|E + s (r, t)E + s (r , t ) + E + s (r , t )E + s (r, t)| 2 ) . ( 139
)
At coincidence (r, t) = (r , t ), both amplitudes are the same and give a contribution 4|E + s (r, t)E + s (r, t)| 2 . However, when (r, t) and (r , t ) differ, a phase appears between the two contributions and reduces the modulus of the sum. We recover here the bunching effect of light as an interference effect between pairs of complex amplitudes. Note that S = 2 sources are enough to create this interference.

If we want to observe the opposite effect, a reduction in the rate at coincidence, we need a single-photon source such as a single two-level atom described by a Hilbert space C 2 excited by a pumping laser field. A natural analogy is that of a gun (the atom) that needs to be reloaded (by the laser pump) before firing a new bullet (the photon). We shall see that the need to reload implies anti-bunching for the fired photons.

What we want to compute is once again the second-order coherence between two detections at the same photodetection location r, but delayed by a time τ , and abusively denoted by G

(2) ρ (r, r, τ ). The whole state space of the source, the field and the detector is then C 2 ⊗ H ∞ ⊗ C 2 . An accurate model of the situation would require a proper description of the interaction between the field, the source and the detector but we only give some elements of the general computation; see [START_REF] Cohen-Tannoudji | Emission et détection de rayonnement : approches semi-classiques et approches quantiques[END_REF] for details. The new element here is the interaction of the field and the source and the main technical problem is to relate the field E + (t) emitted by the source to some of its physical properties. We ask the reader to admit that

E + (t) ∝ σ -(t), where σ -(t) = U † (t)(σ -(0) ⊗ I ⊗ I)U (t).
Physically, this expression means that when the source is excited and then relaxes to its ground state, it emits light. The second-order coherence of the field can then be written as a function of the time evolution of a source operator as

G (2) ρ (r, r, τ ) ∝ σ + (0)σ + (τ )σ -(τ )σ -(0) ρ ∝ σ + (0)U † (τ )σ + (0)U (τ )U † (τ )σ -(0)U (τ )σ -(0) ρ . ( 140
)
Naturally, we assume that the initial state of the atom is excited and the field is in the vacuum so that ρ = |e e| ⊗ |0 0|. With σ -(0) = |g e| we end up with the very simple form:

G (2) ρ (r, r, τ ) ∝ Tr |g g| ⊗ |0 0| U † (τ ) |e e| U (τ ) ∝ e|U (τ ) |g g| ⊗ |0 0| U † (τ )|e . ( 141
)
Two comments are in order. First, we could have also assumed a more general form for the atomic density matrix which would amount to multiply the above expression by the probability p(e) for the atom to be initially excited. Its appearance is natural since the atom has to be excited initially. Second, the physical interpretation of the second-order coherence ( 141) is as follows. Starting from a factorized state in which the atom has already emitted a photon and is thus in its ground state, the system evolves until time τ where we compute the probability for the atom to be excited again. Equation ( 141) is the conditional probability that, knowing that the atom is in its ground state initially, it gets excited at time τ . The oscillations at larger delay correspond to so-called Rabi oscillations between the two energy levels; as the pump power is increased these oscillations increase in frequency. The solid line is a fit to a model. Figures taken from [START_REF] Flagg | Resonantly driven coherent oscillations in a solid-state quantum emitter[END_REF].

From there, predicting photon anti-bunching at coincidence is straightforward: near coincidence τ = 0, the evolution operator is U (0) = 1. Since |e and |g are orthogonal, this implies

G (2) ρ (r, r, 0) = 0 . ( 142
)
In other words, at coincidence, the probability to have a second emission is zero. This result is in fact physically quite intuitive if we think at the single-source level. Indeed, after the atom emits one photon, it is surely in its ground state and cannot emit a second photon right away. It has to be re-excited first by the laser source before emitting the second photon. A typical experimental signal is shown in (3) where we clearly see the coherence signal going to zero at coincidence. The oscillation pattern is also nicely understood as Rabi oscillations between the two atomic levels [START_REF] Degiovanni | Physique quantique, information et calcul: des concepts aux applications[END_REF]. Nowadays, this anti-bunching is used as signature of a good single photon source, along with a second major interference effect called Hong-Ou Mandel; see e.g. [START_REF] Degiovanni | Physique quantique, information et calcul: des concepts aux applications[END_REF].

Electrodetection and fermionic coherences

While photonic optics experiments have reached a high level of sophistication [START_REF] Haroche | Exploring the quantum: atoms, cavities, and photons[END_REF]Raimond, 2006, Grynberg et al., 2010], it is only recently that electronic coherences have started to be probed, thanks to experimental advances in the manipulation, control and measurement of small quantum systems in condensed matter and atomic physics. Formally, the main difference with bosons comes from the anti-commutation of the associated ladder operators, resulting in the Pauli exclusion principle and in determinants appearing in Wick's Theorem 3. Moreover, while we limit ourselves to free (i.e., non-interacting) fermions in this section, in realistic physical models fermions are usually subject to the effects of interaction. These two aspects drastically change the physics compared to photons. Photons are noninteracting objects with a truly empty reference state, the vacuum. Electrons, on the contrary, are subject to Coulomb interaction and their ground state, called the Fermi sea and representing a metal at equilibrium, is full of fermions: their vacuum is not empty. This Section follows the lines of Section 4.2, explaining how similar computations are affected by the fermionic character of the particles. In Section 5.1, we remind what a fermionic field is and discuss how to model the detection of electrons. In Section 5.2, we discuss the coherence functions of arbitrary orders, and use the second-order one in Section 5.3 to justify the default anti-bunching character of fermions. In Section 5.4, we finally discuss the difficulties in recovering classical currents in the way coherent states do for bosonic fields.

Modeling electrodetection events

We proceed as for photodetection by first modeling the field and then the detection, so as to be able to write correlation functions.

Modeling the fermionic field

A fermionic quantum field ψ(r, t) is a collection of operators indexed by space and time and acting on a Fock space H ∞ Fermions of fermions; see Section 3.4. The field models an indefinite number of fermionic particles, like electrons. Similarly to a bosonic field, the state of a fermionic system is built by acting on a reference state with operators that create or annihilate fermions in a given state of a basis of the single-particle Hilbert space. When it is clear that we speak of fermions,35 these ladder operators are usually denoted as c † and c, and indexed by the label of the single-particle state. Unlike for bosons, the creation and annihilation operators satisfy the canonical anticommutation relations (48).

In theory, the reference state can be taken to be the vacuum, meaning the complete absence of any particle, as we did for bosons. But when modeling electronic experiments, physicists often have to take another reference state, labeled as |F , and called the Fermi sea. The Fermi sea is meant to represent, for instance, the metal used in experiments, which is itself full of electrons. In the simplest case, corresponding to temperature zero, the Fermi sea is a state with one fermion in every energy level of the considered system, up to some reference level called the Fermi level. For instance, in a system built using a single-particle Hamiltonian with eigenpairs (ϕ n , ε n ) n∈N , so that c † n creates a particle in mode ϕ n , the Fermi sea is

|F = n:εn<ε F c † n |0 ,
with the Fermi level ε F defined by the experimental setting. Like the vacuum, the Fermi sea can also be built as the limit of a Gaussian state when β → ∞ and the chemical potential ζ is fixed as a function of ε F ; see the computations in Section 6.

Remark 8. Depending on the situation to model, more physical parameters can be introduced to describe the field in addition to the spacetime coordinates. For instance, the spin of the fermionic excitation could contribute, or more fermionic fields can be used to describe different channels, as in the description of the quantum Hall effect [START_REF] Ezawa | Quantum Hall effects: Field theoretical approach and related topics[END_REF].

First-order electronic coherence

As we did for photons, we consider a two-level detector with Hilbert space C 2 probing a fermionic field ψ(r) living in the Hilbert space H ∞ Fermions . 36 The dynamics of the detector is described here by a fermionic operator d that satisfies the anticommutation rule {d, d † } = 0. Their interaction, again assumed to be weak, is modeled by Hamiltonian (in the interaction picture)

H I (t) = d † (t)ψ(r, t) + d(t)ψ † (r, t) , ( 143
)
Note that we place ourselves in the interaction picture, so that the time dependence in (143) results from the evolution of the free part of the Hamiltonian. The first term of the interaction Hamiltonian describes the absorption process of an electron that excites the detector while the second term describes an electron being emitted by the detector. With the same assumptions as in the photodetection problem, we can compute the excitation probability of the detector as a function of an electronic correlation function, naturally called in this context the first-order electronic coherence function, and defined as

G (1e) ρ (r t , r t) = ψ † (r , t )ψ(r, t) ρ . ( 144
)
Remark 9. In concrete condensed matter systems, we argued above that the natural ground state is in fact a Fermi sea |F , contrary to the photonic case where the natural ground state is the true vacuum |0 . In such a context, it is then also interesting to consider the de-excitation probability of the detector, sending back the electron into the system, or equivalently creating a hole in the Fermi sea. This transition probability is controlled by the first-order hole coherence function

G (1h) ρ (r t , r t) = ψ(r , t )ψ † (r, t) ρ . ( 145
)
Note however that the two types of coherence functions are not independent. Indeed, thanks to the canonical anticommutation relations, they satisfy the equal time relation

G (1e) ρ (r t, r t) + G (1h) ρ (r t, r t) = δ(r -r).
In the following, we work with one of the two coherence functions, and we pick G Example 11 (Free fermions at non-zero temperature T ). The simplest example we can consider is a set of free fermions. In particular, we consider the free field Hamiltonian H FF = k ν k c † k c k , and further assume that the energy level ν k = ν k only depends on k = k . We prepare the state in the grand-canonical ensemble (60), which we denote here as ρ = ρ β,ζ . Note that ρ commutes with H FF , so that Schrödinger's equation yields that ρ does not change in time: it is an equilibrium state. In particular, the coherence function G (1e) ρ (r t , rt) depends only on t -t . We arbitrarily set t = t , and write G (1e) ρ (r , r). Introducing the Fourier transform of the annihilation operator

ψ(r) = R 3 c k e i(k•r) dk/(2π) 3 ,
we can write the first electronic coherence function as

G (1e) ρ (r , r) = ψ † (r )ψ(r) ρ = R 6 c † k c k ρ e i(r•k-r •k ) dk (2π) 3 dk (2π) 3 . ( 146
) Now note that 37 c † k c k ρ = f β,ζ (k)δ(k -k ) with k = k where f β,ζ is the Fermi-Dirac distribution f β,ζ (k) = 1 e β(ν k -ζ) + 1 . ( 147
)
We obtain

G (1e) ρ (r , r) = R 3 f β,ζ (k) e i(r-r )•k dk (2π) 3 . ( 148
)
A practically relevant case is the zero-temperature limit, where the grand-canonical ensemble state is replaced by the Fermi sea itself |F F |. The Fermi-Dirac distribution then simplifies to the indicator 1 B(k F ) of the centered ball B(k F ) of radius k F known as the Fermi momentum, which is an input of the model38 . The first-order electronic coherence function is then

G (1e) |F (r , r) = 1 π 2 r -r k F 0 k sin(k r -r ) dk = 1 π 2 r -r sin(k F r -r ) r -r 2 - k F r -r cos(k F r -r ) . ( 149
)
Example 12 (Chiral free fermions at non-zero temperature T ). A second example appears naturally when discussing electronic coherences. This time, we consider a simple model of a chiral quantum wire, that is, electrons moving in one dimension and only in a given direction (e.g. to the right). The Hilbert space is similar to Example 11, C 2 × H ∞ Fermions , with H = L 2 (R) as a single-particle Hilbert space. Rather than describing the Hamiltonian, a typical shortcut to describe the model is to directly discuss the solutions to Schrödinger's equation, or, as we do here for the field operators, how operators evolve in the Heisenberg picture.

Compared to Example 11, we are not at equilibrium, and thus need to keep track of the time dependence of the field ψ(x, t). The Fourier variable conjugated to the time t is denoted by ω, just like the momentum k is conjugated to the position r. Being in one dimension in this example, the vectors r and k are in fact scalar quantities that we will denote x and k respectively. To further constrain the model, we posit an equation of motion for the field, i.e., a relation between its time-and space derivatives39 . The non-zero solutions of this equation are described by a relation between all the Fourier variables, called a dispersion relation. In our quantum wire model, the dispersion relation is assumed to be linear, meaning that the eigenvalues of the free field Hamiltonian are given by ν k = v F |k|, with v F the so-called Fermi velocity of the electron in the quantum wire. Finally, we again assume that the electrons coming into the wire come from a reservoir of particles, whose state is the grand-canonical ensemble ρ = ρ β,ζ at temperature β = 1/k B T and chemical potential ζ.

Now that we have all the additional elements for this model, we can proceed to compute the first-order coherence function. Assuming that the field lives in a large box, from the Fourier decomposition ψ(x, t) = k c k e i(kx-ν k t) , the first-order coherence expands as

G (1e) ρ (x t , x t) = ψ † (x , t )ψ(x, t) ρ = k,k c † k c k ρ e i(kx-k x -(ν k t-ν k t )) . ( 150
)
The average value of c † k c k on the grand canonical state ρ can be shown to be

c † k c k = δ(k -k )f β,ζ (k).
In short, the first-order electronic coherence is simply a Fourier transform of the Fermi-Dirac distribution. Doing the change of variable to ω, we have As it is traditionally done in this field, we will use a rescaled time variable t -x/v F as t. This simplification between time and space is very special to the one-dimensional case thanks to the ballistic propagation of the electron. We then write the simpler expression

G (1e) ρ (x t , x t) = R 1 e β( ω-ζ) + 1 e i(x-v F t-(x -v F t )) ω v F dω 2πv F . ( 151
G (1e) ρ (t , t) = R 1 e β( ω-ζ) + 1 e i(t -t)ω dω 2πv F . ( 152
)
As such, the integral is divergent when t = t and has to be regularized. This divergence has a physical origin: the linear dispersion relation posits a Fermi sea with an infinite number of fermions with negative energy. In practice however, the system has a finite size and the Fermi sea has a finite bandwidth. The strategy is then to introduce a small parameter ε > 0, which corresponds to the inverse of that bandwidth, and compute the integral with the substitutiont -t → t -t + iε. With these modifications, the residue theorem gives

G (1e) F (t , t) = i 2πv F τ th e -i ζ (t-t ) sinh t-t +iε τ th , ( 153
)
where τ th = β/π is called the thermal coherence time. A particular case to have in mind is the zero-temperature limit of the coherence function.

F (t , t) = i 2πv F e -i ζ (t-t ) t -t + iε . ( 155
)

Correlation between electrodetection events

The first-order coherence gives information about the system at the one-particle level. If we want to probe two-particle or higher levels of information like the correlations between electrodetection events, we have to study the coherence functions of higher order. They are defined similarly as the photonic coherences of Section 4. For instance, the second-order electronic and hole coherence functions are defined as

G (2e) ρ (r 1 t 1 , r 2 t 2 , r 2 t 2 , r 1 t 1 ) = ψ † (r 1 , t 1 )ψ † (r 2 , t 2 )ψ(r 2 , t 2 )ψ(r 1 , t 1 ) ρ (156a) G (2h) ρ (r 1 t 1 , r 2 t 2 , r 2 t 2 , r 1 t 1 ) = ψ(r 1 , t 1 )ψ(r 2 , t 2 )ψ † (r 2 , t 2 )ψ † (r 1 , t 1 ) ρ . ( 156b
)
40 The +iε notation is a shorthand notation for the more rigorous distribution

lim ε→0 1 x + iε = P 1 x -iπδ(x) . ( 154
)
This straightforwardly generalizes to higher orders. Using the shorthand notation x = (r, t) for the spatio-temporal coordinates, the n-th-order electronic and hole coherences are

G (ne) ρ (x 1 , • • • , x n , x n , • • • , x 1 ) = ψ † (x 1 ) • • • ψ † (x n )ψ(x n ) • • • ψ(x 1 ) ρ (157a) G (nh) ρ (x 1 , • • • , x n , x n , • • • , x 1 ) = ψ(x 1 ) • • • ψ(x n )ψ † (x n ) • • • ψ † (x 1 ) ρ . (157b)
Example 13 (Continuation of Example 11). The grand canonical ensemble ρ = ρ β,ζ for free fermions is Gaussian, so that Theorem 3 (Wick's) applies and we obtain a 2 × 2 determinant

G (2e) ρ (r 1 , r 2 , r 2 , r 2 ) = ψ † (r 1 )ψ † (r 2 )ψ(r 2 )ψ(r 1 ) ρ (158a) = Wick G (1e) ρ (r 1 , r 1 )G (1e) ρ (r 1 , r 1 ) -G (1e) ρ (r 1 , r 2 )G (1e) ρ (r 2 , r 1 ) . ( 158b
)
The minus sign is of course a consequence of the fermionic character of the electronic field.

In fact, we can already see at this level the repulsion between different fermions. Indeed, if we look for instance at the diagonal part r 1 = r 1 and r 2 = r 2 , in the zero temperature limit, we have

G (2e) |F (r 1 , r 2 , r 2 , r 1 ) ∝ 1 -g (1e) |F ( r 1 -r 2 ) 2 , ( 159
)
where g (1e)

|F ( r 1 -r 2 ) is the first-order electronic coherence function (149) normalized to unity. We clearly see that when r 1 = r 2 , the second-order coherence function is exactly zero: two fermions cannot be at the same place at the same time. This is another avatar of the general anti-bunching effect of fermionic systems.

Example 14 (Second-order coherence of N excitations). To gain some intuitive understanding of the information content of the second-order coherence function, let's continue Example 12 of a quantum wire. Remember that in this simple one-dimensional model, time and space coordinates are identified thanks to the linear dispersion relation and that we are looking for the coherence functions at the fixed position of the detector (position that we omit in the arguments to simplify the notations). Consider a state containing two electrons above the true vacuum |0 (not the Fermi sea) in two orthogonal wavefunctions

|ϕ 1 , ϕ 2 = ψ † [ϕ 1 ]ψ † [ϕ 2 ] |0
, where ψ[ϕ] = ψ(t)ϕ(t)dt. Note again that wavefunctions, strictly speaking functions of x and t, are here functions of x -v F t only because of the linear dispersion relation. Since we further consider a fixed position here (that of the detector), we simply write t → ψ(t) and t → ϕ(t).

Since the true vacuum is a Gaussian state, we can apply Wick's Theorem 3 to the second-order coherence function

G (2e) |ϕ1,ϕ2 (t 1 , t 2 , t 2 , t 1 ) = R 4 0|ψ(y 1 )ψ(y 2 )ψ † (t 1 )ψ † (t 2 )ψ(t 2 )ψ(t 1 )ψ † (y 2 )ψ † (y 1 )|0 ϕ * 1 (y 1 )ϕ * 2 (y 2 )ϕ 2 (y 2 )ϕ 1 (y 1 ) dy 1 dy 1 dy 2 dy 2 = ϕ * 1 (t 1 )ϕ * 2 (t 2 )ϕ 1 (t 1 )ϕ 2 (t 2 ) + ϕ * 1 (t 2 )ϕ * 2 (t 1 )ϕ 1 (t 2 )ϕ 2 (t 1 ) -ϕ * 1 (t 2 )ϕ * 2 (t 1 )ϕ 1 (t 1 )ϕ 2 (t 2 ) -ϕ * 1 (t 1 )ϕ * 2 (t 2 )ϕ 1 (t 2 )ϕ 2 (t 1 ) . ( 160
)
We see that the form of this coherence contains the expected anti-symmetries coming with fermionic statistics. It can be written even more explicitly by introducing the antisymmmetrized wave-function

Φ 12 (t 1 , t 2 ) = det ϕ 1 (t 1 ) ϕ 2 (t 1 ) ϕ 1 (t 2 ) ϕ 2 (t 2 ) . ( 161
)
The second-order electronic coherence function then reads showing the structure of the electronic coherence functions of order up to N in the presence of one up to N fermions in some given wavepackets.

G (2e) |ϕ1,ϕ2 (t 1 , t 2 , t 2 , t 1 ) = Φ * 12 (t 1 , t 2 )Φ 12 (t 1 , t 2 ) . ( 162
) ϕ e1 • • • ϕ e k • • • ϕ e N G (1e) ϕ * e1 ϕ e1 • • • k p=1 ϕ * ep ϕ ep • • • N p=1 ϕ * ep ϕ ep . . . . . . . . . . . . G (ke) 0 Φ * 1,...,k Φ 1,...,k i1<•••<i k Φ * i1,...,i k Φ i1,...,i k . . . . . . . . . . . . G (N e) 0 • • • 0 • • • Φ * 1,...,N Φ 1,...,N Figure
The second-order coherence function is thus essentially the full many-body wavefunction of a two-particle state.

When more than two fermions are present, the second-order coherence function only keeps track of two-particle terms. To see this, consider a state containing N electrons above the vacuum in mutually orthogonal wavepackets

|ψ N = ψ † [ϕ 1 ] . . . ψ † [ϕ N ] |0 . A straightforward, but cumbersome, application of Wick's theorem gives G (2e) |ψ N (t 1 , t 2 , t 2 , t 1 ) = k<l Φ * kl (t 1 , t 2 )Φ kl (t 1 , t 2 ) . ( 163
)
Informally, the second-order coherence function contains only the "two-particle physics".

Example 15 (kth-order coherence of N excitations). Continuing the quantum wire example, we can understand its k-particle physics by computing the kth-order coherence function. Consider again N fermions prepared in mutually orthogonal wavepackets

|ψ N = ψ † [ϕ 1 ] . . . ψ † [ϕ N ] |0
. By Wick's theorem, the kth-order coherence function is

G (ke) |ψ N (t 1 , • • • , t k , t 1 , • • • , t k ) = i1<•••<i k Φ * i1,...,i k (t 1 , • • • , t k )Φ i1,...,i k (t 1 , • • • , t k ) , ( 164
) with Φ i1,...,i k (t 1 , • • • , t k ) = det    ϕ i1 (t 1 ) • • • ϕ i k (t 1 ) . . . . . . . . . ϕ i1 (t k ) • • • ϕ i k (t k )    . ( 165
)
In particular, the N -th coherence function is

G (N e) |ψ N (t 1 , • • • , t N , t 1 , • • • , t N ) = Φ * 1,...,N (t 1 , • • • , t N )Φ 1,...,N (t 1 , • • • , t N ) . ( 166
)
As expected, the N -th coherence function is essentially the full many-body wavefunction of a quantum state with N excitations in mutually orthogonal wavefunctions.

To see the correlation functions of determinantal point processes appear, one can use an ideal broad-band detector as in Section 5. In particular, imagine that we could prepare N electrons in mutually orthogonal wavepackets ϕ 1 , . . . , ϕ N , and set up an ideal detector as we did for photons in Section 4, resulting in a point process with correlation functions given by the diagonal of G

(ne) |ψ N , ρ n (t 1 , . . . , t n ) = G (ne) |ψ N (t 1 , • • • , t n , t 1 , • • • , t n ), (167)
= |Φ 1,...,N (t 1 , • • • , t N )| 2 (168) = det [K(t i , t j )] n i,j=1 , (169)
where

K(t, s) = N i=1 ϕ ei (t)ϕ * ei (s). (170)
We recognize a DPP with (projection) kernel K, as introduced in Section 2.4.

Remark 10.

To illustrate what reference computations different scientists may have in mind, note that, following a physics viewpoint, we computed in Examples 14 and 15 coherence functions applying Wick's theorem, one after the other, starting from the lower ones.

Probabilists and statisticians are more used to derive correlation functions starting from the definition of the point process. Had we started by defining the DPP with kernel K in (170), and asked for the diagonal of, say, the second coherence function in ( 14), it would have been enough to notice, from the definition (184) of the correlation functions of a DPP, that

ρ 2 (t 1 , t 2 ) = det K(t 1 , t 1 ) K(t 1 , t 2 ) K(t 2 , t 1 ) K(t 2 , t 2 ) . ( 171
) By definition (170) of K, ρ 2 (t 1 , t 2 ) = det ϕ 1 (t 1 ) . . . ϕ N (t 1 ) ϕ 1 (t 2 ) . . . ϕ N (t N )    ϕ * 1 (t 1 ) ϕ * 1 (t 2 ) . . . . . . ϕ * N (t 1 ) ϕ * N (t 2 )    . ( 172
)
The diagonal version of (163) then results from the Cauchy-Binet formula, which allows rewriting ρ 2 (t 1 , t 2 ) as a sum of 2 × 2 determinants; see e.g. [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF].

Electron anti-bunching

For photons, we saw that HBT-type experiments can reveal both a bunching effect (for classical light beams) and an anti-bunching effect (for quantum beams), depending on the source. A similar experimental setup as the one pictured for photons in Section 4 can be imagined for fermionic excitations. We defer the discussions on how to properly do this experimentally with fermionic atoms or single electronic excitations to Part II, but we can already foresee from either DPP constructions like the one with kernel (170), or directly from the coherence functions of Section 5.2 that fermions will exhibit anti-bunching detection times. Unlike bosons, this is a direct consequence of the anti-symmetric statistic of fermions, and does not qualitatively change when introducing source models.

If the arrival times follow a DPP with an Hermitian kernel like (170), we saw indeed in Section 2.4 that samples tend to spread regularly, and form less clusters than a Poisson point process. Alternately, this statistical anti-bunching effect can be directly seen from the coherence funtions. Indeed, as for photons, the HBT signal is controlled by the second-order coherence function at the position r of the detector,

G (2e) ρ (r, t 2 , t 1 ) = ψ † (r, t 1 )ψ † (r, t 2 )ψ(r, t 2 )ψ(r, t 1 ) ρ . ( 173
)
At coincidence τ = t 2 -t 1 = 0, because the CAR (48) impose ψ(r, t 1 ) 2 = 0, we obtain

G (2e) ρ (r, t, t) = 0 . (174) If G (2e) ρ
is smooth, then we expect few coincidences at small time differences. As a final remark, we insist that bunching or anti-bunching of quantum excitations is not a signature of the statistics of the elementary excitations. While it is a consequence of the statistics for fermions, bosons also anti-bunch when we add to the model a single-excitation source. Hence, it is better to think of anti-bunching as a signature of the quantum nature of the fundamental excitation.

Recovering a classical current is not as easy as for bosons

For photons, coherent states allowed for relating the quantum formalism to classical fields; see Section 4. Describing a classical electric current or, more generally, a classical theory of fermionic fields, is more difficult. In this section, we examine reasons why it is not possible to build a set of states having all the properties that make bosonic coherent states handy. First, even considering a single mode, bosonic coherent states are superpositions of arbitrary numbers of bosons, which allows for having very small fluctuations around the average boson number. Indeed the number of photons in the mode is a Poisson variable with arbitrarily large parameter, and thus arbitrarily small relative variance. Such small fluctuations are not achievable with fermions, because it is impossible to have more than one fermion per mode.

A second difference is that the eigenvalue of a tentative fermionic coherent state of the fermionic annihilation operator cannot be a complex number. Indeed, were |α an eigenstate of the fermionic annihilation operator a, then the CAR (48) would imply that its eigenvalue α satifies α 2 = 0. Moreover, in the many-mode case, a two-mode coherent state |α 1 , α 2 would satisfy the anti-commutation relation a 41 The eigenvalues of a fermionic coherent state thus cannot be complex numbers. 42A more subtle difference is that fermionic fields cannot have an observable average amplitude,43 unlike bosonic fields prepared in a coherent state.44 This impossibility is fundamental and is related to superselection rules [START_REF] Wick | The intrinsic parity of elementary particles[END_REF]. For concreteness, consider a general fermionic state |ϕ = α |0 + β |1 with |1 = ψ † |0 , where ψ † is any fermionic creation operator. Now, physicists usually impose symmetries on their models, which implies that some states should yield the same measurements. In particular, one would like to model the fact that rotating the physical system around any axis by an angle of 2π should not modify the law of measurements. For fermions, the so-called spin-statistics theorem of relativistic quantum field theory tells us that this kind of rotations should be modelled by the action of the half-integer representations of the group SU (2). For instance, in the representation of order 1/2, a rotation of angle 2π corresponds to multiplying the state by e iπN , where N is the number operator. Applied to |ϕ = α |0 + β |1 , this yields |ϕ := e iπN |ϕ = α |0 -β |1 , where the minus sign comes from the fact that 1 is an odd number. However, since such a rotation should not correspond to any change of the system, we impose that |ϕ = |ϕ up to a phase, and consequently either α or β is zero. Thus |ϕ is in fact not a superposition, implying that ψ |ϕ = 0. More generally, the parity superselection rule imposes that a given state can only consist of superpositions of states with the same parity of the number of fermions. Therefore, a single fermionic annihilation operator can only have a vanishing average. 45For all these reasons, defining a classical regime with fermionic fields is an open question. This actually raises deep conceptual questions when fermions are involved. For instance, it is not clear what is meant by a classical electronic current if we start from a pure quantum description of the electronic current.

1 a 2 |α 1 , α 2 = α 1 α 2 |α 1 , α 2 = -a 2 a 1 |α 1 , α 2 = -α 2 α 1 |α 1 , α 2 .

Wrapping up and open questions

To conclude this Part I, we first show how to use the formalism of Sections 2 to 5 to write a generic permanental or determinantal point process as the detection process of a system of free bosons or fermions, respectively. Second, we propose a list of questions, from simple ones to open problems, that are raised by the connections between point processes and physical measurements in quantum field theory. Third, we announce the sections that are to appear in Part II of this document.

From a point process to free particles

In Sections 4 and 5, we obtained some point processes, including permanental and determinantal point processes, by modeling the detection of physical particles. A natural question is whether all permanental and determinantal point processes arise in this form. We now show that, at least formally and under weak assumptions on the kernel, the answer is yes. More precisely, let µ be a Borel measure on a complete metric space X, and (ϕ i ) an orthonormal family in L 2 (µ). Consider the kernel

K(x, y) := i λ i ϕ i (x)ϕ i (y), (175)
where, for all i, λ i ∈ [0, 1] if η = -1, and λ i ∈ R + if η = +1. We now build a Fock space, a quantum state, and a measurement model that together lead to a permanental or determinantal point process with K as its kernel. Let η = ±1 depending on whether one wants a permanental or a determinantal point process.

Building a Fock space. Build the Fock space like in Section 3.4, using the basis (ϕ i ) as single-particle basis and the symmetrization property dictated by η. Consider the operator

H = i ν i a † i a i ,
where (ν i ) is left as a free parameter for now, and will later be chosen in relation to the spectrum of K. We think of H as a Hamiltonian.

A Gaussian density matrix. Let ζ > 0, and define the so-called grand canonical ensemble as in (60) by the Gaussian density matrix

ρ = e -β(H-ζ i a † i ai) Z η = 1 Z η e -β i (νi-ζ)a † i ai .
We explicitly write the dependence of Z to η, as bosons and fermions lead to different normalization constants. More precisely, since the Fock states |n = |n 1 , n 2 , . . . are eigenvectors of all the number operators a † i a i in the exponential, the normalization constant is

Z η = Tr e -β i (νi-ζ)a † i ai = n n|e -β i (νi-ζ)a † i ai |n = n e -β i (νi-ζ)ni . ( 176
)
For bosons, we sum over all sequences of integers with only p non-zero components, and this for all p ≥ 1. In particular,

Z 1 = p∈N n∈N e -β(νp-ζ)n = p∈N 1 1 -e -β(νp-ζ) , ( 177
)
where we implicitly assumed ν p -ζ > 0 for all p, for the geometric sums in (177) to converge. 46 For Z 1 to be finite, we further need a condition on the (ν p ), such as

p log 1 e β(νp-ζ) -1 < ∞. (178)
For fermions, we only sum in (176) over all sequences in {0, 1} N that have p non-zeros components, for all sizes p ≥ 1, so that

Z -1 = p∈N n∈{0,1} e -β(νp-ζ)n = p∈N 1 + e -β(νp-ζ) . ( 179
)
Guaranteeing convergence is easier than for bosons, e.g. assuming that p (ν p -ζ) converges. Now that we have simple expressions for Z η , η = ±1, we can use the partition function trick 47 to compute the expected number of particles in mode i,

a † i a i ρ = 1 Z η Tr a † i a i ρ (180a) = n n i e -β j (νj -ζ)nj n e -β j (νj -ζ)nj (180b) = - 1 β ∂ ∂ν i log Z η (ν). ( 180c
)
Using ( 177) and ( 179) in turn, this yields

a † i a i ρ = 1 e β(νi-ζ) -η , ( 181
)
which is the so-called Fermi-Dirac distribution 48 for η = -1 and the Bose-Einstein distribution for η = 1. Note that along the lines of Lemma 2, one can show that a † i a j ρ = 0 if i = j.

A measurement basis. Now, consider the alternative basis |x , x ∈ X, with annihilation operator ψ(x) = i x, ϕ i a i ; see (52). We think of |x as the position basis, so that we write x, ϕ i = ϕ i (x), and we actually mean a generalized basis if necessary, e.g., if X = R d ; see Section 3.1. In the basis |x , the state ρ rewrites ρ = 1 Z exp -β i (ν i -ζ) ϕ i (x)ϕ i (y)ψ † (x)ψ(y)dµ(x)dµ(y) .

Wick's theorem applies. By Theorem 3, see also Example 5, for all n ≥ 1, we obtain the coherence functions

ψ † (x 1 )ψ(x 1 ) . . . ψ † (x n )ψ(x n ) = f η (K),
where K = ψ(x i ) † ψ(x j ) ρ 1≤i,j≤n and f η is the determinant if η = -1, and the permanent if η = 1. If the framework is consistent, an ideal detection experiment that measures all particles in ρ in the basis |x should thus correspond to a permanental/determinantal point process on X, with kernel

ψ † (x)ψ(y) ρ = i,j ϕ i (x)ϕ j (y) a † i a j ρ . ( 182
)
47 Statisticians might say the score function trick. 48 The word distribution is used here in a loose sense.

This is precisely K in (175), provided that we choose ν i so that

λ i = (e β(νi-ζ) -η) -1 . ( 183
)
This is achieved by setting

β(ν i -ζ) = log 1 + ηλ i λ i .
Note that when η = 1, one can check that the corresponding permanental point process exists by an adaptation of Macchi's Cox process construction of Section 2.3. When η = -1, the determinantal point process exists by the Macchi-Soshikov theorem; see Section 2.4.

Obtaining any kernel. We conclude with several comments. First, inverting (181) allows us to obtain any spectrum (λ i ) with λ i ∈ (0, 1) for fermions, and λ i > 0 for bosons. With the right assumptions, projection kernels for DPPs can also be obtained, but as a limit as β → +∞, i.e. when the temperature 1/β goes to zero. The corresponding projection operator is onto the span of the ϕ i s for which ν i < ζ. In other words, the chemical potential controls the rank of the limiting projection operator. Note that for DPPs, uniform convergence of the kernel on compact subsets of X implies convergence of the point process in a natural sense; see e.g. [Anderson et al., 2010, Section 4.2.8]. Another way to formally obtain projection DPPs is as a ground state; see Section 3.3.2. Second, the construction in this section, with a state ρ and measuring correlators, is as close as one can hope for a Cox process-like decomposition of determinantal point processes. Loosely speaking, there is a Gaussian object above determinantal point processes with self-adjoint kernels, but it is a Gaussian density matrix, i.e., an operator, not a functional process.49 Third, while we restricted here to Hermitian kernels, thus corresponding to an observable quantity, the Hermitian assumption does not play a major role in the mathematical framework. Furthermore, non-Hermitian free fermions may not be so removed from physics, see e.g. [START_REF] Ashida | Non-hermitian physics[END_REF] for a recent survey and [START_REF] Guo | Entanglement entropy of non-hermitian free fermions[END_REF] for a formal example.

On physical realizability. Finally, we have seen in Section 5 examples of experimentally realizable DPPs, in the sense that the corresponding idealized experimental setups can, in principle, be built in a physics laboratory. We will see in Part II concrete examples of DPP samples obtained from lab experiments. Yet, the generic construction given in Section 6.1 remains rather theoretical, as it is not clear how to build an experimental setup leading to any given DPP, that is, any choice of underlying space (X, µ) and any valid kernel K. One might be able to treat X = R d for d up to 3 by associating mathematical coordinates to space-coordinates. Going slightly beyond 3 might be doable using either time or additional physical degrees of freedom ("pseudo-dimensions") as the next dimensions, but it is unclear to physicists where fermions in arbitrary dimension would be needed as a model. Outside the problem of dimensionality, fixing the spectrum of the Hamiltonian to an arbitrary sequence can also be experimentally difficult.

Going further: from simple to open questions

Many natural questions come to mind when one looks back at the correspondence between permanental point processes and free bosons, or between DPPs and free fermions.

A quantum state is more than a point process

A point process is the result of observing a state ρ with a detector, but the state itself contains a lot more information. For instance, the correlation functions that we identified in Sections 4 and 5 are only the diagonals of the corresponding coherence functions. Relatedly, depending on the measurement basis, a single state ρ can yield many point processes.

As a concrete example, consider for simplicity fermions living in X = {1, . . . , N }, µ the counting measure on X, and the Hamiltonian

H = N i=1 ν i a † (ϕ i )a(ϕ i )
acting on the corresponding Fock space. By (182), measuring in an orthonormal basis (v i ) of C N , with V = ( ϕ i |v j ) i,j , we obtain the DPP with kernel matrix

K = k λ k ϕ k |v i v j |ϕ k 1≤i,j≤N = Vdiag(λ k )V † ,
where (λ i ) is given by ( 183). All Hermitian kernels with spectrum (λ k ) can thus be obtained from the same Hamiltonian, just changing the measurement basis.

At one extreme, measuring in the original basis used to build the Hamiltonian, i.e. taking V = I, leads to a diagonal kernel matrix. This means that the corresponding point process is a set of independent Bernoulli samples for fermions, and independent geometric samples for bosons [START_REF] Hough | Determinantal processes and independence[END_REF]]. Correlation appears when the observation basis (v k ) and the creation basis (ϕ k ) do not match. For instance, let v k = ϕ k for k = 1, . . . , N -2, but apply a special unitary transformation to the last two vectors. That is, for α, β ∈ C such that |α| 2 + |β| 2 = 1, complete the observation basis with

v N -1 = αϕ N -1 -βϕ N and v N = βϕ N -1 + αϕ N .
The resulting permanental or determinantal point process, call it γ α,β , has a block diagonal kernel

K α,β =        λ 1 . . . λ N -2 λ N -1 |α| 2 + λ N |β| 2 (λ N -1 -λ N )αβ (λ N -1 -λ N )αβ λ N -1 |β| 2 + λ N |α| 2        .
Note that since the diagonal of K α,β varies with (α, β), the first correlation function of the point process is also altered. For instance, the marginal probability that item N -1 belongs to the point process is λ N -1 |α| 2 +λ N |β| 2 . If λ N -1 ≥ λ N , then as soon as β = 0, the marginal probability of item N -1 occurring decreases. Note that for η = -1, though, the average number of points in γ α,β remains constant, as the trace of K α,β does not depend on (α, β). Meanwhile, the second correlation function ρ 2,θ (N -1, N ) is the permanent/determinant of the trailing block of K α,β , namely

ρ α,β 2 (N -1, N ) = |α| 2 |β| 2 λ 2 N -1 + λ 2 N + η(λ N -1 -λ N ) 2 + |α| 4 + |β| 4 λ N -1 λ N .
For the sake of illustration, we focus now on η = -1, in which case

ρ α,β 2 (N -1, N ) = λ N -1 λ N |α| 2 + |β| 2 2 = λ N -1 λ N .
The probability of co-occurrence of items N -1 and N thus does not depend on α and β, while the product of their marginal probabilities of occurrence does change: we have thus introduced correlation between the events {N -1 ∈ γ α,β } and {N ∈ γ α,β }. Because γ α,β is a DPP with Hermitian kernel, this correlation is nonpositive.

To go further, an interesting question would be to relate the different ways to quantify repulsiveness in a DPP [START_REF] Biscio | Quantifying repulsiveness of determinantal point processes[END_REF]Lavancier, 2016, Møller and[START_REF] Møller | Couplings for determinantal point processes and their reduced Palm distributions with a view to quantifying repulsiveness[END_REF] to properties of the bases used in the construction of free fermions.

A DPP from bosons in a non-Gaussian state

Another example that illustrates the subtleties of associating correlation functions to states is that (non-Gaussian) density matrices corresponding to interacting (i.e., non-free) bosons can give rise to a DPP for some observables.

One example is the Tonks-Girardeau gas [START_REF] Girardeau | Relationship between systems of impenetrable bosons and fermions in one dimension[END_REF], which consists of impenetrable (statisticians would say hardcore) bosons in X = R. The very strong interactions prevent two bosons to be at the same position, which is reminiscent of the Pauli exclusion principle for fermions. In particular, one considers the usual bosonic field operators a(x), a † (y), satisfying the canonical commutation relations as long as x = y, but the hardcore constraint is enforced by requiring that a(x) 2 = a † (x) 2 = 0 and [a(x), a † (x)] + = 1. These additional constraints are the results of adding a large penalty term to the Hamiltonian at coincidence. Now, there actually is a transformation, called the Jordan-Wigner transformation, that maps these hardcore bosons onto free fermions [START_REF] Lieb | Two soluble models of an antiferromagnetic chain[END_REF]. Concretely, one can write a(x) = S(x)ψ(x) and a † (x) = ψ † (x)S † (x), where ψ(x), ψ † (x) obey the standard canonical anti-commutation relations, and S(x) is a unitary operator called a string operator, which allows for preserving the commutation relations of the bosonic operators. 50By construction, all observables that are built from number operators51 are the same as for free fermions, since a † (x)a(x) = ψ † (x)ψ(x), and the correlation functions of an ideal detection experiment are those of a DPP with kernel ψ † (x)ψ(y) . However, the coherence function a † (x)a(y) = ψ † (x)S † (x)S(y)ψ(y) is very different from that of the free fermions. For instance, the largest eigenvalue of the operator with kernel a † (x)a(y) evaluated in a pure state typically scales like the square root of the number of particles; see e.g. [START_REF] Forrester | Finite one-dimensional impenetrable bose systems: Occupation numbers[END_REF] for investigations on the ground state properties of a Tonks gas in a harmonic trap.

Interacting field theories and point processes

The relationship between, on one side, permanental and determinantal point processes, and, on the other side, bosonic and fermionic quantum field theory, is based on free models, i.e., with a quadratic Hamiltonian in the fields. However, in physical models, Hamiltonians are rarely quadratic, but take into account interactions between the particles. For a given Hamiltonian, success is achieved when it is possible to approximate the coherence functions of the field, which correspond to concrete experimental measurements. The main issue is that Wick's theorem does not apply.

There are interacting Hamiltonians that still yield closed-form correlation functions for the underlying ideal detection experiments. Such systems are usually called integrable by physicists. For instance, the Calogero-Sutherland model is a one-dimensional system of interacting fermions, for which ideal detection leads to well-known point processes, called β-ensembles in random matrix theory; see e.g. [Forrester, 2010, Chapter 11] as well as [START_REF] Stéphan | Free fermions at the edge of interacting systems[END_REF][START_REF] Smith | Full counting statistics for interacting trapped fermions[END_REF].

Yet, most interacting systems are not integrable. In Sections 4 and 5, we have used simple perturbation-theoretic arguments like Assumption 5 to work with the interaction Hamiltonian. There is now a significant and sophisticated toolbox to approximate coherence functions in the presence of interaction. Key tools include mean-field methods [Goldenfeld, 2018, Chapter 3], Feynman diagrams [Folland, 2008, Chapter 6] and renormalisation [Folland, 2008, Chapter 7]. At a high level, it would be interesting to investigate what these methods say about point processes. In words, one of the lessons of quantum field theory is that it is possible, even in the presence of interaction, to get "close" to a free situation in some regimes. Does this imply approximation results for the sophisticated point processes behind models with interaction? Conversely, can we use the mathematical technology of DPPs to perform non-trivial calculations for many-body problems?

Moreover, a central concept in quantum field theory is that of universality, i.e., the fact that many different interacting models behave, in a certain regime usually qualified as the low-energy regime, essentially in the same way. Put differently, the detailed form of the interactions are irrelevant to understand this low-energy regime, only the dimension of the space and the symmetries of the model usually matter. The archetypal phenomenon illustrating this is phase transitions [START_REF] Goldenfeld | Lectures on phase transitions and the renormalization group[END_REF]. Does this universality connect with similar results for point processes, such as those obtained for the eigenvalues of random matrices [START_REF] Anderson | An introduction to random matrices[END_REF]?

Constructive arguments for point processes

The fermionic system "above" a DPP can help address fundamental questions on DPPs, about their invariance or their construction. For instance, a single DPP corresponds to many kernels. In particular, for any f , the kernel x, y → f (x) f (y) K(x, y) yields the same correlation functions (184) as K. Without strong assumptions on the kernel, it has proven difficult to find all the transformations of a kernel that leaves a DPP invariant. It would be interesting if the fermionic framework helped us to understand these invariances; see e.g. recent partial results by [START_REF] Olshanski | Determinantal point processes and fermion quasifree states[END_REF].

As another example, generalizing permanental and determinantal point processes, one can define α-DPPs as having correlation functions

ρ k (x 1 , . . . , x k ) = det α K(x i , x j ) k i,j=1 , k ≥ 1, (184)
where, for an n × n matrix A = ((a ij )), det α (A) := σ∈Sn α n-ν(σ) n i=1 a iσ(i) . [START_REF] Shirai | Random point fields associated with certain Fredholm determinants I: fermion, Poisson and boson point processes[END_REF] have studied the existence of α-DPPs for α ∈ [-1, 1], where α = -1 corresponds to a DPP, and +1 to a permanental point process. Some α-DPPs appear naturally when marginalizing a projection DPP with a separable kernel over coordinates [START_REF] Mazoyer | Projections of determinantal point processes[END_REF]. It is a natural question whether α-DPPs correspond to any physical system of particles. As a partial affirmative answer, [START_REF] Cunden | Free fermions and α-determinantal processes[END_REF] give a limit procedure to construct certain α-DPPs out of fermionic processes.

In the same vein, can we build point processes from particles with more exotic commutation rules for their ladder operators? A natural physical example is anyons [Ezawa, 2008, Chapter 8], for which a phase factor appears for each transposition.

The links of fermions to Pfaffian point processes, another generalization of DPPs are also a promising research direction [START_REF] Koshida | Pfaffian point processes from free fermion algebras: Perfectness and conditional measures[END_REF], as well as the point processes behind quasi-free states [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF][START_REF] Lytvynov | Fermion and boson random point processes as particle distributions of infinite free Fermi and Bose gases of finite density[END_REF][START_REF] Lytvynov | On the correlation measure of a family of commuting hermitian operators with applications to particle densities of the quasi-free representations of the car and ccr[END_REF][START_REF] Olshanski | Determinantal point processes and fermion quasifree states[END_REF], namely (possibly non-Gaussian) states to which Wick's theorem still applies.

A teaser for Part II

This manuscript is intended to become Part I of a monograph. In Part II, we shall present selected topics at the intersection of point processes and quantum optics, using the vocabulary of Part I. We will describe landmark experimental measurements of HBT signals, with both photons and (bosonic and fermionic) atoms. We will discuss the application of determinantal point processes to the study of non-interacting trapped fermions in statistical physics. We will show how the formalism of quantum field theory can help to prove fundamental results on point processes appearing in combinatorics. Finally, we will discuss electronic quantum optics and its interactions with signal processing.

Figure 1

 1 (a) shows a few Poisson samples where µ is the Lebesgue measure on X = R and λ is constant.Example 1 (Deterministic electric fields yield Poisson point processes). Consider the photodetection setup of Example 0. Let us further assume that the field is quasi-monochromatic, i.e., that the modulus |FE| of the Fourier transform of the function E : R → R concentrates around ω and -ω for a single value of the frequency ω > 0. Note that the symmetry of |FE| is a consequence of E being real-valued. Denote by E + the analytic signal of the electromagnetic field E, i.e.,

Figure 1 :

 1 Figure 1: (a, b, c) The intersection with W = [0, 1] of four samples of each of the three types of point processes introduced in Section 2. All three point processes are translationinvariant and are scaled to have the same expected number of points falling in W . (d, e, f) The pair correlation function (3) of the three point processes, and histograms of pairwise distances from 100 independent samples.

Figure 2 :

 2 Figure 2: A sample E of the zero-mean Gaussian process with Lorentz kernel (12), and the modulus |E + | of the corresponding analytic signal.

Proof.

  The proof is simply based on the fact, to be shown below, that [N, O nc,na ] = (n cn a )O nc,na . Then it follows that (n c -n a ) O nc,na = [N, O nc,na ] = Tr (ρN O nc,na -ρO nc,na N ) =0 by the cyclicity of the trace and [ρ, N ] = 0. To prove [N, O nc,na ] = (n c -n a )O nc,na , we proceed by recurrence. It is true for max(n c , n a ) ≤ 1 by the commutation relations of the creation and annihilation operators with the number operator (51). Now assume that it is true up to some n c , n a , and write O nc,na+1 = O mc,ma a i O pc,pa with m c + p c = n c and m a

  1) contractions of order N . Theorem 3 (Wick's theorem). Consider the Gaussian density matrix ρ given by (61) and N = 2k linear combinations b 1 , . . . , b N of creation and annihilation operators, in any basis. Then b 1 . .

Figure 3 :

 3 Figure 3: Experimental measurement of the second-order correlation signal of light coming from a single quantum dot showing anti-bunching at coincidence (zero time delay). The traces show different pump laser powers, all of which are multiples of the lowest power P 0 .The oscillations at larger delay correspond to so-called Rabi oscillations between the two energy levels; as the pump power is increased these oscillations increase in frequency. The solid line is a fit to a model. Figures taken from[START_REF] Flagg | Resonantly driven coherent oscillations in a solid-state quantum emitter[END_REF].

Figure 4 :

 4 Figure4: Representation of a simple model of a quantum wire where fermions are sent from a reservoir at equilibrium (grand-canonical ensemble) into a one-dimensional wire where they can travel and are probed by an electro-detection device. Detection probabilities are given by electronic coherence functions.

G

  It can in fact be computed directly as the Fourier transform (in the sense of distributions) of the zero-temperature Fermi-Dirac energy density distribution f 0,ζ (ω) = Θ(ζ -ω), where Θ is the Heavyside function. It reads 40

  Figure 5: Tableshowingthe structure of the electronic coherence functions of order up to N in the presence of one up to N fermions in some given wavepackets.

3 Elements of quantum field theory

  3.1 The mathematical framework of quantum theory . . . . . . . . . . . . . . . 3.2 Two fundamental one-particle systems . . . . . . . . . . . . . . . . . . . . . 3.2.1 The qubit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 The harmonic oscillator . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Modelling a finite number of particles . . . . . . . . . . . . . . . . . . . . . 3.3.1 Subsystems, bosons and fermions . . . . . . . . . . . . . . . . . . . . 3.3.2 Example N -particle states . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Occupation number representation . . . . . . . . . . . . . . . . . . . 3.4 Modeling an indefinite number of particles . . . . . . . . . . . . . . . . . . . 3.4.1 Fock spaces for bosons and fermions . . . . . . . . . . . . . . . . . . 3.4.2 Creation and annihilation operators . . . . . . . . . . . . . . . . . . 3.4.3 Modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.4 Field operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.5 Observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Wick's theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Gaussian density matrices . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Where permanents and determinants appear . . . . . . . . . . . . . 3.6 Bosonic coherent states model classical fields . . . . . . . . . . . . . . . . . 3.6.1 Definition and properties . . . . . . . . . . . . . . . . . . . . . . . . 3.6.2 The Husimi distribution and time-frequency analysis . . . . . . . . .

Mixed states and Born's rule. In

  practice, most experimental devices are not able to repeatedly produce a given state |ψ ∈ H , but rather a noisy version of it. To model this noise, first associate to each state |ψ ∈ H the projector |ψ ψ| : |h → ψ|h |ψ , which we also abusively call a state. Now consider a linear operator on H defined as

		, Chapters
	VII & VIII].
	ρ = E |ψ ∼p |ψ ψ| ,	(24)

  these are valid correlation functions, namely of the point process formed by k times t 1 , . . . , t k drawn independepently with probability density function proportional to t → |ϕ + (t)| 2 .

Together with a 1974 conference communication, that appeared later as[START_REF] Macchi | The fermion process -a model of stochastic point process with repulsive points[END_REF], and earlier work by Jean Ginibre; see the preface to the recent re-print[Macchi,

2017] of Odile Macchi's thesis for a broader history.2 including some of us prior to this work!

https://dpp-fermions.sciencesconf.org/

https://indico.in2p3.fr/event/25182/

In later sections, the field will be modeled by a collection of operators, and linear polarization will thus correspond to a different mathematical assumption.

Overloading variable names is common in physics, and we shall follow this convention when possible without confusion.

Stricto senu, we are defining here simple point processes, i.e., such that the samples never contain a given point of X more than once. Since all point processes in this work are simple, we take this shortcut and identify a (simple) point process with a random configuration.

Physicists speak of an ensemble average, and denote it using angle brackets • . We stick with the symbol E for classical expectations, and we reserve brackets only for Hilbert space inner products and quantum averages.

Physicists might prefer to call it density of particles, and reserve intensity for other physical quantities. To avoid confusion, we will write first correlation function in full.

Without this assumption, a Poisson point process would not necessarily be a simple point process.

Statisticians call this parameter function the intensity of the process. We refrain from using this naming convention, to avoid confusion with any physical intensity.

However, E + is complex-valued, i.e., the vector in (11) is a complex multivariate Gaussian vector.

This theorem is also known to physicists as Wick's theorem, by analogy with a similar theorem for the quantum average of products of certain operators; see later in Section 3.5.

For brevity, we will stay informal on that point and ignore the subtleties coming from considering a projective Hilbert space rather than H ; see[Folland, 2008, Chapter 3].

In other words, a mixed state can be used to describe epistemic uncertainty on the preparation process, while draws from P A,ψ in (22) represent aleatoric uncertainty. Unlike epistemic uncertainty, aleatoric uncertainty cannot be reduced by better knowledge of the system.

The use of the word matrix for an operator can be confusing at first glance, but it is standard here.

Physicists talk of the expectation value of A in state ρ.

In Section 6, we shall come back to the idealized free fermion construction of this section, introducing a temperature parameter to obtain non-projection DPPs.

Note that like the number operator later on, the ladder operators are defined only on the dense subset of finite linear combinations of Fock states; see[START_REF] Folland | Quantum Field Theory: A tourist guide for mathematicians[END_REF] for discussions on their extension.

We also note that these ladder operators were inspired by those introduced by Dirac to solve the harmonic oscillator; see Section 3.2.2. However, one should keep in mind that they do not act on the same

Arguably, one already has this isomorphism in mind when writing |n 1 , n 2 , . . . for a Fock state.

The name phase space comes from the connection through (37) to the classical phase space with position-momentum coordinates.

Often called correlation functions in physics, but we refrain from using correlation here to avoid confusion with the concept of correlation function of a point process; see Section 2.

Physicists will recognize here an instance of the more traditional form H I (t) = -d(t) • E(r, t), where we assumed the dipolar moment d(t) to be aligned with the polarization u of the field E(r, t) = E(r, t)u.

The normalization constant N k is reinserted into the coherent state parameter in this derivation to simplify the notation.

Strictly speaking, we should allow α to depend on time, but the argument below does not change.

usually the coherence function varies over some timescale, so large time means t sufficiently large compared to that timescale

a † and a are usually kept for photons, or general particles.

Actual experimental settings are more intricate to describe than for photons exciting a two-level atom; Part II will contain examples.

This computation is standard in physics. At this stage, we take it for granted; but see Equation (181) for an explicit derivation.

Equivalently, we could have specified the chemical potential ζ to a certain value ζ F , called the Fermi chemical potential. Moreover, k F is related to the Fermi velocity v F or the Fermi energy ν F by the free Hamiltonian so these too can be used as inputs of the model.

Again, this is a shortcut to avoid describing the Hamiltonian and Schrödinger's equation.

This result is independent on whether one assumes that α i commutes or anti-commutes with a j .

It is however possible to define fermionic coherent states in terms of anti-commuting variables (elements of an exterior algebra) by the action of the analogue of a displacement operator[Berezin, 1966, Chapter 1]. While they are very useful objects to define and use (especially to write path integrals), they do not possess the nice set of properties of bosonic coherent states that allow for deriving coherence functions.

Just like a for photons, ψ is not Hermitian, and thus not an observable, stricto sensu. Yet, for a, the average value corresponds to the intensity of the field, and is thus experimentally accessible as a statistical average. It thus makes sense to wonder whether the average value of ψ can be non-zero.

This property can actually be used as a definition of bosonic coherent states[Cohen-Tannoudji et al., 2019a].

This generalises beyond physically realistic fermions of half-integer spins, and applies to non-relativistic quantum mechanics, see for instance[START_REF] Johansson | Comment on "reasonable fermionic quantum information theories require relativity[END_REF][START_REF] Szalay | Fermionic systems for quantum information people[END_REF].

Note that for a given p, the limit νp → ζ corresponds to the onset of Bose-Einstein condensation, where a macroscopic number of bosons start to occupy the state p.

Alternatively, the formulation in terms of operators can be replaced by a "Gaussian functional integral" over complex commuting (anti-commuting) fields for permanents (determinants). The connection with Gaussian measures is then more direct, at the price of subtle mathematical difficulties. Furthermore, in the case of fermions, this implies generalising the notion of integration to variables that anti-commute, see for instance[Berezin, 1966, Chapter 1].

In particular, it satisfies [S(x), S(y)] = 0 and S(x)ψ(y) + sign(x -y)ψ(y)S(x) = 0.

Usually rather called density operators when we use a generalized basis like here.

Acknowledgments

We thank all participants to the Lille and Lyon workshops, and all participants to the Lille workgroup on point processes and applications; the discussions held in these circles helped a lot in the construction of this document. We hope that the document, in return, participates to developing this cross-disciplinary research field. In particular, RB thanks all physicists with whom he has interacted over the last years for their enthusiasm and their patience, explaining basic physical concepts over and over until we converged to a description that spoke to all of us.

We acknowledge support from ERC grant Blackjack (ERC-2019-STG-851866), ANR AI chair Baccarat (ANR-20-CHIA-0002), ANR grant cl1shot reloaded (ANR-14-CE32-0017), ANR grant QuSig4QuSense (ANR-21-CE47-0012), the Joint Research Project SEQUOIA (17FUN04) within the European Metrology Programme for Innovation and Research (EMPIR), ANR grants QRITiC I-SITE ULNE and ANR-16-IDEX-0004 ULNE, ANR grant COSQUA (ANR-20-CE47-0001-01), ANR grant MENTA (ANR-22-QUA2-0008-01), ANR Grant RaMaTraF (ANR-17-CE30-0027-01), ANR grant Dimers (ANR-18-CE40-0033), ANR grant Combiné (ANR-19-CE48-0011), Labex PALM (ANR-10-LABX-0039-PALM), and Labex MILYON.