Bridging nano-optics and condensed matter formalisms in a unified description of inelastic scattering of relativistic electron beams

In the last decades, the blossoming of experimental breakthroughs in the domain of electron energy loss spectroscopy (EELS) has triggered a variety of theoretical developments. Those have to deal with completely different situations, from atomically resolved phonon mapping to electron circular dichroism passing by surface plasmon mapping. All of them rely on very different physical approximations and have not yet been reconciled, despite early attempts to do so. As an effort in that direction, we report on the development of a scalar relativistic quantum electrodynamic (QED) approach of the inelastic scattering of fast electrons. This theory can be adapted to describe all modern EELS experiments, and under the relevant approximations, can be reduced to any of the last EELS theories. In that aim, we present in this paper the state of the art and the basics of scalar relativistic QED relevant to the electron inelastic scattering. We then give a clear relation between the two once antagonist descriptions of the EELS, the retarded green Dyadic, usually applied to describe photonic excitations and the quasi-static mixed dynamic form factor (MDFF), more adapted to describe core electronic excitations of material. We then use this theory to establish two important EELS-related equations. The first one relates the spatially resolved EELS to the imaginary part of the photon propagator and the incoming and outgoing electron beam wavefunction, synthesizing the most common theories developed for analyzing spatially resolved EELS experiments. The second one shows that the evolution of the electron beam density matrix is proportional to the mutual coherence tensor, proving that quite universally, the electromagnetic correlations in the target are imprinted in the coherence properties of the probing electron beam.


I. INTRODUCTION
Electron energy loss spectroscopy (EELS) performed in a (scanning) transmission electron microscope ((S)TEM) analyzes the energy loss of fast electrons (with energy ranging typically from 30 to 300 keV) after their interaction with a target sample. It allows to probe a wide range of charge excitations in solids such as phonons [1], plasmons [2], excitons [3], and core electron-hole excitations [4] over a wide range of energies, typically from 40 meV to thousands of eV. Moreover, the versatility of the electron optical set-ups allows to achieve either high spatial resolution [5] (better than half an angström) or high momentum resolution [6,7] (smaller than µm −1 ) or to probe directly the symmetry of the excitations [8]. This variety of experimental configurations is illustrated in Fig. 1. A similar diversity is found in the theoretical descriptions of EELS experiments depending on how one treats the following aspects: • The energy range of the probed excitations, usually dispatched in two classes, low-losses (less than tens of an eV) or core-losses (more than a hundred of eVs).
• The classical or quantum character of the beam electrons, corresponding to a description as either point charges or wave functions.
• The classical or quantum character of the target sample.
• The time-dependency of the electron wave function in the quantum formalism.
• The geometry of the experiment, especially whether the scattering events are analyzed in real or reciprocal space • The retardation in the electron-sample interaction and in the fields propagating within the sample.
• The overlap of the electron beam with the sample, e.g., the importance of bulk versus surface effects.
Since the general description (including all of the above cases) to EELS seems hardly viable, a plethora of different theories adapted to different sets of parameters have been developed. For example, the low-loss characterization of electromagnetic surface excitations such as surface plasmons is well interpretable assuming a classical character for both the electron beam and the tar- allow to extract information contained in this final electron state. The gist of our work is to derive the general kinetic equation connecting the initial to the final electron beam state through an interaction kernel embedding all the physical details (classical, relativistic or quantum) of the target.
get [9][10][11]. In such models the EEL cross-section is proportional to a well-known classical electromagnetic quantity, the (retarded) Green's dyad of the system [12]. The Green's dyad description is heavily used in nano-optics, necessitating an accurate description of electromagnetic field propagation within complex geometries of dielectric media. On the opposite side, core-electron excitations are usually described in a Fermi's Golden rule approach rooted in a quantum treatment of the inelastic scattering processes between the beam electron and the target degrees of freedom [13]. This approach can be generalized by employing the mixed dynamic form factor (MDFF) formalism introduced by Kohl and Rose [14] which typically neglects retardation as well as many-body effects in the target beyond the mean field level, when applied to core losses. Both the Green's dyad and the MDFF descriptions have been extended, thereby approaching each other. E.g., by using a fluctuation-dissipation approach Green's dyad has been extended to take into account the quantum character of the electron [11]. Many-body effects, on the other hand, have been absorbed in the MDFF to describe low-loss volume excitations, such as plasmons (on the random phase approximation level) [15,16] or excitons (employing Bethe-Salpeter approximations) [17] in condensed matter systems.
The general relation between the two approaches, however, has still not been established, which is unsatisfactory not only from an intellectual point of view but also concerning the interpretation of EEL spectra in different experimental settings. Indeed, it raises unnecessary conceptual walls between alternative descriptions of the electron beam, of the target excitations, and of the electron to target interaction. For example, it has been demonstrated [18] that the so-called magic angle, at which the electrons are most likely deflected after excitation of specific core-loss transitions in anisotropic materials, strongly depends on the retarded character of the electron to target interaction, which had been considered as negligible in core-loss investigations before. Also, the description of the effect of interferences between inelastically scattered electrons, otherwise speaking, coherence effects, has been long discussed in terms of MDFF in the framework of inelastic holography of bulk plasmons [19] or EELS atomic resolution mapping [20][21][22], while coherence effects in surface plasmon excitations relied on Green's Dyad approaches [23]. Even more, EMCD [24] and phase-shaped EELS applied to plasmonics [8,23] relies physically on the same ingredients -electron phase manipulation to mimick the action of an X or optical photon -yet are described within totally different frameworks.
It is therefore not surprising that early on, researchers have sought for a comprehensive description of EELS in the (S)TEM, or at least for a bridge between different approximations. For example, Ritchie and Howie [25] could explain how the interferences of inelastically scattered electrons are washed out by integrating over all scattering angles rendering the quantum and classical descriptions of the electrons essentially equivalent for most of the experimentally relevant low-loss STEM-EELS cases. In order to model dynamical scattering effects in diffraction, Dudarev, Peng and Whelan [26] drew a direct link between the density matrix of the electron beam and the MDFF in the quasi-static case. Later on, Schattschneider et al. [15] applied a similar approach to relate EELS to the MDFF, and to establish a clear link between reciprocal and real-space representations of the electrontarget interaction. Schattschneider and Lichte [19] later used the MDFF formalism to properly describe coherence effects between electrons in the quasi-static limit, following the pioneering work of Kohl and Rose [14]. Later, García de Abajo proposed a fully relativistic description of low-loss EELS experiments, where the quantum nature of the probe electrons could be taken into account [11,23]. Nevertheless, no universal description of EELS in a (S)TEM has yet been provided, which includes that the relation between the different approximations remained somewhat in the dark.
The present work is an effort towards this goal. By extending and bridging several key works [11,14,15,26] we provide a synthetic description of EELS in a (S)TEM. This description is valid for any sort of classical or quantum description of the electrons, arbitrary equilibrium description of the target object, any sort of excitations (low-loss and core-loss, surface and bulk) and using retarded or non-retarded approaches alike. Incidentally, we are formally establishing the link between the Dyad and the MDDF approach, extending the later to the retarded case, therefore widening the applicability of our work to other spectroscopies not necessarily involving electrons. Working out such a theory is challenging because both quantum and relativistic effects need to be taken into account. As a consequence, the problem of computing the complete beam electron-target interaction cannot be solved in a closed form and different approximation schemes have to be employed, notably diagrammatic perturbation techniques. Accordingly, we won't consider effects connected to the finite temporal length of electron wavepackets [27][28][29][30]. In other words, the wavefunctions encountered throughout this paper do not represent quantum wavepackets but rather electron beams in a steady-state illumination, in strict analogy with photon wave optics [31]. With this in mind, we can synthesise the electron energy-loss processes studied in this paper with a diagrammatic perturbation language roughly schematized in Fig. 2.

FIG. 2:
Diagram representations of the electron energy loss process: a beam electron described by a wavefunction ψ(r, t) interacts with a target represented by a (many-body) wave function ξ(r, t) by exchanging virtual photons. The target wave function and photon propagators are generally dressed (renormalized) by interacting with the various degrees of freedoms of the target (indicated by double lines). Processes involving the emission of photons (represented by gray lines), such as Bremsstrahlung, are not considered in this paper.
Here, the inital and final beam electron state is represented by the wave functions ψ i,f (r, t) and the target by the wave functions ξ i,f (r, t). (Virtual) photons, indicated by wiggly lines, of arbitrary numbers and orderings are exchanged in the course of interaction. In analysing the diagram, we first note that relativistic effects can emerge from both the probe and the target [32]. Indeed, beam acceleration voltages in a TEM typically range from 80 to 300kV leading to electron velocities v between 0.5 and 0.78c, and variations of the Lorentz factor γ = 1 between 1.16 to 1.59 [33]. Therefore, TEM electrons are relativistic, which translates into (A) a modification of the kinetic properties of the electron such as renormalized masses [34], (B) retardation effects in the electromagnetic interaction [35], (C) Cherenkov losses [6,36,37]), and (D) sizeable current interaction effects. While (A) and (D) directly scale with the velocity of the electron beam, the retardation effects in the interaction cannot be neglected anymore when the length scale L of the charge density fluctuations associated with an excitation of energy ω in the target become important ωL/c > 1 (see [32] and references therein). This situation typically occurs in plasmonics which leads to a frequency red-shift, a loss of spatial coherence or even to mode splitting [38]. The large energy of the beam electrons, however, also renders spin-orbit coupling effects in the scattering process itself negligible [39,40], hence a full quantum relativistic modeling of the electron beam is not required. Indeed, approaches to EELS or electron diffraction based on the Dirac equations have been developed [34,[41][42][43][44] and give results comparable to what the Klein-Gordon equations do [26,[45][46][47].
In what follows we deliberately focus on EELS including all possible electron energy-loss mechanisms. However, we do not explicitely compute (secondary) scattering events involving the emission of photons (indicated gray in Fig. 2), should it be due to Cerenkov, Bremsstrahlung or cathodoluminescence for example, in order to keep the length of the paper at bay (indeed the perturbation technique used throughout is well suited to also describe these events and may be easily extended). Most importantly, we will restrict the inelastic interaction to the first order Born level, which, however, does not exclude to consider multiple stacked first order Born events as in Sec. VI. The basis for our considerations is the standard diagram perturbation technique as discussed in, e.g., [48]. Accordingly, in order to properly take into account many body effects in the target (e.g., inelastic interaction with the electron gas at the Fermi level) and partial coherence in the beam (e.g., as produced via inelastic interaction) a generalization of the wave diagram 2 to a density operator diagram (i.e., in the language of second quantization) techniques is in order (Fig. 3). The corresponding diagram is obtained by taking the tensor product of the wave diagram (as indicated) and identifying the initial and final states of the target on both sides (which follows from tracing over all target degrees of freedom).
Here, the gray loop generally contains all possible diagrams up to infinite order, i.e., the exchange of (dressed) virtual photon lines in the target indicated in Fig. 2 the same on the left and right hand side of (a), the target side of the tensor product can be connected forming a loop, i.e. a 4-point correlation (response) function or (generalized) two particle Green's function (b). All ingredients of the resulting diagram are subject to various approximations (c), all of which treated in the paper.
is referred to as 4-point correlation function, two particle function, polarization propagator or generalized Green's function in the literature; and represents the fundamental object encapsulating the physics of the excitations of the target. It is well suited to encompass many-body interaction effects (indicated gray) in the target, which particularly affect the low-loss regime. In this article we will not elaborate on the various sophisticated strategies, which have been developed to compute the 4-point correlation function (typically relying on some infinite series expansions in the interaction). We will rather discuss the various descriptions of the beam electron, the (virtual) photon exchange and the ramifications of relativity.
Indeed, all relevant descriptions of EELS may be identified with certain approximations to the general diagram 2. In particular neglecting many-body interactions (beyond the mean field level) in the target corresponds to the MDFF used in core-loss computations, whereas the Green's dyad is a semiclassical version of the photon propagators including the correlation function. Accordingly, Fig. 3 is well suited to provide an overview of what is treated in this work. In Sec. II we give a more detailed account of the state of the art of EELS in terms of experimental setups and pertinent theoretical descriptions, corresponding to various combinations of approximations indicated in Fig. 3. In Sec. III, we introduce the various conventions and notations used throughout this paper. We also recall some results on free-space photon propagators and their expression in different gauges. In Sec. IV we review electron scattering within the framework of scalar relativistic quantum electrodynamic (QED). The main conclusion of this part is that in addition to the possibly retarded Coulomb interaction (i.e., a charge-charge coupling term) also a magnetic interaction (i.e., a current-current interaction) needs to be considered because of the relativistic velocities of the beam electrons. In Sec. V we focus on the 2-particle function (see Figure 3). For completeness and pedagogy, we first re-derive the expression for the MDFF and connect it to the screened interaction. This way, we draw a parallel between the quasi-static formalism of nano-optics with the one of core-loss spectroscopy. We then focus on the retarded case, where we heavily rely on the fact that the interaction can be modeled in terms of correlation functions and photon propagators (see Figure 3). Using a pedestrian approach of the QED, we calculate the exact photon propagator in the presence of a polarizable material, taking special care of the gauge choices. In complete analogy to what has been done with the MDFF in the quasi-static case, we then connect this photonic kernel to the charge and current density correlation functions of the scatterer. This makes it possible to introduce a fully retarded definition of the MDFF as the spatial Fourier transform of the 4-susceptibility of the target. By doing that, we introduce relations that bridge the world of solid-state physics and optics and that can be applied well beyond the case of electron energy loss spectroscopy. In Sec. VI we then consider the electron probe part in detail. For reasons exposed earlier, we model the electron propagation using the Schrödinger equation with a semi-relativistic correction i.e. mass renormalization [47]. Following the seminal demonstration of Dudarev and collaborators [26], we employ the Bethe-Salpeter equation in the ladder approximation in order to calculate the electron density matrix propagation equation. We then consider separately the quasi-static and the retarded interaction, which give the kinetic equation of Dudarev et al. and its retarded counterpart, respectively. Sec. VII is dedicated to the contextualization of our developments and proposal of different applications. Particularly, by taking the appropriate limits, we demonstrate that our equations encompass all previously obtained results.
As important applications of the formalism developed in this paper, we provide three important formulas. For the sake of pedagogy, we present them here, but they will be reproduced, derived and discussed in length in this manuscript. The first is a concise and universal formula for interpreting the most commonly performed EELS experiments and reads: This formula relates directly, within the framework of Fermi's golden rule (first order Born approximation), the loss probability per unit time Γ (i.e., the EEL spectrum) to the incoming electron wavefunction (ψ i ) at points x and x , the related scattered electron wavefunctions ψ f , and the correlation function of the targetĈ(x, x ). The later can be calculated for the relevant charge excitations such as bulk or surface phonons, plasmons, coreelectrons, etc., in either the quasistatic limit or retarded case, as discussed previously (see Fig. (3)). It is therefore the most general and synthetic formula for computing EELS. The second formula reads: where r = (r ⊥ , z) respectively denotes the transverse (orthogonal to the electron beam trajectory) and longitudinal (along the electron beam trajectory) coordinates. Here, the mutual coherence tensor [49,50] Λ zz (r ⊥ , r ⊥ , q z , −q z , ω) relates the transition from initial to final density matrices, ρ i (r ⊥ , r ⊥ ) and ρ f (r ⊥ , r ⊥ ), of the probe electron by an energy loss process of energy ω and characteristic momentum q z = ω/v. The zz indices hind at the tensorial character of Λ, which we can typically neglect in EELS because of the small scattering angles of the fast beam electrons around the zdirection (paraxial approximation). This generalization of the quasi-static work of Kohl and Rose [14] shows that quite universally, the electromagnetic correlations of the target are imprinted in the coherence properties of the electron beam, and makes straightforward the interpretation of holographic or phase-dependent experiments regardless of how the target is described. The third formula reads: and provides a definition of the relativistic MDFF S µν as the imaginary part of the 4-susceptibility χ µν which can be expressed as a 4-current correlation function in the framework of the Kubo linear response theory. This expression constitutes the root of the relations bridging the world of condensed matter physics and optics. Finally, as guide for the reader, we summarized on table I the main quantities and equations derived in this paper as well as their relation between them.
For brevity, Hartree atomic units ( = e = m e = 1) and Lorentz-Heaviside units for the Maxwell equations will be used from now on, unless otherwise specified. The 3-vectors are labeled by roman letters and written in standard font as x ≡ x a = (x 1 , x 2 , x 3 ). The 4-vectors are labeled by greek letters and written in roman font as

II. STATE OF THE ART
Most of quantum relativistic scattering theories for TEM electrons have been developed for diffraction and holography [26,34,41,42,45,46], i.e. elastic scattering, and only a few deal with electron energy loss spectroscopy [43,44,47], i.e. inelastic scattering. In this section, we will review the principal results and theories of the literature. For analyzing characteristic energy-losses of electrons in a TEM (or EELS) the detector is typically placed in the energy-dispersive plane of an energy filter in the far field of the sample, which allows to discriminate the scattered electrons according to their energy loss and scattering angle. Alternative modes excite the magnetic coils of the TEM such to image the achromatic plane of the filter or its reciprocal on the detector, which can be used to record an energy-filtered diffraction pattern or image. In the following we will focus on the conventional far-field setup, which allows probing the dispersion of characteristic excitations such as plasmons or corelosses, typically yielding the most interesting analysis of the solids electronic structure (e.g., including electronenergy-loss linear dichroism (ELD), electron-energy-loss magnetic chiral dichroism (EMCD) [51]).
The most widespread theory of EELS of optical excitations [52] is built upon a semiclassical description of the scattering process founded in a point-like description of the beam electrons. This notably implicates the initial and final electron states not being (stationary) energy eigenstates, which distinguishes this semiclassical approximation from quantum perturbation theory, where initial and final states are energy eigenstates. We recall that the most important result of the semiclassical formalism is that the retarded electron energyloss probability Γ R appearing in the overall energy loss ∆E = ∞ 0 ωΓ R (ω)dω reads: G 0 is the screened electric Green dyadic of the classical Maxwell equations defined through [11,53,54]: and the free Green's dyad: Here is the frequency dependent dielectric function, which is a local quantity in this classical setting (e.g., given by Drude-Lorentz theory). With the screened retarded Green's dyad the induced electric field is given by an external current source through Considering that the scattering predominantly occurs into very small angles around z−direction and that the energy loss in the low-loss regime is small compared to the total energy of the electron, we can write both the initial and final current as J (r, ω) = ρ 0 (r ⊥ , −q z )e iqzz . Here, q z = ω/v and ρ 0 (r ⊥ , −q z ) denotes the Fourier transform of the density at t = 0 along z−direction and we assume that it is a real quantity (i.e. symmetric along the z−axis). With that the loss probability (5) can be further simplified to This equation exhibits the familiar structure of the first order perturbative approach (depicted in Fig. 3) in that we can distinguish between the (initial and final) electron states and the (generalized) Green's function. Indeed, this structure is frequently encountered throughout regardless of the actual derivation details (e.g., 1st order Born approximation, linear response, semiclassical approximation), since all of them employ the notation of a weak interaction at a certain stage. The quasistatic, non-retarded limit of the Green's dyad is obtained by solving the longitudinal part of the full Maxwellian response, Eq. (7), which corresponds to the solution of the first Maxwell equation. In 1987, Echenique and collaborators proposed a quantum version of the quasistatic transition rate Γ QS i→f between initial and final beam electron state using a selfenergy formalism [55]: which corresponds to Fermi's golden rule as shown by García de Abajo [11]. Note that ψ i and ψ f are energy eigenstates of the electron probe with respective energies i and f , which is one difference between the quantum and the classical treatment. W is the screened Green's function for the electrostatic potential defined as: where φ ind (r , ω) is the scalar potential induced at r by a density of charges n ext (r, ω) located at r. The imaginary part of it, Im{W}, corresponds to the spectral density if the latter is a real quantity in the energy representation, which is typically the case (exceptions will be noted). The Green's function W or its imaginary part contain all the quantum mechanical information about excitations inside the target including valence as well as core excitations. It can thus be applied both for describing low-loss [56][57][58] and core-loss EELS.
Using linear response theory, García de Abajo proposed [59] an extension of (10) to the retarded regime [11]: This equation makes use of the paraxial approximation and has been used in several works [23,60] in order to calculate the dichroism in the interaction between a vortex electron state and a (geometrically) chiral plasmonic nano-particle. We will show below how to generalize this equation making use of quantum propagators.
Although the above loss-probability and transition rate formulas are remarkably elegant and intuitive, they do not provide information on the propagation of the wavefunction in the microscope. However, a proper description of a phase-shaped EELS experiment requires the precise description of the illumination and detection systems. Moreover, information on the coherence of the electron beam, which plays a crucial role in holography, is not explicitly present in these expressions. They are therefore not sufficient to model such experiments.
In order to describe these situations, another fundamental object has to be considered. This is the density matrix operator of the beam electrons, which in an energy eigenbasis {|ψ n } reads: where p n are the occupation probabilities associated to each state vector n. Inserting the completeness relation r |r r| = 1, we obtain the fundamental tool for the description of wave optical experiments: the energydependent density matrix. It is defined as [15,61] (see also Eq. (101)): This quantity is particularly rich in terms of information; for example I = ρ(r, r) gives the intensity at coordinate r in position space (typically identified with (conjugated) image planes in the TEM). Even more importantly, the out-of-diagonal elements measure the mutual coherence of the electron field between positions r and r [16]. In other words, non-zero out-of-diagonal terms entail electron interferences in the image plane.
In 1993, Dudarev, Peng and Whelan [26] demonstrated (the different assumptions leading to this formula will be reviewed in details in this paper) that, in the quasi-static limit, the inelastic scattering of high energy electrons by a polarizable material can be described by the so-called kinetic equation: where F denotes the Fourier transform, ρ i and ρ f are the density matrices of the electron probe before and after the interaction and S(k, k , ω) is the so-called mixed dynamic form factor (MDFF) [14] defined as : where k is a wave-vector, n is the electron density operator and |f is an eigenbasis of the target electronic many-body state. Upon comparision with the quasistatic transition rate (10) we note that the MDFF is the Fourier space (and many-body) version of the spatial integral in said expression. It contains all the information on the correlation of the electronic charge density of the scatterer [16].
S is a hermitian tensor, hence may be decomposed into a real-valued symmetric tensor containig information about the non-chiral (i.e., non-magnetic) transitions and an imaginary antisymmetric tensor (represesented by the vector S) containing the information about chiral (i.e., magnetic) transitions, reading [62] in the dipole approximation (to be detailed further below). Note that the chiral transitions introduce an imaginary component here, which marks one of the few cases, where the spectral density is not a real object.
Remarkably, equation (15) shows that correlation are imprinted in the density matrix (and hence mutual coherence) of the beam during the scattering process. It leads to a fundamental principle of electron holography: generating interferences in order to trace back to the electronic correlations in the target. Indeed, EMCD experiments may be designed such to isolate the second term, hence allowing to characterize magnetic materials in the TEM. It would be rather seducing to use such a formalism in the case of nano-optics and e.g. interpret EELS interference effects on surface plasmons in terms of electric field correlations measurements. The MDFF is particularly adapted to model core-loss spectroscopy as the measured phenomena appear to be quasi-static. However, it gives an incomplete picture of nano-optical phenomena where retardation effects dramatically constrain the coherence properties of the field.

A. Gauge fixing and vacuum photon propagator
The Green's function of equation (A9) is the vacuum photon propagator D µ ν defined by: In the following, we are using Einstein summation convention, with greek letters for covariant indices and latin letters of spatial indices. For readability, all quantities and basic equations have been defined in the Annex. In order to calculate D µ ν one needs to invert the Kernel in (A9). Depending on the gauge, this task can require involved mathematical techniques as it may be singular. In nano-optics, principally three gauges for the electromagnetic field [63] are encountered in the literature: the Coulomb gauge, the (partial) Lorenz gauge and the temporal (or Weyl) gauge -each of them having different specific interests. We will therefore give D µ ν in these cases [64,65] only, but keeping in mind that the vacuum photon propagator can be expressed in arbitrary gauges: • The Coulomb gauge corresponds to the condition: It is of particular interest in standard quantum electrodynamics as it enables a simple quantification of the potentials and leaving the Coulomb interaction in its classical and non-retarded form. In the Coulomb gauge the photon propagator reads [64,65]: • The temporal gauge corresponds to the condition: It is particularly interesting because it drastically facilitates the calculation of conductivity in linear response theory. In the temporal gauge the photon propagator reads [64,65]: • The Lorenz gauge corresponds to the condition: Its main interest is to decouple the motion equation for the four components of the potential. Indeed, in the Lorenz gauge, the propagator reads [65]: B. Mutual coherence tensor, electromagnetic local and cross density of states In section II, we introduced the Green dyadic ↔ G propagating a current source to an electric field as a function of the photon energy ω. From this propagator, as it is commonly done in condensed matter [66], one can define a photonic spectral function (which in the case of optical fields is a tensor), usually called the mutual coherence tensor [49,50] (MCT), as: From the fluctuation-dissipation theorem, the former quantity can be connected to the spectral correlations of the electromagnetic field [67,68]: The mutual coherence tensor is a pillar of nano-optics as it contains all the important information about the optical field. Indeed, from the early work of Agarwal [69], we know that by taking the trace of its diagonal elements (i.e. r = r ), we obtain the electromagnetic density of states (EMLDOS): where, for clarity, we made the summation explicitly appear. One can also consider partial electromagnetic density of states by only taking one of the components e.g. for i ∈ {1, 2, 3}: The EMLDOS is involved in the description of a wide range of phenomena such as the Purcell effect [70] or the Casimir effect. In 2013, Cazé and collaborators introduced [71] the cross density of states (CDOS): which of course allows the definition of a partial CDOS e.g. for i ∈ {1, 2, 3}: From equation (26), it is clear that the CDOS measures the electromagnetic correlations between two points in space at the energy ω. As extensively discussed in [50,71,72], this quantity is particularly relevant to characterize the spatial coherence of optical field and study e.g. the Anderson localization in random media. It is also quite standard to normalize the CDOS by the EMLDOS, thus defining a new quantity: which, depending on the context and for historical reasons, is referred to as complex degree of coherence [49] or mode connectivity [72]. Thanks to the Schwarz's inequality, one can show that 0 ≤ C ≤ 1 where the case C = 0 corresponds to uncorrelated points and C = 1 to the situation of strong connection [72]. Extension to polarization-dependent correlations is straightforward from Λ ij and leads e.g. to the definition of a complex degree of mutual polarization [73].
No matter the plethora of quantities defined in the literature over the last decades, it is crucial to state that they are all contained in the most general mutual coherence tensor Λ ij .

IV. FUNDAMENTALS OF SCALAR RELATIVISTIC QUANTUM ELECTRODYNAMICS
In this section we shortly recapitulate the fundamentals of relativistic QED in order to obtain very general covariant (i.e., fully retarded) expression for the inelastic transition matrix elements and transition rates, which are then further detailed in the following chapter. For the sake of clearness we explicitly note the elementary charge e to highlight the perturbation order of the theory. Moreover, we employ upper and lower indices to indicated covariant and contravariant vectors in this section.
As stated previously we can largely neglect the ramifications of the electron spin in the inelastic scattering process as the spin-orbit coupling is small at the large electron velocities considered here. A scalar relativistic formulation of the scattering process is therefore largely sufficient for our purpose. The fundamentals of a scalar relativistic QED, i.e., the Feynman rules following from a perturbative treatment of the Klein-Gordon field weakly coupled to the electromagnetic field, are derived, e.g., in the book of Greiner [74].
Here, we only repeat the fundamentals of scalar relativistic perturbation theory important for the following computations. First, we restrict ourselves to the first order of the perturbation (i.e., scattering matrix terms linear in the interaction constant e 2 ) as higher order terms such as Vertex corrections can be expected to be very small (their relative strength corresponds to those of the related Lamb shift). This allows us to neglect the quadratic coupling term on the right hand side of (32). The same reasoning allows us to get rid of the quadratic term occurring in the definition of the transition current: which appears in the calculation of the vacuum photon propagator D (Eq. 18). We finish our preparations with noting the scalar product as employed within the framework of scalar relativistic QED and introducing the scattering matrix which describes the transition of one initial beam electron state ψ i to a final state ψ f under an influence of the time evolution operator U (t 2 , t 1 ) containing the actual physics of scattering process. In other words we have for the transition of a state vector (i.e. wave function) and the density operator (i.e., density matrix). The latter is the abstract retarded version of the kinetic Eq. (15) noted previously. With these tools and the corresponding Feynman rules, we may now note the scattering matrix for the first order scattering process depicted in Fig. 4(a) The first line contains the derivatives of the beam electron wave function reminiscent of a symmetric version of Garcia de Abajo's formula Eq. (12). On the second line, the relation to the transition current is highlighted, which provides an ab-initio ratio for to the current-current coupling terms in the linear-response regime (Kubo formalism) discussed further below. Indeed, the transition probability per electron into some final state reads which transforms into (40) upon inserting Eq. (38). Here we introduced a new abbreviation for the correlation of two vector potentials. This quantity represents the central part (i.e., between the two outer vertices) of Fig. 4(b) and hence may be directly related to the four-susceptibility (i.e., the central loop the Feynman graph in Fig. 4(b)) In the last expression we summed over all final states, i.e., ξ = f ξ f i , to take into account every possible final scattering state of the target. The latter expression may be transformed into reciprocal space taking into account the homogeneity of the interaction in time, yielding if employing the photon propagator in Lorentz gauge (24). Indeed, there is a one-to-one relationship between the spectral representation of χ and D and the relativistic generalization of the mixed dynamic form factor where we introduced the 4-current operator j µ (k) in reciprocal space. The relationship can be demonstrated by employing the Dirac identity 1 and hence Note that we again made the assumption that the MDFF is a real quantity here. Indeed, this assumption allows us to to drop the difference between causal Green's functions appearing in the diagrammatic perturbation theory (i.e., Feynman diagrams used in this section) and retarded Green's function in linear response formalism used further below. In the most general case this difference must be treated carefully, leading, e.g., to a more complicated relationship between retarded (and advanced) 4-susceptibility and the MDFF The transition probability for a particular energy loss finally reads In Sec. V we will pick up these strands and further transform this expression to obtain a generalized retarded expression of Eq. (12) noted previously. Additionally, we will show how to derive the celebrated magic angle correction [76][77][78][79] in Appendix C. In this section, we focus on the propagator for the EM field in the presence of a polarizable medium, e.g. a metallic nano-particle. For completeness, we first consider the quasi-static case and show that the screened interaction can be connected to the MDFF, see section V A 1. This way, we connect the standard formalism of EELS to optical quantities.
In section V B, using a Dyson development, we calculate the exact photon propagator in the presence of a polarizable material. In a complete analogy to what has been done with the MDFF, we then connect this photonic kernel to the charge and current density correlation functions of the scatterer.
A. Quasistatic approach: modification of the Coulomb propagator, electron density correlation function We first consider the quasi-static limit c → ∞. In this situation, the calculation of the EM field propagation simply reduces to the resolution of the Poisson equation. Thus, as illustrated on figure 5, we simply need to consider the scalar potential and the electron charge density of the target. For our purpose (i.e., to compute the screened interaction), this Green function needs to be modified.
In vacuum the potential φ ind induced in r by an external charge ρ ext in r is simply given by the Coulomb law. In other terms, the free-space EM propagator is simply given by W 0 (r, r , ω) = 1/|r − r |. This law needs to be modified in the presence of a dielectric medium in order to take into account e.g. the screening effect in the material. Particularly we expect the new propagator W to be energy-dependent as, contrary to the vacuum, the target can be dispersive.
In this section, we will derive the new propagator W and connect it to the mixed dynamic form factor.

Linear response electrostatic susceptibility
In the following, we define the electronic charge density operator n for the target as n(r) = en(r) wheren is the particle number operator for the electrons. We first need to calculate the response of the target of electronic density n(r, t) to an external perturbation φ ext (r, t) . The electronic charge δ n(r, t) ≡ n ind (r, t) = n(r, t) − n(r, t) 0 induced on the target by this electrostatic field can be calculated using the Kubo formula [80,81]: where η → 0 + , t 0 the starting time of the interaction and H(t) is the perturbation Hamiltonian given by: Therefore, one can write: Besides, the linear-response electric susceptibility χ is implicitly defined as: Comparing equations (51) and (52), one can deduce the following expression for χ: We retrieve the well-known linear response susceptibility at equilibrium in the real space. In the spectral domain, the electrostatic susceptibility reads: The latter equation is valid for any temperature as soon as we are at thermal equilibrium. We now take the limit of the latter expression in the limit of null temperature T = 0. This is fully justified when the energy of the electronic excitations are significantly greater than the thermal energy at room temperature k B T ≈ 25 meV. This will be the case in the following developments because e.g. the energy of SPs is typical energy of 1 eV. In equation (54), we can then replace 1 Z n n|.|n exp(−β ω n ) → 0|.|0 and β → 0 which gives: Therefore we see that the latter corresponds to the Fourier transform of the MDFF (16).

Connection between the mixed dynamic form factor and the screened interaction
Now, we are in position to calculate the screened electrostatic propagator W which is formally defined such that W(r, r , ω) = dr 1 dr 2 χ(r 1 , r 2 , ω) |r − r 1 | |r − r 2 | where we used the time-translation invariance of W 0 and χ to simply express the Fourier transform. The latter can be Fourier transformed with respect to r and r which leads to: We used the identity [82] F k 1 r = 4π k 2 where F denotes the Fourier transform. Nevertheless, the screened potential and the electric Green dyadic are related in the quasistatic regime by [83]: Moreover, Fourier transforming the latter gives: Therefore, using (57) we get: We also need to calculate the imaginary part of the screened interaction as it is involved in the definition of the loss probability (10). Therefore, taking the imaginary part of (56) and by using equation (55), we get: Finally using the latter equation combined together with equations (56) and (16), we get: A quick look to the latter formulate clearly indicates that, as expected, the kernel of equations (10) and (15) are the same. Depending on the situation investigated, each of these kernels can be interchangeably used: • When ab initio calculations are required (typically in the case of core-loss spectroscopy), one will preferably use the MDFF as it explicitly displays the quantum mechanical charge density correlations.
• When classical photonic systems are investigated, one will preferably use the screened interaction as it can be simply calculated by e.g. boundary element method [84,85].
Finally, by using the definition of the mutual coherence tensor (25) together with (60), we obtain: This equation shows that, in the quasi-static limit, the CDOS and the EMLDOS are respectively given by the MDFF and the dynamic form factor (DFF, [14]). In other words, in this limit, the electric field correlations (encoded in the MCT) simply reproduce the electronic charge correlations in the target (encoded in the MDFF). This result is of course expected as, in the quasi-static limit, the electric field and the charge density are simply related through: Although intuitive, equations (62) and (63) have, to the best of our knowledge, never been derived. It enables to put on the same level the MDFF formalism for the electronic correlations [15,16,19] and the MCT formalism for the photonic correlations [67,71,86].

B. Retarded approach: Photon propagator and electron four-current correlation function
We now turn to the retarded case where both the scalar φ and the vector potential A need to be considered. In his wonderful review [11], García de Abajo suggested to use the Kubo formalism for the current density to derive a retarded form of the latter equations. However, we could not find such a demonstration in the literature; therefore, in this section, we will follow this suggestion and derive a retarded version of the linear response formalism established in the last section.
The main difficulty in the retarded regime is the choice of the gauge. The developments found in the literature use different choices of gauge depending on the problem, so that a straightforward application is not possible. Some gauge choices are particularly convenient to calculate the EM field in vacuum e.g. the Coulomb gauge. However, these choices may, on the other hand, harden the calculation of the response function of the material. In order to avoid this difficulty while keeping a compact formalism, we will carry, when necessary, the calculation with four-vectors.
In vacuum, the 4-current J ν ext at x will generate a potential A µ ind at x which is related by the vacuum photon propagator D µ ν (x , x), see equation (18). In the same way as in the quasi-static regime, the presence of a polarizable material will modify the EM propagator. The goal of this section is therefore to calculate the exact photon propagator D ν µ in the presence of the target as illustrated on figure 6. However, contrary to the quasi-static case, FIG. 6: Schematics illustrating the problem tackled in this section. The Green function for the electromagnetic field in vacuum is simply the photon propagator which expression depends on the gauge choice. For our purposes (i.e., to compute the screened interaction), this Green function needs to be modified.
we need to take into account both the induced charge and current densities in the medium which are compactly represented by the 4-current density J µ .

Linear response electromagnetic susceptibility
The main interest of the previous quasi-static developments is that, now, the retarded case can be straightforwardly treated by analogy. Particularly, the 4-current δ J ν (r, t) ≡ J ν ind (r, t) = J ν (r, t) − J ν (r, t) 0 induced in the medium by an external perturbation A ν ext is given by the Kubo formula: where the perturbation Hamiltonian is given by: Substituting the latter in the former, we get: we can then define a four-susceptibility χ ν µ as: from which we deduce the linear response foursusceptibility tensor: Note that gauge invariance and four-current conservation imply ∂ µ χ ν µ = ∂ ν χ ν µ = 0 [75], a property, which we will dwell upon in the following. The structure of (69) being exactly analogue to (53), we can immediately deduce the spectral representation of the four-susceptibility at T = 0: From equation (69), one can see that the linear-response four susceptibility has the following structure: where we recall that CÂ ,B denotes the correlator between two quantitiesÂ andB being either/or ρ and j.
The diagonal elements of this tensor are therefore the charge-charge and current-current correlators while the out of diagonal elements correspond to charge-current correlators.
Let's stress an important semantic point. Both susceptibilities (53) and (69) are called retarded as they involved retarded electronic Green functions defined as (A7). Nevertheless, let's keep in mind that, in our case, the retardation needs to be understood in the sense of the EM field, the regime is therefore defined by the value taken for c. To summarize: • In the quasi-static regime (c → ∞), the problem reduces to the Poisson equation and only the scalar potential and charge densities play a role in the response of the system. The light-matter interaction Hamiltonian is then taken to be (50) and the response function of the target is determined, to the first order, by charge-charge correlations.
• In the retarded regime (c finite), both scalar and vector potentials need to be considered and the lightmatter interaction Hamiltonian is then (66).
The problem essentially reduces to a choice of gauge. If one is interested in e.g. the conduction properties of a metal, a suitable choice would be to use the temporal gauge φ = 0 where the electric field is fully determined by the vector potential E = −i(ω/c)A. In this case, the conductivity tensor defined as: can be straightforwardly obtained by the Kubo formula and gives: If one chooses a gauge where both A and φ are non-zero, both the temporal and spatial parts of the Hamiltonian need to be considered and conductivity would include charge densities in its definition.
However, although the temporal gauge seems to simplify the situation on the electronic level, it complicates the expression of the photon propagators. In fact, in our case, where both electrons and photon propagation needs to be taken into account, no gauge seems to give a dramatically simpler solution.

Retarded electric Green dyadic
We are now in position to calculate the propagator for the EM field in presence of the polarizable medium. Thus, we consider the situation described in the introduction of this section: an external source term represented by the four-current J µ ext is positioned at r and we want to calculate the total four-potential A ν tot induced at r . Like we did in the quasi-static case, we apply the Born approximation (see diagram Fig. 6) and get the retarded screened interaction: keeping in mind that there is an implicit summation on the repeated indexes.
In nano-optics, we commonly work with electric and magnetic fields so that the electric Green dyadic ↔ G is one of the most important and fundamental object of the theory. These objects (fields and dyadic) have the advantage of being gauge-independent. However, so far we worked with the potentials (φ, A) as they have simpler transformation laws and symmetries; moreover, they strongly facilitate the connection with the many-particle Kubo formalism. Nevertheless, in this section, we will derive the electric Green dyadic using the results of the previous section in order to obtain formula adapted to discuss nano-optical experiments.
Combining the definition of the Faraday tensor (A11) and the one of the photon propagator (18), we can write: where the prime in ∂ indicates that the derivative is taken with respect to x . Besides, from the explicit form of the Faraday tensor (A12), one can directly deduce that the component E i of the electric field is given by: Therefore, using (75) we get: Now, we decompose the sums over α as: (78) and, Besides, one can write the continuity equation as: Fourier transforming in time the latter gives: which finally gives: Fourier transforming in space equations (78) and (79) and using the latter equation, we get: and, We now Fourier transform (77) and get: We can inject equations (83) and (84) in the latter expression and by substituting the summation index a → j, we get the following: where we used d 4 k = d 3 k dω/c. A Fourier transform with respect to r and r then gives: Besides, for any fields Φ and V the integration by part in R 3 reads: which can be applied to the latter equation in order to get: We can now use the definition (8) in order to identify the screened Green dyadic and get: To the best of our knowledge, this equation has never been derived so far. We can also Fourier transform it back with respect to r and r and get:

Reciprocity theorem and symmetry properties of the Green dyadic
We now impose the following condition: where we define the operator S exchanging the indexes i ↔ j and coordinates r ↔ r . This reciprocity condition [11] states that a current in r creating a EM field in r is equivalent to a current in r creating a EM field in r.
In some particular situations (e.g. chiral meta-materials, moving media or topological materials), the latter condition is no longer true [87]. In this work, we only consider reciprocal media.
Therefore, in a reciprocal medium, the Green dyadic imposes (see appendix B for the detailed derivation): where the first term is a charge-charge correlator while the second term is a current-current correlator. Therefore, the Green dyadic can be written, with the susceptibility tensor components expressed in the Lorenz gauge, as: This equation is probably the most important new result of the section as it generalize the Kubo approach derived in the quasi-static case (60) to the retarded regime. Indeed by taking the quasi-static limit c −→ ∞ we obtain: which corresponds to the formula (60) that we derived in the quasi-static formalism with χ = χ 0 0 . Upon insertion of the screened Green's function (94) into the loss probability (5) we finally obtain a fully retarded version the electron energy loss probability An alternative derivation of Γ R is based on inverting (93) using the divergence-free property of the susceptibility [75]. The obtained D is then inserted in (48). Finally the charge component is replaced using the continuity equation (80). Equation (96) may be transformed to Abajo's expression Eq. (12) by replacing the screened Green's dyad Eq. (91) with the semiclassical one and noting that the transition currents read in the paraxial case. Finally, by combining equations (46) and (93), we can directly connect the mutual coherence tensor to the relativistic MDFF as:

C. Concluding remarks
In the nano-optics framework, one usually calculate the Green dyadic ↔ G as this is the sole quantity required to describe the equilibrium properties of the electromagnetic field as demonstrated by Agarwal [86] and recalled in section III B. On the other hand, from a condensed matter physicist point-of-view, the relevant quantity is the mixed dynamic form factor (or the susceptibility) of the material because it encodes all the information on the space-time dependent electronic correlations as demonstrated by Van Hove [88,89].
These two approaches are completely equivalent and based on the fluctuation-dissipation theorem which connects the response of the system (Green dyadic or susceptibility) to correlations of the underlying fields (electromagnetic or electronic correlation). The essence of section V was to explicitly show this equivalence and demonstrate that, to the first order, the two propagators χ and W (or χ ν µ and G j i ) are simply connected by two vacuum photon propagators.

VI. KINETIC EQUATION FOR THE ELECTRON DENSITY MATRIX
Following the logic of diagram (3), we will focus on the electron probe. Lets consider a fast electron described by the wave-function ψ(r, t). We can then define the single electron density matrix as: ρ(r, t, r , t ) = ψ(r, t)ψ * (r , t ) (99) In the following, we will also consider the case of a density matrix invariant by translation in time ρ(r, t, r , t ) = ρ(r, r , t − t ). In this case, the corresponding Fourier transform reads: If the Hamiltonian is time independent then the corresponding wavefunction becomes separable ψ(r, t) = ψ(r)e i t and the density matrix can be written: which correspond to the already defined above spectral one-electron density matrix. The goal of this section is to calculate the kinetic equation for the density matrix i.e. the equation ruling the evolution of the density matrix during elastic and inelastic scattering events. In the quasi-static limit, this equation has been derived for the first time by Dudarev, Peng and Whelan [26] and reads: It has then been introduced to EELS by Schattschneider and collaborators [15] and later applied to various situations such as EMCD [90], core-loss spectroscopy [22] or diffraction [91]. The goal of this section is to adapt this formula to the case of a retarded interaction kernel. As a matter of fact, apart from the final step, the derivation is essentially the same both in the quasi-static and the retarded case. Therefore, this section is organized as follows. In sections VI A, VI B and VI C, we review the seminal demonstration of Dudarev and collaborators with special emphasis on the different approximations made. The result of the derivation is a kinetic equation in the temporal domain valid for any weak interaction V . In section VI D, we use an explicit expression for V and derive the kinetic equation in the spectral domain in both the quasi-static and retarded interactions case.
To do so, we use the result of section V and assume a steady-state of illumination for the electron beam.

A. Schrödinger equation for the electron propagator
Without loss of generality, we consider a fast electron interacting with a target. The Hamiltonian of the total system {target + e − } is then given by: whereĤ e describes the free propagation of the electron, H t encodes the electronic properties of the target only andĤ int gives the interaction between the excitations in the target and the impinging electron. We now separate the interaction potential into its thermodynamical average and a fluctuating part: The thermodynamical average is taken over the ensemble of realizations of target states: employing the same notations we used in section V. This term now encompass elastic scattering and static field contribution. To simplify the notation, we suppose that Ĥ int = 0 which has no consequence on the following derivation. A non vanishing average could be absorbed by modifying the free electron Hamiltonian as: Besides, the time evolution operatorÛ 0 of the free electron as well as the time evolution operator for the total systemT follow the Schrodïnger equation: Equation (107a) can be straightforwardly integrated and gives:Û However, equation (107b) cannot be explicitly solved. To overcome this difficulty, we first define the operatorÛ as: which corresponds to the evolution operator for the interacting electron. Moreover, we define the Heisenberg representation of the fluctuating part of the interaction V as:V Combining equations (107b), (109) and (110), we get the Schrödinger equation for the time evolution operator of the interacting electron: In the next section, we will use a perturbation approach in order to calculate a good approximation of this evolution operator.

B. Dyson equation for the single electron propagator
We first integrate equation (111) in order to obtain the following integral representation: The latter equation can be solved iteratively by writing: The latter equation can be diagrammatically schematized as: (114) Let's now re-arrange the previous integrals by looking at the second term in (113) and for brevity takingÛ 0 = Id, where Id denotes the identity operator. Separating the integral in two parts and changing the integration variable leads to: The integration limit of the integrals can then be all set to t 0 and t if one introduce the proper Heaviside functions: And using the definition of the time ordering operator (A6), we finally obtain: The same trick can be applied to all order but keeping in mind that, the prefactor for the n th order term is (1/n!). It enables to re-write equation (113) as: In order to use the linear response theory derived in section (V), we now calculate the average value of the exact electron propagator Û (t, t 0 ) ≡Û(t, t 0 ): We now use the Isserlis-Wick theorem which states that for a set of Gaussian random variable {X 1 , . . . , X n }, any monomial of these variables satisfies: where Cov denotes the covariance. And since by construction the mean value ofV is zero, we have . Equation (119) then becomes: (121) where each dotted line represents a covariance product V (t i )V (t j ) . We now turn to the main approximation of this development [92][93][94]: we only keep diagrams involving correlations between neighboring vertexes. For example, we therefore neglect terms (c) and (e) in equation (121). This approximation can be interpreted in two equivalent ways: • First, as pointed out in [26], this approximation consists in treating all the successive scattering events as single independent events which correspond to the Born approximation. Following [26,95], to determine its condition of validity, we introduce the typical correlation length r c of the excitations in the target, v the speed of the traveling electron and |V | the order of magnitude of the interaction. Then this approximation holds if: In other terms, the correlation length should be short enough, or the interaction weak enough, for no dynamical effects to appear. Nevertheless, this Born approximation applies to the fluctuating part of the interaction only while the static part is included a priori. Thus, this approximation is rather a distorted-wave Born approximation [26].
• One can also interpret this approximation in a quantum field theory fashion [96] as the dotted lines can be regarded as a particle exchange. In this case, the approximation above consists in not allowing two (or more) simulataneous excitations, which is valid in the weak interaction limit. We exemplify it on the following diagram: Thanks to the approximation made, equation (121) is dramatically simplified and can be factorized as follow: (124) whereΣ is the self-energy of the probe electrons [97] and reads, in a synthetic form: Equation (125) is the starting point of Echenique's et al. formalism [55] that we will review at the end of this section. The Dyson equation (124) can be re-written in its explicit form, in the time domain, as: C. Bi-linear propagator for the single electron density matrix We will construct the propagator of the single-electron density matrix. To do so, we will use (126) to construct an average propagatorK of the exact density-matrix propagatorK. Starting from the exact electron prop-agatorÛ , one can constructK as a tensor product: Injecting the development (118) in the latter development, we obtain: In the following, for brevity reasons, we will omit the ⊗ symbol in the diagrams. As we did for the electron propagator, we now take the average value ofK. Using the Isserlis-Wick theorem, we obtain the following expression forK: (129) At this point, we will make the same approximation as in the last section and neglect all diagrams with several simultaneous excitations, e.g., diagram (d) in (129).
Diagrams like (f) correspond to coherent back-scattering events which are typically sufficiently small to be neglected [26].This approximation is the so-called forward scattering approximation and is standard in electron microscopy.
Within these approximations, the expansion contains only two building blocks: electron self-energy term (c) and mutual correlations (b). We can of course encounter sequences of these blocks like diagram (e). We can then partially re-sum the self-energy terms which leads to: (130) The latter equation formally corresponds to a Bethe-Salpeter equation in the very specific case where the two bound states correspond to ψ and ψ † and within the socalled ladder approximation. This equation can be resummed and reads: Or in its explicit form: K(r, t; r 0 , t 0 |r , t ; r 0 , t 0 ) =Û(r, t; r 0 , t 0 )Û † (r , t ; r 0 , t 0 ) + t t0 dt 1 dt 1 t t0 dr 1 dr 1Û (r, t; r 1 , t 1 )Û † (r , t ; r 1 , t 1 ) × T {V (r 1 , t 1 )V † (r 1 , t 1 )} K (r 1 , t 1 ; r 0 , t 0 |r 1 , t 1 ; r 0 , t 0 ) (132)

D. The kinetic equation for the single electron density matrix
We are now in position to derive the master equation describing the propagation of the single electron density matrix, i.e., the so-called kinetic equation. Thus let's consider an incident density matrix ρ i (r 0 , t 0 ; r 0 , t 0 ) and propagate it to the point (r, t; r , t ). Taking into account the interaction with the target and using the approximations detailed earlier, the final density matrix ρ f (r, t; r , t ) satisfies: ρ f (r, t; r , t ) = dt 0 dt 0 dr 0 dr 0K (r, t; r 0 , t 0 |r , t ; r 0 , t 0 ) plugging (132) in the latter, we finally get: ρ f (r, t; r , t ) = ρ 0 (r, t; r , t ) + dt 1 dt 1 dr 1 dr 1 U(r, t; r 1 , t 1 )U † (r , t ; r 1 , t 1 )Ĉ(r 1 , t 1 , r 1 , t 1 ) × ρ i (r 1 , t 1 ; r 1 , t 1 ) where the correlation function reads C(r 1 , t 1 , r 1 , t 1 ) = T {V (r 1 , t 1 )V † (r 1 , t 1 )} . Equation (134) is the kinetic equation in the temporal domain where the interaction Hamiltonian is not yet specified. Let's highlight that at this point, the latter equation is very general and can be applied to model, e.g., time-resolved spectroscopy experiments.
We now suppose that the electron beam is in a steady-state of illumination which is the case in the standard EELS experiment we are describing here. In this case, the density matrix only depends on the time difference. We now Fourier transform equation (134) with respect to t and t therefore taking the limits of the integrals over t 1 and t 1 to be ±∞. We therefore obtain: Changing the integration variables leads to: which can be re-written as: And we recognize a convolution product with respect to t 1 − t 1 . Noting −ω the convolution variable, we finally get: ρ f (r, r , E) = ρ 0 (r, r , E) + dr 1 dr 1 U(r, r 1 , E) We now need to calculate the Fourier transform of the correlation function. We will distinguish the quasi-static from the retarded case and denote the corresponding correlation functions withĈ QS andĈ R , respectively.

First case: Quasistatic interaction kernel
The quasi-static interactionV QS between the electron and the target is given by the Coulomb interaction: whereρ is the charge density operator for the target. Therefore, C QS reads: Writing explicitly the time ordering operator, we get: Writing τ = t 1 − t 1 , the Fourier transform reads: which gives: where P denotes the Cauchy principal value. and using the fact that for any complex number z ∈ C we have z + z * = 2Re(z), we finally get: From equation (61), one can see that Therefore, plugging it in (138), we finally obtain: which, thanks to equation (62), can also be identified to the result on Dudarev's paper (15).

Second case: Retarded interaction kernel
The retarded interactionV R between the electron and the target is given by the minimal coupling Hamiltonian: where A µ is the 4-potential associated with the excitations in the target and p µ is the electron 4-momentum operator. Note furthermore that we neglected again the (diamagnetic) e 2 A 2 -term, which is of higher order in the perturbation expansion. Moreover, within the linear response theory, the photon propagator can also be connected to the 4-potential correlation function which gives [65]: (149) where D is again the screened propagator of the EM field (taking into account the polarizability of the medium) which has been introduced in sections IV and V. By strict analogy with formula (53) and (55), we obtain: The Fourier transform of C R can be done in the exact same way as in the quasi-static case and leads to: Thus, plugging it in (138), we finally obtain: ρ f (r, r , E) = ρ 0 (r, r , E) + dr 1 dr 1 U(r, r 1 , E) U † (r r 1 , E) dω Im −D ν µ (r 1 , r 1 , ω) ∂ µ ∂ ν ρ i (r 1 , r 1 , E + ω) (152) Using equation (46), one can express the latter equation in terms of relativistic MDFF: we can now expand the sums over µ and ν with respect to the spatial and temporal coordinates as we did in equa-tions (78) and therefore obtain four terms respectively involving D 0 0 , D i 0 , D 0 j and D i j . As we developed in section V B 3 and appendix B, to the price of the hypothesis that the medium is reciprocal, one can neglect the antisymmetric terms D i 0 and D 0 j . This being so, we finally get: ρ f (r, r , E) = ρ 0 (r, r , E) + dr 1 dr 1 U(r, r 1 , E) U † (r r 1 , E) dω Im −D 0 0 (r 1 , r 1 , ω) ∂ 0 ∂ 0 ρ i (r 1 , r 1 , E + ω) The first term is a Coulomb term analogue to the second term in (146) while the second one is the retarded part of the interaction.
We now move to the temporal gauge φ = 0 where, as we explained in section III A and detailed in [64], the temporal part of the vacuum photon propagator cancels D 0 0 = D 0 j = D i 0 = 0. Using the Dyson developments (74) and the expression of the Green dyadic (90), equation (152) can be directly reduced to: All the quantities involved in the latter equation being gauge-independent, expression (155) must be valid in the general case of arbitrary gauge. Equations (146), (152) and (155) are the essential results of this section. Before concluding, we will apply them to the case of electron energy loss spectroscopy.

VII. SINGLE SCATTERING APPROXIMATION: APPLICATION TO ELECTRON ENERGY LOSS EXPERIMENTS
We now apply the previous result to the specific case of inelastic energy loss spectroscopy. We will therefore make furthers approximations: 1. The first term of the right hand side of equation (138) describes the elastic part of the interaction. As we are going to discuss EELS experiments in the following, we will not consider this term.
2. As done by Schattschneider, Nelhiebel and Jouffrey [15], we will consider a monochromatic electron, of energy 0 and density matrix ρ i , interacting a single time with the sample. It enables us to replace U by the free space electron Green's functions U 0 .

3.
As we are now interested in energy-resolved quantity, we remove the integral over ω.
Under these assumptions equation (138) reads: where we intentionally do not specify the operatorĈ in order to not lose generality as both the quasi-static (Eq. VI D 1) and retarded (Eq. 152) interactions can be used indifferently.

A. Electron energy loss probability
From equation (156), one can deduce the wave-optical EELS probability (10) and (12). To do so, we first decompose the final density matrix as (14): while the initial electron can be considered as a monochromatic pure state [15] i.e.: We multiply each side of equation (156) by ψ * n (r)ψ n (r ), which leads to: We now perform an integral over r and r which leads to: drdr ρ(r, r , f )ψ * n (r)ψ n (r ) = dx dx Since the Green function U 0 is symmetric with respect to the positions x and r, we have by definition of the electron propagator: Thus, we get: Coming back to the definition of the density operator (13), one can write: Therefore, one can write: dr dr ρ(r, r , f )ψ * n (r)ψ n (r ) = drdr m p m r|ψ m ψ m |r r |ψ n ψ n |r δ( n − f ) Using dr |r r| = Id, we get: dr dr ρ(r, r , f )ψ * n (r)ψ n (r ) = dr m p m r|ψ m ψ m |ψ n ψ n |r δ( n − f ) which leads to: dr dr ρ(r, r , f )ψ * n (r)ψ n (r ) = dr p n r|ψ n ψ n |r δ( n − f ) Replacing the latest equation in (162), we obtain: Summing over n we finally obtain: Finally, observing that ρ(r, r, f ) is the probability of finding an electron at r with the energy f , one can directly identify the integral as the total electron energy loss probability Γ(ω) and get: ReplacingĈ by either its quasi-static or the retarded form, one respectively obtain equations (10) and (12).

B. Application to coherence measurements of optical fields
In the following, we will note p f and p i respectively the wave-vectors of the final and initial electrons. The subscript z will denote the component of vectors parallel to the propagation axis while the subscript ⊥ denotes the plane perpendicular to z. The vector k correspond to the conjugate variable of r therefore indexing the reciprocal space. First of all, let's calculate the Fourier transform of equation (156) in the plane ⊥: where for brevity we omitted the energy in the argument of the density matrices. The free particle Green function reads [15]: Therefore its Fourier transform is given by [15,26]: and we moreover have The latter inserted in equation (170) gives: Now the Fourier transforms with respect to r z and r z become trivial and give: We can integrate over k z and k z as they are not observed experimentally [15]: where we used p f,z ≈ mv/ . We now consider the case of the retarded interaction and make use of the paraxial approximation but in a slightly different formulation: where L denotes the interaction length between the probe electron and the sample. Moreover the incident electron kinetic energy being principally contained in its zcomponent, one can write [11]: Plugging equations (176), (177) and the retarded form of K in (175), one gets: The integration over x z and x z gives: Using the definition of the MCT (25), one can then conclude that: which simply reads: Finally, one can come back in the real space and deduce the rather elegant formula: As we discussed earlier, Agarwal demonstrated [86], using the fluctuation-dissipation theorem, that the MCT is proportional to the electromagnetic correlation function. Thus, equation (182) shows that, when an electron is scattered by an optical field, the electromagnetic correlations are imprinted in the coherence properties of the electron beam. Producing electronic interferences thus constitutes a measurement of these correlations.
We can now connect equation (182) to the standard theory of electron holography. Indeed, during an inelastic interaction and for small scattering angles, the final and initial density matrices can be connected by the relation [61]: Here, T (r ⊥ , r ⊥ , ω) denotes a general tensor which only depends on the scatterer and the energy loss ω and is usually referred to as the mutual object transparency (MOT). In 1985, Kohl and Rose demonstrated that, in the quasistatic limit, the MOT corresponds to the MDFF [14] but, so far no equivalent relation has been established for the retarded regime. Remarkably, equation (182) constitutes the extension to the retarded case of their results and rigorously demonstrates that in this case, the MOT corresponds to the Mutual coherence tensor.
The formalism recalled or developed here is the building block of inelastic electron holography. Such an experiment can be schematized in three steps: 1. We prepare an initial electron state which density matrix ρ i (r ⊥ , r ⊥ , E) corresponds to a pure state. In standard off-axis electron holography, it simply corresponds to a plane-wave but, with modern phase-shaping techniques, it could corresponds to e.g. a vortex with a pure OAM.

2.
The initial electron state is scattered by the sample to a set of final states. After an energy loss ω and for small scattering angles, the final density matrix is given by ρ f (E − ω) = T (ω) ρ i (E) where the mutual object transparency corresponds: (1) to electronic charge correlation in the quasi-static regime or (2) to photon correlation in the retarded regime. In other words, the scattering event imprints the signature of the correlations in the target onto the beam density matrix. The final density matrix does not correspond to a pure state anymore but rather to mixed electron states i.e. a partially coherent electron beam [61].
The off-diagonal elements of the density matrix, which modulus gives the mutual coherence of the field [16], encodes the correlations in the scatterer. 3. We produce interference in order to retrieve these off-diagonal elements and therefore obtain information on the electronic or photonic correlations in the target.
Our formalism thus paves the road toward electron holography of optical field as preliminary investigated in the case of surface plasmon in e.g. [98].

VIII. CONCLUSION AND PERSPECTIVES
In this work, we have established a fully retarded formalism of fast electron inelastic scattering which can be used to described any type of TEM spectroscopy experiments such as low-loss and core-loss EELS, inelastic holography or energy-filtered 4D-STEM. Our work is built upon general response tensors including both quantum and relativistic aspects. Also, we made no assumptions on the peculiar details of the sample under consideration. Consequently, our formalism can be applied to a large set of systems and combined with any numerical methods from ab-initio to classical electrodynamics simulations. The core of our work relies on the introduction of a relativistic extension of the celebrated mixed dynamic form factor as the Fourier transform of the 4-susceptibility. By connecting this new quantity to the mutual coherence tensor -the central object of the theory of optical fields -we have drawn a formal and rigorous connection between the condensed-matter and nano-optical approaches, thus encompassing all the existing theoretical developments for EELS in an unique and general framework. Then, taking a careful consideration of most of the possible approximations, we have demonstrated that all of the approaches generally employed in the literature can be deduced from our formalism. Beyond this effort of synthesis and unification, we believe that our work paves the road toward new experiments and interpretations. First of all, as thoroughly discussed in section VII B, by introducing a relativistic form of the kinetic equation, our work enables to take into account retardation effects in holographic experiments. This constitutes the key ingredient to design and model new experiments enabling the measurement of the cross density of states [71,72] directly at the nano-scale, following an original idea of García de Abajo [11]. Moreover, our work now enables the application of all the powerful tools developed for electron holography [61,99] to the nano-photonic domain. Particularly, recent developments in differential phase contrast or ptychography for plasmonics should be described with this language. Secondly, our work bridging nano-optical and condensed-matter formalisms, we foresee a mutual and beneficial exchange of concepts between these two approaches. On the one hand, introducing the retarded screened interaction to solid state systems, one could interpret core-loss spectroscopy experiments in terms of X-ray photon exchange. This would give a natural interpretation of the already known mathematical close identity between EELS and inelastic X-ray scattering [100]. Such a comparison would be exactly similar to the standard analogy between EELS and optical extinction experiments on surface plasmons [101]. On the other hand, employing the relativistic MDFF to describe lowloss spectroscopy experiments directly paves the way toward the ab-initio modelling of EELS experiments on nano-optical systems, with potential and far-reaching applications in, e.g., quantum plasmonics [102] or excitonplasmon [103,104] and phonon-plasmon coupling physics [103,105]. Also, this should ease bridging the antagonist descriptions of EELS phonon spectroscopy experiments, either quantum [106][107][108][109] or classical, in the quasi-static [110,111] or retarded [112] regimes, or attempts to mix the two [113]. Besides, by putting the relativistic MDFF at the core of our approach, we enable the extensive use of ab-initio methods for the modelling of EELS experiment, in the same vein as the pioneering works employing density functional theory [114,115]. Finally, the introduction of the 4-current correlation function enables to employ the time-dependent current-density-functional theory [116] to model TEM spectroscopy experiments, which has never been done so far to the best of our knowledge. Such an approach would be particularly well suited to model e.g. EMCD experiments or relativistic effects in core-loss EELS. Under this convention, the raising or lowering of a spatial index changes the sign of a tensor; raising or lowering the temporal index leaves the sign unchanged. Unless otherwise specified, we will always use the implicit Einstein summation on repeated indices: x µ x µ ≡ 4 µ,ν=0 g µ,ν x µ x ν = c 2 tt − x.x (A1) The Fourier transform in M 4 is defined as: where the 4-wavevector is defined as k µ = (ω/c, k). We define the 4-gradient as: We can therefore define the 4-impulsion operator: and in presence of a EM field, one needs to perform the minimal substitution p µ → p µ − qA µ , q being the charge of the particle. Moreover, the 4-current associated with a wavefunction ψ as:

Correlators and Green functions
The time ordering operator T between two fields A(x) and B(y) is defined as: T {A(r, t)B(r , t )} = θ(t − t )A(r, t)B(r , t ) ± θ(t − t)B(r , t )A(r, t) where a + sign applies for bosons and a − sign for fermions. For a scalar field A, one can also define two different Green functions: • The retarded Green function: G R (r, r , t, t ) = −iθ(t − t ) [A(r, t), A(r , t )] ± 0 (A7) • The causal Green function: G C (r, r , t, t ) = −i T {A(r, t)A(r , t )} 0 (A8) In each case, . represents the statistical average value at thermal equilibrium and [, ] ± represents the fermion anti-correlator (resp. boson correlator).

Lagrangian form of the Maxwell equations
The four-potential defined as A ν = (φ/c, A) and the fourcurrent defined as J ν = (cρ, j) are connected by the equation of motion for the EM field: where A µ is defined up to a scalar gauge function Λ: The anti-symmetric Faraday tensor F µν is defined as: which explicitly reads: For any anti-symmetric tensor T , we also introduce the Hodge dual as: where αβµν is the Levi-Civita pseudotensor defined as: , for an even permutation of (0, 1, 2, 3) −1, for an odd permutation of (0, 1, 2, 3) 0, otherwise (A14) The Maxwell equations then read: The last equation can be derived from the Lagrange equation applied to the standard EM Lagrangian density defined as L: The first term concerns only the EM field while the second is the field-source interaction.
recall the definition of the retarded screened interaction (74): D β α (r , r, ω) = dr 1 dr 2 D β µ (r , r 2 ) χ µ ν (r 2 , r 1 , ω) D ν α (r 1 , r) (B3) Let's consider the first term of equation (90) and examine its symmetry. The relation (B2) leads to: D 0 0 (r, r ) = dr 1 dr 2 D 0 0 (r, r 2 ) χ 0 0 (r 2 , r 1 ) D 0 0 (r 1 , r ) (B4) The vacuum photon propagators can be straightforwardly reversed as D 0 0 (r, r 2 ) = D 0 0 (r 2 , r) and D 0 0 (r 1 , r ) = D 0 0 (r , r 1 ). From equation (70), the electron part can then be written: χ 0 0 (r 2 , r 1 , ω) = 2 n 0| J 0 (r 2 ) |n n| J 0 (r 1 ) |0 where Θ(ω) = 1 ω+iη . One can then see that χ 0 0 (r 2 , r 1 ) = (χ 0 0 (r 1 , r 2 )) † because the lowering and raising of 0 indexes won't bring any sign changes. Finally, we notice that ∇ j ∇ i = ∇ i ∇ j because the raising of i and the lowering of j will give both a minus sign. Indeed, ∇ j = δ j i g ij ∇ j = −∇ j because g jj = −1 by definition of the metric we have chosen. Thus, we finally have: The same arguments leads to [117]: We finally need to look at the last part of G i.e.: M i j (r , r) ≡ ∂ j D i 0 (r , r) + ∂ i D 0 j (r , r) (B8) We thus calculate: S M i j (r , r) = ∇ i D j 0 (r, r ) + ∇ j D 0 i (r, r ) (B9) = −∇ i D j 0 (r, r ) − ∇ j D 0 i (r, r )(B10) Moreover, the first photon propagator reads: D j 0 (r, r ) = dr 1 dr 2 D 0 0 (r, r 2 ) χ 0 a (r 2 , r 1 ) D j a (r 1 , r ) (B11) where we used (B2). We can again reverse the vacuum photon propagators which are obviously symmetric. However, the susceptibility term is antisymmetric. Indeed, it corresponds to a charge-current correlator and lowering the time part will keep the sign unchanged, while raising the spatial part will give a minus sign. Therefore: The above approximations are generally valid for TEM-EELS (i.e., employing fast electrons and considering energy losses well below the beam electron energy). The following two approximations to S more specifically apply to the core-loss regime. We first consider non-relativistic target atomic wave functions ξ, which are not subject to spin-orbit coupling or other relativistic corrections.
In that particular case we can approximate the timecomponent of the (relativistic) target transition current in the general expression (C2) with the electron rest mass i| j 0 (k) |f ≈ 2m i| e ikr |f (C5) Moreover, we can employ the (non-relativistic) commutator p = im Ĥ , r to replace the p z operator in the z− component of the transition current yielding i| j 3 (k) |f = imω i| ze ikr |f (C6) Last but not least we use the dipole approximation to simplify e iqr ≈ 1 + iqr and ze iqr ≈ z, i.e., only linear terms in spatial coordinates are kept in the strongly localized integrals of the atomic wave functions. Inserting all these approximations into (C1) and (C2) we obtain with S R (k, k , ω) = S QS (k ⊥ , k ⊥ , k z − k i ω/E i , k z − k i ω/E i , ω) (C8) Accordingly, we obtained a simplified approximation for the retarded loss probability in the core loss regime, which basically consists of rescaling the arguments in the non-retarded expression. Notably, the correction factor for the z-momentum (i.e., scattering angle) in the MDFF reads k i ω E i = γm e v i ω γm e = v i ω (C9) (= v i ω/c 2 in SI units) which is exactly the correction of the z−momentum obtained in Refs. [76][77][78]. The above MDFF (eventually including further approximations) is frequently used in EELS computations of core losses.