Ultrafast dynamics and decoherence of quasiparticles in surface bands: Preasymptotic decay and dephasing of quasiparticle states

P. Lazić,1 V. M. Silkin,2,3 E. V. Chulkov,2,4 P. M. Echenique,2,4 and B. Gumhalter3 1Rudjer Bošković Institute, HR-10000 Zagreb, Croatia 2Donostia International Physics Center, E-20018 San Sebastian, Spain 3Institute of Physics, HR-10000 Zagreb, Croatia 4Departamento de Física de Materiales and Centro Mixto CSIC-UPV/EHU, UPV/EHU, Apdo. 1072, 20080 San Sebastián/Donostia, Basque Country, Spain Received 21 February 2007; revised manuscript received 28 May 2007; published 20 July 2007


I. INTRODUCTION
Recent applications of high resolution electron spectroscopies to the investigations of surfaces have enabled important new insights into the electronic structure and dynamics of probed systems. [1][2][3][4][5][6][7][8] Studies of the electronic properties of confined systems are usually based on electron excitations from or injections into the system, either in the one-step processes as in the direct photoemission ͑PE͒ or inverse photoemission ͑IPE͒ spectroscopy, or in the two step laser pumpprobe induced transitions as in the two-photon photoemission ͑2PPE͒ spectroscopy. Common to all these experiments is a nonadiabatic creation of quasiparticles, electrons in the unoccupied and holes in the occupied states, whose subsequent motion is subjected to the final state ͑in PE and IPE͒ or intermediate state interactions ͑in 2PPE͒ with the remainder of the system. This gives rise to the processes of decoherence and decay of primary excitation 9,10 which affect the spectra of excited quasiparticle͑s͒. Experiments carried out on the systems whose typical relaxation times are much shorter than the duration of measurement provide information on the asymptotic steady state relaxation of excited systems in their passage towards thermodynamic equilibrium. This type of relaxation is commonly described in terms of the rate constants that control the decay and dephasing of quasiparticle states. 7,8 The sources of decay and dephasing of quasiparticles in surface bands are their interactions with the degrees of freedom that constitute the heatbath of the system 10 ͑excitations of electronic charge density, spin density, vibrations, etc.͒ as well as with the localized scatterers 11 ͑impurities, defects, steps, etc.͒. The manifestations of these processes in the measured spectra enable the assessments of dynamical properties of the system. Of special interest in this context are the systems which support unoccupied and occupied quasi-twodimensional ͑Q2D͒ surface bands because that makes possible the studies of dynamics of surface localized electrons and holes in interaction with the same heatbath by employing complementary experimental techniques. A paradigmatic example of this type of structure is the Cu͑111͒ surface with a partly occupied surface state band and an unoccupied image potential state band that were among the first Q2D electron bands extensively studied by high resolution PE 12 and time resolved 2PPE spectroscopy, 3,13 respectively.
New developments and applications of time resolved spectroscopies have pushed the limits of detection of ultrafast phenomena towards the subfemtosecond scale. In this regime the act of measurement may proceed on the time scale comparable to or shorter than that of relaxation processes. Recent spectroscopic measurements utilizing ultrashort laser pulses 14 and novel applications of x-ray techniques 15 have demonstrated the possibility of probing the early evolution of excited quasiparticles with high temporal resolution. However, since the early quasiparticle evolution may considerably differ from the asymptotic Markovian behavior described by the rate constants, 9,16 reliable interpretations of the results of this kind of measurements should be based on nonasymptotic treatments of the dynamics of excited systems. Nonasymptotic treatments enable pinpointing the intervals in which the descriptions of quasiparticle propagation in terms of few rate constants cease to be valid and thereby provide greatly improved physical insights into the evolution of intermediate states in 2PPE and the final states in PE and IPE from surface electronic bands.
The modeling and analyses of electron dynamics in time resolved spectroscopies of surfaces have been carried out dominantly in the framework of optical Bloch equations parametrized in terms of rate constants that were fitted to reproduce the measured data. 17,18 However, the phenomenologically introduced rate constants do not carry information on the early quasiparticle evolution and hence the validity of this approach is restricted to the asymptotic regime in which the time scale͑s͒ of perturbation of the system by external probe͑s͒ is separated from the time scale of relaxation processes that affect spectral shapes. 6,18 To overcome this shortcoming of the phenomenological approaches we develop in the present paper a microscopic description of presymptotic evolution of electrons and holes subsequent to their promotion into the surface localized states. The formalism is applied to illustrate and interpret ultrafast quasiparticle dynamics in Q2D bands on the benchmark Cu͑111͒ surface in the course of spectroscopic measurements. Extension to other systems that exhibit surface electronic bands is straightforward as it only requires modifications in the input of unperturbed band structure.
The present work is focused on the description of decay and dephasing of quasiparticles caused by the interaction with the heatbath of electronic excitations. We model this interaction within the linear response formalism developed earlier to study dynamical electronic properties of metallic surfaces. 7,8 In Sec. II we first outline a general method for calculating the self-consistent ͑SC͒ response of interacting electrons confined in a thick slab which enables the identification of the bulk and surface localized electronic states in the system and of the ensuing spectrum of single pair and collective excitations. The method is applied to calculate the electronic response of a Cu͑111͒ surface and to demonstrate its equivalent boson representation needed for tractable descriptions of quasiparticle interactions with the heatbath. In Sec. III we define the survival amplitudes and probabilities which contain the desired information on ultrafast dynamics of quasiparticles promoted into surface bands. Upon establishing the relation between electron and hole propagators in the real time domain and quasiparticle survival amplitudes the latter are calculated by combining the description of preasymptotic quasiparticle dynamics elaborated in Ref. 19 ͑hereafter to be referred to as I͒ and the formulation of bosonized electronic response developed in Sec. II. Analyses of the thus calculated survival probabilities of hot quasiparticles excited in surface bands on Cu͑111͒ enable the identification and assessment of electronic relaxation processes characteristic of the final states in PE and IPE and of the intermediate states in 2PPE spectroscopy of surfaces. In concluding Sec. IV we summarize the main results of application of the developed formalism to description of quasiparticle dynamics in surface bands and discuss its merits in the context of complementing and refining the earlier theoretical treatments of ultrafast dynamics in which the decoherence and decay of quasiparticles were introduced either on a semiphenomenological basis 13,17,18 or were restricted to asymptotic regimes. 7,8 To facilitate applications of the developed description in interpretations of time resolved measurements we formulate a simple interpolation scheme for fitting the non-Markovian ultrafast decay of quasiparticles that could be used for analyzing the results of high resolution spectroscopic studies.
In Appendix A we present a detailed derivation of the electronic response function for a thick atomic slab required in the calculations of survival probabilities in Sec. III. In Appendix B we use this formulation to derive a quantum analog of the classical image potential needed in the discussion of final state screening effects in Sec. II B. Preliminary results of this work have already been published in Ref. 16.

A. Electronic response of a metallic slab
Dynamical interaction of a probe charge with the charge density fluctuations in a polarizable medium can be conveniently described within the formalism of dielectric response. Dielectric properties of planar metallic surfaces have been extensively studied in the past four decades, mainly through the descriptions of electronic response by the linear densitydensity response function ͑r , rЈ , t͒ for interacting electron gas homogeneous in the direction parallel to the surface. 20,21 In our calculation of the properties of response functions pertinent to metallic surfaces that support Q2D electronic bands we start from the definition ͑r,rЈ,t͒ = − i͗0͉͓ ͑r,t͒ † ͑r,t͒, ͑rЈ,0͒ † ͑rЈ,0͔͉͒0͑͘t͒, ͑1͒ where ͓..., ...͔ denotes a commutator, † ͑r , t͒ and ͑r , t͒ are the electron field creation and annihilation operators in the Heisenberg picture, respectively, with r = ͑ , z͒ and rЈ = ͑Ј , zЈ͒ denoting the coordinates of the field with ͑Ј͒ parallel and z ͑zЈ͒ perpendicular to the surface plane, and ͉0͘ represents the ground state of the system. In the calculations of response functions corresponding to electron systems with translational invariance along the surface we shall represent the various propagators and interactions by their twodimensional ͑2D͒ spatial and temporal Fourier transforms ͑FT͒ following the notation of Ref. 22. Summations over 2D wave vectors parallel to the surface are performed according to ͚ Q → ͑L /2͒ 2 ͐ d 2 Q and inverse Fourier transforms are obtained as f͑z , zЈ ; ͒ = ͑1/L 2 ͚͒ Q exp͑iQ͒f͑z , zЈ ; Q͒ = ͐ exp͑iQ͒f͑z , zЈ ; Q͒d 2 Q / ͑2͒ 2 , where L plays the role of box quantization length in x and y directions. With these conventions the 2D spatial and temporal FT of the response function ͑1͒ reads ͑hereafter ប =1͒: ͑z,zЈ;Q,͒ = ͵ d 2 e −iQ ͵ dte it ͑z,zЈ; ,t͒, ͑2͒ where = ͑ − Ј͒ and Q is a 2D wave vector parallel to the surface. The dimension of the thus defined ͑z , zЈ ; Q , ͒ is ͑length͒ −4 ϫ ͑energy͒ −1 .
In this section we construct ͑z , zЈ ; Q , ͒ corresponding to interacting electron gas at zero temperature that is confined within a thick atomic slab translationally invariant in the surface plane because this geometry has proved convenient in the assignments and numerical treatments of elec-tronic states localized at atomically flat surfaces. [23][24][25] We start from 0 ͑z , zЈ ; Q , ͒ corresponding to a noninteracting electron gas in the same geometry to obtain self-consistent ͑z , zЈ ; Q , ͒ in the random phase approximation ͑RPA͒ and then make use of the spectral properties of response functions to calculate the excitation spectrum of interacting electrons in the slab.
In the framework of SC RPA the response function ͑z , zЈ ; Q , ͒ is obtained by solving the integral equation ͑z,zЈ;Q,͒ = 0 ͑z,zЈ;Q,͒ The ingredients of this equation are 0 ͑z , zЈ ; Q , ͒ calculated below and the 2D Fourier transform V͑z , zЈ ; Q͒ of the bare Coulomb interaction ͓see Eqs. ͑A4͒ and ͑A5͔͒. The noninteracting electrons are described by the single particle wave functions ͉͗r͉K,n͘ = K,n ͑r͒ = 1 ͱ L 2 e iK n ͑z͒,

͑4͒
where K and n are the 2D wave vector and the quantum number describing the particle motion parallel and perpendicular to the surface, respectively. The single particle orbitals n ͑z͒ are normalized to the length perpendicular to the surface so that ͗KЈ , nЈ ͉ K , n͘ = ␦ K,K Ј ␦ n,n Ј , and the corresponding energies ⑀ n are solutions of the 1D Schrödinger equation Here V MP ͑z͒ is the effective model potential described in Ref. 26 that is constructed so as to follow the asymptotic form of the semiclassical image potential for an electron placed at distance z outside the metal surface with e denoting the electron charge and z 0 the effective position of the image plane, respectively ͓note this difference in selection of V MP ͑z͒ with respect to Refs. 20 and 21͔. For slabs of finite thickness analogous asymptotic behavior of the potential is also introduced outside the other surface ͑see below͒. The potential V MP ͑z͒ reproduces the key properties of the surface electronic band structure which in the case of ͑111͒ surfaces of noble metals are the presence of a band gap at the center of the 2D Brillouin zone and the existence of surface localized Shockley and image potential states as the initial unperturbed state features ͑hereafter to be referred to as SS and IS states, respectively͒. Upon expanding the particle field operators in Eq. ͑1͒ in terms of unperturbed wave functions ͑4͒ and electron creation and annihilation operators c K,n † and c K,n for the corresponding eigenstates, respectively, viz. ͑r͒ = ͑,z͒ = ͚ K,n e iK n ͑z͒ ͱ L 2 c K,n ,

͑7͒
we find that the 2D FT of the response function for noninteracting electrons takes the form 0 ͑z,zЈ;Q,͒ = 2 ͚ n,nЈ n ͑z͒ n Ј ͑z͒ n ͑zЈ͒ n Ј ͑zЈ͒ where the prefactor 2 stands for spin, f K,n is the Fermi distribution, ␦ is a positive infinitesimal, and denotes the energy of a particle in the eigenstate ͉K , n͘ in the nth band of the slab. The thus obtained 0 ͑z , zЈQ , ͒ has the required dimension ͑length͒ −4 ϫ ͑energy͒ −1 and is independent of the box quantization length L.
In the following, we shall take n ͑z͒ calculated in the geometry corresponding to a thick atomic slab with the lefthand side ͑lhs͒ and right-hand side ͑rhs͒ crystal edges ͑sur-faces͒ located at z L Ͻ 0 and z R = 0, respectively, and a large enough number of atomic layers in between so as to retrieve both the bulk and surface electronic properties of the crystal. We shall also assume that the electron density vanishes at a large distance z s from either crystal edge which enables us to expand the one-dimensional wave functions n ͑z͒ in a Fourier series of the form 24 where n =1,2,... takes the role of a band index, the distance d is given by and N and d 0 are the number of atomic layers and the interlayer spacing, respectively. This yields z L =−͑N −1͒d 0 . Due to the symmetry of the slab model potential that enters Eq. ͑5͒, the eigenfunctions n ͑z͒ are either even ͑with c n,l − =0͒ or odd ͑with c n,l + =0͒ with respect to the mirror symmetry plane z =−͑N −1͒d 0 / 2. In the case of the hereafter studied ͑111͒ surface of Cu we take N = 31 and z s =20d 0 with d 0 = 3.943 a.u. to fix d in Eq. ͑11͒, and l max in Eq. ͑10͒ to correspond to an energy of 150 eV.
In the limit N → ϱ the lhs crystal edge z L → −ϱ and the pairs of wave functions 2k−1 and 2k ͑k =1,2, ...͒ become degenerate. The limit of a very thick slab is appropriate to modeling of a semi-infinite crystal with the surface at z R =0 and in such situations of large N it may turn out more convenient to introduce an equivalent basis set of orthonormal wave functions that are localized on the left ͑L͒ or the right ͑R͒ surface of the slab. These wave functions are obtained by applying the transformation

͑12͒
For the above choice of the slab parameters the lhs and rhs crystal edge SS wave functions are given by Eq. ͑12͒ with k = 16, i.e., they are formed by taking the linear combinations of 31 in which the coefficients n 1 ,n 2 0,± ͑Q , ͒ are calculated following the procedure described in Ref. 24. The SC response function ͑z , zЈ ; Q , ͒ of interacting electrons can be expanded in the same type of series and expressed in terms of the coefficients n 1 ,n 2 ± ͑Q , ͒ which are determined by substituting Eq. ͑13͒ into integral equation ͑3͒ and solving the ensuing matrix equation in which the coefficients V n Љ ,n ٞ ͑Q͒ are obtained by integrating the product of 2D FT of the bare Coulomb interaction V͑z , zЈ ; Q͒ and the eigenfunctions of the Fourier series for 0 and .
In the present study of electron dynamics at surfaces we shall investigate how the excited quasiparticles are affected by the interactions with charge density fluctuations in the system that are described by . To this end we develop a method for renormalization of quasiparticles based on the equivalent boson representation that proves particularly convenient in the studies of quasiparticle dynamics in the real time domain. In doing so we start from the above outlined linear response formalism as a suitable framework for the calculation of quasiparticle self-energies at the level of GW approximation ͑GWA͒. In this scheme the short range exchange and correlation effects associated with the probe quasiparticle and the screening electrons are omitted which is justified in common practice by the large mutual cancellation between these effects in the final expressions for selfenergies ͑see discussion in Sec. 2.1.1 of Ref. 8͒. As demonstrated in Appendix A, resorting to this scheme is equivalent to studying the quasiparticle intra-͑nЈ = n͒ and inter-band transitions ͑nЈ n͒ induced by the interactions with bosonized excitations of wave vector Q and energy of the electronic density in the system. These excitations are modeled by a boson propagator D that has a Lehmann representation D n,n Ј ;n Ј ,n ͑Q,͒ = ͵ 0 ϱ dЈS n,n Ј ;n Ј ,n ͑Q,Ј͒ whose spectrum is calculated from ͑z 1 , z 2 ; Q , ͒ for interacting electrons as S n,n Ј ;n Ј ,n ͑Q,͒ = ͵ dz 1 ͵ dz 2 f n,n Ј ͑z 1 ,Q͒ ϫ ͫ − 1 Im ͑z 1 ,z 2 ;Q,͒ ͬ f n Ј ,n ͑z 2 ,Q͒.

͑16͒
Here the generalized oscillator strengths f n,n Ј ͑z 1 , Q͒ and f n Ј ,n ͑z 2 , Q͒ are defined in Eq. ͑A6͒ and the calculation of coefficients determining the expansion of expression on the rhs of Eq. ͑16͒ is demonstrated in Appendix A. This finally yields the SC RPA electronic excitation spectrum S n,n Ј ;n Ј ,n ͑Q , ͒ in the form given by Eq. ͑A12͒, with the dimension ͑length͒ −2 ϫ ͑energy͒ −1 . Note also that the role of S n,n Ј ;n Ј ,n ͑Q , ͒ in the phase space of quantum numbers ͑Q , n , nЈ , ͒ is analogous to the role played by the imaginary part of Lindhard's function in the case of a 3D noninteracting electron gas. The evolution of a wave packet comprising bosonized electronic excitations that build up the spectrum ͑16͒ is described by D n,n Ј ;n Ј ,n ͑Q , t͒ which is obtained as the inverse temporal Fourier transform of Eq. ͑15͒. One of the most important features of the thus defined SC response function for interacting electrons and its boson representation D is that they correctly reproduce the static limit of a retarded screened interaction between the probe charges and metal surfaces in the asymptotic form of a classical image potential ͑cf. Appendix B͒. This signifies that the selfconsistent linear response formalism satisfies the limit of perfect screening at surfaces which imposes important sum rules on the spectrum of surface electronic excitations. [27][28][29][30] In Sec. III we exploit the above-developed slab approach in the calculation of electron dynamics in Q2D bands localized at metallic surfaces. In order to apply the thick slab model to this situation one needs to express the spectrum ͑16͒ in the ͑L , R͒-basis spanned by the wave functions ͑12͒, i.e., compute it for the initial states localized at one surface of the slab and for the various combinations of final states partaking in intra-and inter-band transitions, viz. IS→ IS, IS→ SS, SS→ SS, and similar ones involving the bulk final states. This requires the representation of the corresponding oscillator strengths in the ͑L , R͒-basis which are derived at the end of Appendix A.
Quite generally, interacting electron gas in a thick slab may exhibit incoherent excitations ͑singleand multi-pair ex-citations͒ and coherent excitations ͑collective modes͒. In the SC RPA formalism described above the incoherent excitations are single electron-hole pairs and coherent excitations include bulk plasmons ͑BP͒, surface plasmons ͑SP͒, and multipole plasmons ͑MP͒ which all contribute to the spectral function S n,n Ј ;n Ј ,n ͑Q , ͒. Figure 2 shows the computed intensity of the n = nЈ = IS intraband component of the excitation spectrum S IS,IS ͑Q , ͒ obtained for the IS-band localized on one side of a 31 layer thick Cu͑111͒ slab, over the phase space of excitation energies and wave vectors Q relevant to our calculations. The contributions to the intensity come from the collective mode which disperses along a parabolic curve starting at the point ͑Q =0, = 7.6 eV͒ and dominates the spectrum, another distinct mode with higher energy whose dispersion curve starts at = 11.6 eV, and the electron-hole ͑e-h͒ quasicontinuum with nonvanishing intensity in the region encompassed by the parabolas = Q 2 /2m n ± Qv F,n and 0 Յ Յ Qv F,n − Q 2 /2m n ͑v F,n is the Fermi velocity in the nth band͒. Here the assignment of the various collective modes are additionally complicated by the effective electron masses m n different from the free electron value due to which the magnitudes of surface plasmon frequency s and bulk plasmon frequency p need not be equal to those of the equivalent jellium model of the slab characterized by the corresponding free electron density parameter r s . To assist assignments and estimate the relative contributions of collective excitations to the surface excitation spectra we show in Fig. 3 the calculated shapes of the spectrum S IS,IS ͑Q , ͒ for small values of wave vectors from the interval 0 Ͻ Q Յ 0.09 a.u. in which the e-h continuum component that is linear in Q is still small. While the maxima developing from the peak centered at ϳ7.6 eV for small Q can be clearly identified with the monopole surface plasmon excitations in the slab, the assignment of the mode giving rise to the peaks at ϳ11.6 eV is less straightforward. Their very weak dispersion and the enhancement of intensity with the increase of Q point towards the multipole surface plasmon character of this mode 31,32 rather than of the bulk plasmon. This assignment is additionally corroborated by inspection of the spectra S n,n Ј ;n Ј ,n ͑Q , ͒ involving IS→ bulk band transitions and intrabulk band transitions ͑not shown here͒ which exhibit a third peak at ϳ12.5 eV that can be assigned the bulk plasmon character. Further elaboration of this issue requires analyses going beyond the scope of the present work in which the various S n,n Ј ;n Ј ,n ͑Q , ͒ enter as numerical inputs.
The n = nЈ = SS intraband component S SS,SS ͑Q , ͒, shown in Fig. 4 exhibits similar features. However, owing to the larger oscillator strengths f SS,SS Ͼ f IS,IS the overall magnitude of S SS,SS ͑Q , ͒ is larger than S IS,IS ͑Q , ͒ by a factor ϳ3. In contrast to this, the interband n =IS,nЈ = SS surface excitation spectrum S IS,SS ͑Q , ͒ shown in Fig  = 0 due to the orthogonality of wave functions n ͑z͒, and for Q Ͼ 0 exhibits the magnitude that is smaller than S IS,IS ͑Q , ͒ by a factor ϳ45 because f IS,SS f IS,IS .

B. Image potential as an adiabatic final state effect
The effective potential V MP ͑z͒ employed in Eq. ͑5͒ to calculate the single-electron wave functions in the slab has been taken to incorporate the asymptotic form of the semiclassical image potential outside the slab, as given for the rhs surface by Eq. ͑6͒. As a consequence, a series of states localized at either surface is produced by a combined effect of the long ranged image potential tail and the surface projected bulk band gap. However, an external probe charge e that is nonadiabatically ͑suddenly͒ promoted to a distance z outside the metal surface does not initially feel the image potential that would arise due to instantaneous formation of the image charge density in the surface region. Rather, the quantum analog of classical image potential appears as a result of the retarded interaction of the external charge with the induced screening charge described by D n,n Ј ;n Ј ,n ͑Q , t͒ that, according to the results presented in Sec. II A and Appendix A, builds up in a finite time interval, typically of the order of inverse surface plasmon frequency s −1 for which the spectra ͑16͒ exhibit maxima ͑cf. Figs. 2-5͒. In other words, the full semiclassical form of the image potential is reached only after the external charge has interacted sufficiently long with the electronic density excitations in the metal in order to develop a fully relaxed screening charge density. Thus in PE spectroscopy the creation of a positive hole charge in an initially occupied state at the surface is followed by the formation of negative screening charge density which gives rise to a final state image potential energy shift of the photoinduced hole. On the other hand, in IPE and the first step of 2PPE the image potential appears as a relaxation effect arising from the positive charge density induced in the metal in response to the electron promotion into an empty state ͑for the discussion of formation of image potential shifted energy levels see Refs. 33 and 34, and the Introduction of Ref. 36͒. Hence a natural choice of a complete set of states used in the linear response description of dynamical screening at surfaces may be the set of eigenstates of the one-body potential Ṽ MP ͑z͒ that does not initially incorporate image potential effects. However, as a time-dependent perturbation approach starting from this basis of initial states may converge very slowly or require summations of infinite series of the scattering contributions, resorting to different schemes proves necessary. A similar problem arises also in connection with energy level renormalizations in the field theoretical formulations of scattering 37 and we shall adopt the same method of solution in the present treatment.
We proceed by assuming that the energy spectrum of the unperturbed Hamiltonian which describes the probe electron moving inside and near the boundary of a semi-infinite metal is solved in the same fashion as implemented in the derivation of the electronic response function in Sec. II A. The effective interaction U crys ͑z͒ of the probe particle with the periodic structure of the crystal is complemented with the adiabatic ͑static͒ surface correction which is nonvanishing in the interval z Ͼ z 0 outside the crystal surface located at z = 0. Here f 1 ͑z − z 0 ͒ is a kind of Tang-Toennies attenuating function 38 which guaranties a smooth match of the asymptotic form of the semiclassical image potential v͑z͒ given by Eq. ͑6͒ with U crys ͑z͒ at z = z 0 , and which saturates at unity for z → ϱ. Thus the effective potential acting on the probe electron is the same as the one appearing in Eq. ͑5͒ and given by V MP ͑z͒ = U crys ͑z͒ + U͑z͒. In the case of slabs of large but finite thickness an analogous correction is introduced also at the opposite surface. This enables construction of a Hamiltonian H 0 for noninteracting electrons in the second quantization representation in terms of the electron creation and annihilation operators introduced in Eq. ͑7͒ and the one-electron energies Eq. ͑9͒ as Here the eigenstates ͉K , n͘ and eigenenergies E K,n corresponding to H 0 incorporate the final state effects of image potential on the noninteracting electrons as generated by the Schrödinger equation ͑5͒.
However, when considering the effects of dynamic response of the electronic density to the introduction of a single electron ͑hole͒ into the system it is necessary to take into account that U͑z͒ is equal to the final state total energy correction or the energy relaxation shift which originates from the static limit of the dynamic many-body response of the metal to the presence of external charge ͑cf. Appendix A͒. Hence in order to avoid overcounting the same final state potential should be subtracted from the effects produced by the dynamical many-body interaction V of the electron ͑hole͒ with the electronic response. This can be achieved by making use of the procedure outlined in Sec. 5 of Ref. 37 where H bath is the unperturbed Hamiltonian of bosonized electronic charge density fluctuations constituting the heatbath. The interaction VЈ, obtained from V by subtracting the final static image effects, then reads where is the coupling constant and V K,n,K Ј ,n Ј is the matrix element of the afore-mentioned dynamic interaction of the excited quasiparticle ͑electron or hole͒ with the bosonized charge density fluctuations described by the response function ͑or equivalently, by the boson propagator D͒, and U n,n Ј = ͗ n ͑z͉͒U͑z͉͒ n Ј ͑z͒͘ with U͑z͒ ϰ 2 given by Eq. ͑17͒. By introducing the creation and annihilation operators a Q, † and a Q, , respectively, for bosonized charge density modes ͑Q , ͒, we find that V K,n,K Ј ,n Ј ϰ ͑a Q, † + a −Q, ͒ which is offdiagonal in the number of boson excitation quanta ͑cf. Refs. 35 and 36, and Fig. 14͒. Now, in a full quantum calculation of renormalization of quasiparticle energy in the adiabatic limit we start from definitions ͑18͒-͑20͒ and find that the first order Rayleigh-Schrödinger perturbation correction in VЈ for the energy level ⑀ n yields the shift −U n,n ϰ 2 , whereas V K,n,K Ј ,n Ј gives no contribution since by construction it is off-diagonal in the quantum numbers ͑Q , ͒. However, V K,n,K Ј ,n gives a contribution to the second order perturbation correction which by construction of the image potential shift within the linear response is ϰ 2 and cancels out the first order correction −U n,n in the limit of a quasiclassical particle ͑see Appendix B͒. On the other hand, the diagonal U n,n gives no contribution to the second order perturbation correction ͑neither to any higher order one͒, whereas the offdiagonal U n,n Ј produces a contribution ϰ 4 , likewise the fourth order corrections in V K,n,K Ј ,n Ј . All such terms that are ϰ 4 represent higher order corrections to the image potential and should be discarded in the quasiclassical limit in order to preserve consistency with the image potential terms included in Eq. ͑18͒ that are ϰ 2 . Hence the present choice of H 0 and VЈ eliminates overcounting of the final state image potential energy shifts up to O͑ 4 ͒ and preserves the consistency of perturbation treatment. In the adiabatic limit of quasiparticle motion this cancellation holds to all orders in . Thereby the final state image potential energy shifts are treated at the level of GW approximation.
We shall pursue the study of dynamical effects induced by the perturbation VЈ consistent with the above-discussed treatment of static image potential shifts arising from the same perturbation. This means that we should consider only those dynamical processes induced by the perturbation V which lead to the results equivalent to the renormalization of quasiparticle propagators through the response function ͑3͒. In the energy representation this leads to the quasiparticle selfenergy renormalization up to the order 2 , which was demonstrated in Appendix A to be equivalent to the description of quasiparticle interactions with bosonized charge density excitations represented by the propagator D n,n Ј ;n Ј ,n ͑Q , Ј͒ ͓cf. expression ͑A8͒ and Fig. 14͑b͔͒. This provides a rationale for resorting to the bosonization approximation in the treatment of ultrafast phenomena at surfaces which will be developed in the next section. Clearly, the main advantages of such an approach are a good representation of the selfconsistent screened interparticle interactions and a justification of the ad hoc inclusion of the image potential U͑z͒ in H 0 as a quasiparticle final state energy relaxation shift originating from the second order perturbation treatment of the same interactions.

III. DYNAMICS OF QUASIPARTICLES INJECTED INTO SURFACE BANDS
Quantum mechanical transition amplitudes describing the processes of PE, IPE, and 2PPE from surfaces can be expressed in terms of the propagators of excited quasiparticles ͑electrons and holes͒ whose dynamics is renormalized by the interactions with the environment ͑heatbath͒ and among each other ͑cf. Refs. 40-47͒. Effects of the various interactions involving the heatbath of the system can be conveniently illustrated on the example of the pump-probe pulse induced 2PPE from SS-bands on ͑111͒ surfaces of noble metals because these events are affected by the decoherence and decay processes associated with the evolution of both the electron and the hole photoexcited in surface bands. In the representation of the above-described bosonized interaction of quasiparticles with the heatbath and among themselves the amplitude of the pump and probe pulse induced 2PPE transition in a two band model is schematically illustrated in Fig. 6. Here we assume that prior to the pump photon absorption the system is in a neutral ground state ͉0͘ in which the motion of electrons is described by H 0 given in Eq. ͑18͒, and ͗0 ͉ V ͉ 0͘ = 0 by construction. Annihilation of a pump photon gives rise to creation of an electron in one of the unoccupied states above the Fermi energy and a hole in the SS-band, i.e., to a formation of an intermediate electron-hole ͑e-h͒ pair in an eigenstate of H 0 . Now, as the injection of a quasiparticle in an excited state of the system effectively switches on the interaction VЈ ͓cf. Eq. ͑20͔͒, the latter will act so as to destroy the primary coherence of the photoexcited e-h pair in the course of time ͑for a more detailed discussion see Sec. I of Ref. 35͒. Consequently, the photons from the time delayed second pulse, which eject electrons from the intermediate states into the states above the vacuum level, act as a probe for the decohered states of intermediate e-h pairs. Note also that since we are interested in the response of the system to pulsed photon fields, the relevant quantities to study are the populations of final photoelectron states ͉f͘ above the vacuum level rather than the steady state photoelectron currents as is the case in standard formulation of photoemission in terms of current-current correlation functions. 41,42,47 The required population of a final state ͉f͘ is obtained from the absolute square of the sum of all 2PPE transition amplitudes of the kind shown in Fig. 6, and this is equivalent to the population formulated in the density matrix approach utilizing the Keldysh nonequilibrium Green's function method developed in Ref. 46.
Quasiparticle interactions with the bosonized heatbath give rise to two types of renormalizations of the transition amplitudes by the propagator D denoted by the wiggly lines in Fig. 6. The heatbath-mediated interparticle interactions illustrated by the leftmost wiggly line give rise to processes described by interband vertex corrections, whereas the interactions in which the quasiparticles emit real and virtual quanta of the heatbath excitations give rise to processes described by the dressing of single quasiparticle propagators. Earlier estimates of the effects of the particle-heatbath interactions on the transition amplitudes showed that in the ultrafast regime major contributions to the decay and decoherence of intermediate states in 2PPE arise from single quasiparticle interactions with the heatbath ͑see Figs. 1 and 2 in Ref. 9͒. In the language of renormalization of transition amplitudes by heatbath excitations, these are the processes arising from the self-energy-type of renormalization of the propagators of excited quasiparticles throughout the inter-val͑s͒ in which the interaction VЈ is effective. Hence in the following we shall concentrate on the role which single quasiparticle decoherence and decay processes play in the various stages of evolution of photoexcited states partaking in the PE, IPE, and 2PPE events.
We assess the ultrafast dynamics of electrons and holes promoted in surface bands from the survival probability L K,n ͑t͒ which yields information on the evolution of a quasiparticle upon its injection into a 2D momentum eigenstate ͉K , n͘ within the nth surface band. Creation of a charged quasiparticle in a state ͉K , n͘ of the initially neutral system switches on the interaction VЈ which causes the decoherence and decay of the quasiparticle state in the course of time. For electron ͑hole͒ promotion into an empty ͑occupied͒ band the survival probability is defined as L K,n ͑t͒ = ͉͉͗⌿ K,n 0 ͑t͉͒⌿ K,n ͑t͉͒͘ 2 , ͑21͒ where ⌿ K,n 0 ͑t͒ and ⌿ K,n ͑t͒ are the unperturbed and perturbed wave functions of the system subsequent to the quasiparticle injection into the nth band at t = 0, whose temporal evolutions are governed by H 0 + H bath and full H defined in Eqs. ͑18͒ and ͑19͒, respectively. Assuming these initial conditions we may write and ͉⌿ K,n ͑t͒͘ = e −iHt c K,n † ͉0͘. ͑23͒ This enables expressing the survival amplitude as ͉͗⌿ K,n 0 ͑t͉͒⌿ K,n ͑t͒͘ = ͗0͉c K,n e iH 0 t e −iHt c K,n † ͉0͑͘t͒ = G 0 * ͑K,n,t͒G͑K,n,t͒, ͑24͒ where G 0 ͑K , n , t͒ and G͑K , n , t͒ are the unperturbed and perturbed Green's functions of a single electron or a hole in the state ͉K , n͘ of H 0 , respectively. Hence expressions ͑21͒-͑23͒ define the quantum mechanical probability that a quasiparticle, prepared at t = 0 in the eigenstate ͉K , n͘ of H 0 and subjected to a perturbation VЈ arising from the coupling to charge density excitations in the system, retains its identity and be recovered in the same initial state at a later instant t Ͼ 0. Expression ͑24͒ establishes the sought for relation between the survival amplitudes and quasiparticle propagators constituting the PE, IPE, and 2PPE spectra and yields. In the case of single electron injection into an empty band the propagator G͑K , n , t͒ can be expressed as 19,48,49 G͑K,n,t͒ = G 0 ͑K,n,t͒exp͓C͑K,n,t͔͒, ͑25͒ where G 0 ͑K,n,t͒ = − i exp͑− iE K,n t͒͑t͒ ͑ 26͒ is the unperturbed propagator of an electron of effective mass m n , charge e, and energy E K,n given by Eq. ͑9͒. Hence is obtained as a sum of all cumulants C l ͑K , n , t͒ generated by the interaction VЈ defined in Eq. ͑20͒. Once the boson representation of VЈ is established the cumulants C l ͑K , n , t͒ can be readily calculated in terms of the propagators G 0 ͑K , n , t͒ and D n,n Ј ;n Ј ,n ͑Q , t͒ following the method described in I. Employing this procedure we were able to estimate that for the coupling strengths typical of quasiparticles in Q2D surface bands the cumulant series ͑28͒ is represented with high accuracy by the sum of first and second order cumulants. 19 Hence in the following we shall restrict our calculations to the derivation of quasiparticle propagators in this approximation. This approach fulfills the unitarity condition and enables a systematic treatment of intra-and inter-band quasiparticle transitions on equivalent footing. For the interaction ͑20͒ the first order cumulant is trivially given by C 1 ͑K,n,t͒ = − i͑− U n,n ͒t, ͑29͒ and since being a purely imaginary function it can give rise only to a level shift ϰ 2 that was discussed in Sec. II B. The relevant information on quasiparticle evolution is contained in the second order cumulant which is a complex function of t. Using the method presented in Sec. 3.1 of Ref. 19 for the calculation of second order cumulants involving bosonized interactions, we find that C 2 ͑K , n , t͒ can be expressed in terms of the density of substrate electronic excitations S n,n Ј ;n Ј ,n ͑Q , ͒ derived in Sec. II A. Following that procedure we obtain for the real and imaginary parts of the second order cumulant 16 Re C 2 ͑K,n,t͒ = − 1 Im C 2 ͑K,n,t͒ = − 1 Here V Q =2e 2 / Q and the index nЈ denotes unoccupied parts of the bands allowed in intra-and inter-band electron transitions so that E K+Q,n Ј Ն E F where E F is the Fermi energy. The contribution of U n,n Ј to the second cumulant is of the order O͑ 4 ͒ and consistent with our earlier approximation for the image potential effects it will be neglected in comparison with Eqs. ͑30͒ and ͑31͒ which are ϰ 2 . Reformulation of Eqs. ͑30͒ and ͑31͒ to describe the dynamics of a hole in an occupied Q2D surface band is carried out on noting that the hole excitation energies are negative, i.e., in the case of a hole excitation the quasiparticle energy differences appearing in Eqs. ͑30͒ and ͑31͒ are replaced by E K,n hole − E K+Q,n Ј hole = ͑K + Q͒ 2 /2m n Ј − ͑K͒ 2 /2m n + ⑀ n Ј − ⑀ n , where ͑K , n͒ and ͑K + Q , nЈ͒ now range over the occupied states of the respective bands. In the following the restrictions on the final electron and hole quantum numbers in summations over Q and nЈ on the rhs of Eqs. ͑30͒ and ͑31͒ and expressions ͑33͒-͑36͒ deriving thereof are implicit. It is easily verified that the obtained cumulants are dimensionless quantities independent of the quantization length L. At this stage a word of caution is required concerning the use of cumulant expansion in the calculation of survival probabilities leading to Eqs. ͑30͒ and ͑31͒. The representation of single particle propagators in the form ͑25͒ is possible provided the initial quasiparticle energy E K,n is sufficiently far from the Fermi level of the system so as that the assumption of a single quasiparticle propagating in only one time direction and interacting with bosons may be applied. A reasonable requirement for the validity of this approximation is that in the course of interaction the magnitude of the mean energy transfer to the heatbath be smaller than ͉E K,n − E F ͉. This implies that interactions with bosons should not cause a drop of the quasiparticle to the vicinity of E F where the retarded propagators ͑25͒ and ͑26͒ must be replaced by the causal ones that allow for propagation in both time directions. As was found in Ref. 9 this condition is safely fulfilled in the case of electrons in IS bands on Cu͑111͒, whereas in the case of SS holes it holds for the states in the lower half of the occupied part of the band.
Quite generally, Re C 2 ͑K , n , t͒ describes the decay and Im C 2 ͑K , n , t͒ describes the energy renormalization and dephasing of the initial quasiparticle state ͉K , n͘. The early evolution of Re C 2 ͑K , n , t͒ in the ultrashort interval set by the Heisenberg uncertainty follows the universal behavior of the form −t 2 /2 n 2 + O͑t 4 ͒. This gives rise to a ballistic or "Zeno-like" initial decoherence of the survival probability 50,51 L K,n ͑t → 0͒ → exp͓− t 2 / n 2 + O͑t 4 ͔͒, ͑32͒ where n −2 = 1 In the same temporal interval the energy renormalization and dephasing processes described by the linear and nonlinear terms on the rhs of Eq. ͑31͒, respectively, cancel out each other so that Im C 2 ͑K , n , t → 0͒ = O͑t 3 ͒ and only the energy shift described by the first cumulant ͑29͒ survives. As the time grows larger the nonlinear component in t on the rhs of Eq. ͑31͒ builds up very fast and then saturates at a finite value, whereby it gives rise only to a pure phase change. In this regime the linear in t component on the rhs of Eq. ͑31͒ that describes renormalization of the quasiparticle energy through the processes quadratic in V is by construction cancelled out by the first cumulant ͑29͒ up to the order O͑ 4 ͒. In other words, the effect of interaction ͑20͒ on the imaginary phase of the renormalized electron propagator is to produce an early transient which saturates already in a subfemtosecond time interval, but it does not give rise to a contribution to final state energy shift up to O͑ 4 ͒.
In the opposite limit of long times several scenarios concerning the temporal dependence of C 2 ͑K , n , t͒ are possible, depending on the variation of excitation spectrum S n,n Ј ;n Ј ,n ͑Q , ͒ across the resonant limit E K,n = E K+Q,n Ј + of the integrands in Eqs. ͑30͒ and ͑31͒. A nonvanishing and smoothly varying S n,n Ј ;n Ј ,n ͑Q , ͒ across the resonant limit gives rise to the behavior 16,19 C 2 ͑K,n,t ϳ ⌫ K,n −1 ͒ տ − ͑⌫ K,n /2 + i⌳ K,n ͒t − w K,n . ͑34͒ Here the decay rate or inverse lifetime that has the appearance of Fermi's golden rule ͑FGR͒, arises from the adiabatic limit ␦͑E K,n − E K+Q,n Ј − ͒t of the time dependent factor in the integrand on the rhs of Eq. ͑30͒. The energy relaxation shift which arises from the linear in t component in the integrand on the rhs of Eq. ͑31͒ has the form of a Rayleigh-Schrödinger correction to unperturbed energy that is ϰ 2 and builds up on the extremely short time scale s ϳ s −1 ϳ 0.2 fs. 36,52 As is shown in Appendix B, in the quasiclassical limit of adiabatic motion of the quasiparticle outside the surface this correction yields asymptotically the image potential shift ͑6͒ that is cancelled out by the contribution from the first order cumulant ͑29͒. The time-independent term w K,n is an off-resonant correction to the first two terms on the rhs of Eq. ͑34͒ and measures the nonadiabaticity of excitation processes induced by the interaction VЈ that is switched on with the quasiparticle injection into a state ͉K , n͘. In the opposite situation of S n,n Ј ;n Ј ,n ͑Q , ͒ varying discontinuously across the resonant limit or in the case of quasiparticle transitions to the threshold, here given by the lower bound of the final energy E K+Q,n Ј = E F , the asymptotic form of C 2 ͑K , n , t͒ exhibits a more complicated behavior leading to a nonexponential decay of the survival probability in the limit t → ϱ. [53][54][55] Hence further insight in the quasiparticle dynamics, both on the ultrashort ͑preasymptotic͒ and long ͑asymptotic͒ time scale requires the evaluation of integrals on the rhs of expressions ͑30͒ and ͑31͒ for concrete forms of S n,n Ј ;n Ј ,n ͑Q , ͒.
In the following we illustrate and quantify the abovediscussed different stages of ultrafast dynamics of quasiparticles for the benchmark surface Cu͑111͒. The desired information is deduced from the survival probabilities ͑21͒ which in the present approach are fully determined by the behavior of Re C 2 ͑K , n , t͒ because the imaginary parts of the cumulants cancel out in the absolute square of Eq. ͑27͒. Now, Re C 2 ͑K , n , t͒ is determined by the ͐dQ ͐ d K,Q ͐ d-integral over the product of unbounded time-independent factor W n,n Ј ;n Ј ,n ͑K,Q,͒ = Q ͑2͒ 2 ͉V Q ͉ 2 S n,n Ј ;n Ј ,n ͑Q,͒ ͑E K,n − E K+Q,n Ј − ͒ 2 ,

͑37͒
which becomes singular for E K,n − E K+Q,n Ј = , and a bounded time dependent factor ͑1 − cos͓͑E K,n − E K+Q,n Ј − ͒t͔͒, which may quench this singularity. Reexpressing Eq. ͑37͒ in the basis ͑12͒ we show in Fig. 7 a three-dimensional plot of W IS,IS ͑K , Q , ͒ corresponding to intraband electron transitions in the IS-band over the allowed segments of the ͑Q , ͒ plane for the initial quasiparticle wave vector K =0. Its structure is dominated by the maxima arising from the plasmon peaks in S IS,IS , and the singularities produced by the zeros of the denominator in Eq. ͑37͒. We find that the region around the infrared singularity at ͑Q =0, =0͒ gives the dominant contribution to Re C 2 ͑K =0,IS,t͒ in the femtosecond domain and thereby to decoherence of the quasiparticle. However, it does not contribute to its exponential decay because for K = 0 only the transitions into lower bands can give a contribution to ⌫ K,IS , as is evident from the argument of the ␦-function on the rhs of Eq. ͑35͒. For K Ͼ 0 the singularity moves away and its weight is redistributed over a larger segment of the ͑Q , ͒-phase space. 56 This effective increase of the phase space for quasiparticle scattering by substrate excitations gives rise to an onset of intraband contributions and thereby to enhancement of the total decay rate ⌫ K,IS . A notably different situation is encountered in the case of interband electron transitions from the IS band into the unoccupied part of the SS band on the same side of the slab that is described by W IS,SS ͑K , Q , ͒ and whose behavior is depicted in Fig. 8. Here due to the Pauli exclusion the maximum deexcitation or electron recoil energy is equal to ⌬ IS−SS = ͓⑀ IS − ⑀ SS − ͑K F SS ͒ 2 /2m SS ͔ for Q = K F SS , which then diminishes quadratically with the increase of Q. Hence W IS,SS ͑K , Q , ͒ is zero for Q Ͻ K F SS and its maxima for Q Ͼ K F SS that rise from the e-h continuum are associated with the zeros of the denominator in Eq. ͑37͒, as is clearly seen in Fig. 8. Figure 9 shows a three-dimensional plot of W SS,SS ͑K , Q , ͒ corresponding to a hole created at the bottom of the SS band. Here the dominant maxima in the plot are of the same origin as in Fig. 7, however, the whole structure is cut off at the Fermi wave vector in the SS band, Q = K F SS , which is the upper bound for the SS-hole momentum recoil. Since the largest phase space for intraband hole decay is for initial states near the band bottom ͑K =0͒ from which the holes can decay towards the Fermi level, the contribution of W SS,SS ͑K , Q , ͒ to the Q integral in Eq. ͑30͒ would generally diminish with the increase of K.
The final step in determination of Re C 2 ͑K , n , t͒ is the phase space integration of the multidimensional integrands on the rhs of Eq. ͑30͒ that appear in the form of products of W n,n Ј ͑K , Q , ͒ and the corresponding factors ͑1 − cos͓͑E K,n − E K+Q,n Ј − ͒t͔͒, in which the latter introduce a parametric time dependence. Animations of the intensity evolution of the integrands on the rhs of Eq. ͑30͒ as functions of Q ͑horizontal axis͒ and ͑vertical axis͒ for K = 0.01, which yield the temporal dependence of Re C 2 ͑K ,IS,t͒ and Re C 2 ͑K ,SS,t͒ in the interval 0 Ͻ t Ͻ 10 fs, are available online ͑see Ref.

57͒.
Temporal behavior of the survival probability ͑21͒ is now readily calculated for the initial band indices and values of K using the above-obtained results. Figure 10 shows L K,IS ͑t͒ together with its major contributions for the case of an electron created in the IS band on the Cu͑111͒ surface with the initial wave vector near the band bottom, K = 0.01 a.u. The early universal ballistic decay of the survival probability ͑32͒ is superseded by a superposition of oscillations arising from the nonadiabatic ͑off-resonant͒ excitations of surface plasmons in the slab and a gradual buildup of the w K,IS -modified FGR decay ͑34͒ arising from the resonant excitation of e-h pairs. Due to the off-the-energy-shell character of plasmon excitations their amplitude diminishes with the diminution of Heisenberg uncertainty for transfer of energy as t grows larger. All these features signify non-Markovian dynamics in the early evolution of quasiparticles that is inaccessible to asymptotic and adiabatic theories. We find that for K = 0.01 the interband IS→ SS transitions contribute about 39% to the total decay rate ⌫ K,IS , and the remaining 61% arise from interband transitions into the bulk bands, in agreement with the earlier calculations. 7 The contribution from intraband IS→ IS electron transitions to the resonant decay rate ⌫ K,IS is insignificant at this value of K as only very small initial kinetic energy is available for excitation of other electrons in   10. ͑Color online͒ Intra-and inter-band components and total survival probability L K,IS ͑t͒ for an electron promoted into the first image state band on Cu͑111͒ with initial state wave vector K = 0.01 a.u. ͑see legend for assignments͒. Asymptotic decay given by the interpolation fit of Eq. ͑42͒ and "bare" FGR exponential decay law are shown by thin dash-double dotted and dash-dotted curves, respectively.
the system. On the other hand, their role is dominant in offresonant excitation of surface plasmons and the buildup of nonadiabatic correction w K,IS = 0.15. Qualitatively similar trends in the behavior of survival probability are retained for initial K = K F SS /2ϳ 0.05 a.u., as illustrated in Fig. 11. Figure 12 shows various contributions to the survival probability ͑21͒ for a hole created in the SS band on Cu͑111͒ with the initial wave vector K = 0.01. Here the intraband SS → SS transitions contribute the major part ͑ϳ70% ͒ to the total decay rate ⌫ K,SS , and the remainder is due to interband transitions into the bulk bands, again in accord with earlier calculations. 7,8 Intraband transitions also give dominant contribution ͑ϳ60% ͒ to the short time behavior of L K,SS ͑t͒ reflected in the magnitude of the total nonadiabatic correction w K,SS = 0.21. A similar structure of the survival probability is recovered for a SS hole with initial K = K F SS /2ϳ 0.05 a.u., as shown in Fig. 13.

IV. DISCUSSION
The results for survival probabilities shown in Figs. 10-13 enable the identification of three distinct regimes of ultrafast dynamics of electrons and holes subsequent to their creation in the IS-and SS-bands on the Cu͑111͒ surface, respectively. The early ballistic or Zeno regime ͑0 Ͻ t Ͻ 1 fs͒ is followed by preasymptotic non-Markovian evolution with superimposed off-resonant excitation of surface plasmons and resonant excitation of e-h pairs. This structure persists up to t ϳ 10 fs and only past that time the offresonant plasmon excitations die out and the steady state evolution governed by the modified FGR decay ͑34͒ takes over. However, even long past that time the "bare" FGR decay commonly used in asymptotic descriptions of quasiparticle evolution is not yet reached. This signifies that manifestations of quasiparticle dynamics which become important on the ultrashort scale of time resolved experiments can be assessed only within preasymptotic treatments of electronic relaxation processes. A non-Markovian evolution of the survival probability decaying faster than Eq. ͑38͒ will manifest itself through the modifications of relative intensities of direct and indirect transitions in 2PPE and of the peak shapes in PE, IPE, and 2PPE spectra from surface bands relative to the ones deduced from the asymptotic approaches 18,46,47 utilizing Eq. ͑38͒. These features should be taken into account particularly in the interpretations of time resolved measurements. Time resolved 2PPE spectroscopy has enabled an unprecedented insight into the hot carrier dynamics in the bulk and surface electronic bands. Identifications and qualitative descriptions of the 2PPE signal from metal surfaces have been frequently carried out by resorting to optical Bloch equations ͑OBE͒ for modeling of the population dynamics of a restricted number ͑usually two or three͒ of the levels partaking in 2PPE. 17,18 Application of the OBE to interpretation of 2PPE from surface bands has enabled the assignments of features in the measured 2PPE spectra as well as the estimates of rate constants that control state populations and coherences in the asymptotic regime. However, the OBE employed in these analyses are based on two assumptions that render the whole approach highly phenomenological. ͑i͒ Surface bands partaking in 2PPE processes are modeled by dis-   crete energy levels whose populations are governed by the rate constants which give rise to a Markovian decay of quasiparticles of the form ͑38͒ during each stage of the experiment, and ͑ii͒ these rate constants are represented by a sum of the rate constants describing the decay ͑in that context interlevel transitions͒ and the dephasing ͑intralevel elastic scattering͒ of quasiparticles. To this end, two types of rate constants ⌫ n,n Ј that model the decay and dephasing processes are introduced in the OBE and varied in order to fit the experimental data. However, the validity of such descriptions is by construction restricted to the regime of separate time scales of perturbation and relaxation of the system 6,18 and to the use of asymptotic representation of quasiparticle decay in the form ͑38͒. Hence the analyses of quasiparticle dynamics on the femtosecond scale require a treatment that is free from these limitations. In contrast to the phenomenological approaches embodied in the OBE, the theoretical description of excited quasiparticle propagation developed in Secs. II and III treats both the intra-and inter-band transitions during all stages of quasiparticle evolution on equivalent footing. Thereby it avoids a division of irreversible quasiparticle relaxation into separate decay and dephasing processes that is justified only in the Markov approximation. However, conditions for experimental observation of the different stages of ultrafast quasiparticle dynamics discussed above may be technically very demanding and therefore critical. Whereas it would certainly be desirable to detect non-Markovian oscillatory transients in the preasymptotic evolution of excited quasiparticles, it may be more feasible to observe a coarse grained ͑CGR͒ behavior of transitions from the early ballistic to the intermediate regime in which the quasiparticle decay is described by modified FGR law ͑MFGR͒. Representing C͑K , n , t͒ in the early interval by 58 C CGR ͑K,n,t͒ = − t 2 2 K,n 2 , ͑39͒ and in the intermediate regime by one can obtain a temporally coarse grained interpolation of the relevant real part of C͑K , n , t͒ in the femtosecond regime in a Padé approximantlike form: This leads to the following three-parameter interpolation fit for the survival probability: L K,n fit ͑t͒ = exp ͫ − ͑t 2 / K,n 2 ͒͑⌫ K,n t + 2w K,n ͒ ͑t 2 / K,n 2 + ⌫ K,n t + 2w K,n ͒ ͬ . ͑42͒ Such fits of the computed survival probabilities are also shown in Figs. 10-13. Interpolations based on Eq. ͑42͒ should be simple enough to allow also the fitting of quasiparticle decay estimated from the time resolved experiments. Hence it may be of considerable interest to exploit expres-sion ͑42͒ in the analyses of experimental data in order to demonstrate limitations of the phenomenological form of quasiparticle decay ͑38͒ used so far in the interpretations of 2PPE amplitudes and yields. 17,18,44,46,47 In summary, the microscopic model of ultrafast hot carrier dynamics combined with self-consistent calculations of electronic response developed and described in Secs. II and III provide unified treatments of quasiparticle evolution in surface bands in which the various mechanisms that control ultrafast phenomena can be traced back to fundamental quantum processes rather than being described by phenomenological parameters. The need for such descriptions was recognized in the first fully quantum-mechanical formulations of 2PPE from surfaces 44,46 but concrete realizations were hindered by the lack of adequate microscopic descriptions of ultrafast relaxation dynamics. Hence further progress in the interpretations of PE, IPE, and 2PPE spectra from surfaces should be sought in terms of the above elucidated quasiparticle evolution with which the corresponding approximate descriptions based on asymptotic and phenomenological rate constants need be correlated and compared. The work on realization of this task is in progress.

APPENDIX A: BOSONIZATION OF THE INTERACTION OF QUASIPARTICLES WITH THE RESPONSE OF A DIELECTRIC SURFACE
The problem of mapping the interaction of a probe electron with the charge density fluctuations in a metallic slab to an equivalent electron-boson interaction can be conveniently illustrated on the example of calculations of electron selfenergies. We start from the expression for the self-energy ⌺ n ͑K , ͒ of an electron in the nth surface band at the level of GW approximation regarding the final state image potential energy shifts and short range exchange-correlation effects. 7,8 Exploiting the translational invariance of the system in the x and y directions parallel to the surface we obtain the expression for the probe electron self-energy in the form ⌺ n ͑K,͒ = ͵ dz ͵ dzЈ n ͑z͒⌺͑z,zЈ;K,͒ n ͑zЈ͒, ͑A1͒ which represents a correction to the unperturbed energy given by expression ͑9͒. By introducing the surface response function defined by Eqs. ͑1͒ and ͑2͒ expression ͑A1͒ can be written as ͓cf. Fig. 14͑a͔͒ ⌺ n ͑K,͒ = 1 ϫ n Ј ͑z͒V͑z 1 ,z,Q͒͑z 1 ,z 2 ;Q,Ј͒V͑z 2 ,zЈ,Q͒ is the unperturbed retarded Green's function of the probe electron in the nth unoccupied band, and the z-component of the 2D Fourier transform of the bare Coulomb interaction e 2 / ͉r 1 − r ͉ = e 2 / ͱ͑ 1 − ͒ 2 + ͑z 1 − z͒ 2 is according to definitions of Fourier transforms in Sec. II A given by that has the dimension ͑energy͒ ϫ ͑length͒ 2 . Next we introduce generalized oscillator strengths f n,n Ј ͑Q,z 1 ͒ = ͵ dz n ͑z͒e −Q͉z−z 1 ͉ n Ј ͑z͒, ͑A6͒ which enable us to represent ⌺ n ͑K , ͒ in the form ⌺ n ͑K,͒ = 1 1 ͵ dz 2 f n,n Ј ͑z 1 ,Q͒ ϫ͑z 1 ,z 2 ;Q,Ј͒f n Ј ,n ͑z 2 ,Q͒G n Ј ͑K − Q, − Ј͒.

͑A7͒
This expression, which has the dimension of energy and is independent of the box quantization length L, can be visualized as the self-energy renormalization of the single electron Green's function G n ͑K , ͒ through a boson propagator D n,n Ј ;n Ј ,n ͑Q , Ј͒ describing the projection of electronic excitations of the system of wave vector Q and energy Ј onto the nth band ͑cf. Fig. 14͒, viz.

͑A15͒
in which the values of A, B, C, and D are given by where ␦ p,q are Kronecker symbols and L 1 = In the case of electronic excitation spectra ͑16͒ projected onto the band states localized on one side of a thick slab, viz. S IS,IS ͑Q , ͒, S SS,SS ͑Q , ͒, S IS,SS ͑Q , ͒, etc., we must first compute the oscillator strengths describing the transitions between the states described by wave functions ͑12͒, viz.
where ␣ = L or R, which are symmetric with respect to permutation of indices of wave functions in the integrand. In the limit of a thick slab, z L → ϱ, the wave functions k ͑L͒ ͑z͒ and where ␣ = L and ␤ = R or vice versa, describe electronic transitions from the states on one side of the slab to the states on the other side. Hence the oscillator strengths in the ͑L , R͒-basis defined by wave functions ͑12͒ are readily obtained from the oscillator strengths in the n-basis spanned by wave functions ͑10͒. For a 31 layer thick slab discussed in Sec. II A we get for oscillator strengths corresponding to one surface of the slab: f IS,IS = 1 2 ͑f 33,33 + f 34,34 ͒, f SS,SS = 1 2 ͑f 31,31 + f 32,32 ͒, f IS,SS = 1 2 ͑f 33,31 + f 34,32 ͒, and analogously so for other intra-and inter-band transitions. Finally, the electronic excitation spectrum of a thick slab that is projected onto the states described by the wave functions ͑12͒ is obtained by substituting the oscillator strengths ͑A18͒ and ͑A19͒ in expressions ͑16͒ and ͑A11͒.

APPENDIX B: SEMICLASSICAL IMAGE POTENTIAL ENERGY SHIFTS FROM THE BOSONIZED RESPONSE
In this appendix we present a simple demonstration of consistency of the use of a bosonized interaction of particles with charge density excitations ͑cf. Appendix A͒ in the construction of dynamic interaction VЈ ͓Eq. ͑20͔͒ from which the final state image potential shifts have been subtracted due to their inclusion in H 0 . We consider a semiclassical limit of the interaction of an external point charge e located at r a = ͑ , z a Ͼ 0͒ with the semi-infinite interacting electron gas occupying the halfspace z Ͻ 0. This represents a semiclassical analog of the classical image potential in which the level shift ͑36͒ appearing in the imaginary part of the second cumulant ͑31͒ should exactly cancel out the contribution −U n,n from the first cumulant ͑29͒. To demonstrate this property in the semiclassical limit we model the z component of the point charge density distribution with a ͑z͒ = ͉ n=a ͑z͉͒ 2 = ␦͑z − z a ͒ where z a Ͼ 0 is assumed to lie outside the tail of the electronic charge density extending across the surface. The diagonal matrix element U n,n calculated using Eq. ͑6͒ and the semiclassical density distribution a ͑z͒ reads U a,a = − e 2 4͑z a − z 0 ͒ . ͑B1͒ On the other hand, the semiclassical limit of the energy level relaxation shift ͑36͒ is obtained by setting E K,n = E K+Q,n Ј = E a in the denominator of the integrand and taking S a,a;a,a ͑Q , ͒ calculated from Eq. ͑A9͒ with nЈ = n = a and for the above point charge density distribution, viz.
To obtain S a,a;a,a ͑Q , ͒ from D a,a;a,a ͑Q , ͒ calculated with the semiclassical a ͑z͒ we observe that in this case z a Ͼ z 1 and z a Ͼ z 2 in expressions ͑A6͒ for the oscillator strengths entering expression ͑A9͒. This property enables factorization of the integration coordinates, yielding D a,a;a,a ͑Q,͒ = e −2Qz a ͵ dz 1 ͵ dz 2 e Qz 1 ͑z 1 ,z 2 ;Q,͒e Qz 2 .

͑B3͒
Recalling the definition of the surface response function to external charges [27][28][29] R͑Q,͒ = 2e 2 Q ͵ dz 1 ͵ dz 2 e Qz 1 ͑z 1 ,z 2 ;Q,͒e Qz 2 = 1 − ͑Q,͒ 1 + ͑Q,͒ , ͑B4͒ with ͑Q , ͒ denoting the surface dielectric function 20 for laterally homogeneous electron gas, we may write D a,a;a,a ͑Q,͒ = e −2Qz a Q 2e 2 R͑Q,͒. ͑B5͒ Now, the leading contributions to Eq. ͑B5͒ are obtained from the first two terms in the small Q-expansion of R͑Q , ͒. These terms were derived by Feibelman 32 who showed that the series expansion for R͑Q , ͒ assumes the form where ͑͒ is the long-wavelength limit of the surface dielectric function that is equal to the bulk dielectric function in the same limit. 20 The definitions of d ʈ ͑͒, which is the distribution of the currents induced parallel to the surface by external tangential electric field, and of d Ќ ͑͒, which is the centroid of the induced surface charge due to external electric field normal to the surface, are given in Refs. 32 and 59. The required S a,a:a,a ͑Q , ͒ can now be calculated from Eqs. ͑B5͒ and ͑B6͒ in the two pole approximation appropriate to the long wavelength response of simple metal surfaces. 59,60 This is effectuated by using the Drude expression ͑͒ =1 − p 2 / 2 for bulk dielectric function which enables representing expression ͑B6͒ near the plasmon poles in the form that is valid to the order O͑Q 2 ͒. Here the dispersion of monopole surface plasmon frequency for small Q is given by where z B is the position of the positive jellium background edge, and d Ќ ͑͒ is approximated by a multipole plasmon interpolation ansatz 58 This renders a model expression for R͑Q , ͒ that exhibits a correct static ͑ → 0͒ and extreme dynamical ͑ → ϱ ͒ limits required in the calculation of ⌳ a , the poles at monopole and multipole plasmon frequencies Q and m , respectively ͑cf. Fig. 3͒, but misses out on the effects of finite lifetime of plasmon excitations caused by their decay into electron-hole pairs which are not accounted for in Eq. ͑B8͒. Substituting Eq. ͑B7͒ in Eq. ͑B5͒, taking the imaginary part −͑1/͒Im D a,a:a,a ͑Q , ͒ to obtain the corresponding S a,a:a,a ͑Q , ͒, and using the latter in Eq. ͑B2͒ yields after integration over and summation over Q the energy shift ⌳ a sc = − e 2 4z a ͩ1+ z 0 z a ͪ+O͑z a −3 ͒ Ӎ − e 2 4͑z a − z 0 ͒ ͑B9͒ that cancels out the shift ͑B1͒ in the sum of the first and second order cumulants, Eqs. ͑29͒ and ͑31͒, calculated in the semiclassical limit. Hence the bosonized surface response yields the correct semiclassical image potential introduced in V MP ͑z͒ which, in turn, enables a construction of the dynamical quasiparticle-surface interaction ͑20͒ that avoids overcounting of the final state image potential energy shifts in the calculations of quasiparticle evolution.