# One-step theory of two-photon photoemission

###### Abstract

A theoretical frame for two-photon photoemission is derived from the general theory of pump-probe photoemission, assuming that not only the probe but also the pump pulse is sufficiently weak. This allows us to use a perturbative approach to compute the lesser Green function within the Keldysh formalism. Two-photon photoemission spectroscopy is a widely used analytical tool to study non-equilibrium phenomena in solid materials. Our theoretical approach aims at a material-specific, realistic and quantitative description of the time-dependent spectrum based on a picture of effectively independent electrons as described by the local-density approximation in band-structure theory. To this end we follow Pendry’s one-step theory of the photoemission process as close as possible and heavily make use of concepts of multiple-scattering theory, such as the representation of the final state by a time-reversed low-energy electron diffraction state. The formalism is fully relativistic and allows for a quantitative calculation of the time-dependent photocurrent for moderately correlated systems like simple metals or more complex compounds like topological insulators. An application to the Ag(100) surface is discussed in detail.

###### pacs:

78.47.D,78.47.J,79.60.-i## I Introduction

Different pump-probe photoemission experiments have been developed in recent years as powerful techniques to study the non-equilibrium dynamics of the electronic degrees of freedom in condensed-matter systems. LLM+05 (); CMU+07 (); PFB+08 (); BPW10 (); RHW+11 (); GDP+11 (); GPM+11 (); VTT+12 (); CDF+12 (); SYA+12 (); RVB+12 (); AGJH12 (); WHS+12 (); DAB+13 (); MBV+13 () A widely used variant of time-resolved photoemission experiments is given by two-photon photoemission (2PPE). AZ09 (); WSPD10 (); Fau12 (); RGK+14 (); BMA15 (); KRGH16 () A prominent application of 2PPE is to study ultra-fast demagnetization processes in magnetic solids. CAH+06 (); LLB+07 (); PSG+08 (); SPD+10 () Besides the general interest in such non-equilibrium phenomena, a variety of spectroscopic issues like dichroic phenomena, SPA+08 (); MPL+08 (); SGV+16 () spin-dependent life-times of electronic states, WSMD07 (); GRM+07 (); PSWD10 () quantum beats Sch07 (); MSS+11 () or delay-dependent spectral-line width Sch07 (); MSS+11 () are of great scientific interest.

While many theoretical investigations on time-dependent correlation effects in solids had been performed in recent years, FKP09 (); SL13 (); SKM+13 (); USL14 (); KSM+14 (); FNF14 (); SCK+15 (); Bon16 () their main focus has been on the detailed understanding of the electronic non-equilibrium dynamics of the solid and not so much on the photoemission process itself. Many-body calculations SL13 (); Bon16 () are often performed for simple systems and certain simplified model Hamiltonians, such as Hubbard-type models, to allow for an application of certain techniques such as dynamical mean-field theory GKKR96 (); KSH+06 (); Hel07 () or time-dependent density functional theory. KDE+15 ()

For a quantitative description of experiments, however, it is of utmost importance to additionally account for the effects of transition-matrix elements. This is essential, for example, to address 2PPE or pump-probe experiments performed with linearly or circularly polarized light. PSG+08 (); SPD+10 () The corresponding polarization-dependent effects visible in the measured spectra SPA+08 (); PSG+08 (); SPD+10 () are encoded in the matrix elements. Furthermore, the explicit consideration of the solid surface is inevitable since many pump-probe experiments are realized by pumping into surface states WSPD10 (); SHR+98 (); Fau12 (); NF14 (); NOH+14 (); GMTH14 () which serve as intermediate states. A comprehensive and quantitative theory of time-resolved photoemission, including 2PPE, must therefore account for surface-related, final-state and transition-matrix effects in a quantum-mechanically coherent and consistent way in addition to pure modeling of the time-dependent electronic structure.

Our time-dependent one-step theory introduced recently BRP+15 () represents the first step toward such a theory. From the very beginning, it is based on the formalism of Green functions on the Keldysh contour in the complex time plane. Kel65 () This provides the natural interface with many-body theory and thereby in principle allows to incorporate time-dependent correlation effects beyond a simple mean-field level or beyond the local-density approximation (LDA) of band-structure theory. For the time-independent equilibrium state, it is already possible to account for many-body correlations by combining the density-functional with dynamical mean-field theory (DFT+DMFT), MCP+05 () as has been demonstrated successfully for various correlated systems in the last years. BME+06 (); SFB+09 (); BMM+10 (); MBE13 () Moreover, the approach does quantitatively incorporate all matrix-element, surface-related and final-state effects and allows us to directly compare calculated spectra with corresponding experimental data.

The goal of the present paper is to apply the general formalism to the setup of 2PPE experiments. In this case we can benefit from the fact that the initial pump pulse is weak and can be treated perturbatively. This greatly simplifies the evaluation of the theory and provides us with a numerically tractable approach if, in addition, many-body correlations beyond the LDA can be disregarded. As a proof of principle, we consider 2PPE from the Ag(100) surface. The theory, however, is in principle applicable to a wide range of materials and problems. Examples comprise ultrafast demagnetization processes, which are often studied for moderately correlated systems. WSMD07 (); PSG+08 (); MSS+11 (); GMTH14 () On the same level of accuracy our theoretical approach works for Rashba systems or topological insulators which are recently of high scientific interest for spintronics applications. CAH+06 (); NF14 (); RGK+14 (); KRGH16 ()

## Ii One-step description of two-photon photoemission

Angle-resolved photoemission from a system in thermal equilibrium is conventionally treated by Pendry’s one-step model Pen74 (); Pen76 (); HPT80 () which describes the photoemission as a single coherent quantum process, including final-state multiple-scattering effects in the bulk and at the surface potential, dipole selection rules and effects of the transition-matrix elements in general by describing the photoemission final state as a time-reversed low-energy-electron-diffraction (LEED) state. Coulomb-interaction effects are accounted for on a mean-field-like level via the one-particle Green function in the local-density approximation of band-structure theory. HK64 () Adopting the sudden approximation implies that the description of the final state is disentangled from the multiple-scattering theory for the initial state. The one-step model has successfully been applied to a wide range of problems and spans photocurrent calculations ranging from a few eV to more than 10 keV GPU+11 (); GMU+12 (); MBE13a (); BMK+14 () at finite temperatures and from arbitrarily ordered BMM+13 () and disordered systems. BMM+10 () Strong electron correlations can be accounted for in addition via an improved many-body modeling of the initial-state Green function. BME+06 (); SFB+09 ()

There are already a few steps toward a general theory of time-resolved photoemission by Freericks et al., FKP09 (); MDF10 (); SKM+13 () and Eckstein et al. EK08 (); RFE16 () followed by work from other groups. Ing11 (); USL14 () Moreover, a first realistic description of two-photon photoemission has been worked out. UG07 () The main complication for a numerically efficient computation of a time-resolved pump-probe photoemission spectrum consists in the determination of the lesser Green function which depends on two independent time variables. This adds to the necessity to consider realistic geometries, e.g., a semi-infinite stack of atomic layers, to realistically model the surface region and to incorporate realistic electronic potentials typically obtained from band-structure methods like the Korringa-Kohn-Rostoker (KKR) method. Kor47 (); EKM11 ()

To address two-photon photoemission experiments, we here start from our recently proposed ab initio theory of pump-probe photoemissionBRP+15 () where the system is assumed to be exposed to a strong light pulse, described by a light-matter interaction term added to the Hamiltonian . This drives the system’s state out of equilibrium. The pump pulse is assumed to have a finite duration. After the pump electronic relaxation processes set in on a femtosecond time scale followed by slower relaxation mechanisms involving lattice degrees of freedom. The latter are disregarded here. Under these conditions the following expression for the transition probability can be derived: FKP09 (); BRP+15 ()

(1) | |||||

Here, refers to the quantum numbers of the photoelectron, is the envelope function for the laser probe pulse and the photoemission matrix element. For the final state we have adopted the sudden approximation and assumed that the Coulomb interaction of the high-energy photoelectron with the low-energy part of the system can be neglected. The time-dependent correlation function in Eq. (1) can be identified with the lesser Green function which depends on two time variables and is a matrix in the orbital indices and . It is given by an equilibrium expectation value but is time inhomogeneous as the Heisenberg time dependence is governed by the full and explicitly time-dependent Hamiltonian . Accordingly, the central problem consists in computing the lesser Green function which describes the temporal evolution of the electronic degrees of freedom after the pump pulse.

The calculations simplify substantially when considering a system of effectively independent electrons. In this case, the lesser Green function can be written as: BRP+15 ()

(2) | |||||

where denotes the Fermi distribution function and is a time just before the perturbation representing the pump pulse is switched on. and are the energy-dependent retarded and advanced equilibrium Green functions of the unperturbed system, i.e. before the time . The retarded Green function for the perturbed system must be obtained from the following integral equation

(3) |

where in its operator representation still may be an arbitrarily strong perturbation. The corresponding advanced Green function is simply given by .

To make use of these expressions we turn to the real-space representation and to a fully relativistic four-component formulation. Eq. (1) then reads as:

where is the component of the wave vector parallel to the surface, and is the kinetic energy of the photoelectron. Here one may use for the lesser Green function the expansion

(5) |

where denotes the relativistic spin-angular functions Bra96 () with the spin-orbit () and the magnetic () quantum numbers combined to =(). Furthermore, is the single-particle-like final state of the photoelectron as usual in the form of a time-reversed spin-polarized LEED state. The perturbation describing the probe pulse, which is assumed as weak, is given by:

(6) |

where is the electronic charge and using the dipole approximation denotes the spatially constant amplitude of the electromagnetic vector potential corresponding to the radiation field and its polarisation . The three components of the vector are defined as the tensor product for and denote the 22 Pauli spin matrices.

With this one can finally write the photocurrent within the one-step model as:

where represents the high-energy wave field Bra96 () and where the matrix element is defined by:

Summations run over the contributions coming from atomic layers I, atomic cells n in layer I and sites q within cell n.

In dealing with the perturbation (see Eq. (6)) it is in practise often more convenient to use the so-called gradient-V form. In this case, the corresponding single electron potential at () may depend in the most general case on the electronic spin. Explicit expressions for the corresponding matrix elements split into an angular and radial part which can be found elsewhere. Bra96 ()

To evaluate entering Eq. (LABEL:eq:mat.el) by use of Eq. (2) together with the Dyson equation (3), the real-space representation of the retarded Green function is needed. This is obtained by Fourier transformation from the energy-dependent retarded Green function which in turn may be evaluated in a direct way by means of the KKR or multiple-scattering formalism. MCVP89 (); EBKM16 ()

Here we assume that is in the atomic cell () while is in (). The first term in Eq. (II) is the single-site contribution to the Green function made up of the regular () and irregular () solutions to the single site Dirac equation for site () where the index n can be dropped because of the two-dimensional periodicity in a layer . The sign ”” distinguishes left-hand-side solutions for the Dirac equation from the standard right-hand-side ones. EBKM16 () The second term in Eq. (II) represents the back-scattering term which accounts for all scattering events between sites () and () in a self-consistent way. EBKM16 () To calculate the so-called structural Green function occurring in that term several techniques are available. EKM11 (); MCVP89 () Finally, the energy-dependent factor represents essentially the relativistic momentum. EBKM16 ()

In the following an application of the scheme to 2PPE spectroscopy will be presented. Accordingly, we now assume that both, the pump and the probe pulse , are weak in intensity. As mentioned before, this situation quantitatively describes the scenario of two-photon photoemission spectroscopy. Consequently, Eq. (3) can be solved perturbatively in first-order approximation by replacing by on the right side of Eq. (3). This leaves us with a simple integral expression for while is available from Eq. (II).

As a first step, we calculate the atomic-like contribution using the single-site part of Eq. (II) only. Substituting the Fourier transform of the single-scattering contribution into Eq. (3), this site-diagonal term is obtained as:

with the four matrix element functions

(11) | |||||

(12) | |||||

(13) | |||||

(14) |

Here the real space representation for the pump pulse has been split in analogy to Eq. (6) for the probe pulse and is the critical radius of the sphere bounding the atomic cell and use has been made of the fact that the system has two-dimensional periodicity. The symbols and in Eq. (II) are defined as and . In Eq. (11)-(14) and serve as dummy variables for the corresponding integration boundaries. Furthermore, for and the constraint must hold, otherwise they are zero. In addition to this atomic-like part, the retarded Green function leads to two mixed contributions between single-scattering and multiple-scattering events, as well as a double multiple-scattering contribution. First we present the two mixed contributions. The first one is given by:

(15) | |||||

Here the structural Green function accounts for all multiple-scattering events for the propagation from site to site .

For the second mixed contribution we find:

(16) | |||||

As the last bulk-like contribution the multiple-to-multiple-scattering part is obtained as

The total double-time dependent bulk-like radial part of the retarded Green function follows by summing the four different contributions presented above:

It remains to calculate the surface contribution of the time-resolved photocurrent. Especially in the case where the image-potential states GBB+93 () serve as intermediate states in the two-photon photoemission process, it is essential to describe all surface-related electronic states including the image states in a realistic way. In the static version of the one-step model this is usually realized by employing a Rundgren-Malmström barrier potential. MR80 () This approach can be generalized to the time-dependent case. We start with an appropriate formulation of the retarded free-electron Green function

(19) |

where denotes a two-dimensional reciprocal lattice vector and the corresponding wave vector for energy . The Fourier transform is

(20) | |||||

The scattering properties of the surface potential can be expressed in terms of the barrier scattering matrix . With the Kronecker delta, is a diagonal matrix, i.e., corrugation effects in the surface potential are neglected. MR80 () This matrix is represented as an additional surface layer in the formalism and is typically located in front of the first atomic layer at the distance , which defines the so-called image plane. MR80 () Given this matrix, plane-wave amplitudes and can be defined, where the first amplitude is emitted by the surface barrier and the second one is emitted by the semi-infinite stack of atomic layers. HPT80 (); Bra96 () The corresponding amplitude emitted from the semi-infinite bulk may be denoted by , and represents its reflected counterpart. These four amplitudes satisfy the following linear system of equations:

(21) |

where denote the free electron Green functions which propagate the wave field between the first atomic layer and the surface layer. The amplitude can also be calculated by standard multiple-scattering techniques. The result is:

This plane wave amplitude appears directly proportional to the gradient of the surface potential . Here, denotes the bulk-reflection matrix and represents the wave function of the intermediate state in the barrier region. The amplitude is also known from standard multiple-scattering techniques applied to the calculation of the initial-state wave function between the different layers of the semi-infinite bulk. HPT80 () Consequently, the two amplitudes and can be calculated from the system of linear equations which is defined above. For we find:

and results in:

(24) |

Having these coefficients calculated, the wave field of the intermediate state at the image plane can be expressed within a plane-wave representation:

The corresponding expansion in spherical harmonics gives:

(26) |

with the spherical coefficients

(27) | |||||

For an explicit calculation of the retarded surface Green function the radial part of the spherical wave field is needed. Unfortunately, the surface potential is a function of the Cartesian coordinate only, and so is the corresponding wave function. This means that the wave function is not directly available. Nevertheless, as a good approximation we define:

(28) |

where represents the spherical wave function which belongs to an empty-sphere potential of the first vacuum layer. This approximation works well because in self-consistent TB-SPRKKR EKM11 () electronic-structure calculation a set of vacuum layers is typically used, which is located on top of the first atomic layer. Due to the charge transfer from the first atomic layer into the empty sphere potentials of the first and second vacuum layers, these empty-sphere potentials represent in a reasonable approximation the polynomial region of the surface potential. Within the local spin-density formalism, the surface retarded Green function is then given by

(29) |

where the -function has been omitted since the condition is fulfilled. Inserting the expression for in Eq. (3), the retarded surface Green function in first-order approximation reads

with

(31) |

and where

(32) |

and

(33) |

The radial part of the surface matrix element is defined as:

The total double-time-dependent radial part of the retarded Green function follows by summing the five contributions discussed above:

(35) | |||||

Having calculated numerically, can be obtained via Eq. (2). Finally, with Eqs. (II) and (LABEL:eq:mat.el) one may compute the time-dependent 2PPE signal.

## Iii 2PPE from Ag(100)

As a first application the theory is applied to the (100) surface of Ag, i.e., to a prototypical simple paramagnetic metal. We compute the lesser Green function and the 2PPE spectrum within multiple-scattering theory in a fully relativistic way by using the Munich SPRKKR program package in its tight-binding version.SPR-KKR6.3 () The spherically symmetric potential was obtained within atomic-sphere approximation (ASA) and the corresponding single-site wave functions serve as input quantities for the calculation of the lesser Green function. The latter is obtained in two steps: First we determine the retarded Green function from the respective Dyson equation (3) where we treat the perturbation to first order in the pump-pulse strength. Second, the lesser Green function is calculated from Eq. (2). The solutions of the conditional equations for the Green functions are obtained by numerical matrix operations, where the expansion into spherical harmonics includes orbital quantum numbers up to . The two radial coordinates and are restricted to the ASA sphere. With respect to the dynamical degrees of freedom, the equations are Fourier-transformed from time to energy space. To this end we choose an equidistant mesh for the two time variables and with a time step of fs. The energy-dependent retarded KKR Green function is calculated for a complex energy, with an energy-dependent imaginary part in eV, to account for damping effects due to inelastic scattering events. Therewith, the life-time broadenings of the first and second images states, meV and meV, are accounted for quantitatively. SFFS92 () Concerning the numerical effort, this also helps to reduce the number of energy grid points necessary for a numerically accurate Fourier transformation. Converged results are obtained for an energy window of about 3.0 eV around the Fermi level .

Our calculations refer to the energy-resolved operational mode of an 2PPE experiment as schematically shown in Fig. 1. Pic07 () In this mode the time delay between the pump and probe pulses is fixed, and a kinetic-energy spectrum of the excited electrons is recorded. For a given kinetic energy , the energy of the intermediate state , with respect to the Fermi level , is obtained by from the photon energy of the probe pulse and the work function . The initial-state energy is