{'fraction_non_alphanumeric': 0.042275508929686526, 'fraction_numerical': 0.010412763831647309, 'mean_word_length': 4.456303236797274, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Single and double electron emission: combination of projection operator and nonequilibrium Green's function approaches

This work provides a unified theoretical treatment of the single and correlated double-electron emission from a general electronic system. Using Feshbach projection method, the states of interest are selected by the projection operator; the Feshbach-Schur map determines the effective Hamiltonian and the optical potential for the emitted electrons. On the other hand, the nonequilibrium Green's functions method is demonstrated to be a complementary approach and an explicit correspondence between both methods is established. For a self-contained exposition some results on single electron emission are re-derived using both formalisms. New insights and results are obtained for the correlated electron-pair emission: This includes the effective two-electron Hamiltonian, the explicit form of the Feshbach self-energy in terms of the many-body self-energies, and the diagrammatic expansion of the two-particle current. As an illustration of the diagrammatic technique the process of the two-particle emission assisted by the excitation of plasmons is explicitly worked out.


I. INTRODUCTION
Scattering experiments deliver the most detailed information on the structure of matter. For instance, the fully resolved spectra of an electron emitted from an electronic system upon photon or particle impact encode the spin and momentumresolved spectral properties of the sample [1][2][3][4][5]. For direct information on the two-particle properties the detection of a correlated electron pair is necessary which is usually performed in a one-photon double-electron emission [4] or in a swift particle-impact double-electron emission experiment [6]. Calculations of the electron emission spectra from atomic and molecular systems [1,4,[7][8][9] as well as from condensed matter [1][2][3] are done routinely. The underlying theories and techniques differ, however. The issue addressed here concerns the formulation of a unified and numerically accessible theoretical framework of single and double photoelectron emission (SPE and DPE) from finite and extended electronic systems. A method of choice for this purpose is the nonequilibrium Green's functions (NEGF) approach [10][11][12][13]. In full generality the response function describing electron emission is more involved than the optical response which is related to timeordered particle-hole (p*h) Green's function (GF) for which well established approximations exist. Even for a single electron emission the response function can only be defined on the Keldysh contour and after performing the calculations, the times are projected on the real observable times. The second complication is that for a fixed energy and momentum of the detected electron the sample maybe left in an excited state. A typical example is the plasmon satellites in core-level photoemission [14]. There, the target is left with one excited plasmon [15]. The conservation of energy and momentum allows to focus on, e. g., the no-loss current. The response function is then determined by the product of two vertex functions and three single-particle Green's functions [16]. If an approxima- * yaroslav.pavlyukh@physik.uni-halle.de tion is made for one of the constituents, it has to be taken over consistently to the others. The notion of a conserving approximation is rooted in this requirement.
First theories of electron emission were empirical: E.g. for surfaces, following Berglund and Spicer [17] the photoemission is regarded as a three stages process: excitation, transport to the surface (during this stage the particle may loose energy), and the transformation into a scattering (detector) state. In 1970 Gerald D. Mahan wrote "we have not yet been able to derive a simple, time-ordered, correlation function which would serve as the starting point for a closed-loop type of calculation. That is, we have not yet found a "Kubo formula for photoemission."" [18]. Shortly thereafter Schaich and Ashcroft [19] and Langreth [20] employed a time-ordered formalism for the response function, and Caroli et al. [21] introduced the nowadays standard NEGF formulation. The well-known Fermi Golden rule expression for the photocurrent the Auger effect [25][26][27][28]. In both cases the system is left in an excited state that relaxes subsequently either due many-body effects or results in the emission of a secondary electron. It should be noted, however, that the borderline in such a classification is blurred: one can consider the Auger effect as a two-step process, in which the decay is treated independently from the primary ionization or as the no-loss double photoemission [29]. The former point of view yields a description of the Auger effect in terms of an equilibrium two-hole Green's function [26,30,31]. The goal here is to generalize the nonequilibrium approach as to treat single and double electron emission. We will mostly discuss processes related to the absorption of one photon. Particle impact is discussed only in the optical limit as specified in the Appendix A. In particular, this work provides a detailed discussion of DPE, a process that was experimentally realized for various systems [4,32]. For a self-contained presentation we start by defining observables and introducing basic formulas solely based on the time-dependent perturbation theory and the assumption of adiabatic switching of the light-matter interaction (Sec. II A). Already on this level one can reformulate these expressions in the Fermi golden rule form and demonstrate how the sudden approximation can be used to reduce the many-body to two-body description (Sec. II B). Such reduction, however, neglects the energy loss of an emitted electron on its way to detector. These extrinsic losses are treated by means of the projection operator technique (Sec. III). For single photoemission (SPE) this approach was established in works of Almbladh [16], Bardyszewski and Hedin [33], Fujikawa and Hedin [34], Hedin, Michiels and Inglesfield [35], and for DPE by Brand and Cederbaum [36]. The notion of the optical potential is central to this approach. While the case of elastic scattering was considered in a classical work of Bell und Squires [37], the inelastic case, which is especially relevant for photoemission, is more involved and has a long history with a recent progress due to Cederbaum [38,39]. In Sec. IV we closely follow the derivation of Almbladh and extend the theory to the two-electron case. There are important differences as compared to the singleelectron emission. Under some assumptions DPE is only possible for interacting systems [40]. We demonstrate that the vertex function is the source of this electronic correlation effect. Finally, we corroborate our findings by performing a diagrammatic expansion of the derived DPE response function in terms of Green's function on the Keldysh contour (Sec. V). We consistently use atomic units.

II. THE TWO-ELECTRON CURRENT
For DPE from atomic and molecular systems [41,42] a variety of very successful techniques, based on a full numerical solution or using approximate correlated scattering states of the few-body Schrödinger equation, were put forward. The wave-function-based methods and, consecutively, the scattering approach are less suitable for extended degenerate fermionic systems. Such DPE experiments were first performed for Cu(001) and Ni(001) crystals [32] and mean-while for a variety of other samples. Here comes the response formalism into play: the expectation values of products of the creation and annihilation operators are computed over the ground state of a (many-body) system, and perturbative expansions are evaluated with the help of Wick's theorem. If the studied process can be regarded as a multi-step event, then the rate equations are often a very efficient tool. They can be derived either from the density matrix or from the NEGF formalisms using some additional assumptions. For instance the generalized Kadanoff-Baym Ansatz has been used to derive the quantum master equations starting from NEGF approach to describe the transport in molecular systems [43].
Here we present a self-contained derivation of the twoparticle current starting from the time-dependent perturbation theory. The resulting formula (Eq. (12)) is, however, less useful for practical applications because it requires (generally unknown) many-body states. One has either a choice to completely neglect the target-ejected particles interaction which still might be relevant for higher energies (Sec. II B), or, as will be demonstrated in the next section (III) to properly reduce the formulations as to work with effective residual interactions (i.e. optical potentials).
A. Basic definitions a. Hamiltonian: A system of interacting fermions is considered that has the Hamiltonian where the field operatorψ (ψ † ) with argument x ≡ (r, σ) annihilates (creates) a fermion in position r with spin σ. Needed below is the anti-symmetrized interaction One may wish also to change the basis for the representation of creation and annihilation operators viâ where the sum runs over a complete set of one-particle states and we consistely skip· · · on c i and c † i . To study photoemission we need to further classify the states according to their geometric character. A state will be called bound (φ i ∈ B) if for any > 0 there is a compact set B ⊂ R 3 such that for all times t the state remains in B: χ B c e itĤ φ i < , where B c is the complement of B, χ B c denotes the corresponding characteristic function. Analogically for the scattering states (φ k ∈ C) we adopt the following definition: they are the vectors for which lim T →∞

2T
T −T χ B e itĤ φ k dt = 0 for all compact sets B ⊂ R 3 , i.e. they leave any bounded region. It is clear that B ⊥ C and according to the RAGE theorem [44] all the states from the discrete (point) spectrum are bound, whereas the continuum states (absolutely continuous and singularly continuous) are the scattering states. Thus, parallels between the geometric and the spectral classification allows us to use continuum and scattering, and point and bound terms interchangeably, although for the purpose of the present work the geometric classification is preferred. Finally we note that if our theory is to be applied to solids the use of localized Wannier functions [45] is preferred, at least for systems where their existence can be proved [46].
We will use the letters (abcd) for general orbitals, (i jnm) for bound orbitals and bold-face letters for continuum states. In these notations: b. Initial state preparation: The above Hamiltonian determines the quantum state of the target (wave-function |Ψ 0 with corresponding energy E 0 ) in the remote past (t = −∞). When the system is perturbed by the interaction with external fields it evolves to a new state. As a typical mechanism we consider here the light-matter interaction In this expressionV(t) is adiabatically turned on allowing to introduce a typical interaction time ∼ (2η) −1 . The form (6) permits generalizations: In Appendix A we consider the process of impact ionization caused a charged projectile particle (e.g. an electron) impinging on the target system. At high energy the projectile can be regarded as distinguishable from electrons of the system. This allows to average the projectiletarget interaction over the projectile's states and write the perturbation in essentially the same form as in Eq. (6), i.e. as a single-particle operator. From the first-order time-dependent perturbation theory we obtain the approximate eigenstate |Ψ (+) of the full Hamilto-nianĤ +V(t) at time t = 0: Readers will immediately notice parallels of Eq. (7) with the scattering theory where the Møller operatorsΩ (±) convert an eigenstate ofĤ (the Hamiltionian of the target system) at t = ∓∞, into an eigenstate ofĤ +V(0) (the full Hamiltonian) |Ψ (±) α =Ω (±) |Ψ α at time t = 0 (cf. Eqs. (14.66) of Joachain [47]). The scattering theory is required when electromagnetic fields are quantized. For classical fields Eq. (7) follows from the first order expansion (in∆) of the Møller op-eratorΩ (+) . To emphasize the similarity we denote the state given by Eq. (7) as the scattering state. In what follows we omit the tilde which we used to denote its approximate character.
c. Observables: Assuming we know the quantum state of the target at t = 0 some observables can be computed. Since we are interested in photoemission these are the expectation values of the current operators. The safe way to introduce them is to use the continuity equation which is gaugeinvariant. The one-electron current J k is defined as the number of electrons N k with a given momentum k outside the target divided by the effective interaction time (2η) −1 . There is a detailed discussion [16] on why electrons in the sample give a negligible contribution to the current. Same arguments are valid for the two electron case. Thus, we analogically define the two-electron current as In the expression above (and all subsequent derivations) we do not explicitly spell out the spin quantum numbers. The dependence on the spin can be recovered by substituting the continuum quantum numbers like k by kσ (likewise for bound indices). The second term excludes the one-electron current in the case when two momenta are equal. Eq. (8) gives access to the differential cross-section through the following relation: where I/ω is the photon flux density [48]. For the velocity gauge∆ = 1 , where A 0 is the amplitude of the vector potential. Similar expressions can be given for the length gauge.
The average in Eq. (8) is performed over the perturbed state (7): where we used the usual anti-commutation relations for the fermionic operators. The current is quadratic in∆ or linear in the number of absorbed photons. The first order in∆ gives the linear conductivity current and is of no interest here. [21] To derive the Fermi golden rule for DPE we insert a complete set of the (N − 2)-particle states and use the scattering theory to evaluate matrix elements of the type: We will generally use lower indices to distinguish quantum states and upper indices to indicate the charge of the system or the nature of the state (±), i.e., incoming or outgoing. For a scattering process with the following energy balance the Møller operatorΩ (−) translates a wave-function in the remote future into a incoming (they are sometimes called inverted LEED states [35]) scattering state at t = 0: Following Almbladh [16] we obtain: resulting in the Fermi golden rule for DPE for an adiabatic switching ofV(t): This is essentially an exact equation if strong field effects are neglected, i.e. if the first-order perturbation theory in field strength is adequate. Now we discuss some common approximations. In the sudden approximation the Møller operator is set to the identity operator and it follows |Ψ (−) leading, e.g., to Eq. (1) of Napitu and Berakdar [49]. The sudden approximation is broadly used to interpret the single photoemission. However, it is easy to construct an example when it completely fails: Consider photoemission from a system surrounded by a impenetrable potential barrier. Irrespective of the photon energy there will be zero current in the detector. Thus, it is extrinsic losses [35] that are missing in the sudden approximation.

B. Sudden approximation
In the sudden approximation for SPE it is possible to reduce the many-body description to a single-particle picture which also allows to approximately treat the Møller operator and accommodate extrinsic losses. The central object in such an approach are the Dyson orbitals [50]. The hole Dyson orbital is defined as an overlap of (N − 1) many-particle state with the N-particle initial state: A rather extensive review of such overlap operators as well as the proof on the last "dressed in the fancy outfit of the occupation number formalism" identity can be found in Ref. [51]. Practical approaches for their computation are overviewed in Refs. [52,53]. By introducing a similar two-hole Dyson orbital: and neglecting the Møller operator we obtain for the twoparticle current (12): where |k 1 k 2 is asymptotic two-particle state, i.e. antisymmetrized product of two plane-waves. The two-hole orbital is anti-symmetric with respect to the interchange of particle coordinates and in general has norm ≤ 1. To derive (15) it is instructive to consider first a corresponding matrix element for SPE: and it is time to make another very important assumption: It is not valid in general, however, one can use the same arguments as Almbladh (see discussion around his Eq. (11)) to demonstrate that it gives a vanishing contribution. For homogeneous electron gas this is even a generally valid statement. Besides allowing to compute the matrix elements the assumption (16) also justifies why terms resulting from the secondorder perturbation theory give vanishing contributions to the current.
In this way (see Appendix D) For DPE we analogically analyze the matrix element entering Eq. (11) and neglect terms with two holes at momenta k 1 and k 2 (i.e. c k 2 c k 1 Ψ 0 ≈ 0) as compared to the terms with only one hole (Appendix D). Notice that for SPE we neglected one hole term as compared to zero hole contribution (cf. Eq. (16)). It is obvious that the sudden approximation is only valid for large momenta k 1,2 and it is indifferent to the state in which the system is left in (the final double ionized state can be an excited state). Thus, it is desirable to generate improved approximations to Eq. (12) by rewriting it in the two-particle form, but with an improved final state (such as Eq. (4) of Fominykh et al. [54] or Eq. (2) of Fominykh et al. [55]).

III. EXTRINSIC EFFECTS
A many-body target interacts with light such that certain number of electrons is emitted. Here, the fundamental question is whether it is legitimate to describe the process in such a way that only quantum numbers of ejected particles are considered and remaining degrees of freedom are traced out, i.e. put into some effective interactions. The projection operator formalism is a general method to treat this kind of problems. In this section we introduce the basic concepts of this theory and demonstrate the reader that a deep connection with the nonequilibrium Green's function formalism exist. We conclude this rather mathematical section by considering two examples. Based on these examples the Fermi golden rule is derived in the subsequent section.

A. Nonequilibrium Green's functions
In the Keldysh formalism [13] the field operators evolve on the time-loop contour C shown in Fig. 1. Operators on the minus-branch are ordered chronologically while operators on the plus-branch are ordered anti-chronologically. Letting z 1 and z 2 be two contour-times, the Green's function G(x 1 z 1 , x 2 z 2 ) can be divided into different components G αβ (x 1 t 1 , x 2 t 2 ) depending on the branch α, β = +/− to which z 1 and z 2 belong. As before, x i , denote a composite coordinate comprising space and spin variables. For α = β = − we have the time-ordered Green's function In this expression the average . . . is taken over a given density matrixρ and T is the time-ordering operator. The subscript "H" attached to a general operatorÔ signifies that that operator is in the Heisenberg picturê whereÛ(t 1 , t 2 ) is the time-evolution operator and t 0 is an arbitrary initial time. Reversing the time arrow the G −− is converted into the anti-time-ordered Green's function whereT orders the operators anti-chronologically. Finally, choosing z 1 and z 2 on different branches we have The last two components are equivalently written as G −+ = G < (lesser Green's function) and G +− = G > (greater Green's function), and describe the propagation of an added hole (G < ) or particle (G > ) in the medium. It is often convenient in addition to time ordered and antiordered functions to introduce the retarded and advanced components: In order to find their representation in frequency space we multiply the retarded GF by e −ηt with η → 0+ in order to enforce the convergence and compute the Fourier integral: The forward branch is denoted with a "−" label while the backward branch is denoted by a "+" label.
Let further introduce (for general z ∈ C) the particle-type and hole-type GF by From Eqs. (23) follows Finally, let us present the equation of motion (EOM) for the retarded GF in the form: The two-particle Green's functions are much more diverse. However, we will only need those containing creation operators with the same time argument and the same holds for annihilation operators. To specify the relative order of creation (or annihilation) operators infinitesimally small times are added. Because such Green's functions depend on two times only, the same nomenclature as in the single-particle case can be used. Thus, we define They are the constituents of the retarded and advanced twoparticle Green's functions: For the retarded function the following equation of motion can be derived:

B. Two projection operators
In the previous section we have seen that relevant types of Green's functions can be written in the form of a resolvent To be more specific about the state over which the averaging is performed we select from all possible states of the target and emitted particles the relevant ones for the effect of interest by employing projection operators. In the following we consistenly skip· · · when writing these operators and use 1 to denote the identity operator. Hence P + Q = 1 are two complementary projection operators with the idempotence (P 2 = P, Q 2 = Q) as their defining property and the basis formula for computing resolvents whereĤ P = PĤP,Ĥ Q = QĤQ, and the self-energy operator is defined as:Σ The map F p :Ĥ →Ĥ p +Σ P (E) is called the Feshbach-Schur map, it relates the eigenvalue problem on the full Hilbert space and to that on its subspace. We summarize relevant matrix identities in Appendix C. Due to the presence of the bath HamiltonianĤ Q in Eq. (27) this definition cannot be used for practical computation of the self-energy. Fortunately, a connection with the many-body perturbation theory (MBPT) exists [56,57]. If, for example, starting from the N-particle wheren(r) is the hole-density of ionized state α, i.e. n(r) ≡ Ψ + α ψ (r)ψ † (r) Ψ + α , the eigenvalue problem on the Psubspace (C3) ( Ψ + α ψ(r)(Ĥ P +Σ P (E) − EÎ P )P Ψ 0 = 0) is the Lipmann-Schwinger equation for the hole Dyson orbital (13). Notice thatĤ P contains the electrostatic and exchange part of self-energy, whereasΣ P (E) → 0 for E → ±∞. Similarly, in 1959 Bell and Squires [37] considered a one-body potential for the scattering of a particle incident on a complex (manybody) target. They demonstrated that this optical potential is exactly given by the sum of all proper linked diagrams, i.e. many-body self-energy in the time-ordered formulation. In fact, their Eq. (7) directly corresponds to Eq. (C3) when P is a projection yielding a particle Dyson orbital.
In order to study single and double photoemission we introduce two special projection operators. The main goal of this section is to establish an equivalence between the abstractly defined self-energy (Eq. (27)) and the self-energy of the manybody perturbation theory. We consider the expression appearing in the first line of Eq. (26) i. e. resolvents of the type We will demonstrate that the formalism of nonequilibrium Green's functions is easily paralleled with the Feshbach projection algebra (FPA). The basic relation for the subsequent derivations are the operator identities We will show below that witĥ the operator identity (28) has a structure of the Dyson equation for certain Green's functions. For SPE we consider the projection operator where the sum runs over scattering states. It is common to select these single-particle states |ϕ k to be eigenfunctions of some reference Hamiltonian with proper boundary conditions. We request that |Ψ + α is a completely bound remainder of the ionization event and does not emit a second electron at a later stage (Auger electrons are a typical example for these kind of processes). There are many equivalent ways to impose this restriction, for instance we will assume i.e., implying |Ψ + α is a vacuum state for photoelectrons. From the assumption follows the indempotency (P 2 α = P α , see Appendix D for proof) and, thus, P α represents a true projection operator. The application of P α restricts the possible processes which might occur upon excitation to the definite emission of one photoelectron, whereas the ionized system is left in a (possibly excited) bound state |Ψ + α . From the assumption Eq. (31) follows another restriction: where the first term is equal to zero because each bound state (i) is necessarily given by a square integrable function (converse is not true). In the following we will use another consequence of assumptions Eqs. (31,32): lim The projection operator for DPE we define as Here, |Ψ 2+ β is the doubly-ionized reference state, to which two photoelectrons with continuum quantum numbers p and p are added. We can easily show the indempotency of the projection operator (35) if we require, similar to Eq. (31), With the definition of the SPE projection operator P α in Eq. (30), we find where we applied the definition of the particle-type GF Eq. (23a). Note that the GF is defined for a particular subspace spanned by the operator P α and should therefore always be understood as the GF associated to |Ψ + α . For brevity, however, we omit labelling GF by α.
Using these notations the operator identity reads With the help of our assumption Eq. (31) we can now remove the sum by applying Ψ + α |c p from the left and c † q |Ψ + α from the right as Eq. (31) As stated above, we can think of |Ψ + α as a vacuum state for free particles (cf. Eq. (31)). The hole-type GF is identically zero. Therefore, (37) we realize its equivalence to Eq. (24). In other words, by applying the FPA we can derive EOM for the retarded Green's function.
e. Effective Hamiltonian: In Eq. (28a)Â −1 plays the role of the reference Green's function. Correspondingly, PĤP is the effective Hamiltonian. Using the standard anticommutation algebra and the assumption (31), we find i.e. it consists of the total energy of the ionized system and the Hartree-Fock Hamiltonian for continuum states. The latter is computed with the density matrix of the target: Letĥ be an operator acting on the subspace of continuum states with matrix elements given by Eq. (38). Its resolvent relates to the reference retarded GF as g R pq (ω) = g (p) f. Self-energy and the Dyson equation: The second correlator in the EOM (37) amounts to With Eq. (29) inserted into the identity Eq. (28a) we apply P α from left and right, use the same trick to multiply with suitable states from left and right, and find With z = E + α + ω + iη Eq. (41) has a structure of a Dyson equation for the retarded Green's function in the subspace of continuum states: The second sum runs over the full set of orbitals (bound and continuum). This is the most general form and without additional analysis it cannot be reduced to the Dyson equation with the self-energy from the projection formalism (cf. Eq. (27)). Let us compare Eq. (41) and Eq. (42). At first we notice that Eq. (39) defines the reference Hamiltonian only on the subspace of scattering states. We might extend the definition and request, for instance, that all the basis functions (bound and scattering) are the eigenstates of the reference Hamiltonian. This impliest pq = ε p δ pq andt nq = 0. Thus, mean-field terms of the Hartree-Fock Hamiltonian are then cancelled by the frequency independent part of the last correlator in Eq. (41). In the case when the reference Hamiltonian is not diagonal in the chosen basis the embedding self-energy terms additionally appear. In the simplest case (no interaction), they can be written as Σ em pq (z) = mn t pn g nm (z)t nq . Let us now assume that the single-particle basis is such that no embedding selfenergy appear. What would be the diagrammatic structure of (a) we determine the Green's function in this sector (G bc ) and substitute in Eq. (42): where for brevity the subscripts b and c denote the bound and the continuum sectors. Expression in square brakets (Eq. (44)) can now be compared with the self-energy from the projection formalism (27). Notice, that the reference Green's function was assumed to be diagonal, i.e. g b ≡ g bb and g bc = 0. g. Dominant scattering mechanisms: Let us recapitulate what led us to Eq. (44). We have chosen a projection operator in the form (30). This specifies the state of a system after the photoionization as containing one photoelectron in the scattering state plus the bound ionized target. Next, we obtained an effective Hamiltonian (38) acting on the P subspace and used it to define the reference Green's function (40). We want to understand what is the diagrammatic content of the Feshbach self-energy (27). It is not possible to use this equation directly because it involves the effective Hamiltonian on the complementary Q-subspace. However, it is possible to use another matrix identity (28a) and to formulate the Dyson equation for the full Green's function in the P subspace (41) avoiding the use of the QĤQ resolvent. This equation can be put in a direct correspondence with the Dyson equation for the retarded GF from the many-body perturbation theory. The difference between them is the domain where the selfenergies are defined: the Feshbach self-energy operates on the continuum sector only, whereas many-body perturbation theory does not impose such a restriction. By writing another Dyson equation (43) in the bound-continuum sector we can finally obtain the Dyson equation with an effective self-energy in the continuum-continuum sector. This self-energy is an exact counterpart of the Feshbach self-energy (27). To the best of our knowledge it is the first explicit example of such correspondence. Critical for our derivation was the choice of the single-particle basis. We have demonstrated that it is the projection operator that determines the effective Hamiltonian, and if the basis is such that the Hamiltonian is diagonal the embedding self-energy vanishes and one arrives at Eq. (44). No further assumptions have been made and Eq. (44) is so far exact.
Let us analyze the meaning of different terms of the photoelectron self-energy (Fig. 2). As discussed in details by Bardyszewski and Hedin [33], Almbladh [16] and Fujikawa and Hedin [34] scattering states vanish in the sample (damped) represent the real photoelectron states more precisely. One can derive explicitly the residual interaction that they experience. The reasoning is easier to perform in real space where the Coulomb interaction depends on two coordinates only (cf. Eq. (2)) as opposite to the Coulomb matrix elements which are four index quantities. Since the scattering states are damped in the sample, there are only two nonvanishing Green's functions G vv and G VV operating exclusively in the inner (v), outer (V) spaces, respectively. The Green's function starting in the sample and ending outside of it (G Vv ) and the reverse (G vV ) vanish. We can rewrite Eq. (44) in these new notations, however, it is not even necessary as it amounts to mere replacement b → v and c → V. What has changed is the interaction lines in the diagrammatic expansion of the self-energy. They can connect v and V domains and generate therefore nonzero contributions. It is easy to see, however, that the second self-energy term vanishes: a diagrammatic expansion of Σ vV necessarily contains at least one g vV line which is zero according to our assumption. Thus, only Σ VV needs to be analyzed. By explicitly forbidding the particle exchange with the sample we arrived exactly at the case of elastic electron scattering considered in the seminal paper of Bell and Squires [37]. We will see below that the structure of Σ VV is quite general and appears in the diagrammatic consideration of other processes, remarkably, in the parquet diagram treatment of the Fermi edge singularities [58]. There, however, a similar diagrammatic expansion arises due to the specific choice of the interaction between the deep hole (labeled by m) and the conduction electrons:Ĥ 1 = kk V kk c † k c k c m c † m . In contrast to their work, what induces a special structure of diagrams for Σ VV is not a specific form of the interaction matrix elements, but rather the absence of the off-diagonal blocks in g. It is easy to construct the electron self-energy fulfilling these restrictions: it consists of one open photoelectron line (depicted as solid line on Fig. 2 and a number of closed bound electron loops (depicted as dashed lines). Because of the restriction (33) there are no photoelectron loops.
The topic of the present section is quite extensive and such an aspect as the Lehmann representation of the Green's functions mentioned here was completely left out of our discussion. This is, however, very relevant for the treatment of finite systems, with important recent progress, e.g., [59].

D. Example of DPE
h. Equation of motion: The derivation for the twoparticle case goes along the same lines. We insert the definition of the projection operator (Eq. 35) in the identity , and as for SPE compute the matrix elements of the whole expression. The final results read as The prefactor 1/4 originating from the product of two projection operators is cancelled because of the symmetries of the particle-particle GF and of the second term on the right-hand side of Eq. (45): Inserting z = E 2+ β + ω + iη shows the equivalence of Eq. (45) to the equation of motion (25).
i. Effective two-particle Hamiltonian: Analogically to the SPE case we consider the Feshbach-projected Hamiltonian in the subspace defined by P β and describing two electrons including their interaction and their mean-field interaction with the ionized system: where the last term can be expressed as follows The first correlator in the square brackets evaluates in terms of the density matrix with respect to |Ψ 2+ β with bound state indices to: Here we have written it in terms of the matrix elements of the anti-symmetrized Coulomb interaction (2) V abcd ≡ abcd − abdc . Similarly, the second correlator is obtained from this expression by the index exchange p ↔ p . The effective twoparticle Hamiltonian (47) is so expressible as a Hartree-Fock Hamiltonian (39) for two independent electrons plus the interaction (Fig. 3): j. Kernel and Dyson equation: We return to the matrix identity (28a) and insert the splitting (29) with P = P β (Eq. (35)): and define the reference two-particle GF Invoking again the symmetries (46), which also hold for the reference GF, and applying the same states from left and right, we obtain nn qq (z) . (51) It is instructive to divide the kernel entering the equation of motion (second line of Eq. (45)) or the Dyson equation (second line of Eq. (51)) into the terms containing higher correlation functions and those expressible in terms of two-particle GFs: The latter gives rise to the particle-particle embedding selfenergy. We can now formally introduce the correlated frequency-dependent and the static kernels: The static part is exactly cancelled by the density-dependent part of the effective Hamiltonian: The embedding self-energy originates from the kernel as well as from the effective Hamiltonian (49):  × K em kk nn + K c kk nn (ω) G R nn qq (ω) . (55) Eq. (55) has a form of the Dyson equation for the two-particle Green's function, however, the reference GF g R pp qq (ω) is not given as a product of fully-interacting single-particle GFs, but rather is the full two-particle GF -the resolvent of the effective Hamiltonian (47) which includes the full electronelectron repulsion and the mean-field contribution from the ionized system.

IV. FERMI GOLDEN RULE
A. Single photoemission SPE was treated by several authors. We recapitulate the main points. The total observed current is proportional to the expectation value of the electron number operatorN k = c † k c k . Out of all possible final states of the target we discard all unbound states, i.e. c k |Ψ + α = 0 and choose only those relevant for a specific experiment. Let λ α be a corresponding distribution function. For instance when the target is left in the ground state we can set λ 0 = 1 and λ α = 0 for all excited states. Modified particle number operator for this process reads: The same expression can be obtained from the Langreth approach starting from the Wigner distribution function [20].
Let now the SPE current be the expectation value of this operator We only consider the case where we neglect the off-diagonal term in Eq. (26) and de-fineΣ (±) P (ω) =Σ P (ω ± iη). We omit the subscript α where it does not cause a confusion. A simple calculation leads to the modified matrix element Using the same assumption for the computation of the matrix element of∆, Ψ + α |c p∆ |Ψ 0 = p|∆|φ α and the definition of the Green's function on the P α subspace: we obtain for the current where ε α = E 0 − E + α . As shown in Appendix B we can express the particle Green's functions in terms of Møller operators This finally leads to the current A standard definition of the spectral function entails tô Therefore, we can recast the expression for the current in a more familiar response form where the tilde denotes a spectral function with restrictions imposed by the weighting factors λ α and µ is the chemical potential, or in the Fermi golden rule form: The major distinction from other approaches is that both, initial and final states are dependent on the final state of the target α. Formally, |χ (−) k,α is the incoming scattering state of an electron in the optical potential of the ionized target in the state |Ψ + α . Notice that the current has been obtained using the approximation (57). Exact calculation leads to the appearance of the vertex functions that describe a screening of the optical field by the electrons of the target [16]. We will stop on this point when treating DPE process.

B. Double photoemission
The total observed current is given in terms of the expectation value of the electron number operatorsN k 1 k 2 =N k 1N k 2 − δ k 1 ,k 2N k 1 , viz. Eq. (8). Out of all possible final states of the target we discard all unbound states, i.e. c k |Ψ 2+ β = 0 and introduce weights λ β selecting the relevant ones. The modified observable reads: This allows us to improve upon Eq. (15): Using assumption (57) Eq. (62) can be written in the Fermi golden rule form with a modified matrix element Using the matrix elements of∆, Ψ 2+ β |c q c p∆ |Ψ 0 = pq|∆|φ (2) β (cf. Eq. (D8)), and the properties of the two-particle Green's functions (Appendix B) we finally obtain for Eq. (10) is the negative of second ionization potential, |ψ (−) k 1 k 2 ,β is the incoming damped two-electron scattering state in the optical potential of doubly ionized target andÂ (2) (ζ) is the two-particle spectral function, which can be written in terms of two-hole Dyson orbitals: with ε (2) β = E 0 − E 2+ β . Notice that the current has been obtained using the approximation (57). Exact calculation leads to the appearance of the vertex functions resulting from Q β∆ |Ψ 0 and describing a screening of the optical field by the electrons of the target [16].
In the valence shell the DPE mechanism is typically due to ground state electron correlation, i.e. due to the correlated two-particle spectral function entering (64). In contrast, when core electrons are involved a dominant mechanism for DPE is due to the final state relaxation (so called shake-off). Multiple stages are then described by introducing corresponding projection operators for each intermediate stage. In the following, we focus on the diagrammatic approach because it allows us to treat all these effects on equal footing.

V. DIAGRAMMATIC APPROACH
Treatment of the off-diagonal part of the Hamiltonian resolvent is the main difficulty of the Feshbach projection algebra. It is even more aggravated in the two-particle case. The diagrammatic technique provides a natural and practical solution to this problem.

A. Derivation
Eq. (10) when transformed to the time domain gives rise to the following ground state correlator: where the field operators are in the Heisenberg representation and t, t ∈ (−∞, 0] are physical times. For clarity, we omitted the indices in the notation of the correlator. It can be evaluated diagrammatically by adiabatically switching on the interaction in the remote past, i. e.Ĥ δ =Ĥ 0 + e −δ|t|Ĥ 1 . Now the average is performed over the noninteracting ground state |Φ 0 and the times t − 2 ≺ t + 1 lie on forward, backward branches of Keldysh contour γ (Fig. 1), respectively: T here is the usual contour ordering operator [13] with the order relation ≺.Ĥ δ is such that it is equal to the Hamiltonian of noninteracting system H 0 in the remote past and is identical toĤ at t = 0. Notice that it is different from adiabatic switching on of the electromagnetic field in Eq. (7). |Φ 0 is the ground state ofĤ 0 . Using Wick's theorem we can contract the product of field operators in order to express the correlator in terms of products of single-particle Green's functions. Zeroth order obviously yields four fermionic lines. However, if we use the same assumption as in Sec. IV B any zeroth order diagram vanishes. This is easy to understand by comparing with SPE case. There, no-zero contributions are coming from the following contraction: This is the only combination that results in greater GFs when one of the arguments is a scattering state (and is compatible with (33)). In particular, the above contraction equals to In DPE two creation operators with continuum state indices need to be contracted with two annihilation operators on the positive track. However, there is only one such operator. Hence, 0th order in interaction is zero. The argument that excludes the first order diagram is slightly different and is based on the fact that bare interaction is instantaneous, i.e. corresponding time-arguments necessarily lie on the same, positive or negative, track. Second order nonvanishing contributions contain products of two Coulomb interaction operators (e.g. at contour times t + andt − ) and already a familiar product of six operators as in Eq. (67). From all possible contractions (they yield eight fermionic lines) we have to exclude many terms. Some of them immediately vanish because of the assumption (33) for noninteracting GF. Others, represent the Hartree-Fock renormalization of two fermionic lines and likewise vanish because of the same assumption for the full fermionic propagators, Fig. 4 (a). Then, there are diagrams ( Fig. 4 (b)) containing isolated islands of pluses and minuses which also vanish because otherwise the two-particle current cannot be written in the Fermi Golden rule form [60,61]. Finally, there are only four (times two for exchange) nonzero diagrams. Two of them are depicted at Fig. 4 (c,d).
It is clear now how more general diagrams for the twoelectron current can be constructed: i) One replaces all bare fermionic propagators and interaction lines with the dressed ones; ii) Each pair of parallel fermionic lines are replaced by the corresponding two-particle propagator, Fig. 5 (a). In doing so one obtains, in principle, diagrams given by Fig. (1b) of Fominykh et al. [54] with a small correction that zeroth and the first-order two-particle GF should be excluded from the vertical track; iii) Next class of the diagrams are those that describe the screening of the optical field, Fig. 5 (b); iv) Processes involving intrinsic or extrinsic losses are given by the diagrams with interaction lines connecting points on different tracks, i. e. "+ −", "+ 0", "− 0". They cannot be obtained by the renormalization of fermionic or bosonic propagators, one such example shown at Fig. 5 (c) reveals a process with extrinsic losses.
Finally, we give a description of a general diagram for a photoemission process. Examining SPE and DPE diagrams we see that all of them are constructed from the common ancestor: the density-density response function χ < ≡ χ −+ having a form of two islands with time arguments belonging to either forward or backward tracks of the Keldysh contour. Now we introduce detectors (shown as black squares at Fig. 5 (d) measuring J k 1 ,k 2 . As explained before i) the lesser GF with one of the indices being a continuum state vanishes because of the assumptions (31,32); and ii) observation is made at the rightmost point of the contour (i. e. at t − = t + = 0 in our notations), thus, each detector measuring particle numbers N k i is connected to two greater GF. In view of this, the detectors "lie" on the fermionic lines flowing from the "−" (forward track) to "+" (backward track) islands. Each response function constructed in this way has an important property that it can be represented in the Fermi Golden rule form, such construction obviously generalizes to an arbitrary number (n) of emitted particles. Simple counting shows that these processes are of at least 2(n − 1) order in the Coulomb interaction.
The diagram in Fig. 5 (d) is a generic one describing all the DPE processes including the ones with losses such as shown at Fig. 5 (c). One can go a step further and give a prescription for classes of lossless diagrams. A detailed analysis of this particular situation is possible and will be done elsewhere. Here, we mention without a derivation that such diagrams can be split into the scattering part (the two-particle propagators can be written in terms of the scattering states |ψ (−) k 1 k 2 ,β , cf. Eq. (63)) and the spectral part (containing the two-particle spectral function, Eq. (65)).

B. Example of plasmon assisted DPE
As an example we consider the processes depicted in Fig. 6. The diagrams show a very common situation where a primary electron excited by the laser pulse is loosing its energy on the way to the detector by exciting a secondary electron. There could be either bare or screened Coulomb interaction between the two electrons. In the latter case some resonant phenomena related to the excitation of e. g. plasmon are expected. The SPE case (Figs. 6 (a,b)) is identical to the process of secondary electron excitation considered by Caroli et al. [21]. All DPE processes covered by the diagram at Fig. 6 (c) form a subset of the SPE process. The only difference between the two scenarios is whether primary, secondary or both electrons are observed in the detector. It is obvious that one reduces the DPE diagram to the SPE ones by integration over the energy and momentum of the secondary, or primary electrons, respectively.
Since we do not take into account the interaction between the two emitted electrons (as given, for, instance by two Γblocks at Fig. 5 (a) one can express the final result for the current as a matrix element over the direct product of two singleparticle scattering states. This is typically a good approximation for the case when two electrons have different energies (momenta), or for approximately equal k 1 and k 2 in the case of larger energies [48].
To work this out consider a part of the DPE diagram that contains a product of two GFs involving the external momentum k. Introducing the Fourier representations for each of the GFs , expressing the interacting GF as a product of the Møller operator and the free-particle Green's function (see Appendix B) we obtain expressions similar to Eqs. (60). Thus, in the time domain the product of two interacting singleparticle GFs reduces to a simple propagator computed on the scattering states with incoming boundary conditions: As an exercise let us evaluate the diagram at Fig. 7 (a) describing the SPE process with extrinsic plasmon losses. The current is given by the following expression in the time domain: Representing the lesser Green's function on the vertical track in terms of the electron spectral function (normalized as b µ −∞ dζ 2π A bb (ζ) = N, N is the number of electrons in the system) and the greater component of the screened interaction in terms of the plasmon spectral function representing time-ordered G −− x c (τ , t ) and anti-time-ordered G ++ ax (t, τ) as Fourier integrals and using expression (68) we obtain: Now the limits can be taken making use of an identity discovered by C. O. Almbladh [16] (see Appendix E). It transforms the product of four fractions in the equation above into the product of three δ-functions (2π) 3 δ(ω 1 − ω − ζ)δ(ω 2 − ω − ζ)δ(ξ + k − ω + ζ), and after the frequency integration we obtain The two-particle current is obtained along the same lines using the energy flow as shown on Fig. 7 (b). Similarly to the previous case, the limits η → 0, δ → 0 yield a product (of five) δ-function which were subsequently used to perform three frequency integrations here (see Appendix E). All the quantities in Eqs. (73,74) can be expressed in terms of the spectral functions. We can, for instance, start with a general expresion for the time-ordered function in terms of functions on the Keldysh contour: where in the first equation τ ≡ t − − t − is equal to the timedifference on the forward branch of the contour, and τ ≡ t + −t + is equal to the time-difference on the backward branch of the contour in the second equation. After the Fourier transform f (ω) = ∞ −∞ dτ e iωτf (τ), we havê The fluctuation-dissipation theorem at zero temperature allows to express the lesser and greater propagators in terms of the corresponding spectral functions (Kubo-Martin-Schwinger (KMS) conditions [10]): The screened interaction obeys KMS conditions for bosonic propagators: with the symmetry property for the spectral functionB(−ω) = −B(ω) (follows e. g. from the fact thatŴ R (t, t ) is a real function or, more precisely, a Hermitian matrix). We have already used these equations (cf. Eqs. (70,71)) to express SPE current in terms of spectral functions. Using Eq. (76) we can write the spectral representation of the fermionic propagator where µ is the Fermi energy. The anti-time-ordered GF is obtained similarly G ++ (ω) = −[G −− (ω)] † . The screened interaction is expressed as an integral over the positive frequencies: Let us consider plasmon-mediated DPE. This process is of relevance for metallic and large molecular systems. Since plasmon is a long wavelength or small momentum electronic excitation it is useful to go from the abstract basis to momentum representation and write W −− (k, ω) in a short form as where ω p (k) is the plasmon dispersion, ω p ≡ ω p (0) is the classical plasmon frequency, and v k = 4π k 2 is the matrix element of Coulomb interaction. It is clear that in this form the plasmon peak completely exhausts the f -sum rule. Such plasmon pole approximation for the screened interaction is broadly used in the electronic structure calculation when full-fledged calculations are not feasible. Similarly, it can be used to simplify Eq. (74).

C. Numerical results
Let us make some simplifications. Usually it is a good approximation to start with the mean-field Green's functions where n a is the occupation number of the state a and n a ≡ 1 − n a . After straightforward, but tedious calculation the frequency integrations in Eq. (74) can be performed (for technical reasons it is better to start from the time rather then frequency expression, and it can be obtained by directly transcribing the diagram at Fig. 7 (b) using standard rules) yielding the following expression for the two-particle current: with the following matrix elements Notice that it is not necessary to separately treat the bare Coulomb interaction, it can be recovered as ω p → ∞ limit as explained in [62]. Let us compare Eq. (79) with the general result obtained using the Feshbach projection formalism (64). For the meanfield approximation (78) the two-particle spectral function is diagonal and is given by the convolution of two single-particle spectral densities: The energy conservation for the whole process, which is given by the δ-function in the numerator of (79), is expressed in terms of the two-particle spectral function A (2) (ε k 1 + ε k 2 − ω), (cf. Eq. (65)). The denominator of the first line reflects the resonant character of the considered two-step process. From the resonance conditions (zeroes of the denominator) we see that the double photoemission is enhanced when a and c are continuum states and therefore we denote them as k a and k c . We replace the scattering states |χ (−) k 1 and |χ (−) k 2 entering the matrix elements (80) by the plane-waves and perform the integration yielding f q k a k = δ(k−k a −q). Combining all together we obtain the following concise expression for the plasmonassisted DPE process: We have seen that the plane-wave approximation for the scattering states (i. e. the Møller operator is given by the identity operator) results in a great simplification for the two-particle current: it is given by a sum over two bound states (they correspond to two lesser propagators in the diagrammatic representation of this process) and by the two momentum integrals corresponding to the propagators of the secondary electron. In contrast, in the full-fledged calculations based on Eq. (79) the momenta of the secondary electron and the emitted electrons are not rigidly related. Therefore, in general two additional momentum integrations are required. This will be the subject of a forthcoming publication where this formalism is applied to a large molecular system. The DPE process described by Eq. (82) is suited to probe the plasmon dispersion and damping. First, let us look at the classical plasmon that carries vanishing momentum and otherwise is strongly damped. This leads us to consider the case k a ≈ k c ≈ k 1 , and ε d − ε k 2 = ω p is the condition for the plasmon resonance. In this case the second line reduces to | k 2 |d | 2 /ω 2 p , and is clearly off-resonance. The situation greatly changes if we allow for the plasmon to carry finite momentum q c and consider a large momentum of the secondary electron k a ≈ k c ≈ k 1 > √ ω p . For simplicity take a symmetric situation when both screened interaction lines carry approximately the same energy and momentum and denote K ≈ 1 2 (k a + k 1 ) ≈ 1 2 (k c + k 1 ) and q ≈ k a − k 1 ≈ k c − k 1 . In this case one achieves the resonant enhancement when Thus for colinear k a , k c and k 1 the probability for the plasmon-assisted emission of the secondary electron is enhanced when K reaches the value of ω p /q c .
In order to illustrate the features arising due to the plasmonassisted process in an experiment, we computed the current for a simple model system. To be concrete, we consider the basic jellium model for the C 60 molecule (treated as spherically symmetric) [63,64], which is known for its pronounced (dipolar) plasmon resonance at ω p ∼ 22 eV. Inserting a smoothed box-like potential as approximation to the Kohn-Sham potential, we solved the Schrödinger equation for the 120 orbitals required (240 electrons in total). This procedure yields the single-particle energies ε d associated to the orbitals φ d (r), from which we can compute all quantities in Eq. (82). Because of the spherical symmetry, we can separate the radial and the angular dependence, that is φ d (r) = u d (r) r Y d m d (r) (Y m (r) are the spherical harmonics) and only solve the radial Schrödinger equation. For the optical matrix elements, we choose the length gauge and assume a linear polarization along the z axis (∆ = z). Since we are not interested in the absolute scale a prefactor proportional to the field strength will not be included. The matrix elements ∆ kb attain the form where j denotes the spherical Bessel function. The coefficients C m b m b are obtained from the standard Clebsch-Gordan algebra [65,66]. Similarly, the Fourier-transformed orbitals k|d = φ d (k) can be expressed in terms of the Bessel transformation: Next we transform the summation over k a and k c into integrations and substitute them by the integration over the momentum transfer vectors q a,c = k 1 − k a,c . At this stage, no further simplification can be made, such that the six-dimensional integral has to be evaluated. However, it is reasonable to consider q a,c as small, since the plasmon branch enters the particle-hole continuum for growing momentum, where it is strongly damped. Hence, we introduce the momentum cutoff q max and assume k 1 , k 2 q max . Thus, we approximate ∆ k a,c b = ∆ k 1 −q a,c b ≈ ∆ k 1 b and φ d (k 1 +k 2 −k a,c ) = φ d (k 2 +q a,c ) ≈ φ d (k 2 ). Furthermore, we integrate over the spherical angles of k 1 and k 2 , keeping only the dependence on their magnitude. Thus, the two-electron current can be written as where  (d) Equal energy sharing (ε k 1 = ε k 2 ) for the current and trace of the two-particle spectral density (shaded curve).
Note that we inserted the imaginary shift iΓ in the energy argument accounting for a finite width (lifetime in the time domain) of the plasmon resonance (which is assumed dispersionless for simplicity).
In an experiment, the distinction between primary (k 1 ) and secondary electron (k 2 ) is, of course, not possible. For this reason, the photo-current needs to be symmetrized (let us denote it by J sym ). Representing the J sym as a function of ε k 1 and ε k 2 yields the typical energy-sharing diagrams (Fig. 8). Spectral properties of the system (dominated by A (2) (ε)) display themselves along the main diagonal, as only the sum ε k 1 + ε k 2 enters. Dominant scattering events mediated by the (screened) interaction on the other hand are visible along lines ε k 1 = const. (or ε k 2 = const.). As Eq. (83), indicates the two-particle current contains contributions from (i) the bare Coulomb (two interacting lines in Fig. 7 (b) are not screened), (ii) plasmonic scattering (both lines are screened), and the interference terms. (i), (ii) and the total contribution is shown at panels (a), (b) and (c) of Fig. 8, respectively. For vanishing Γ the current is dominated by sharp plasmonic resonances. For finite damping parameter such as used for the present simulations (Γ = 0.1, we use a realistic value as in Ref. [67]) the interference terms are important: we still have a large plasmonic contribution (viz. Fig. 7 (b)), however, the bare Coulomb contributes with the opposite sign. Therefore, in total current the large peak at ε k 2 ≈ 0.15 becomes less pronounced and ad-ditional peaks at higher energies (e.g. at ε k 2 ≈ 0.5) appear. The whole spectral width of the signal is limited by the twoparticle spectral function shown at Fig. 7 (d) as a shade curve.

VI. CONCLUSIONS
There is a large number of theoretical works devoted to the interaction of light and matter which involves the emission of one or more electrons. This contribution is meant to expose parallels between the single and the double electron photoemission in a formal way. We started by defining corresponding observables and deriving expressions for one-and two-particle currents based on the first-order time-dependent perturbation theory. These expressions are suitable if exact formulas in terms of many-body states are required. In order to obtain computationally useful expressions many-body effects should also be accounted for in a perturbative fashion. Thus, in the first part of the manuscript we applied the projection operator formalism. Starting from the explicit form of the projection operators dividing the whole Hilbert space of the system into that of the emitted electron(s) and the target we derived the effective one-and two-particle Hamiltonian, discussed integral equations for the Green's functions describing emitted particles and demonstrated a close connection of this formalism to the nonequilibrium Green's function theory. For the latter, one can easily derive the diagrammatic expansions for one-and two-particle currents starting from the timedependent perturbation theory and using the adiabatic switching of the electron-electron interaction. Hence, we have electromagnetic field switched on at the remote past (as e ηt ) and independently adiabatically switched on the interaction such that the total Hamiltonian takes a formĤ δ =Ĥ 0 + e −δ|t|Ĥ 1 . We analyzed in details the diagrammatic structure of one-and two-particle currents. It is surprisingly simple: one starts with the density-density response function χ < which necessarily contains two blocks associated with the forward ("−") and backward ("+") parts of the Keldysh contour. Requesting that one or two lines flowing from "−" to "+" blocks are associated with scattering states (with momenta k i ) one obtains exactly the diagrams for SPE and DPE currents showing the close connection between these types of light-matter interaction. It is not difficult to generalize this approach to an arbitrary number of particles. Finally, we presented a detailed analysis of the plasmon-assisted DPE and showed that if one of the emitted particles is unobserved, its diagrammatic representation reduces to the one describing external losses in the SPE process considered by Caroli et al. [21]. Plasmon pole approximation was employed to derive computationally manageable expressions. We illustrated the distinct features to be expected in an experiment by analyzing the simple and yet realistic jellium model for the C 60 molecule. This will extended used in the forthcoming paper devoted to the ab-initio treatment of this large molecular system.

ACKNOWLEDGMENTS
This work is supported by the DFG grants No. SFB 762 and No. PA 1698/1-1.

Appendix A: Particle-impact ionization
Under some circumstances the formalism developed in the main text can be extended to other mechanisms of ionization, e.g., particle-impact ionization. The basic requirement we impose is the distinguishability of the projectile from the target electrons. This applies also for a projectile electron if the impact energy is high and the small momentum transfer is small (optical limit).
The target we describe by the Hamiltonian Eq. (4). The Coulomb-interaction between the projectile (with charge Z) and the sample readŝ is the annihilation (creation) operator of the projectile states |ν . These states can be chosen as the eigenstates of the projectile Hamiltonianĥ p with energy ε ν .
Assuming that the projectile initially possesses the momentum k i , we can construct the asymptotic state prior to the interaction (that is, at t = −∞) as the product state |Ψ 0,k i = |k i ⊗ |Ψ 0 .
V is switched on reaching its full strength at t = 0. Assuming that this strength is still way smaller than the internal interaction within the sample, we can apply the first order perturbation theory (i.e., the first Born approximation in the projectile-target interaction [47]). Denoting the full Hamiltonian byĤ +ĥ p one may write The projectile has a well defined final momentum k f . In analogy to Sec. IV A we introduce the particle number operator N k → P fNk P f with P f = |k f k f | projecting only onto the projectile space. N k acts on the system's states only (including the ejected electrons upon particle impact). Evaluating then the current as in Sec. II A and approximating the projectile states by plane waves r|k = e ik·r yields where q = k i − k f is the momentum transfer, andV eff (q) is the effective single-particle operator acting on the target, explicitlŷ In this optical limitV eff (q) = 4πZ q 2 e iq·r .
acts similar to the light-matter interaction∆; the transferred energy (or energy loss) ε k i − ε k f resembles the photon energy.
and Sigal [71]. It is further possible to write the inverse of the matrix M explicitly [72]: This identity is natural to apply to compute resolvents. For instance, Eq. (26) is given the first line of Eq. (C4). This formula can also be found in Almbladh as Eq. (19) [16].