Nonlinear response in the cumulant expansion for core-level photoemission

Most currently used approximations for the one-particle Green’s function G in the framework of many-body perturbation theory, such as Hedin’s GW approximation or the cumulant GW + C approach, are based on a linear-response approximation for the screened interaction W . The extent to which such a hypothesis is valid and ways to go beyond have been explored only very little. Here we show how to derive a cumulant Green’s function beyond linear response from the equation of motion of the Green’s function in a functional derivative formulation. The results can be written in a compact form, which opens the possibility to calculate the corrections in a ﬁrst-principles framework using time-dependent density functional theory. In order to illustrate the potential importance of the corrections, numerical results are presented for a model system with a core level and two valence orbitals.


I. INTRODUCTION
Core-level x-ray photoemission spectroscopy (XPS) is a sensitive probe of correlation properties in condensed matter [1].At high photon energies, where extrinsic and interference effects due to scattering of the outgoing photoelectron can be neglected, the XPS photocurrent is directly related to the spectral function associated with the core-hole Green's function.The importance of manybody effects for spectra related to core excitations has long been recognized, and the corresponding research has a rich history, going back at least to the 1920's [2][3][4].Of particular interest are edge singularities and asymmetric line shapes.These could be explained by the coupling between the core hole and the valence electrons using studies based on model Hamiltonians, e.g., in the seminal works of Mahan, Nozières and De Dominicis, and coworkers [1,[5][6][7][8][9].Several exact and approximate results were obtained, for example assuming the interaction potential is separable [5].Another commonly used assumption is the absence of interaction between the valence electrons themselves, which implies that their scattering from the core hole potential leads to excitation of independent electron-hole pairs.In the same framework, the core-hole problem was treated by solving two coupled Bethe-Salpeter equations for core and valence electrons to lowest order in a parquet approximation [6], to which self-consistency in self-energy and vertex were added in [7].In the subsequent paper of the same series [8], an effective one-body approach was proposed, based on the calculation of the transient response to the sudden creation of a core hole.This was complemented by Langreth [10], with a more compact derivation, and in [11], by an alternative approach based on a finite number of electrons in a box, which validated the approximation * of a separable potential and added numerical illustrations.
The picture of excitations created by the sudden appearance of the core hole is also naturally reflected in fermion-boson coupling Hamiltonians, where the fermion refers to a deep core orbital, and the bosons are electronhole excitations, plasmons [12,13], or phonons [14][15][16][17].Remarkably, for a single fermion with linear coupling to a dispersing boson, the model can be solved exactly [18][19][20].The fermionic spectral function then consists of a quasi-particle peak followed by a Poisson series of satellites.For increasing coupling strength the quasi-particle loses increasing weight to the satellites, and the envelope of the satellite spectrum becomes a Gaussian in shape.For this model, this exact solution is equivalent to the spectral function of the second-order (in the coupling constant) cumulant Green's function [1,15].For more than one fermion energy and/or for higher-order coupling to bosons, the cumulant solution is not exact [15]; nevertheless since it contains the essential physics of the system responding by bosonic excitations to the the creation of an additional electron or hole, it is often a good approximation that is widely used, in particular for the core hole problem [21][22][23][24][25][26][27].
Of course the electron-boson picture is just another way to look at the old problem of many-body effects in core spectroscopy, but it highlights three questions that one might ask, namely: i) How well does a Hamiltonian with linear electron-boson coupling describe the real problem?; ii) Which excitations should be contained in the boson?; and iii) Is the second-order cumulant solution good enough for the real application, or if not, how can one go beyond?Heretofore these questions have not found a definite answer, not least because the domain of application of the cumulant approach in condensed matter is very wide, including for example coupling to phonons [16,28], or valence electron spectroscopy [29][30][31][32][33][34][35][36][37][38][39][40].
Analogous questions arise in an at first sight, different framework which is the calculation of Green's func-tions from a Dyson equation, where interaction effects are contained in the self-energy, the integral kernel of the equation [41].This is the framework now commonly adopted in first principles calculations.Currently the most widely used approximation is Hedin's GW approximation (GWA) [42], where the self-energy is approximated as the product of the one-electron Green's function G and the screened Coulomb interaction W to first order.Indeed, one can view the GW approximation as an approximate solution of the electron-boson problem [1], where the bosonic excitations are the excitations contained in W . Also in this case, since the coupling is only linear, the boson contains only a selection of excitations (spin flip excitations, for example, are neglected), and the problem is solved only approximately.While the GWA has been very successful for the calculation of quasi-particle energies, it typically gives a poor description of satellite spectra [43].This shortcoming can be traced to the third of the above questions -the approximation used in the solution of the Hamiltonian problem -since this approximation leads to problems such as the appearance of a spurious plasmaron solution.This can be seen in the GW solution of the exactly solvable model [1].Changing from the GWA to a cumulant Green's function leads to a significant improvement without much additional effort since the same boson appears through W but the electron-boson problem is solved more accurately.The GW plus cumulant (GW+C) approach for the oneelectron Green's function [29,32] has recently become a popular first principles method to describe spectral functions, including satellite series.Here, the satellites are mostly due to plasmons [36], the dominant excitations in W .
The success of the GW and GW+C approximations may seem surprising in view of the relative simplicity of these approximations being linear in G and W .While it is physically reasonable, e.g., in many simple semiconductors, that plasmonic excitations should represent the dominant boson, it is less obvious that the physics is sufficiently described by linear response.This is especially evident in systems with few electrons, or when the removal or addition orbital is very localized, where one might expect that response terms of higher order should play a role.These terms are not contained in the electron-boson model with linear coupling, because the boson itself is fixed and does not respond to the excitation; analogously, they are not contained in GW nor GW+C because W is calculated on the ground state, in presence of the original N electrons.The fact that non-linear screening should be present and should show up in a cumulant solution has been pointed out early on by Mahan [5], who calculated the leading correction to the standard second-order (in the coupling constant) linear-response cumulant solution.In that work the valence electrons that respond to a core hole excitation are described by an independent-electron picture, and even so, going to yet higher orders turned out to be too complicated.While this pioneering work, as well as the seminal solution of Nozières and De Dominicis [8], trace a way to go towards the inclusion of non-linear screening effects, these model approaches are not directly transferable to first principles calculations, for several reasons: i) because of the approximations involved from the very beginning on the interaction potential; ii) because of the absence of interaction between valence electrons, which would lead to a poor description of screening with the absence of plasmons; and iii) because it is not clear how an expansion that is orderby-order concerning the response functions (i.e., linear response, second order response, etc.) would converge.Progress in several directions has been made, such as a better description of screening [44,45], higher-order cumulant solutions of the electron-boson model, such as in [15,46,47] the derivation of higher order correlation functions from the Mahan-Nozières-DeDominicis (MND) or electron-boson Hamiltonians [48], or more recent progress with the description of non-linear electron-boson couplings in Ref. [49] and in the equation-of-motion, coupledcluster method [50].However, there is still a gap to be filled concerning the ab initio calculation of spectral functions.
The present work aims at developing a robust, firstprinciples derivation of non-linear screening effects in a cumulant Green's function.While we do not claim to invent new physics here, our derivation has the advantage of being compact, and situated within the framework now commonly used in the ab initio community.Our approach does not depend on a contact or separability approximation on the interaction potential.Moreover, the approach highlights the underlying physics, which is crucial if one wishes to understand when non-linear effects are important, and where should be the limits of their applicability.Finally, it allows us to propose a way to put the equations into practice, by combining the many-body perturbation formalism with time-dependent density functional theory (TDDFT) [51].In this way, an order-by-order expansion of the response is avoided, and the problem of calculating the effects of screening on the Green's function is separated from that of the calculation of the screening itself.Contrary to numerous previous works on core spectroscopy, we are not so much interested in asymmetry of line shapes but rather, in the satellites, for which the standard second-order linear response cumulant approach is now a well established first principles approach and convenient starting point.
The paper is organized as follows: The background concerning the formalism and the GW approximation is contained in Sec.II.Next, we derive the core-hole cumulant Green's function in its various approximations in Sec.III, and we analyze the results in Sec.IV.Numerical results for an illustrative model are given in Sec.V, and a short conclusion is given in Sec.VI.

II. BACKGROUND
Beginning with its equation of motion following the approach of Martin and Schwinger [52], the one-body Green's function can be described by a functional differential equation which is often referred to as the Kadanoff-Baym equation [53].The equation was initially derived for temperature dependent, non-equilibrium quantum systems.However, it can also be used to create the diagrams describing equilibrium and/or zero-temperature physics [43], as we will do here.In this formalism, the fully interacting propagator G is given by 13)G( 33 + )G( 12) 13) δG( 12) where G 0 is the non-interacting propagator in the presence of u, v(12 the Coulomb interaction, and u(3) is a local, time-dependent external perturbing potential that simulates interaction effects due to the propagation of particles; it will be taken to zero at the end of the calculation.Here and below we employ the usual notation of an integer for a set of space, spin and time variables (1 → (r 1 , σ 1 , t 1 )), and bars for variables that are integrated over: f ( 1)g( 1) ≡ d1f (1)g (1).All quantities are functionals of u.The classical Hartree term (the second term on the right) depends on the interacting density as given by the diagonal part of the propagator n(1) = −iG(11 + ), where 1 The last term contains the functional derivative of G, which accounts for the exchange and correlation effects.This term turns Eq. ( 1) into a first-order nonlinear functional differential equation with respect to u.
Introducing the total classical potential u H , one obtains a set of coupled equations for u H and G, In extended systems, screening plays an important role.For this reason, is is often convenient to rewrite the functional derivative using the chain rule with respect to the classical potential, where Note that here the screened potential W is a functional of the external potential, such that the equation remains exact.
Often, the dependence of W on u is neglected.The solution of this linear-response version of ( 5) in a simple model has been discussed in [54].In general, however, even in the linear response approximation, the equation cannot be solved exactly.Therefore, the functional derivative on the right side is usually approximated such that the limit u → 0 can be taken directly.One of the most widely used approximations for the Green's function is Hedin's GW approach [42], which is obtained by approximating the functional derivative, δG(12) δu H (4) ≈ G( 14)G( 42).
Then the u → 0 limit can be taken, and one has a Dyson equation for G, with The GW approximation (GWA) has been very successful for the calculation of quasi-particle energies.However, a major shortcoming is its poor treatment of the satellite part of electron addition or removal spectra [29,32,36,43].Calculations based on a cumulant Green's function which yields a much much better satellite spectrum have been found to be an advantageous alternative.The cumulant approach, which avoids the approximation in Eq.( 7) that leads to the GW self-energy, will be discussed in the following section.

III. DERIVATION OF THE CORE-HOLE CUMULANT
In order to describe the photoemssion spectra from core states, we need the projection of the Green's function on the core orbital G cc .Its associated spectral function simulates the core photoemission spectrum, especially at high photon energies, where the photoelectron and associated effects of extrinsic losses can be ignored.Thus here we make the approximation that the core hole is decoupled from all other orbitals, except for the screening of the interaction, which is due to the valence electrons.The decoupling approximation is physically reasonable for the case of localized electrons, such as a deep-core state, which has little overlap with the valence electrons.For this case, we may suppose that the interacting and the Hartree Green's functions G and G H do not have matrix elements linking the core and a valence state, so Eq. ( 4) written in a basis of single-particle orbitals reads where the matrix elements of the Coulomb interaction are defined as In principle, this equation is understood to be on a contour, since in presence of the time-dependent external potential the system is out of equilibrium [55].As an additional approximation, we only retain the contributions that correspond to the propagation of a hole, i.e. t 2 > t3 > t 1 .This allows us to make the ansatz that the interacting propagator is proportional to G H , i.e.
With the identity δG δU = δG H δU F + G H δF δU , we obtain As in the derivation of the GWA, we can now use the chain rule with the total classical potential, and we again suppose that G H has no off-diagonal elements linking core and valence single-particle states.This leads to with the screened interaction W c ≡ W cccc , where .
At this stage W c (t 1 t 4 ; u) still depends on the potential u, which has not yet been set to zero.Also note that ( 14) suggests a response of the density; however, one has to keep in mind that here we are not in a retarded framework, but only keep parts corresponding to the time ordering defined above.With the constraint t 2 > t 1 imposed on the time orderings, G H cc = ie −iε H c (t1−t2) factorizes, and the factors cancel.This yields .
Next, we define F in cumulant form and take the derivative of Eq. ( 15) with respect to t 1 .This yields a differential equation for the cumulant function C With the boundary condition C(t 2 t 2 ) = 0, the integral equation for C is given by The last term of this expression is commonly neglected, so one can directly set the external potential to zero.This yields the widely used cumulant in the linear response approximation For the core hole, together with Eq. ( 14) the expression can be interpreted as the integral over the variation of the Coulomb potential due the valence density at time τ , induced in linear response of the system by the creation of a hole at time τ .The response is causal with τ > τ , and the process is integrated over the interval (t 1 , t 2 ).Implicit in the evaluation of the cumulant C 0 is that only the positive frequency components of the Fourier transform W c (ω), corresponding to lossy excitations, are present in the linear response approximation [18].
To go beyond linear response, one can iterate (18).The lowest order corrections stems from the derivative of C 0 , taking into account that δW c /δu = 0.This causes a second order response function to appear, in a contribution C 1 given by This term contains the variation of the response of the valence density due to the density change caused by the linear response to the creation of the core hole.Higher order corrections involve higher order derivatives of W c , and therefore higher order non-linear response functions.They should be evaluated from the density-density correlation functions in the ground state, i.e. before the removal of the core-electron.The general form of the solution is thus given by where the m th contribution C m is obtained recursively using the relation with C 0 given by Eq. ( 19) in terms of W c .The recursive formulation suggests that the series might in practice be truncated at a given order n.Once all orders of the correction have been found, the limit of u → 0 is applied.
The full recursive solution implied by Eq. ( 22) is the first main result of this paper.It may be seen as an extension of Mahan's approach in [5], which contains the first nonlinear response contribution.The main difference is a more general formulation such that i) a separable potential is not needed needed for the derivation: ii) the result is valid to infinite order in the response; and most importantly, iii) the valence electrons are in principle fully interacting.As we will see below, the fact that the result is formulated in general terms of response functions rather than explicit sums over transitions, allows us to benefit from the use of TDDFT for an efficient approximation to the non-linear response.
IV. ANALYSIS

A. Effective interaction
The cumulant C 0 (tt ) in Eq. ( 19) and consequently the higher order terms, are double integrals of a two-time function over time.We can therefore express the cumulant in terms of a new function w(τ τ ) ≡ m=0 w m (τ τ ), where w 0 (τ τ ) = W c (τ τ ) and w m stands for the order m correction in W c .Each cumulant term is related to w through the expression A recursive relation, similar to Eq. ( 22), holds between the different orders of w, where τ > τ , and τ the time that refers to the variations of the density.The interaction −v cckl couples the core level from one side, to the variations of the external potential taken with valence levels on the other side.The negative sign can be understood in terms of the sign of a core-hole charge and w m can also be viewed in terms of the orders of an expansion to this core-hole potential.Therefore the full matrix w(τ τ ) = n m w m (τ τ ) plays the role of an effective interaction, accounting for all orders of the density variations due to the propagation of a core-hole.

B. Induced density variations
In order to enable further interpretation and practical use, it is convenient to expand Eq. ( 24) for m > 0, With w 0 (τ τ ) = W c (τ τ ) and ( 14) the cumulant is given in terms of the induced valence density variations as where the factor 1/(m + 1)! is due to the extension of the integration domain, and repeated indices are summed over.
The first term yields the exchange correction to the energy of the core level.The remainder can be summed: it is an expansion of the Coulomb potential acting on the core hole, created by the change in density ∆n at time τ , which is in turn due to the potential −v cc (r)θ(t − t 1 created by the switching on of a core hole at time t 1 .This potential is then integrated in τ over the time of propagation of the core hole.The compact result can be written as However, as pointed out earlier, ∆ñ is not equal to the causal response of the density.In particular, in order to extend the integration domains and obtain this compact result, the symmetry of the time-ordered W has been used.In order to use the result, we therefore have to make an approximation.One simple possibility is to use the real time-dependent induced density in Eq. ( 27).However, this approach produces negative spectral weight on the high energy side of the quasiparticle (see Fig. 1).As an alternative approximation here we restrict ∆ñ ij (τ ) to include only the positive frequency components of the induced density, as the case in linear response.This result has a physical interpretation which moreover has a significant practical advantage: since the core orbital is kept fixed the perturbation is known and its effect can be calculated directly.For example, the TDDFT in real space and time can be used, as proposed in [56] for the case of linear response.It is then not necessary to calculate the rather clumsy higher order response functions, which thereby overcomes the issues discussed in [5].
When the core hole is suddenly switched on, the system may react violently.However, an interacting system will eventually reach a new equilibrium, given by the final state of the system with a static core hole.In this limit ∆ñ ij (τ ) becomes independent of τ .The simplest approximation is therefore a shift of the core level due to the fact that the valence density should be calculated in presence of the core hole.Since this approximation completely neglects dynamical effects, it cannot lead to satellites in the spectral function.In order to do better, one could apply the same reasoning starting from the next order, which means in practice, to evaluate explicitly the lowest order cumulant with the valence density calculated in presence of the core hole.Such an intuitive approach, which has been used successfully (see e.g., [57][58][59]), may be justified by our derivation.

C. Self-energy
Let us now compare the above result with what one would obtain in the context of a Dyson equation.For this case, we start again from the exact expression in Eq. ((5)) and then using δG/δu H = −G(δG −1 /δu H )G, the exact Dyson equation becomes which has a self-energy given by Again taking the matrix element in the core orbital, and using the relation This self-energy has the same structure as the GWA, but instead of a linear response W , one has an effective interaction If the contribution in the last line is neglected, i.e., if one neglects the variation of the core hole Green's function, w eff = w is the effective interaction (24) that appears in the cumulant.To lowest order neither the neglected term nor the functional derivative of w eff contribute.This is also the reason why the GW self-energy and the lowest order cumulant have the same effective interaction, namely the linear response screened Coulomb interaction W .For the higher orders, however, it is important that the effective interactions are different, since w has to be used in the cumulant and w eff in the Dyson equation, while both are exact in principle within the approximations made here.For example, again making use of the product property of the core hole Green's function, the last term gives rise to a contribution , which is neglected in the GWA, but already included in the cumulant Green's function through the lowest order cumulant function.
Finally, the form Σ xc = iGw eff is reminiscent the Tmatrix formalism [43,60].Such T-matrix approximations can be derived by supposing that the self-energy consists of a Green's function and an effective interaction which, inserted into Eq.( 29), allows one to derive a Dyson equation for the effective interaction.However, for our present core-hole problem we have assumed that the core Green's function does not vary whereas derivatives of the interaction are taken into account.Instead, to obtain the scattering diagrams in the T-matrix approximations the Green's function is varied and the effective interaction is kept fixed [43,60].This is more important in a situation of low density (partial filling), which is quite different from the core hole problem.

V. MODEL CALCULATIONS
As a concrete illustration of the approach presented here, we introduce a simple 3-state model, similar to that used by Lee, Gunnarsson and Hedin in [61], which has also been used to treat charge-transfer satellites in x-ray spectra [62].This model system with two electrons, a core electron and a valence electron propagating in two atomic levels a, b, is described by the Hamiltonian where 0 are atomic energies evaluated in the presence of the core-electron, and U is the potential from the corehole nh = 1 − nc coupling only to one of the two levels, a.The valence levels belong to different atoms and this justifies the approximation to consider weak coupling to the level b.The hybridization between the levels a and b is represented by the interaction parameter t.This model aids a physical understanding of core-photoemission in molecules with flat valence bands, where charge transfer excitations between different atoms modify the spectrum.The initial state, where n h = 0, is described by the two-particle state |ψ i 0 = sin φ|a |c + cos φ|b |c that mixes core and valence levels, with tan 2φ = 2t/ and = 0 a − 0 b .The energy of the initial state is given by The final states where n h = 1 have only one particle, and are given by the single-particle wavefunctions |Ψ f 1 = cos θ|a − sin θ|b and |Ψ f 2 = sin θ|a + cos θ|b , with tan 2θ = 2t/( + U ) and In the model the core-hole potential couples only to the level a and therefore the interaction with the time-dependent occupation of the level a will appear.
In order to obtain the cumulant solution for this model, we apply Eq. ( 26) where n 0 a is the occupation of the state a in the groundstate, and the first term, linear in U replaces the bare Coulomb interaction seen in Eq. ( 27) for the model Hamiltonian, and causes an overall quasiparticle shift.Note that there are also further quasiparticle shifts that come from the higher order terms in Eq. (33).] Here, the density is given in terms of the timedependent wavefunction, n a (t) = | a|ψ(t) | 2 , where |ψ(t) is the state of the system at time t after the appearance of the core-hole, and is initially equal to the ground state valence wavefunction, |ψ(0) = |ψ i 0 .This is the only component that contributes to the non-linear cumulant solution due to the fact that the core-level couples only to the atomic level a.As pointed out above, the density variations are linear and non-linear response functions, evaluated in the initial ground state, i.e. without a core hole.Finally we apply Eq. ( 27) for the cumulant, which yields where the time-dependent potential due to the core hole is v a (t) = 0 for t < t 1 and v a (t) = U for t > t 1 .We can now compare the results obtained from the two lowest orders of Eq. ( 33) and the full Eq.(34).27) (blue), the same but with only positive frequencies included in the density ∆ñ (red), the linear response approximation (green), and from including the first non-linear (U 3 ) correction term (purple).Note that the results from the time-dependent density are in good agreement for the positions of the quasi-particle and satellite peaks for U = 1 and U = 7.Even at U = 20 the excitation energy (difference between quasiparticle and satellite) matches that of the exact result, although the result from the real density produces small negative satellites on the high energy side of the quasiparticle.In contrast, the linear response and U 3 corrected approximations vastly under-estimate the satellite energy and give unphysical results, such as large boson-like satellite progression and negative spectral weight.
Results for the core-level spectral function A c (ω) = − 1 π Im G(ω) vs. ω/∆ , where ∆ = f 1 − f 2 is the valence excitation energy in the presence of the core-hole, are shown in Fig. (1).The position of the peaks reflects the energy of the transition of the valence electron from the initial (with the core electron) to the final (with the corehole) state, while the height of the peaks corresponds to the probability amplitude of the transition.Only a nonnegligible value of the parameter t allows for the corepotential to affect the transition energies between the initial and the final state, since otherwise no screening can happen.From top to bottom, the curves show results for core-hole strength U = 1 (top), U = 7 (middle), and U = 20 (bottom), corresponding to weak, intermediate, and strong coupling, all with the parameters t = 3, and = 1.The top set of curves represent the weak coupling limit, which can be seen by the lack of any visible satellite, and by the agreement of all curves, which should be nearly identical in the linear response regime.At intermediate coupling (U = 7, middle), the various approximations now give appreciably different results.The linear response approximation (green) underestimates the splitting between the satellite and quasiparticle positions by nearly a factor of two, underestimates the size of the quasiparticle peak, and produces a second satellite indicating a boson-like progression, as expected.The lowest order (U 3 ) corrected result (purple) gives a small correction to the quasiparticle weight and position, but does little to correct the splitting between the quasiparticle and satellite.The issues with the linear response and U 3 corrected approximations become even more apparent in the strong coupling regime (U = 20), where the linear response vastly underestimates quasiparticle weight and quasiparticle-satellite splitting, and produces a long progression of satellites, initially with increasing weight and overestimates the quasiparticle peak by a large amount.This is unphysical, since a single valence electron cannot produce multiple excitations.Also the inclusion of the U 3 correction produces spurious negative spectral weight.In contrast, the results of Eq. ( 27) with the real-time density used for ∆ñ(t) produces good agreement for the quasiparticle and satellite peak positions for U = 7, and even reproduces the satellite splitting at U = 20, although the satellite weight is underestimated in both cases, and the approximation produces negative spectral weight above the quasiparticle peak.To correct this unphysical be-havior, we have used a form for ∆ñ(t) with the negative frequencies filtered out.This produces a Landau type form for the cumulant, similar to that of the linear response.Results for this approximation (red) show little difference in comparison to those obtained with the real density, apart from the lack of any negative spectral weight, showing that this approximation is a reasonable method for obtaining a physical result.

VI. CONCLUSIONS
We have demonstrated that the Kadanoff-Baym functional differential equation is a convenient starting point to derive the form of the cumulant Green's function beyond the linear response approximation.For a single level that can be considered as decoupled from the rest of the system, such as a localized deep core level, the result can be formulated in a compact way, which highlights the essential physics: i.e., the sudden switching-on of a core hole perturbs the valence density, and the subsequent time-integral of the change in density leads to a quasi-particle correction and to satellites in the spectral function.The approach is tested on a simple 3-level model system similar to that of Lee, Gunnarsson and Hedin.The numerical results suggest that for molecular systems with strong core-hole effects non-linear effects can significantly improve the calculated photoemission spectra and generally need to be taken into account.We suggest that coupling this non-linear cumulant approach with real-time TDDFT can be a promising way to include the non-linear effects in ab initio calculations of the photoemission spectra for such systems.

FIG. 1 .
FIG. 1. Normalized spectral function A(ω) = (−1/π)Im Gc(ω) vs frequency ω for model parameters = 1eV and t = 3eV .The coupling U = 1 simulates the weak coupling to the core-hole (top), while U = 7 simulates intermediate coupling regime (middle), and U = 20 strong coupling.Results are shown for the exact spectral function (black), the real time-dependent density following Eq.(27) (blue), the same but with only positive frequencies included in the density ∆ñ (red), the linear response approximation (green), and from including the first non-linear (U 3 ) correction term (purple).Note that the results from the time-dependent density are in good agreement for the positions of the quasi-particle and satellite peaks for U = 1 and U = 7.Even at U = 20 the excitation energy (difference between quasiparticle and satellite) matches that of the exact result, although the result from the real density produces small negative satellites on the high energy side of the quasiparticle.In contrast, the linear response and U 3 corrected approximations vastly under-estimate the satellite energy and give unphysical results, such as large boson-like satellite progression and negative spectral weight.