text:
# Response theory for time-resolved second-harmonic generation and two-photon photoemission

## I Introduction

During the last decade time-resolved spectroscopy of condensed-matter systems has become a very active area of experimental research Ste92 ; HMKB97 ; Sim98 ; KSR98 ; Gud99 ; Schoe ; Aesch1 ; Aesch2 ; Wolf1 ; Hoef97 ; Cao97 ; Petek ; Wolf2 ; Leh99 ; Kno00 ; GWH00 ; SBW02 ; ONP02 ; Gerb ; Ae.cluster ; Petekrev ; Fausterrev . This is mainly due to the progress in experimental technique, in particular the ability to create ultra-short laser pulses with a duration of the order of a few femtoseconds femto . Since this is similar to the relaxation times of excited electrons and collective excitations in solids, these experiments allow to study non-equilibrium physics, e.g., the time evolution of excited electrons before and during thermalization. Of particular interest are non-linear techniques such as time-resolved sum-frequency generation (SFG) and two-photon photoemission (2PPE), which are sensitive to excited electron states Shen . A theoretical understanding of these processes is crucial. Petek and Ogawa Petekrev  noted in 1997 that a theory for time-resolved 2PPE is still lacking, and, despite the efforts of many theorists, much remains to be done. The situation for SFG is similar. The construction of such a theory is a formidable task—the main problems are (a) the desription of the time-dependent response and (b) the treatment of the surface. Our main concern is with the first point. A simplified description of the surface using Fresnel factors has been employed successfully to describe SFG from metals SMD87 ; HBB94 ; Luce97 ; LuceLDA . A detailed discussion of boundary conditions at the surface, focusing on the nonlinear optical response of magnetic systems, can be found in Ref. AtK02 .

In the present paper we discuss the electronic processes taking place during time-resolved SFG (in particular second-harmonic generation, SHG) and 2PPE and derive the dependence of the SFG light intensity and the 2PPE photoelectron yield on the time dependence of the exciting laser field. We show that most effects observed for time-resolved 2PPE appear similarly for SFG, such as their dependence on energy relaxation, dephasing, and detuning of intermediate states. Other examples are the enhancement of the response due to collective excitations and the sensitivity regarding the ultra-fast spin-dependent relaxation. We develop a unified time-dependent response theory for SFG and 2PPE, starting from the self-consistent field approach of Ehrenreich and Cohen EC ; HB , which can be applied to specific materials described by their band structure, relaxation rates, and dipole matrix elements. For illustration, we apply the theory to a generic tight-binding model for a metal to study interference effects in both pump-probe single-color SFG and 2PPE and their dependence on relaxation rates and detuning. We exhibit the strong similarities between both methods.

In SFG Ste92 ; HMKB97 ; Sim98 ; KSR98 ; Gud99  electrons are excited by absorbing two photons and they subsequently emit a single photon at the sum frequency. In Fig. 1 we illustrate the type of process yielding SFG. For simplicity we talk about SFG in the following, although *difference*-frequency generation is automatically included in our theory. Time-resolved measurements Ste92 ; HMKB97 ; Sim98 ; KSR98 ; Gud99  usually employ the pump-probe technique, where two laser pulses of the same (single-color) or different (two-color) frequency are applied with a time delay $`\mathrm{\Delta }T`$ between them. This time delay controls the time between the two absorptions and thus the relaxation dynamics of the electron in the intermediate state $`|2`$ is crucial, see Fig. 1. SFG is strongly surface-sensitive, since the SFG response of the bulk of an inversion-symmetric crystal vanishes in the dipole approximation. The inversion symmetry can also be broken by nanostructures. The most important case of SFG is second-harmonic generation (SHG), where the electrons are excited by approximately monochromatic light of frequency $`\omega `$ and light of frequency $`2\omega `$ is detected. Note, in the case of ultra-short laser pulses the spectrum is necessarily broadened and a full treatment of SFG is required even for these single-color experiments. Also note that a single laser pulse, depending on its duration and shape, involves time-delayed absorptions.

Time-resolved 2PPE experiments of metal surfaces Schoe ; Aesch1 ; Aesch2 ; Wolf1 ; Hoef97 ; Cao97 ; Petek ; Wolf2 ; Leh99 ; Kno00 ; GWH00 ; SBW02 ; ONP02  as well as of clusters Gerb ; Ae.cluster  employing the pump-probe technique have been performed more often than time-resolved SFG. Reviews can be found in Refs. Petekrev  and Fausterrev . Figure 2 shows a sketch of the processes yielding 2PPE. An electron is excited above the vacuum energy $`E_{\mathrm{vac}}`$ due to the absorption of two photons. The interplay between the relaxation of the electrons in intermediate states and the time between the two absorptions will determine the resulting photoelectron current. The probability of electrons above the vacuum level actually leaving the solid is also crucial. The limited mean free path of the electrons makes photoemission surface-sensitive, but in general less than in the case of SFG. In both SFG and 2PPE interference effects Ste92 ; Sim98 ; Aesch2 ; Cao97 ; Petek ; Wolf2  appear, which our theory allows to study. Of course, these interference effects are expected to depend on the pulse shapes.

The response theory presented here goes beyond previous theoretical treatments of ultra-fast processes rem.slow  in SFG and 2PPE in metals, which mainly fall into four classes: (a) density functional theory and approaches based thereon Lie89 ; Ull97 ; LuceLDA ; Lie99 ; Koh99 ; WDS ; Camp , (b) rate equations Kno00 ; Rol1 ; Rol2 , (c) optical Bloch equations Wolf1 ; Petek ; Loudon , and (d) perturbative methods HB ; HBB94 ; PPK00 ; Sha00 . At first, density functional theory has been applied in the time-dependent local-density approximation for jellium models Lie89 ; Ull97 ; Lie99 ; Koh99 . In the jellium approximation one ignores the potential of the ion cores and, consequently, any band-structure effects. Thus this approach is not suitable if single bands or surface states or quantum-well states in thin films are important. On the other hand, collective excitations are usually described rather well Lieplas . Going beyond the jellium model, Luce and Bennemann have employed the local density approximation to calculate dipole matrix elements as they enter also in our approach LuceLDA . Additionally taking excited states into account within the GW approximation, Schöne et al. WDS  have calculated electronic lifetimes. Hole dynamics have also been studied with density functional methods Camp .

However, one would like to gain more general physical insight than the numerical results can provide. To this end one may consider rate equations for the occupation of excited states, e.g., the Boltzmann equation Kno00 ; Rol1 ; Rol2 . This approach allows to incorporate important effects such as secondary electrons due to relaxation from higher-energy states and to Auger processes as well as transport into the bulk Kno00 ; Rol1 ; Rol2 . However, rate equations neglect the electric polarization of the electron gas, its dephasing, and any quantum-mechanical interference effects, resulting from the superposition of the laser field and the induced fields. To include these effects one has to solve the equation of motion for the entire density matrix $`\rho `$, not only for its diagonal components, i.e., the occupations. This can be done in response theory. Its simplest form yields the optical Bloch equations: The system is modelled by a small number of levels and the von Neumann equation of motion (master equation) for the density matrix is integrated numerically Loudon ; Wolf1 ; Petek . However, this approach is limited to a small number of levels so that a realistic band structure cannot be described. Furthermore, many-particle effects like collective excitations are not included.

On the other hand, the response theory presented here does include the band structure and collective excitations. It generalizes the theory of Hübner and Bennemann HB  to SFG due to incident light of arbitrary time dependence and spectrum. The previous theory HB  has been used successfully for SHG from metal surfaces, thin films, quantum wells, and metallic monolayers due to continuous-wave, monochromatic light HB ; HBB94 ; Luce97 ; LuceLDA ; NOLI ; Luce96 ; Luce98 ; AnH02 . However, the dependence of SHG on the pulse shape and the effect of energy relaxation and dephasing were not discussed. We also derive the response expressions for time-dependent 2PPE within the same framework. Since our theory is explicitly formulated for continuous bands, it can also serve as a basis for the discussion of the averaging effects due to bands of finite width discussed in a more heuristic framework using optical Bloch equations for discrete levels in Ref. Weida . Since the full time or frequency dependence is included, effects of frequency broadening of short pulses and of finite frequency resolution of the detector (for SFG) Weida  are easily studied.

Our theory employs a generalized self-consistent-field approach EC ; HB , which is equivalent to the random-phase approximation (RPA) Dasgupta ; Lind54 ; HG70 ; Ch98 ; Mahan . We employ the electric-dipole approximation, which is valid for small wave vector $`𝐪`$ of the electromagnetic field and has been used successfully to describe SHG from metal surfaces HBB94 ; Luce97 ; LuceLDA ; NOLI ; Luce96 ; Luce98 ; AnH02 ; Shala . This is reasonable, since the skin depth, which is the length scale of field changes, is about one order of magnitude larger than the lattice constant. One has to take care in interpreting SFG experiments for inversion-symmetric crystals, since the surface contribution only dominates over higher multipole bulk contributions for surfaces of low symmetry HBB94 ; HB . Similar in spirit to our response theory, Ueba Ueba  has studied continuous-wave 2PPE from metal surfaces, Pedersen et alPPK00  have considered continous-wave SHG from metallic quantum wells, and Shahbazyan and Perakis Sha00  have developed a time-dependent, but linear response theory for metallic nanoparticles.

It is important to understand that at the surface of a metal, in thin films, and in nanostructures the light couples to collective plasma excitations. The field within the metal is of course not purely transverse Barton ; Board . Its transverse and longitudinal components couple with the conduction electrons to form plasmon-polaritons and plasmons, respectively Board . The (longitudinal) plasmon modes only decouple from the applied field for a structureless jellium model of the solid Barton ; Board . However, we consider a more realistic model that incorporates the crystal structure. Also, we will see that the induced nonlinear polarization couples to (longitudinal) plasmon modes.

On general grounds one may expect that the discussion of the intimate relationship between 2PPE and SFG also helps to understand the dependence of 2PPE on light polarization. It has been shown that the light-polarization dependence of SFG is important for the analysis of the electronic structure and magnetism Oxf .

The organization of the remainder of this paper is as follows: We first summarize the response theory for SFG and 2PPE in Sec. II. This lays the ground for our discussion in Sec. III of time-dependent SFG and 2PPE. Details of the response theory are given in appendices A and B.

## II Response theory

### II.1 Sum-frequency generation

We first outline the response theory for SFG. We consider a semi-infinite solid with single-particle states $`|𝐤_{}l`$ with energies $`E_{𝐤_{}l}`$ described by the momentum $`𝐤_{}`$ parallel to the surface, which is assumed to be perpendicular to the $`z`$ direction, and a set of additional quantum numbers $`l`$. For bulk states, which may be affected by the surface but are not localized close to it, the composite band index is $`l=(k_z,\nu ,\sigma )`$, where $`k_z`$ is the $`z`$ momentum component *in the bulk*, $`\nu `$ is a band index, and $`\sigma `$ is the spin quantum number. $`k_z`$ has a continuous spectrum. On the other hand, for states localized at the surface, $`l`$ is discrete. Examples are image-potential states, adsorbate states, quantum-well states in a thin overlayer, and proper surface states.

Part of the electron-electron interaction is included by the self-consistent-field approximation or RPA EC ; Dasgupta . The remaining electron-electron scattering is approximately taken into account by inserting phenomenological relaxation rates PS00  into the single-electron Green functions and by shifting the band energies $`E_{𝐤_{}l}`$ Louisell . We assume that $`E_{𝐤_{}l}`$ are quasiparticle energies containing these shifts. Note, the electron-phonon interaction only becomes relevant on longer time scales and is not considered here rem.slow . Also, *intra*band contributions to the response are not considered for simplicity, which is reasonable at optical frequencies.

The electrons are coupled to the effective electric field $`𝐄`$ within the solid through a dipolar interaction term (for simplicity we assume that the dipole coupling dominates). The optically induced polarization $`𝐏`$ within the solid is expanded in orders of the electric field $`𝐄`$. The linear response is given by

$$P_i^{(1)}(𝐪,t)=\frac{1}{2\pi }\underset{q_z^{}}{}𝑑t^{}\chi _{ij}(𝐪,𝐪^{};tt^{})E_j(𝐪^{},t^{}),$$
(1)
where $`\chi _{ij}`$ is the linear susceptibility, $`𝐪=(𝐪_{},q_z)`$, and $`𝐪_{}=𝐪_{}^{}`$ due to conservation of momentum parallel to the surface. Summation over repeated indices is always implied. The non-conservation of $`q_z`$ is explicitly taken into account.

We assume throughout that the photon momentum $`𝐪`$ is small compared to the dimensions of the Brillouin zone and that the band energies, relaxation rates, and transition matrix elements change slowly with momentum so that the difference between the parallel crystal momentum of an electron before and after the interaction, $`𝐤_{}`$ and $`𝐤_{}^{}`$, respectively, can be ignored. If we further neglect the frequency dependence of the transition matrix elements the self-consistent-field approach gives the time-dependent linear susceptibility

$`\chi _{ij}(𝐪,𝐪^{};tt^{})={\displaystyle \frac{e^2}{v}}{\displaystyle \frac{2\pi i}{\mathrm{}}}\mathrm{\Theta }(tt^{})`$ (2)
$`\times {\displaystyle \underset{𝐤_{}}{}}{\displaystyle \underset{ll^{}}{}}D_{𝐤_{}l^{};𝐤_{}l}^i(q_z)D_{𝐤_{}l;𝐤_{}l^{}}^j(q_z^{})`$
$`\times [f(E_{𝐤_{}l^{}})f(E_{𝐤_{}l})]\mathrm{exp}\left[i{\displaystyle \frac{E_{𝐤_{}l^{}}E_{𝐤_{}l}}{\mathrm{}}}(tt^{})\right]`$
$`\times \mathrm{exp}[\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}l^{}}(tt^{})],`$
where $`v`$ is the volume of the system. Note, the last two factors explicitly describe the oscillations and decay of the linear induced polarization. In the dipole approximation the transition matrix elements are

$$𝐃_{𝐤_{}l;𝐤_{}l^{}}(q_z)𝐤_{}l|𝐫|𝐤_{}l^{}.$$
(3)
The matrix elements are given without approximations in App. A. The linear susceptibility is represented by the usual Feynman diagram shown in Fig. 3.

The finite lifetime of electrons due to their interaction enters Eq. (2) through the dephasing rates $`\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}^{}l^{}}`$, which describe the decay of the superposition of states $`|𝐤_{}l`$ and $`|𝐤_{}^{}l^{}`$ and thus of the polarization. The change of the occupation of states is described by the energy relaxation rates $`\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}l}\tau _{𝐤_{}l}^1`$, where $`\tau _{𝐤_{}l}`$ are the lifetimes. $`\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}l}`$ is the rate of spontaneous transitions out of the state $`|𝐤_{}l`$. Since the depopulation of the states $`|𝐤_{}l`$ or $`|𝐤_{}^{}l^{}`$ certainly leads to the destruction of the polarization, the dephasing rates can be expressed in terms of the lifetimes as Louisell

$$\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}^{}l^{}}=\frac{\tau _{𝐤_{}l}^1+\tau _{𝐤_{}^{}l^{}}^1}{2}+\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}^{}l^{}}^{\mathrm{ph}},$$
(4)
where $`\mathrm{\Gamma }^{\mathrm{ph}}`$ describes additional dephasing.

The induced second-order polarization is given by

$$P_i^{(2)}(t)=\frac{1}{(2\pi )^2}𝑑t_1𝑑t_2\chi _{ijk}^{(2)}(tt_1,t_1t_2)E_j(t_1)E_k(t_2).$$
(5)
The second-order susceptibility $`\chi ^{(2)}`$ depends only on two time differences due to homogeneity in time. Obviously, $`|t_1t_2|`$ is the time interval between the two absorptions. For a single laser pulse this interval is controlled by the pulse width. For two pulses we expect a contribution for $`|t_1t_2|`$ of the order of the pump-probe delay time $`\mathrm{\Delta }T`$. Note, the light polarization is characterized by the components $`E_j`$.

To express the electric field $`𝐄`$ within the solid in terms of the applied external light field $`𝐄_{\mathrm{las}}`$ and similarly the electric field $`𝐄_{\mathrm{out}}`$ of the outgoing light in terms of the polarization $`𝐏`$ one should employ Fresnel formulae, which are also of importance for the coupling of the light to collective excitations, as we discuss below. We do not present the Fresnel formulae here, since they can be found in the literature SMD87 ; HBB94 . See also Refs. Born75 ; Board  for effective Fresnel factors for systems of several layers, such as the important case of a coupling prism separated from the metal by a thin layer of air or vacuum Otto . Of course, it would be of interest to repeat Fresnel’s analysis for SFG, in particular deducing phase shifts etc.

Equation (5) is the basis for time-dependent SFG. Clearly the pulse shape of the applied light described by $`𝐄(t)`$ affects the induced polarization $`𝐏^{(2)}(t)`$. Note, for simple pulse shapes (Gaussian, Lorentzian, rectangular) it is possible to evaluate the integrals in Eq. (5) further. The light polarization dependence is controlled by the symmetries of the tensor $`\chi _{ijk}^{(2)}`$. The symmetries of $`\chi ^{(2)}`$ for magnetic and nonmagnetic crystals under monochromatic light have been discussed in Ref. Pan . They are determined by the symmetry operations that leave the particular surface invariant. These symmetry arguments are unchanged for general time dependence of the applied laser field.

The intensity of SFG light is $`I^{(2)}(t)[E_{\mathrm{out}}^{(2)}(t)]^2[P^{(2)}(t)]^2`$. So far, typical experiments do not resolve the time dependence of the intensity, but measure the time-integrated SFG yield

$$^{(2)}𝑑tI^{(2)}(t)𝑑t\left[E_{\mathrm{out}}^{(2)}(t)\right]^2𝑑t\left[P^{(2)}(t)\right]^2.$$
(6)
For simplicity we here sum over polarization directions.

Time-resolved SFG may be performed by measuring $`^{(2)}(\mathrm{\Delta }T)`$ as a function of the time delay $`\mathrm{\Delta }T`$ between the applied field pulses. Omitting surface effects (Fresnel factors) to emphasize the structure of the results, the yield can be written as

$`^{(2)}(\mathrm{\Delta }T){\displaystyle 𝑑t𝑑t_1𝑑t_2𝑑t_3𝑑t_4\chi _{ijk}^{(2)}(tt_1,t_1t_2)}`$ (7)
$`\times \chi _{ilm}^{(2)}(tt_3,t_3t_4)E_j(t_1)E_k(t_2)E_l(t_3)E_m(t_4),`$
which is of fourth order in the incoming light field and thus of second order in its intensity. As mentioned above, the typical time differences dominating the response are controlled by the delay $`\mathrm{\Delta }T`$, besides the pulse durations.

It is useful to write the second-order polarization $`𝐏^{(2)}`$ also in frequency space,

$$P_i^{(2)}(\omega )=𝑑\omega ^{}\chi _{ijk}^{(2)}(\omega ,\omega ^{})E_j(\omega ^{})E_k(\omega \omega ^{}),$$
(8)
where $`\chi _{ijk}^{(2)}(tt_1,t_1t_2)=𝑑\omega 𝑑\omega ^{}e^{i\omega (tt_2)}e^{i\omega ^{}(t_2t_1)}\chi _{ijk}^{(2)}(\omega ,\omega ^{})`$ or

$$\chi _{ijk}^{(2)}(\omega ,\omega ^{})=\frac{1}{4\pi ^2}𝑑t𝑑t^{}e^{i\omega t}e^{i\omega ^{}t^{}}\chi _{ijk}^{(2)}(t+t^{},t^{}).$$
(9)
Note that we employ the convention of Eq. (8) in Ref. HB  for the Fourier transformation. The frequency representation is better suited to discuss transition *energies*. $`𝐏^{(2)}`$ has components at the sum of two frequencies of the incoming light. Since the Fourier transform of the *real* electric field contains positive and negative frequencies, the difference frequency also appears.

If at the first step we ignore screening effects, then Eqs. (5) and (7) only contain the second-order *irreducible* susceptibility

$`\chi _{\mathrm{irr};ijk}^{(2)}(𝐪,𝐪_1,𝐪_2;tt_1,t_1t_2)={\displaystyle \frac{e^3}{v}}\left({\displaystyle \frac{2\pi i}{\mathrm{}}}\right)^2\mathrm{\Theta }(tt_1)\mathrm{\Theta }(t_1t_2){\displaystyle \underset{𝐤_{}}{}}{\displaystyle \underset{ll^{}l^{\prime \prime }}{}}D_{𝐤_{}l;𝐤_{}l^{\prime \prime }}^i(q_z)(D_{𝐤_{}l^{\prime \prime };𝐤_{}l^{}}^j(q_{1z})`$ (10)
$`\times D_{𝐤_{}l^{};𝐤_{}l}^k(q_{2z})\left[f(E_{𝐤_{}l})f(E_{𝐤_{}l^{}})\right]\mathrm{exp}\left[i{\displaystyle \frac{E_{𝐤_{}l}E_{𝐤_{}l^{}}}{\mathrm{}}}(t_1t_2)\right]\mathrm{exp}[\mathrm{\Gamma }_{𝐤_{}l^{};𝐤_{}l}(t_1t_2)]`$
$`D_{𝐤_{}l^{};𝐤_{}l}^j(q_{1z})D_{𝐤_{}l^{\prime \prime };𝐤_{}l^{}}^k(q_{2z})\left[f(E_{𝐤_{}l^{}})f(E_{𝐤_{}l^{\prime \prime }})\right]\mathrm{exp}\left[i{\displaystyle \frac{E_{𝐤_{}l^{}}E_{𝐤_{}l^{\prime \prime }}}{\mathrm{}}}(t_1t_2)\right]`$
$`\times \mathrm{exp}[\mathrm{\Gamma }_{𝐤_{}l^{\prime \prime };𝐤_{}l^{}}(t_1t_2)]\left)\mathrm{exp}\right[i{\displaystyle \frac{E_{𝐤_{}l}E_{𝐤_{}l^{\prime \prime }}}{\mathrm{}}}(tt_1)]\mathrm{exp}[\mathrm{\Gamma }_{𝐤_{}l^{\prime \prime };𝐤_{}l}(tt_1)],`$
which is derived in App. A. We neglect the photon momenta relative to the crystal momentum. This expression, which forms the basis of our discussion of SFG, goes beyond the one given in Ref. HB  in that it is valid for a time-dependent and spatially varying laser field. Furthermore, it includes the transverse response explicitly. Equation (10) already exhibits the interplay between the time interval $`|t_1t_2|`$ between absorptions, the photon frequencies, the dephasing times, and the transition frequencies.

We now discuss the physics contained in Eq. (10) with the help of Fig. 4. We consider the first of the two terms in Eq. (10). The interpretation of the second term is similar rem.term2 . The step functions incorporate the time ordering $`t_2<t_1<t`$ and thus guarantee causality. The system is in equilibrium until the first absorption at time $`t_2`$ creates a superposition of the two states $`|𝐤_{}l`$ and $`|𝐤_{}l^{}`$, denoted by the wavy line in Fig. 4. This important physics is lost in the interpretation illustrated by Fig. 1. The Fermi functions make sure that one of the states is initially occupied and the other is empty. Let us say state $`|𝐤_{}l`$ is occupied. Since the system is in a superposition of two eigenstates, the polarization oscillates with the frequency $`(E_{𝐤_{}l}E_{𝐤_{}l^{}})/\mathrm{}`$, as follows from the first exponential in the parentheses in Eq. (10). Such superpositions are described by the off-diagonal components of the density matrix rem.linear . The diagonal components denoting the occupation numbers of states are *not* changed by a single absorption. The superposition decays with the dephasing rate $`\mathrm{\Gamma }_{𝐤_{}l^{};𝐤_{}l}`$ associated with this transition, making it clear why the dephasing rates rather than the energy relaxation rates dominate the response. A second absorption at the later time rem.simul  $`t_1`$ changes the state into a superposition of the originally occupied state and the state $`|𝐤_{}l^{\prime \prime }`$ with its own characteristic oscillation frequency $`(E_{𝐤_{}l}E_{𝐤_{}l^{\prime \prime }})/\mathrm{}`$ and dephasing rate. This oscillating polarization can emit a photon at that frequency. After the emission the electron is again in the pure eigenstate $`|𝐤_{}l`$. The nonlinear susceptibility in Eq. (10) contains a sum over many contributions of this type from different momenta and bands Weida . Note, the product of *three* dipole matrix elements appearing in $`\chi _{\mathrm{irr}}^{(2)}`$ is responsible for the surface sensitivity of SFG, since in inversion symmetric crystals the product of dipole matrix elements connecting three states vanishes except when inversion symmetry is explicitly broken, e.g., by the surface.

In frequency space the nonlinear susceptibility is given by

$`\chi _{\mathrm{irr};ijk}^{(2)}(𝐪,𝐪_1,𝐪_2;\omega ,\omega ^{})={\displaystyle \frac{e^3}{v}}{\displaystyle \underset{𝐤_{}}{}}{\displaystyle \underset{ll^{}l^{\prime \prime }}{}}`$ (11)
$`\times {\displaystyle \frac{D_{𝐤_{}l;𝐤_{}l^{\prime \prime }}^i(q_z)}{\mathrm{}\omega +E_{𝐤_{}l}E_{𝐤_{}l^{\prime \prime }}i\mathrm{}\mathrm{\Gamma }_{𝐤_{}l^{\prime \prime };𝐤_{}l}}}`$
$`\times [D_{𝐤_{}l^{\prime \prime };𝐤_{}l^{}}^j(q_{1z})D_{𝐤_{}l^{};𝐤_{}}^k(q_{2z})`$
$`\times {\displaystyle \frac{f(E_{𝐤_{}l})f(E_{𝐤_{}l^{}})}{\mathrm{}\omega +\mathrm{}\omega ^{}+E_{𝐤_{}l}E_{𝐤_{}l^{}}i\mathrm{}\mathrm{\Gamma }_{𝐤_{}l^{};𝐤_{}l}}}`$
$`D_{𝐤_{}l^{};𝐤_{}l}^j(q_{1z})D_{𝐤_{}l^{\prime \prime };𝐤_{}l^{}}^k(q_{2z})`$
$`\times {\displaystyle \frac{f(E_{𝐤_{}l^{}})f(E_{𝐤_{}l^{\prime \prime }})}{\mathrm{}\omega +\mathrm{}\omega ^{}+E_{𝐤_{}l^{}}E_{𝐤_{}l^{\prime \prime }}i\mathrm{}\mathrm{\Gamma }_{𝐤_{}l^{\prime \prime };𝐤_{}l^{}}}}].`$
This shows that the contribution of intermediate (virtual) states falls off with the inverse of the initial-state energy plus the photon energy minus the intermediate state energy, i.e., with the inverse of the *detuning*. This is not related to the lifetime broadening, but is due to Heisenberg’s uncertainty principle, which allows energy conservation to be violated on short time scales. The frequency picture also allows to incorporate a weight factor to account for the frequency resolution of the detector Weida . It is of interest to note that Eq. (11) and Fig. 4 can also describe spin-selective electron excitations due to circularly polarized light. Including the electron spins our response theory and in particular Eq. (11) apply also to magnetic systems.

To prepare the analysis of the effect of collective plasma excitations on the nonlinear optical response we now include the screening of the electric fields. Screening enters in two ways: First, the effective field $`𝐄`$ within the solid is not identical to the external field because of linear screening, which is expressed by the Fresnel formulae SMD87 ; HBB94 ; Luce97 ; LuceLDA ; AtK02  containing the dielectric function $`\epsilon `$, which can be determined in the RPA. Secondly, the second-order polarization $`𝐏^{(2)}`$ of the electron gas, which corresponds to a displacement of charge, leads to an additional electric field Jackson

$$E_i^{(2)}(𝐫)=d^3r^{}\underset{j}{}\left[\frac{3(r_ir_i^{})(r_jr_j^{})}{|𝐫𝐫^{}|^5}\frac{\delta _{ij}}{|𝐫𝐫^{}|^3}\frac{4\pi }{3}\delta _{ij}\delta (𝐫𝐫^{})\right]P_j^{(2)}(𝐫^{}).$$
(12)
Fourier transformation leads to

$$𝐄^{(2)}(𝐤,t)=4\pi \widehat{𝐤}\widehat{𝐤}𝐏^{(2)}(𝐤,t),$$
(13)
where $`\widehat{𝐤}`$ is the unit vector in the direction of $`𝐤`$. Thus only the component of $`𝐏^{(2)}(𝐤,t)`$ parallel to $`𝐤`$, i.e., its longitudinal part, is accompanied by an electric field $`𝐄^{(2)}`$ Jackson , which is also longitudinal. Note that a longitudinal component of the electric field and of the induced polarization generally exists even for a transverse applied external field for lattice models Barton ; Board , see Eq. (11).

Due to the *linear* polarizability of the solid the additional field $`𝐄^{(2)}`$ leads to a polarization contribution of the form $`\chi 𝐄^{(2)}`$. Since the field $`𝐄^{(2)}`$ in Eq. (13) is of second order in the applied field, see Eq. (5), this polarization contribution must be taken into account in $`𝐏^{(2)}`$. Doing this selfconsistently corresponds to the summation of an RPA series HB , as shown in App. A. Then, $`P_i^{(2)}𝑑t_1𝑑t_2\chi _{ijk}^{(2)}E_jE_k`$ where now the nonlinear susceptibility $`\chi _{ijk}^{(2)}`$ obtains an additional factor and is given by

$`\chi _{ijk}^{(2)}(𝐪,𝐪_1,𝐪_2;tt_1,t_1t_2)={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{m}{}}{\displaystyle \underset{\stackrel{~}{𝐪}}{}}{\displaystyle 𝑑\stackrel{~}{t}}`$ (14)
$`\times \epsilon _{\mathrm{long};im}^1(𝐪,\stackrel{~}{𝐪};t\stackrel{~}{t})\chi _{\mathrm{irr};mjk}^{(2)}(\stackrel{~}{𝐪},𝐪_1,𝐪_2;\stackrel{~}{t}t_1,t_1t_2).`$
Here the irreducible susceptibility $`\chi _{\mathrm{irr}}^{(2)}`$ is given by Eq. (10). Since only the longitudinal component of $`𝐏^{(2)}`$ is accompanied by an electric field, the screening factor $`\epsilon _{\mathrm{long}}^1`$ appears only for the longitudinal component. This is expressed by the factor $`\widehat{\stackrel{~}{q}}_m\widehat{\stackrel{~}{q}}_j`$ in the explicit expression $`\epsilon _{\mathrm{long};ij}(𝐪,\stackrel{~}{𝐪};t\stackrel{~}{t})\delta _{ij}+4\pi \chi _{im}(𝐪,\stackrel{~}{𝐪};t\stackrel{~}{t})\widehat{\stackrel{~}{q}}_m\widehat{\stackrel{~}{q}}_j`$. $`\epsilon _{\mathrm{long};im}^1(𝐪,\stackrel{~}{𝐪};t\stackrel{~}{t})`$ is the inverse matrix with respect to the indices $`(i,q_z)`$ and $`(m,\stackrel{~}{q}_z)`$.

This analysis is illustrated by Fig. 5. Figure 5(a) shows the diagram of the second-order susceptibility $`\chi ^{(2)}`$. The square vertex represents the additional factor of $`\epsilon _{\mathrm{long}}^1`$. It is obtained from the Dyson equation in Fig. 5(b). The expression for $`\chi _{\mathrm{irr}}^{(2)}`$ in Eq. (10) is called irreducible since its diagram Fig. 5(a) with the square vertex replaced by a normal one cannot be cut into two by severing a single photon line.

The response theory clarifies how collective plasma excitations affect SFG. They essentially enter in two ways, both of which are controlled by the full (not only longitudinal) dielectric function $`\epsilon `$:

First, the effective electric field is expressed in terms of the external field by means of Fresnel formulae SMD87 ; HBB94 ; Luce97 ; LuceLDA ; AtK02 , which contain contributions of order $`1/\epsilon `$ for small $`\epsilon `$. The dielectric function $`\epsilon `$ becomes small if the frequency of the external field is close to the plasma frequency. This contribution can be interpreted as *field enhancement*. In addition, the outgoing (sum-frequency) electric field $`𝐄_{\mathrm{out}}`$ also contains terms that are enhanced for small $`\epsilon `$ due to the Fresnel factors. This enhancement is most pronounced if the sum frequency is close to the plasma frequency.

Secondly, the longitudinal component of the nonlinear polarization $`𝐏^{(2)}`$ of the electron system is accompanied by an electric field $`𝐄^{(2)}`$ given by Eq. (13). Thus, the factor $`\epsilon _{\mathrm{long}}^1`$ appears in the nonlinear susceptibility in Eq. (14) and thus in $`𝐏^{(2)}`$ HB . This leads to an enhancement of the SFG light due to the longitudinal part of $`𝐏^{(2)}`$ if the sum frequency is close to the plasma frequency.

### II.2 Two-photon photoemission

To demonstrate the similarities between SFG and 2PPE, we continue by summarizing the results of the response theory for 2PPE. We consider the same band structure as for SFG, which is characterized by single-electron energies $`E_{𝐤_{}l}`$. We emphasize that this band structure contains the bulk states with the $`z`$-component $`k_z`$ of $`𝐤`$ included in $`l`$.

The response theory starts from the observation that the photoelectron current $`j(t;𝐤,\sigma )`$ of electrons of momentum $`𝐤`$ and spin $`\sigma `$ is given by the change of occupation of the vacuum state $`|𝐤\sigma ,\mathrm{out}`$ outside of the crystal. However, in practice the time-dependence of $`j`$ is not measured, but only the total photoelectron yield $`𝒩(𝐤,\sigma )=𝑑tj(t;𝐤,\sigma )`$. This is similar to SFG, where only the time-integrated intensity is measured. The response theory directly determines the photoelectron yield $`𝒩`$. To prepare the discussion it is useful to first consider ordinary single-photon photoemission.

Single-photon photoemission: The photoelectron yield is given by

$$𝒩(𝐤,\sigma )=\underset{𝐪}{}𝑑t_1𝑑t_2\eta _{ij}(𝐪;t_1,t_2;𝐤,\sigma )E_i(𝐪,t_1)E_j(𝐪,t_2),$$
(15)
with the response function (see App. B)

$`\eta _{ij}(𝐪;t_1,t_2;𝐤,\sigma )={\displaystyle \frac{e^2}{\mathrm{}^2}}{\displaystyle \frac{\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}}{\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}}}{\displaystyle \underset{\lambda }{}}D_{𝐤\sigma ,\mathrm{in};𝐤_{}\lambda }^i(q_z)`$ (16)
$`\times \mathrm{exp}\left[i{\displaystyle \frac{E_{𝐤_{}\lambda }E_{𝐤\sigma ,\mathrm{in}}}{\mathrm{}}}(t_2t_1)\right]e^{\mathrm{\Gamma }_{𝐤_{}\lambda ;𝐤\sigma ,\mathrm{in}}|t_2t_1|}`$
$`\times f(E_{𝐤_{}\lambda })D_{𝐤_{}\lambda ;𝐤\sigma ,\mathrm{in}}^j(q_z).`$
Here, $`|𝐤\sigma ,\mathrm{in}`$ is a state with momentum $`𝐤`$ and spin $`\sigma `$ inside the crystal but above the vacuum energy. We have again neglected the momentum transferred by the photon. The standard diagrammatic representation of ordinary photoemission is shown in Fig. 6 Feder . The effective field $`𝐄`$ within the solid should again be expressed in terms of the external light field with the help of the proper boundary conditions. The response function $`\eta `$ will play a role when we discuss the various contributions to 2PPE.

We briefly commend on the structure of this expression: The prefactor $`\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}/\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}`$ describes the probability that electrons excited above the vacuum energy actually leave the crystal. Photoemission is often described by a three-step picture BS ; FFW ; WD : First, electrons are excited, then they are transported to the surface, and finally they leave the crystal. In this work we are mainly interested in the first step. The second and third steps are incorporated phenomenologically by effective relaxation rates $`\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}`$, which describe electrons dropping below $`E_{\mathrm{vac}}`$ before they reach the surface, and effective transition rates $`\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}`$ from states above $`E_{\mathrm{vac}}`$ within the solid to free electron states outside of the solid. Note, the yield is proportional to the electric field squared and thus to the intensity of the incoming light.

Two-photon photoemission: The total 2PPE yield consists of the three contributions

$$𝒩^{2\mathrm{P}\mathrm{P}\mathrm{E}}=𝒩_{\mathrm{irr}}^{2\mathrm{P}\mathrm{P}\mathrm{E}}+𝒩_{\mathrm{red},1}^{2\mathrm{P}\mathrm{P}\mathrm{E}}+𝒩_{\mathrm{red},2}^{2\mathrm{P}\mathrm{P}\mathrm{E}}$$
(17)
corresponding to Fig. 7(a), (b), and (c), respectively. The second and third term arise from the nonlinear optical properties of the solid: Close to the surface the effective field $`𝐄`$ leads to a second-order polarization $`𝐏^{(2)}`$, see Eq. (5), which is accompanied by an electric field $`𝐄^{(2)}`$. This field may contribute to photoemission, leading to the processes in Fig. 7(b) and (c). The second-order field $`𝐄^{(2)}`$ is also responsible for SFG accompanying 2PPE. However, this SFG is usually a small effect since the SFG light *intensity* is of *sixth* order in dipole matrix elements $`𝐃`$, see Eqs. (7) and (10), whereas the 2PPE current is of *fourth* order, as is shown below in Eq. (24). This changes if the sum frequency is close to the plasma frequency, in which case SFG is enhanced as discussed at the end of Sec. II.1.

Since the diagram in Fig. 7(a) cannot be cut into two by severing a single photon line, the first term $`𝒩_{\mathrm{irr}}^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ is irreducible, while the other two are reducible. The irreducible contribution in Eq. (17) can be written as

$$𝒩_{\mathrm{irr}}^{2\mathrm{P}\mathrm{P}\mathrm{E}}=𝑑t_1𝑑t_2𝑑t_3𝑑t_4\eta _{ijkl}^{2\mathrm{P}\mathrm{P}\mathrm{E}}(t_1,t_2,t_3,t_4)E_i(t_1)E_j(t_2)E_k(t_3)E_l(t_4),$$
(18)
which is of *fourth* order in the electric field and of *second* order in the incoming intensity. This is already clear from the simple picture in Fig. 2: The occupation of the intermediate state $`|2`$ is proportional to the light intensity. To reach state $`|3`$ above $`E_{\mathrm{vac}}`$ another absorption is required, leading to a total proportionality to the intensity squared.

Obviously, the structure of Eq. (18) is very similar to Eq. (7) for the SFG yield:

$`^{(2)}`$ $``$ $`{\displaystyle 𝑑t𝑑t_1𝑑t_2𝑑t_3𝑑t_4\chi _{ijk}^{(2)}(tt_1,t_1t_2)\chi _{ilm}^{(2)}(tt_3,t_3t_4)}`$ (19)
$`\times E_j(t_1)E_k(t_2)E_l(t_3)E_m(t_4).`$
Hence, we expect similar interference effects in both cases.

The other two contributions to $`𝒩^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ are

$`𝒩_{\mathrm{red},1}^{2\mathrm{P}\mathrm{P}\mathrm{E}}=4\pi {\displaystyle }dt_1dt_2dt_3\eta _{ijk}^{(3)}[P_i^{(2)}(t_1)E_j(t_2)E_k(t_3)`$ (20)
$`+E_i(t_1)P_j^{(2)}(t_2)E_k(t_3)+E_i(t_1)E_j(t_2)P_k^{(2)}(t_3)],`$
and

$$𝒩_{\mathrm{red},2}^{2\mathrm{P}\mathrm{P}\mathrm{E}}=(4\pi )^2𝑑t_1𝑑t_2\eta _{ij}P_i^{(2)}(t_1)P_j^{(2)}(t_2),$$
(21)
where $`\eta `$ is given in Eq. (16). Since the nonlinear susceptibility $`\chi ^{(2)}`$ and hence $`𝐏^{(2)}`$ contains three dipole matrix elements, the reducible contributions to the photoelectron current are of higher order in dipole matrix elements and are thus usually small. However, the longitudinal component of $`𝐏^{(2)}`$ contains a factor $`\epsilon _{\mathrm{long}}^1`$. If the nonlinear polarization is enhanced due to a bulk plasma resonance at the sum frequency, one expects significant contributions from the reducible terms. The response functions $`\eta ^{(3)}`$ and $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ are given in App. B.

We next consider the response functions $`\eta `$, $`\eta ^{(3)}`$, and $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ which determine the yield $`𝒩^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$. The functions $`\eta ^{(3)}`$ and $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ appearing in Eqs. (20) and (21), respectively, are of the same general form as $`\eta `$, Eq. (16), but have more terms resulting from different orders of the time arguments. We now first present the structure of the response expression for the main, irreducible contribution to the 2PPE yield, Eq. (18), and then discuss its physical interpretation. Fully written out, Eq. (18) reads

$`𝒩_{\mathrm{irr}}^{2\mathrm{P}\mathrm{P}\mathrm{E}}(𝐤,\sigma )`$ $`=`$ $`{\displaystyle \underset{𝐪_1𝐪_2𝐪_3}{}}{\displaystyle 𝑑t_1𝑑t_2𝑑t_3𝑑t_4\eta _{ijkl}^{2\mathrm{P}\mathrm{P}\mathrm{E}}(𝐪_1,𝐪_2,𝐪_3;t_1,t_2,t_3,t_4;𝐤,\sigma )E_i(𝐪_1,t_1)E_j(𝐪_2,t_2)}`$ (22)
$`\times E_k(𝐪_3,t_3)E_l(𝐪𝐪_1𝐪_2,t_4).`$
Defining the complex transition energy

$$\mathrm{\Omega }_{𝐤_{}l;𝐤_{}^{}l^{}}\frac{E_{𝐤_{}l}E_{𝐤_{}^{}l^{}}}{\mathrm{}}i\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}^{}l^{}},$$
(23)
we obtain the response function

$`\eta _{ijkl}^{2\mathrm{P}\mathrm{P}\mathrm{E}}(𝐪_1,𝐪_2,𝐪_3;t_1,t_2,t_3,t_4;𝐤,\sigma )={\displaystyle \frac{e^4}{\mathrm{}^4}}{\displaystyle \frac{\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}}{\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}}}{\displaystyle \underset{\lambda _1\lambda _2\lambda _3}{}}D_{𝐤\sigma ,\mathrm{in};𝐤_{}\lambda _1}^i(q_{1z})D_{𝐤_{}\lambda _1;𝐤_{}\lambda _2}^j(q_{2z})`$ (24)
$`\times D_{𝐤_{}\lambda _2;𝐤_{}\lambda _3}^k(q_{3z})D_{𝐤_{}\lambda _3;𝐤\sigma ,\mathrm{in}}^l(q_{1z}q_{2z}q_{3z})`$
$`\times \{\mathrm{\Theta }(t_1t_2)\mathrm{\Theta }(t_2t_3)\mathrm{\Theta }(t_3t_4)e^{i\mathrm{\Omega }_{𝐤_{}\lambda _1;𝐤\sigma ,\mathrm{in}}(t_1t_2)}e^{i\mathrm{\Omega }_{𝐤_{}\lambda _2;𝐤\sigma ,\mathrm{in}}(t_2t_3)}e^{i\mathrm{\Omega }_{𝐤_{}\lambda _3;𝐤\sigma ,\mathrm{in}}(t_3t_4)}[f(E_{𝐤_{}\lambda _3})]`$
$`\mathrm{\Theta }(t_1t_2)\mathrm{\Theta }(t_2t_4)\mathrm{\Theta }(t_4t_3)e^{i\mathrm{\Omega }_{𝐤_{}\lambda _1;𝐤\sigma ,\mathrm{in}}(t_1t_2)}e^{i\mathrm{\Omega }_{𝐤_{}\lambda _2;𝐤\sigma ,\mathrm{in}}(t_2t_4)}e^{i\mathrm{\Omega }_{𝐤_{}\lambda _2;𝐤_{}\lambda _3}(t_4t_3)}\left[f(E_{𝐤_{}\lambda _3})f(E_{𝐤_{}\lambda _2})\right]`$
$`\mathrm{}\}.`$
There are eight terms in the curly braces, which correspond to different temporal orders of interactions with the electric field. Note, the dependence of 2PPE on light polarization is incorporated in the symmetries of the tensor $`\eta _{ijkl}^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$, which depend on the dipole matrix elements $`𝐃`$. Unlike for the nonlinear optical response, these symmetries have not been discussed so far. It would be very interesting to determine the symmetries for surfaces of nonmagnetic and magnetic solids.

Equation (24) forms the basis for our discussion of 2PPE. To clarify the time dependence exhibited in Eq. (24) we discuss the second term, the others are in principle similar but correspond to different orders of the times $`t_i`$. The processes described by this term are illustrated in Fig. 8. The system starts in equilibrium from the state $`|𝐤_{}\lambda _2`$. The first interaction with the electric field takes place at time $`t_3`$ and creates a superposition of the states $`|𝐤_{}\lambda _2`$ and $`|𝐤_{}\lambda _3`$, leading to oscillations at the frequency $`(E_{𝐤_{}\lambda _3}E_{𝐤_{}\lambda _2})/\mathrm{}`$ expressed by the third exponential factor in this term. The Fermi factors ensure that one of the states is initially occupied and that the other one is empty. Let us assume that $`|𝐤_{}\lambda _2`$ is occupied. The second interaction at $`t_4`$ changes the state into a superposition of $`|𝐤_{}\lambda _2`$ and the vacuum state $`|𝐤\sigma ,\mathrm{in}`$, leading to oscillations at the corresponding difference frequency (second exponential factor), and the third interaction at $`t_2`$ creates a superposition of the vacuum state and $`|𝐤_{}\lambda _1`$. After the fourth interaction the electron is in a pure state above $`E_{\mathrm{vac}}`$ and can leave the solid with finite probability. Of course, due to the sum over bands there are usually several contributions of this type. Only if the superpositions decay very rapidly compared to the pure states, a description in terms of rate equations, as suggested by Fig. 2, is applicable Kno00 ; Rol1 ; Rol2 . Also compare the discussion of SFG above, see Fig. 4.

While SFG is only governed by the dephasing rates but not the energy relaxation rates, 2PPE depends on both. This is because in the 2PPE response function $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ the change of occupation of states enters besides the polarization of the electron gas, whereas SFG only depends on the latter.

Note, the 2PPE yield contains four dipole matrix elements. Thus, even for inversion-symmetric crystals parity does not forbid 2PPE from the bulk. However, 2PPE is sensitive to a surface region of a thickness given by the mean free path of electrons above $`E_{\mathrm{vac}}`$. The optical penetration depth is typically significantly larger than the mean free path and thus does not enter here. Equations (16) and (24) also illustrate that 2PPE is sensitive to specific points in the Brillouin zone: The photoelectron momentum $`𝐤`$ measured by momentum-resolved 2PPE is approximately the same as the lattice momentum of the original unperturbed electron and also of the intermediate state due to the small photon momentum. These effects obviously require a theoretical description that considers the $`𝐤`$-dependent states in the solid, like our approach does as opposed to both the random-$`𝐤`$ approximation and Bloch equations. In view of the importance of angle-resolved (ordinary) photoemission spectroscopy (ARPES) for, e.g., cuprate high-$`T_c`$ superconductors, $`𝐤`$-resolved 2PPE is expected to yield interesting results in the future. On the other hand, if one only measures the total number of photoelectrons, the $`𝐤`$-space resolution is lost and 2PPE and SFG give very similar information. It is obvious that 2PPE has the disadvantage of being limited to frequencies $`\omega _1`$, $`\omega _2`$ such that $`E_F+\mathrm{}\omega _1+\mathrm{}\omega _2`$ lies above the vacuum energy, unlike SFG. As stated already the SFG accompanying 2PPE is usually small since additional dipole matrix elements are involved.

The response expressions show that collective plasma excitations affect 2PPE in two ways: First, exactly like for SFG the effective field $`𝐄`$ within the metal differs from the external field due to linear screening and is enhanced close to the plasmon resonance. Secondly, the reducible contributions in Eqs. (20) and (21) depend on the second-order polarization $`𝐏^{(2)}`$, which contains a factor of $`\epsilon _{\mathrm{long}}^1`$, see Eq. (14). $`𝐏^{(2)}`$ is enhanced if the sum frequency is close to the plasma frequency. In 2PPE this enhancement enters only in the reducible contributions in Eqs. (20) and (21).

## III Discussion

The aim of the present section is to discuss and illustrate the results of the response theory for time-resolved SFG and 2PPE. In particular, we consider time-dependent effects on the femtosecond time scale. Our results exhibit the intimate relation between SFG and 2PPE. We can already gain insight by studying the general structure of the response expressions for SFG and 2PPE, for example Eqs. (7) and (18), respectively, independently of the specific approximations made here. For clarity we apply our response theory to a simple model system.

### III.1 Time-dependent effects in SFG and 2PPE

The response expressions of the preceding section are valid for any time dependence of the exciting laser field. The time enters the response expressions for both SFG and 2PPE in two ways, apart from the step functions from causality, cf. Eqs. (10), (16), and (24): The difference between the time arguments of electric fields appears in exponentials oscillating at the transition frequency of the involved electron states and in exponentials decaying with the dephasing rate of the superposition of the two states, and, for 2PPE, also exponentials decaying with the energy relaxation rate of an intermediate state. (See the discussion of Figs. 4 and 8 for the interpretation of SFG and 2PPE in terms of electronic excitations.) The time passing between absorptions can be controlled by the pulse shape of the exciting laser pulses: If the total duration $`T`$ of a pulse of arbitrary shape is much larger than typical relaxation times $`\tau `$ then the yield depends on the probability to absorb two photons within a time interval $`\tau `$, which is independent of $`T`$. On the other hand, for $`T\tau `$ there is almost no relaxation during the pulse. Thus the response theory reproduces the well-known result that $`\tau `$ can only be inferred from SFG or 2PPE experiments if the total pulse duration is $`T\tau `$.

To be more specific, in most experiments two approximately Gaussian pulses are used (pump-probe method) Ste92 ; HMKB97 ; Sim98 ; KSR98 ; Gud99 ; Schoe ; Aesch1 ; Aesch2 ; Wolf1 ; Hoef97 ; Cao97 ; Petek ; Wolf2 ; Leh99 ; Kno00 ; GWH00 ; SBW02 ; ONP02 ; Gerb ; Ae.cluster . If the two pulses are of different mean frequencies $`\omega _1`$ and $`\omega _2`$ (two-color case) and one measures the SFG or 2PPE response at the sum frequency $`\omega _1+\omega _2`$ then it is obvious which photon was absorbed out of which pulse. Then for long time delay $`\mathrm{\Delta }T`$ compared to the single-pulse duration the relaxation rate of intermediate states can be read of directly from the $`\mathrm{\Delta }T`$ dependence of the total yield. In pump-probe experiments with two pulses of the same mean frequency $`\omega `$ (single-color case), photons can be absorbed out of the same or different pulses. However, the contribution with all absorptions out of the same pulse obviously does not depend on $`\mathrm{\Delta }T`$, just leading to a constant background. Note, in all these cases only a typical relaxation time enters, which usually is a weighted average over relaxation times of many states Weida . If only a single relevant intermediate state is present, e.g., for a quantum-well state, or if there are many but of similar relaxation rate, the relaxation time extracted from experiment will be the actual dephasing time of intermediate states. However, if intermediate states with very different dynamical properties are involved, for example if both sp and d bands are relevant, the measured relaxation time does not describe any single excited electron state.

In pump-probe SHG Ste92 ; HMKB97 ; Sim98 ; KSR98 ; Gud99  or pump-probe single-color 2PPE Aesch2 ; Cao97 ; Petek ; Wolf2 ; Kno00  experiments, time-dependent interference effects are especially pronounced. Their origin is the following: The first absorption of a photon of frequency $`\omega `$ sets up an oscillating polarization of the excited electrons. Now the probability of a second absorption depends on the relative phase of the oscillating polarization and the second photon. Since the oscillating polarization is described by the off-diagonal components of the density matrix $`\rho `$, a description in terms of rate equations, which omits these components, is unable to describe interference.

For further illustration of this interference, we show results for SHG and 2PPE for a simple model. Unless stated otherwise, this model consists of three bands. The lowest one is a three-dimensional tight-binding band $`1`$ with band center at rem.param  $`3.33\mathrm{eV}`$ (all energies are measured relative to the Fermi energy) and half width $`3.81\mathrm{eV}`$. The band maximum is at $`𝐤=0`$. The second, rather flat tight-binding band $`2`$ is centered at $`2.29\mathrm{eV}`$ with half width $`0.48\mathrm{eV}`$ and maximum also at $`𝐤=0`$. Finally, there is a free electron band $`3`$ representing electrons above the vacuum energy $`E_{\mathrm{vac}}=4.29\mathrm{eV}`$. There exist points in the Brillouin zone for which the energy differences between bands 2 and 1 as well as between bands 3 and 2 both equal the photon energy of $`\mathrm{}\omega =3.05\mathrm{eV}`$. We assume that the relaxation rates $`\mathrm{\Gamma }_{n_1n_2}`$ only depend on the band indices $`n_1`$, $`n_2`$ but not on the $`𝐤`$ vector (see below). We use the energy relaxation rates $`\mathrm{}\mathrm{\Gamma }_{22}=0.191\mathrm{eV}`$ (corresponding to the lifetime $`\tau _2=3.5\mathrm{fs}`$) and $`\mathrm{}\mathrm{\Gamma }_{33}=0.381\mathrm{eV}`$ ($`\tau _3=1.7\mathrm{fs}`$) and no additional dephasing, i.e., $`\mathrm{\Gamma }_{n_1,n_2}^{\mathrm{ph}}=0`$ in Eq. (4). These short lifetimes are assumed to bring out the time-dependent effects more clearly. The dipole matrix elements are treated as constants.

In the following we use this model to show how time-dependent effects emerge from our response theory. For clarity we neglect the Fresnel formulae, which do not change the results qualitatively. We demonstrate that our theory gives reasonable results for a moderately complicated system. Obviously, it can be applied to a more realistic band structure at the expense of computation time. The boundary conditions (Fresnel factors) are also omitted for simplicity.

In Fig. 9(a) we show the 2PPE photoelectron yield for a particular momentum $`𝐤`$ as a function of the delay time $`\mathrm{\Delta }T`$ between two identical Gaussian pump and probe pulses. A mean photon energy of $`\mathrm{}\omega =3.05\mathrm{eV}`$ is assumed, corresponding to a wave length of about $`\lambda =400\mathrm{nm}`$, and the duration of each pulse is $`10.3\mathrm{fs}`$ (full width at half maximum of the Gaussian envelope of the electric field). The vector $`𝐤`$ is chosen so that the transition energies between the bands match $`\mathrm{}\omega `$. In Fig. 9(b) we show the total SHG photon yield for exactly the same system. Unlike 2PPE, SHG integrates over the whole Brillouin zone. Nevertheless, the overall similarity of Figs. 9(a) and (b) demonstrates the similarity of the response expressions for SHG and 2PPE, compare Eqs. (7) and (18), for example. It means that similar information, e.g., about the relaxation rates, can be obtained from both. The SHG curve is quite similar to the case of flat bands, shown in the inset in Fig. 9(b). This means that only a small region of $`𝐤`$ space contributes. The resulting interference between different $`𝐤`$ points becomes apparent in the tail of the interference pattern, where the main plot in Fig. 9(b) is more irregular and decays faster. This is the averaging effect discussed in Ref. Weida . More precisely it is an *interference* effect between different oscillation frequencies of superpositions of different states.

The 2PPE and SHG interference patterns in Fig. 9 show the well-known $`8:1`$ enhancement of the signal for $`\mathrm{\Delta }T=0`$. This enhancement is due to the yield being of fourth order in the field: For a single pulse the signal would be proportional to $`E^4`$, for two isolated pulses this becomes $`2E^4`$, but for two overlapping pulses the amplitude is doubled, leading to $`(2E)^4=16E^4`$.

For both SHG and 2PPE, the central part of the interference pattern, which corresponds to short delay times $`\mathrm{\Delta }T`$ up to about the single-pulse duration $`T`$, is dominated by the four-field autocorrelation function

$$𝒜_{ijkl}^{(4)}(\mathrm{\Delta }T)𝑑\omega _1𝑑\omega _2𝑑\omega _3E_i(\omega _1)E_j(\omega _2)E_k(\omega _3)E_l(\omega _1\omega _2\omega _3).$$
(25)
This central part stems from the overlap of the two pulses and would be present even for very fast relaxation: Then the response functions $`\chi ^{(2)}`$ and $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ are very sharply peaked in time and thus nearly constant in frequency space, leading to $`^{(2)}𝒜^{(4)}`$ and $`𝒩_{\mathrm{irr}}^{2\mathrm{P}\mathrm{P}\mathrm{E}}𝒜^{(4)}`$ for the SHG and 2PPE yield, respectively, see Eqs. (7) and (18). The autocorrelation signal alone is shown in Fig. 10(b).

In Sec. II we have discussed the response expressions for time-dependent SFG and 2PPE, (10) and (24), respectively. The first interaction creates an oscillating polarization. There is interference if the phase information is still preserved when the second photon is absorbed. This is governed by the dephasing time $`\mathrm{\Gamma }_{21}`$. Thus the interference effects should decay with the time constant $`\mathrm{\Gamma }_{21}^1`$ for large delays $`\mathrm{\Delta }T`$. This is shown in Fig. 10(a) for moderately fast ($`\tau _2=6.9\mathrm{fs}`$) and extremely fast ($`\tau _2=0.86\mathrm{fs}`$) relaxation. For the slower relaxation the tail indeed decays with $`\mathrm{\Gamma }_{21}^1`$ but to observe this one obviously has to look at rather large $`\mathrm{\Delta }T`$ where the interference is already weak. For fast relaxation the curve is nearly indistinguishable from the autocorrelation in Fig. 10(b).

However, there is another crucial origin of the decay of interference: Intermediate states with energies that do not exactly match the energy of the original state plus the photon energy lead to beats at the frequency of the detuning. This effect can be seen from the response expressions. We now discuss this for the case of SFG: For pulses of short duration $`T`$, the times $`t_1`$ and $`t_2`$ in Eq. (10) can be approximated by the pulse centers if we are interested in phenomena at frequencies small compared to $`T^1`$. Then Eq. (10) shows that the polarization of the electron system shortly before the second interaction at $`t_1`$ is proportional to

$$\mathrm{exp}\left[i\frac{E_{𝐤_{}l}E_{𝐤_{}l^{}}}{\mathrm{}}(t_1t_2)\right]e^{\mathrm{\Gamma }_{𝐤_{}l^{};𝐤_{}l}(t_1t_2)}e^{i\omega t_2},$$
(26)
omitting the sum over states. The last factor stems from the electric field of frequency $`\omega `$ describing the first interaction at time $`t_2<t_1`$. The decaying exponential obviously describes the decay of interference with the dephasing rate $`\mathrm{\Gamma }_{𝐤_{}l^{};𝐤_{}l}`$. The oscillating terms are of the form $`\mathrm{exp}[i\delta \omega (t_1t_2)]\mathrm{exp}(i\omega t_1)`$ with $`\delta \omega =(E_{𝐤_{}l^{}}E_{𝐤_{}l})/\mathrm{}\omega `$. Thus we expect slow beats with the detuning frequency $`\delta \omega `$, which lead to an initial decay of the signal on a time scale of $`(\delta \omega )^1`$. There should be a recurring signal at large delay times, but this is in practice suppressed by relaxation. A similar argument can be made for 2PPE using Eq. (24). The effect is clearly seen in Fig. 11 for 2PPE: The width of the pattern is reduced by the detuning. Its tails also become more irregular.

Next, we turn to the effect of the band structure. We first discuss SHG. The SHG yield is determined by a sum over many transitions of different energies and dephasing rates, see Eq. (10). If the decay is governed by dephasing one observes the smallest dephasing rate at large time delays $`\mathrm{\Delta }T`$. However, at intermediate $`\mathrm{\Delta }T`$ one sees an averaged rate. The dephasing rate $`\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}l^{}}`$ for two bands $`l`$ and $`l^{}`$ should usually not change dramatically with $`𝐤`$. On the other hand, the contribution of detuning is necessarily different for transitions with different transition frequencies. Thus, in the interference pattern a continuum of beating frequencies appears. Consequently, the initial decay is governed by an average detuning and later, probably unobservable, recurring signals are strongly reduced by destructive interference of different beating frequencies. The averages are weighted by a factor approximately inversely proportional to the detuning, as seen from Eq. (11). For narrow valence and intermediate state bands the average is restricted to a small effective band width $`W`$. For broader bands but constant relaxation rates throughout each band the dependence of numerical results (not shown) for the SHG yield on the width of the intermediate band turns out to be weak, since in this case all contributing processes are governed by the same relaxation rates rem.SHGband . Hence, if only a single intermediate state or a few states contribute significantly Luce96  or if there is a narrow band with uniform relaxation rates, the Bloch equations should work well. In this case our expressions reduce to a perturbative solution of the Bloch equations.

On the other hand, if there is strong electron-electron scattering at certain $`𝐤`$ vectors, e.g., due to Fermi surface nesting, the rates can be strongly $`𝐤`$ dependent. If many states of different relaxation rates enter the SHG photon yield, then for broad bands the description of SHG using optical Bloch equations with a single intermediate state is not justified. If interference patterns are fitted with results from Bloch equations, there is no simple relation between the extracted relaxation rate and the dephasing rates of the excited electrons.

For 2PPE the situation is quite different, since this method allows to probe specific momenta $`𝐤`$ in the Brillouin zone. Here, the band width is not crucial. There are typically several contributions to the photoelectron yield, since there are several unoccupied bands. The contributions are again weighted with the inverse detuning, but now there is only a small number of states involved for fixed photoelectron momentum $`𝐤`$. Thus a description in terms of a small number of states, e.g., by optical Bloch equations, is valid. However, if one experimentally integrates over $`𝐤`$, 2PPE behaves much like SFG.

As mentioned above, 2PPE is generally accompanied by SFG, although the latter is generally smaller in intensity due to the appearance of additional dipole matrix elements. It might be interesting to perform both SFG and 2PPE experiments on the same sample. The two techniques are complementary in that 2PPE gives information about specific points in the Brillouin zone whereas SFG averages over the whole zone. Furthermore, comparison of Eqs. (10) and (24) shows that while the general form of the expressions for SFG and 2PPE is similar, they do depend on the material parameters in quite different ways. We give three examples: First, 2PPE also depends on the energy relaxation rates (lifetimes) directly, whereas SFG only depends on the dephasing rates. Secondly, SFG crucially depends on the dipole matrix element $`𝐃_{31}`$ of the transition from the excited state above $`E_{\mathrm{vac}}`$ to the original state in the Fermi sea, whereas 2PPE does not. Thus comparison of SFG and 2PPE may proof useful for measuring the dipole matrix elements. Thirdly, SFG generally results from a much thinner surface region than 2PPE and the relaxation rates obtained from 2PPE are more bulk-like, allowing to study the dependence of the rates on the distance from the surface. Finally, we have shown that there is a contribution to 2PPE from SFG light generated within the solid, see Eqs. (20) and (21) as well as Figs. 7(b) and (c). Simultaneous measurement of SFG and 2PPE may allow to detect this interesting effect.

### III.2 Collective plasma excitations

Since SFG and 2PPE may be strongly enhanced by collective plasma excitations, it is useful to discuss them in the framework of the response theory. Our goal is to show plasmon enhancement of SFG and 2PPE *in principle*, even though the plasma frequency is larger than currently accessible laser frequencies in some metals (but not, for example, in silver, many heavy-fermion metals, and the interesting compound $`\mathrm{MgB}_2`$). Note, the plasma frequency is smaller in clusters, for which the same general picture applies.

We have seen in Sec. II that in both SFG and 2PPE field enhancement of the effective electric field $`𝐄`$ is described by the Fresnel formulae SMD87 ; HBB94 ; Luce97 ; LuceLDA ; AtK02 . This mechanism is relevant at the frequency of the exciting laser field and, in the case of SFG, also at the sum frequency for the outgoing SFG light. It corresponds to the coupling of the external light field to *plasmon-polaritons* in the solid Board .

The second important origin of plasmon enhancement is the screening of the nonlinear polarization $`𝐏^{(2)}`$, which is caused by the effective electric field $`𝐄^{(2)}`$ accompanying the longitudinal part of $`𝐏^{(2)}`$ and appears at the sum frequency. Since the electric field $`𝐄^{(2)}`$ is longitudinal, a true *plasmon* excitation is involved. It is important to remember that the exciting light does couple to plasmons in real solids; this coupling is only absent in simple jellium models Barton ; Board .

Note, in the case of pump-probe SFG with two pulses of different mean frequencies $`\omega _1`$ and $`\omega _2`$, one of them and the sum frequency $`\omega _1+\omega _2`$ can be close to the plasma frequency. This so-called double resonance leads to a particularly strong enhancement Lie99 . Surface plasmons lie outside the scope of this paper, since they require a more detailed description of the surface. See Ref. Board  for a discussion.

Motivated by 2PPE experiments on clusters Gerb , we briefly consider the plasmon decay. A plasmon decays into a single particle-hole pair BVV00 . The probability of this decay is determined by the phase space available for the final electron-hole pair. It is only energetically possible if the plasmon dispersion lies within the electron-hole continuum at the plasmon momentum $`𝐪`$, leading to Landau damping. On the other hand, decay into a single electron-hole pair may be possible close to the surface, since translational symmetry is broken and $`q_z`$ is not conserved. If the energy of the electron is higher than the vacuum energy, photoemission may result. Creation of several pairs is possible by subsequent inelastic electron-electron scattering. A plasmon also looses energy through inelastic scattering of the virtual electrons and holes in the loop in Fig. 5(b). This process is governed by the single-particle relaxation rates. The plasmon lifetime is thus shorter than typical lifetimes of the relevant excited electrons.

A plasma mode can be multiply excited. In a recent 2PPE experiment a doubly excited plasma mode of silver nanoparticles decays into a *single* electron-hole pair Gerb ; rem.clusters . What is actually observed is an enhancement of the 2PPE yield when the sum frequency is close to *twice* the plasma frequency. The origin of this effect is that 2PPE with the incident-light frequency close to the plasma frequency is enhanced due to field enhancement. The general process is not specific to clusters but is also relevant for flat surfaces. It would be interesting to look for this effect experimentally.

### III.3 Further remarks

Concerning the range of validity of the second-order response theory we remark the following. We consider first pump-probe SFG with a very long delay time $`\mathrm{\Delta }T`$. The first-order density operator $`\rho ^{(1)}`$ describes the result of the first interaction. It contains a finite polarization (off-diagonal components) but no change of occupation (diagonal components), see Eq. (38). A change of occupation is only obtained from $`\rho ^{(2)}`$ and higher-order contributions, which involve a larger number of dipole matrix elements and are usually small compared to $`\rho ^{(1)}`$. However, the off-diagonal components usually decay faster than the diagonal ones so that for long delay times the higher-order change of occupation can dominate over the second-order polarization and the second-order approximation breaks down. On the other hand, our expressions for 2PPE already include the change of occupation due to the first pulse, since we have directly calculated the photoelectron yield to fourth order. Thus the results should hold even for long delay times. For the case that the polarization has decayed at the time of the second pulse, but the non-equilibrium occupation has not, the resulting limiting form of $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ is given by Eq. (72). It only depends on the energy relaxation rates. This is the case where rate equations are appropriate Kno00 ; Rol1 ; Rol2 . Due to the vanishing polarization there are no interference effects.

There is an alternative and physically appealing description of pump-probe SFG and 2PPE as a two-step process: The first pulse creates a non-equilibrium distribution, which is probed by the second one. We now discuss the validity of calculations based on this picture. In App. A we derive an expression for the linear susceptibility of an electron gas in an arbitrary non-equilibrium state described by the density matrix $`\rho _{\mathrm{neq}}`$, see Eq. (45). If we insert $`\rho ^{(1)}`$ due to the first pulse for $`\rho _{\mathrm{neq}}`$, we obtain a two-step description for $`\chi ^{(2)}`$ and the polarization $`𝐏^{(2)}`$. Omitting the details, we only state that the result is identical to the one obtained directly for the second-order polarization $`𝐏^{(2)}`$, Eq. (5), but with the full susceptibility $`\chi ^{(2)}`$ replaced by its irreducible part $`\chi _{\mathrm{irr}}^{(2)}`$ of Eq. (10). Thus, by assuming two separate interaction processes and treating each in a first-order approximation, we loose the screening of the second-order polarization. This is not justified if the sum frequency lies close to the plasma resonance.

Next we consider a two-step description of 2PPE: The 2PPE photoelectron yield $`𝒩^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ is expressed in terms of an arbitrary non-equilibrium density matrix $`\rho _{\mathrm{neq}}`$ as discussed in App. B. Then the second-order density matrix $`\rho ^{(2)}`$ due to the first pulse is inserted for $`\rho _{\mathrm{neq}}`$. We reobtain the full irreducible fourth-order result $`𝒩_{\mathrm{irr}}^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ of Fig. 7(a), but only part of the reducible contributions, Fig. 7(b), (c): The two-step description neglects contributions of two photons out of different pulses being converted into one SHG photon. These contributions may become important if the sum frequency is close to a plasma resonance. In conclusion, the two-step picture of SFG and 2PPE is valid unless the response at the sum frequency is enhanced by plasmon effects.

Finally, we emphasize that our theory can also describe time-resolved SFG and 2PPE from ferromagnetically ordered systems. Ultimately, the light couples to the (spin) magnetization through spin-orbit coupling, which is incorporated, in principle, in the dipole matrix elements $`𝐃`$ and the band structure. The spin-dependent matrix elements can be calculated by a perturbative expansion in the spin-orbit coupling HB ; NOLI . SFG and 2PPE also depend on magnetic order through the band energies $`E_{𝐤_{}l}`$ and relaxation rates $`\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}^{}l^{}}`$, since $`l`$ also contains a spin index $`\sigma `$. Of particular importance for magnetically-ordered materials is the rotation of the polarization of SHG light relative to incident light (NOLIMOKE) HB ; NOLI ; Luce96 . As mentioned above, the light polarization is controlled by the symmetries of the tensor $`\chi _{ijk}^{(2)}`$, which are known for low-index surfaces Pan .

Compared to NOLIMOKE, 2PPE for magnetic systems has the advantage that in principle one can obtain information on the spin-dependent lifetimes of electrons in specific states $`|𝐤_{}l`$. The dependence of 2PPE on the light polarization for magnetic systems has not been studied so far. In the response theory this dependence is controlled by the symmetries of the tensor $`\eta _{ijkl}^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$, as mentioned in Sec. II.2.

For pump-probe experiments with *long* time delay $`\mathrm{\Delta }T`$ the main contribution to 2PPE comes from the change of occupation brought about by the pump pulse. *Only* in this case the photoelectron yield is proportional to the occupation of the corresponding intermediate states. If in addition the matrix elements and the relaxation rates out of vacuum states depend only weakly on spin, then Eq. (60) shows that the 2PPE yield becomes proportional to the spin-dependent occupation of these intermediate states:

$$𝒩^{2\mathrm{P}\mathrm{P}\mathrm{E}}(\mathrm{\Delta }T;𝐤,\sigma )\rho _{\mathrm{neq};𝐤\nu \sigma ;𝐤\nu \sigma }(\mathrm{\Delta }T)=n_{\mathrm{neq};𝐤\nu \sigma }(\mathrm{\Delta }T),$$
(27)
where $`n_{\mathrm{neq};𝐤\nu \sigma }`$ denotes the non-equilibrium occupation of the state $`|𝐤\nu \sigma `$ after the pump pulse. Then the difference of the spin-up and spin-down 2PPE yield, $`𝒩^{2\mathrm{P}\mathrm{P}\mathrm{E}}(\mathrm{\Delta }T;𝐤,)𝒩^{2\mathrm{P}\mathrm{P}\mathrm{E}}(\mathrm{\Delta }T;𝐤,)`$, is proportional to the difference of the occupations and thus to the transient magnetization of the intermediate states.

Note, circularly polarized light might excite electrons spin-selectively due to angular-momentum conservation. In our response theory these selection rules are incorporated in the dipole matrix elements $`𝐃`$. Conversely, spin-selective excitation will lead to corresponding polarization of the SFG light. The use of circularly polarized light in 2PPE and SFG is of particular interest regarding ferromagnets and transient magnetizations.

### III.4 Conclusions

To summarize, we have presented a unified perturbative response theory for time-resolved SFG and 2PPE. The theory is fully quantum-mechanical and contains the interference effects described by off-diagonal components of the density matrix. It does not rely on any assumption about the time or frequency dependence of the exciting laser pulses. The solid is described by its band structure, relaxation rates, and dipole matrix elements. We have discussed metals but the response theory can be applied to semiconductors and insulators as well, see, e.g., Ref. ShA03 . Since the theory is formulated directly in the time domain, it presents a suitable framework for the discussion of the time-dependent physics of SFG and 2PPE. We have shown that similar information as from 2PPE can be gained from SHG. Of course, 2PPE is sensitive to specific momenta $`𝐤`$ in the Brillouin zone, while SHG in general is not. A simple tight-binding model of a metal has been studied in order to show that the theory gives reasonable numerical results and to illustrate the following effects important for the understanding of SFG and 2PPE.

We have shown how relaxation rates and detuning affect the interference patterns in single-color pump-probe SHG and 2PPE experiments: The lifetime in the intermediate states and their detuning with respect to the photon energy lead to a similar narrowing of the interference patterns. The effect of detuning must be taken into account in order to extract meaningful life times from such experiments. Also, in particular in SHG the measured relaxation rate is a weighted average over the relaxation rates of many excited states. Furthermore, the weights in this average change with the pump-probe delay. Thus different rates govern the decay of the interference pattern depending on the pump-probe delay—the decay is not simply exponential. We have also discussed the range of validity of the optical Bloch equations, applicable if only a few states contribute, and of semiclassical rate equations valid for very long pump-probe delays. Both approaches are limiting cases of our theory.

Finally, we have considered the role played by collective plasma excitations. Plasmon effects in both SFG and 2PPE can only partly be understood in terms of field enhancement at the surface. One also has to take the electric field accompanying a nonlinear polarization of the electron system into account. This effect leads to interesting additional contributions to 2PPE, in which incoming photons are converted into sum-frequency photons which then lead to *ordinary* photoemission. These contributions should be observable if the sum frequency is close to the plasma frequency.

###### Acknowledgements.

We thank Drs. W. Pfeiffer, G. Gerber, G. Bouzerar, and R. Knorren for valuable discussions. Financial support by the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich 290 is gratefully acknowledged.

## Appendix A Response theory for the nonlinear optical response

In this appendix we derive the transverse second-order susceptibility and polarization for arbitrary pulse shapes of the exciting laser field. The resulting expressions allow to calculate the SFG yield for arbitrary pulse shapes, thereby going beyond the results of Hübner and Bennemann HB  for continuous-wave, monochromatic light. The self-consistent-field approach EC  is applied to a solid described by its band structure and relaxation rates. The flat surface is assumed to lie at $`z=0`$ with the solid at $`z<0`$. We neglect the intraband contribution, which is reasonable for optical frequencies.

The single-particle Hamiltonian is $`H=H_0+V`$, where $`H_0`$ describes the unperturbed solid with the normalized eigenstates $`|𝐤_{}l`$ and eigenenergies $`E_{𝐤_{}l}`$. $`𝐤_{}`$ is the crystal momentum parallel to the surface and all other quantum numbers, discrete as well as continuous, are collectively denoted by $`l`$ (see the discussion at the beginning of Sec. II.1).

The time-dependent perturbation is Schiff

$$V(𝐫,t)=\frac{ie\mathrm{}}{mc}𝐀(𝐫,t)\mathbf{}\frac{ie\mathrm{}}{2mc}[\mathbf{}𝐀(𝐫,t)]$$
(28)
with the vector potential $`𝐀`$, which is treated classically. We have made the usual approximation to neglect the quadratic term in $`𝐀`$ and have used a gauge with vanishing scalar potential. We will later need the temporal Fourier transform (using the convention of Ref. HB )

$`V(𝐫,\omega )`$ $`=`$ $`{\displaystyle \frac{dt}{2\pi }e^{i\omega t}V(𝐫,t)}`$ (29)
$`=`$ $`{\displaystyle \frac{e\mathrm{}}{m\omega }}𝐄(𝐫,\omega )\mathbf{}+{\displaystyle \frac{e\mathrm{}}{2m\omega }}[\mathbf{}𝐄(𝐫,\omega )],`$
where we have used $`𝐄=(1/c)𝐀/t`$. Inserting the spatial Fourier transform we get

$$V(𝐫,\omega )=\frac{e\mathrm{}}{m\omega }\underset{𝐪}{}e^{i𝐪𝐫}𝐄(𝐪,\omega )\left(\mathbf{}\frac{i𝐪}{2}\right).$$
(30)
The matrix elements of $`V`$ are

$`𝐤_{}l|V|𝐤_{}+𝐪_{},l^{}`$ $`=`$ $`{\displaystyle \frac{e\mathrm{}}{m\omega }}{\displaystyle \underset{q_z}{}}𝐄(𝐪,\omega )𝐤_{}l|e^{i𝐪𝐫}\left(\mathbf{}i𝐪/2\right)|𝐤_{}+𝐪_{},l^{}`$ (31)
$``$ $`e{\displaystyle \underset{q_z}{}}𝐄(𝐪,\omega )𝐃_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}(q_z,\omega )`$
using momentum conservation and $`𝐪=(𝐪_{},q_z)`$. Here,

$$𝐃_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}(q_z,\omega )=\frac{\mathrm{}}{m\omega }𝐤_{}l|e^{i𝐪𝐫}\left(\mathbf{}i𝐪/2\right)|𝐤_{}+𝐪_{},l^{}.$$
(32)
If the field $`𝐄`$ were purely transverse the term $`i𝐪/2`$ in the parentheses would drop out, but this is not guaranteed close to a surface. It is not our goal to calculate $`𝐃`$ explicitly. We only remark that if one uses the dipole approximation HK90 ; HBB94 ; Luce97 ; LuceLDA ; NOLI ; Luce96 ; Luce98 ; AnH02 ; Shala  $`e^{i𝐪𝐫}1`$, the contribution from $`i𝐪/2`$ vanishes since we neglect the intraband contributions so that $`l^{}l`$, and the remainder gives Klings ; Voon

$`{\displaystyle \frac{\mathrm{}}{m\omega }}𝐤_{}l|\mathbf{}|𝐤_{}+𝐪_{},l^{}={\displaystyle \frac{1}{\mathrm{}\omega }}𝐤_{}l|[H,𝐫]|𝐤_{}+𝐪_{},l^{}`$ (33)
$`={\displaystyle \frac{E_{𝐤_{}l}E_{𝐤_{}+𝐪_{},l^{}}}{\mathrm{}\omega }}𝐤_{}l|𝐫|𝐤_{}+𝐪_{},l^{}.`$
Since the response is dominated by contributions with $`\mathrm{}\omega E_{𝐤_{}+𝐪_{},l^{}}E_{𝐤_{}l}`$ the prefactor can be further approximated by unity if the frequency spectrum of the incoming light is sufficiently narrow. The dipole approximation should be justified since the electric field changes slowly on the scale of the lattice constant (the skin depth is about one order of magnitude larger than the lattice constant). However, our response theory for SFG does not require the dipole approximation to be made.

The time evolution of the density operator $`\rho `$ is described by the master or von Neumann equation Louisell ; Loudon

$$\frac{d}{dt}\rho =\frac{1}{i\mathrm{}}[H,\rho ]+[\rho ].$$
(34)
The functional $`[\rho ]`$ represents relaxation terms made explicit below. Matrix elements of $`\rho `$ are written as $`\rho _{𝐤_{}l;𝐤_{}^{}l^{}}𝐤_{}l|\rho |𝐤_{}^{}l^{}`$. The master equation then reads Loudon

$`{\displaystyle \frac{d}{dt}}\rho _{𝐤_{}l;𝐤_{}^{}l^{}}={\displaystyle \frac{1}{i\mathrm{}}}𝐤_{}l|[H_0+V,\rho ]|𝐤_{}^{}l^{}+\delta _{𝐤_{}𝐤_{}^{}}\delta _{ll^{}}`$ (35)
$`\times {\displaystyle \underset{𝐤_{}^{\prime \prime }l^{\prime \prime }}{}^{}}\gamma _{𝐤_{}l;𝐤_{}^{\prime \prime }l^{\prime \prime }}\rho _{𝐤_{}^{\prime \prime }l^{\prime \prime };𝐤_{}^{\prime \prime }l^{\prime \prime }}\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}^{}l^{}}\rho _{𝐤_{}l;𝐤_{}^{}l^{}}.`$
Here, $`\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}l}\tau _{𝐤_{}l}^1`$ is the inverse lifetime of state $`|𝐤_{}l`$, which arises mainly from inelastic electron-electron scattering. $`\gamma _{𝐤_{}l;𝐤_{}^{}l^{}}`$ gives the rate of spontaneous transitions from state $`|𝐤_{}^{}l^{}`$ to state $`|𝐤_{}l`$. Because of conservation of electron number

$$\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}l}=\underset{𝐤_{}^{}l^{}}{}^{}\gamma _{𝐤_{}^{}l^{};𝐤_{}l}.$$
(36)
Primed sums run over all states except $`|𝐤_{}l`$. $`\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}l}`$ and $`\gamma _{𝐤_{}l;𝐤_{}^{}l^{}}`$ describe energy relaxation, i.e., the change of the diagonal components of $`\rho `$, whereas the dephasing rate $`\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}^{}l^{}}`$ with $`|𝐤_{}l|𝐤_{}^{}l^{}`$ describes relaxation of the off-diagonal components.

To solve the master equation (35) perturbatively, the density operator is expanded in powers of the perturbation $`V`$ as $`\rho =\rho ^{(0)}+\rho ^{(1)}+\rho ^{(2)}+\mathrm{}`$ In thermal equilibrium the unperturbed density matrix is expressed in terms of the Fermi function, $`\rho _{𝐤_{}l;𝐤_{}^{}l^{}}^{(0)}=\delta _{𝐤_{}𝐤_{}^{}}\delta _{ll^{}}f(E_{𝐤_{}l})`$. The temporal Fourier transform of Eq. (35) reads

$`i\omega \rho _{𝐤_{}l;𝐤_{}^{}l^{}}(\omega )=\left({\displaystyle \frac{E_{𝐤_{}l}E_{𝐤_{}^{}l^{}}}{i\mathrm{}}}\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}^{}l^{}}\right)\rho _{𝐤_{}l;𝐤_{}^{}l^{}}(\omega )`$ (37)
$`+{\displaystyle \frac{1}{i\mathrm{}}}{\displaystyle \frac{dt}{2\pi }e^{i\omega t}\underset{𝐤_{}^{\prime \prime }l^{\prime \prime }}{}\left[𝐤_{}l|V(t)|𝐤_{}^{\prime \prime }l^{\prime \prime }\rho _{𝐤_{}^{\prime \prime }l^{\prime \prime };𝐤_{}^{}l^{}}(t)\rho _{𝐤_{}l;𝐤_{}^{\prime \prime }l^{\prime \prime }}(t)𝐤_{}^{\prime \prime }l^{\prime \prime }|V(t)|𝐤_{}^{}l^{}\right]}+\delta _{𝐤_{}𝐤_{}^{}}\delta _{ll^{}}{\displaystyle \underset{𝐤_{}^{\prime \prime }l^{\prime \prime }}{}^{}}\gamma _{𝐤_{}l;𝐤_{}^{\prime \prime }l^{\prime \prime }}\rho _{𝐤_{}^{\prime \prime }l^{\prime \prime };𝐤_{}^{\prime \prime }l^{\prime \prime }}(\omega ).`$
Keeping only terms linear in $`V`$ one obtains

$`\rho _{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(1)}(\omega )=e{\displaystyle \underset{q_z}{}}𝐃_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}(q_z,\omega )𝐄(𝐪,\omega )`$ (38)
$`\times {\displaystyle \frac{f(E_{𝐤_{}+𝐪_{},l^{}})f(E_{𝐤_{}l})}{\mathrm{}\omega +E_{𝐤_{}+𝐪_{},l^{}}E_{𝐤_{}l}+i\mathrm{}\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}}}.`$
Note that the diagonal components vanish: There is no change of occupation to first order.

The polarization is given by the thermal average of $`e𝐃`$, which is the conjugate of the electric field according to Eq. (31),

$`𝐏(𝐪,\omega )`$ $`=`$ $`{\displaystyle \frac{e}{v}}{\displaystyle \underset{𝐤_{}}{}}{\displaystyle \underset{ll^{}}{}}\rho _{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}(\omega )𝐃_{𝐤_{}+𝐪_{},l^{};𝐤_{}l}(q_z,\omega ),`$ (39)
where $`v`$ is the volume. To first order

$$P_i^{(1)}(𝐪_{},q_z,\omega )=\underset{q_z^{}}{}\chi _{ij}^{(1)}(𝐪_{},q_z,q_z^{},\omega )E_j(𝐪_{},q_z^{},\omega ),$$
(40)
with the linear susceptibility of Lindhard form

$`\chi _{ij}^{(1)}(𝐪_{},q_z,q_z^{},\omega )={\displaystyle \frac{e^2}{v}}{\displaystyle \underset{𝐤_{}}{}}{\displaystyle \underset{ll^{}}{}}D_{𝐤_{}+𝐪_{},l^{};𝐤_{}l}^i(q_z,\omega )`$ (41)
$`\times D_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^j(q_z^{},\omega )`$
$`\times {\displaystyle \frac{f(E_{𝐤_{}+𝐪_{},l^{}})f(E_{𝐤_{}l})}{\mathrm{}\omega +E_{𝐤_{}+𝐪_{},l^{}}E_{𝐤_{}l}+i\mathrm{}\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}}},`$
shown diagrammatically in Fig. 3. It takes into account that the $`z`$ component of momentum is not conserved.

If we neglect the frequency dependence of $`𝐃`$ we obtain a result in the time domain. Equation (40) and the Fourier transform of Eq. (41) give

$`P_i^{(1)}(𝐪_{},q_z,t)`$ $`=`$ $`{\displaystyle \underset{q_z^{}}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dt_1}{2\pi }}\chi _{ij}(𝐪_{},q_z,q_z^{},tt_1)`$ (42)
$`\times E_j(𝐪_{},q_z^{},t_1),`$
with

$`\chi _{ij}(𝐪_{},q_z,q_z^{};tt_1)={\displaystyle \frac{e^2}{v}}{\displaystyle \frac{2\pi i}{\mathrm{}}}\mathrm{\Theta }(tt_1){\displaystyle \underset{𝐤_{}}{}}{\displaystyle \underset{ll^{}}{}}D_{𝐤_{}+𝐪_{},l^{};𝐤_{}l}^i(q_z)D_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^j(q_z^{})`$ (43)
$`\times [f(E_{𝐤_{}+𝐪_{},l^{}})f(E_{𝐤_{}l})]\mathrm{exp}\left[i{\displaystyle \frac{E_{𝐤_{}+𝐪_{},l^{}}E_{𝐤_{}l}}{\mathrm{}}}(tt_1)\right]\mathrm{exp}[\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}(tt_1)].`$
The response is thus only non-zero if $`t>t_1`$, which expresses causality.

The above results have been obtained under the assumption that the system is initially in thermal equilibrium. We now consider the non-equilibrium case. Then the unperturbed polarization $`𝐏_{\mathrm{neq}}^{(0)}`$ is, in general, non-vanishing so that the electrons experience an effective field even in the absence of an external perturbation, leading to a master equation that is nonlinear in the non-equilibrium density operator $`\rho _{\mathrm{neq}}^{(0)}`$. We assume, however, that this effect of electron-electron interaction is negligible. Then the linear response of a non-equilibrium system can be written as

$$P_{\mathrm{neq};i}^{(1)}(𝐪_{},q_z,\omega )=\underset{𝐪_{}^{},q_z^{}}{}_{\mathrm{}}^{\mathrm{}}𝑑\omega ^{}\chi _{\mathrm{neq};ij}^{(1)}(𝐪_{},q_z,\omega ;𝐪_{}^{},q_z^{},\omega ^{})E_j(𝐪^{},\omega ^{}),$$
(44)
with

$`\chi _{\mathrm{neq};ij}^{(1)}(𝐪_{},q_z,\omega ;𝐪_{}^{},q_z^{},\omega ^{})={\displaystyle \frac{e^2}{v}}{\displaystyle \underset{𝐤_{}}{}}{\displaystyle \underset{ll^{}l^{\prime \prime }}{}}D_{𝐤_{}+𝐪_{},l^{};𝐤_{}l}^i(q_z,\omega )`$ (45)
$`\times [D_{𝐤_{}l;𝐤_{}+𝐪_{}^{},l^{\prime \prime }}^j(q_z^{},\omega ^{})\rho _{\mathrm{neq};𝐤_{}+𝐪_{}^{},l^{\prime \prime };𝐤_{}+𝐪_{},l^{}}^{(0)}(\omega \omega ^{})`$
$`\rho _{\mathrm{neq};𝐤_{}l;𝐤_{}+𝐪_{}𝐪_{}^{},l^{\prime \prime }}^{(0)}(\omega \omega ^{})D_{𝐤_{}+𝐪_{}𝐪_{}^{},l^{\prime \prime };𝐤_{}+𝐪_{},l^{}}^j(q_z^{},\omega ^{})].`$
This equation gives the linear susceptibility in terms of an arbitrary density operator $`\rho _{\mathrm{neq}}^{(0)}`$.

To return to the response of an equilibrium system, we now consider the second-order contribution. We collect the terms in the master equation (37) that are of second order in the effective electric field,

$`i\omega \rho _{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}(E_{𝐤_{}l}E_{𝐤_{}+𝐪_{},l^{}})\rho _{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(2)}`$ (46)
$`+{\displaystyle \frac{1}{i\mathrm{}}}𝐤_{}l|[V^{(1)},\rho ^{(1)}]|𝐤_{}+𝐪_{},l^{}+{\displaystyle \frac{1}{i\mathrm{}}}𝐤_{}l|[V^{(2)},\rho ^{(0)}]|𝐤_{}+𝐪_{},l^{}`$
$`+\delta _{𝐪_{}0}\delta _{ll^{}}{\displaystyle \underset{𝜿_{}\lambda }{}^{}}\gamma _{𝐤_{}l;𝜿_{}\lambda }\rho _{𝜿_{}\lambda ;𝜿_{}\lambda }^{(2)}\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}\rho _{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(2)}.`$
$`V^{(1)}V`$ is the perturbation by the effective field. Higher-order perturbations $`V^{(n)}`$, $`n2`$, result from the electric field due to the displaced charge calculated at order $`n`$. From the expression for the electric field due to a polarization $`𝐏^{(n)}`$, Eq. (12) Jackson , one obtains $`𝐄^{(n)}(𝐪)=4\pi \widehat{𝐪}\widehat{𝐪}𝐏^{(n)}(𝐪)`$, where $`\widehat{𝐪}𝐪/|𝐪|`$. Importantly, this additional field also has to be taken into account as a perturbation HB . Specifically, the longitudinal part of $`𝐏^{(n)}`$ leads to a perturbation with matrix elements

$$𝐤_{}l|V^{(n)}|𝐤_{}+𝐪_{},l^{}=e\underset{q_z}{}𝐄^{(n)}(𝐪,\omega )𝐃_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}(q_z,\omega ),$$
(47)
see Eq. (31). Note, due to the reduced symmetry at the surface a transverse electric field in general leads to a polarization with a longitudinal component.

Since the field $`𝐄^{(2)}=4\pi \widehat{𝐪}\widehat{𝐪}𝐏^{(2)}`$ is explicitly of *second* order in the applied field, see Eq. (5), it must be taken into account in our calculation of the second-order response. From Eqs. (39) and (46) we then obtain for $`𝐪_{}0`$

$`P_i^{(2)}(𝐪,\omega )`$ $`=`$ $`{\displaystyle \frac{e^2}{v}}{\displaystyle d^2k_{}\underset{ll^{}}{}\frac{D_{𝐤_{}+𝐪_{},l^{};𝐤_{}l}^i(q_z,\omega )}{\mathrm{}\beta +E_{𝐤_{}+𝐪_{},l^{}}E_{𝐤_{}l}+i\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}}}`$ (48)
$`\times {\displaystyle }d^3q^{}{\displaystyle \underset{\lambda }{}}{\displaystyle }d\omega ^{}[D_{𝐤_{}l;𝐤_{}+𝐪_{}^{},\lambda }^j(q_z^{},\omega ^{})\rho _{𝐤_{}+𝐪_{}^{},\lambda ;𝐤_{}+𝐪_{},l^{}}^{(1)}(\omega \omega ^{})`$
$`\rho _{𝐤_{}l;𝐤_{}+𝐪_{}𝐪_{}^{},\lambda }^{(1)}(\omega \omega ^{})D_{𝐤_{}+𝐪_{}𝐪_{}^{},\lambda ;𝐤_{}+𝐪_{},l^{}}^j(q_z^{},\omega ^{})]E_j(𝐪^{},\omega ^{})`$
$`+4\pi \left[f(E_{𝐤_{}+𝐪_{},l^{}})f(E_{𝐤_{}l})\right]{\displaystyle 𝑑q_z^{}D_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^j(q_z^{},\omega )\widehat{q}_j\widehat{q}_kP_k^{(2)}(𝐪^{},\omega )},`$
where in the last term $`𝐪_{}^{}=𝐪_{}`$. Note that the nonlinear polarization $`𝐏^{(2)}`$ appears on both sides of the equation. Solving this equation for $`𝐏^{(2)}`$ we obtain

$$P_i^{(2)}(𝐪,\omega )=\underset{𝐪^{},𝐪^{\prime \prime }}{}_{\mathrm{}}^{\mathrm{}}𝑑\omega ^{}\chi _{ijk}^{(2)}(𝐪,𝐪^{};\omega ,\omega ^{})E_j(𝐪^{},\omega ^{})E_k(𝐪𝐪^{},\omega \omega ^{}),$$
(49)
with the second-order susceptibility

$$\chi _{ijk}^{(2)}(𝐪,𝐪^{};\omega ,\omega ^{})=\underset{m}{}\underset{\kappa _z}{}\epsilon _{\mathrm{long};im}^1(𝐪_{},q_z,\kappa _z,\omega )\chi _{\mathrm{irr};mjk}^{(2)}((𝐪_{},\kappa _z),𝐪^{},𝐪𝐪^{};\omega ,\omega ^{})$$
(50)
and

$$\epsilon _{\mathrm{long};ij}(𝐪_{},q_z,\kappa _z,\omega )\delta _{ij}+4\pi \chi _{im}(𝐪_{},q_z,\kappa _z,\omega )\frac{(𝐪_{},\kappa _z)_k}{|(𝐪_{},\kappa _z)|}\frac{(𝐪_{},\kappa _z)_j}{|(𝐪_{},\kappa _z)|}.$$
(51)
Here, $`\epsilon _{\mathrm{long};ij}^1(𝐪_{},q_z,\kappa _z,\omega )`$ is the *inverse matrix* of $`\epsilon _{\mathrm{long}}`$ with respect to the indices $`(i,q_z)`$ and $`(j,\kappa _z)`$. Equation (51) means that $`\epsilon _{\mathrm{long}}`$ only acts on the longitudinal component and is unity for the transverse ones. Solving for $`𝐏^{(2)}`$ is found to be equivalent to the summation of an RPA series for the electron-electron interaction mediated by the electromagnetic field. In this language the interaction is absorbed into $`\chi _{im}`$ in Eq. (51).

The irreducible susceptibility in Eq. (50) reads

$`\chi _{\mathrm{irr};ijk}^{(2)}(𝐪,𝐪^{},𝐪^{\prime \prime };\omega ,\omega ^{})={\displaystyle \frac{e^3}{v}}{\displaystyle \underset{𝐤_{}}{}}{\displaystyle \underset{ll^{}l^{\prime \prime }}{}}{\displaystyle \frac{D_{𝐤_{}+𝐪_{},l;𝐤_{}l^{\prime \prime }}^i(q_z,\omega )}{\mathrm{}\omega +E_{𝐤_{}+𝐪_{},l}E_{𝐤_{}l^{\prime \prime }}+i\mathrm{}\mathrm{\Gamma }_{𝐤_{}l^{\prime \prime };𝐤_{}+𝐪_{},l}}}`$ (52)
$`\times [D_{𝐤_{}l^{\prime \prime };𝐤_{}+𝐪_{}^{},l^{}}^j(q_z^{},\omega ^{})D_{𝐤_{}+𝐪_{}^{},l^{};𝐤_{}+𝐪_{},l}^k(q_z^{\prime \prime },\omega \omega ^{}){\displaystyle \frac{f(E_{𝐤_{}+𝐪_{},l})f(E_{𝐤_{}+𝐪_{}^{},l^{}})}{\mathrm{}\omega +\mathrm{}\omega ^{}+E_{𝐤_{}+𝐪_{},l}E_{𝐤_{}+𝐪_{}^{},l^{}}+i\mathrm{}\mathrm{\Gamma }_{𝐤_{}+𝐪_{}^{},l^{};𝐤_{}+𝐪_{},l}}}`$
$`D_{𝐤_{}+𝐪_{}𝐪_{}^{},l^{};𝐤_{}+𝐪_{},l}^j(q_z^{},\omega ^{})D_{𝐤_{}l^{\prime \prime };𝐤_{}+𝐪_{}𝐪_{}^{},l^{}}^k(q_z^{\prime \prime },\omega \omega ^{}){\displaystyle \frac{f(E_{𝐤_{}+𝐪_{}𝐪_{}^{},l^{}})f(E_{𝐤_{}l^{\prime \prime }})}{\mathrm{}\omega +\mathrm{}\omega ^{}+E_{𝐤_{}+𝐪_{}𝐪_{}^{},l^{}}E_{𝐤_{}l^{\prime \prime }}+i\mathrm{}\mathrm{\Gamma }_{𝐤_{}l^{\prime \prime };𝐤_{}+𝐪_{}𝐪_{}^{},l^{}}}}],`$
shown diagrammatically in Fig. 5(a). Finally, the time dependence of the polarization $`𝐏^{(2)}`$ is obtained by Fourier transformation of Eq. (49) using Eq. (9), leading to Eq. (10).

## Appendix B Response theory for photoemission

In this appendix we give details of the derivation of the time-integrated photoelectron yield. We also present the analytical expressions for the response functions omitted in Sec. II.2. The starting point is again the master equation (34). The terms of order $`n1`$ can be calculated recursively,

$`{\displaystyle \frac{d}{dt}}\rho _{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(n)}={\displaystyle \frac{1}{i\mathrm{}}}(E_{𝐤_{}l}E_{𝐤_{}+𝐪_{},l^{}})\rho _{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(n)}`$ (53)
$`+{\displaystyle \frac{1}{i\mathrm{}}}{\displaystyle \underset{m=1}{\overset{n}{}}}𝐤_{}l|[V^{(m)},\rho ^{(nm)}]|𝐤_{}+𝐪_{},l^{}`$
$`+\delta _{𝐪_{}0}\delta _{ll^{}}{\displaystyle \underset{𝜿_{}\lambda }{}^{}}\gamma _{𝐤_{}l;𝜿_{}\lambda }\rho _{𝜿_{}\lambda ;𝜿_{}\lambda }^{(n)}`$
$`\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}\rho _{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(n)}.`$
Here, $`V^{(m)}`$ is the perturbation of order $`m`$, see the discussion leading to Eq. (47). Among the states $`|𝐤_{}l`$ etc. appearing in Eq. (53) are states $`|𝐤\sigma ,\mathrm{in}`$ lying in the crystal above the vacuum level. We assume that electrons leaving the crystal are in states $`|𝐤\sigma ,\mathrm{out}`$ and are detected without further interaction and without returning to the solid. Then the only way their occupation can change is through spontaneous transitions out of $`|𝐤\sigma ,\mathrm{in}`$, governed by the rate $`\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}`$. Note, in principle higher vacuum bands appear by shifting the (nearly) free electron dispersion back into the first Brillouin zone. We omit these bands for notational simplicity.

First, we consider the irreducible part $`\rho _{\mathrm{irr}}^{(n)}`$. This is the contribution of only the direct, first-order perturbation $`V^{(1)}=V`$ at every step of the recursion. The resulting equation reads

$`{\displaystyle \frac{d}{dt}}\rho _{\mathrm{irr};𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(n)}(t)={\displaystyle \frac{1}{i\mathrm{}}}(E_{𝐤_{}l}E_{𝐤_{}+𝐪_{},l^{}})\rho _{\mathrm{irr};𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(n)}(t)`$ (54)
$`+{\displaystyle \frac{e}{i\mathrm{}}}{\displaystyle \underset{𝐪^{}\lambda }{}}\left[𝐃_{𝐤_{}l;𝐤_{}+𝐪_{}^{},\lambda }(q_z^{})\rho _{\mathrm{irr};𝐤_{}+𝐪_{}^{},\lambda ;𝐤_{}+𝐪_{},l^{}}^{(n1)}(t)\rho _{\mathrm{irr};𝐤_{}l;𝐤_{}+𝐪_{}𝐪_{}^{},\lambda }^{(n1)}(t)𝐃_{𝐤_{}+𝐪_{}𝐪_{}^{},\lambda ;𝐤_{}+𝐪_{},l^{}}(q_z^{})\right]𝐄(𝐪^{},t)`$
$`+\delta _{𝐪_{}0}\delta _{ll^{}}{\displaystyle \underset{𝜿_{}\lambda }{}^{}}\gamma _{𝐤_{}l;𝜿_{}\lambda }\rho _{\mathrm{irr};𝜿_{}\lambda ;𝜿_{}\lambda }^{(n)}(t)\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}\rho _{\mathrm{irr};𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(n)}(t).`$
We assume $`𝐃`$ to be frequency-independent, see App. A. In Eq. (54) we write the first sum in the form $`_{𝐪^{}\lambda }`$ as a reminder that all components of the external momentum $`𝐪^{}`$ are summed over. Hence, we here exclude $`q_z^{}`$ from $`\lambda `$. The solution for the off-diagonal elements is

$`\rho _{\mathrm{irr};𝐤_{}l;𝐤_{}+𝐪_{},l^{}}^{(n)}(t)={\displaystyle \frac{e}{i\mathrm{}}}{\displaystyle _{\mathrm{}}^t}𝑑t_1{\displaystyle \underset{𝐪^{}\lambda }{}}e^{\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}+𝐪_{},l^{}}(tt_1)}\mathrm{exp}\left[i{\displaystyle \frac{E_{𝐤_{}l}E_{𝐤_{}+𝐪_{},l^{}}}{\mathrm{}}}(tt_1)\right]`$ (55)
$`\times \left[𝐃_{𝐤_{}l;𝐤_{}+𝐪_{}^{},\lambda }(q_z^{})\rho _{\mathrm{irr};𝐤_{}+𝐪_{}^{},\lambda ;𝐤_{}+𝐪_{},l^{}}^{(n1)}(t_1)\rho _{\mathrm{irr};𝐤_{}l;𝐤_{}+𝐪_{}𝐪_{}^{},\lambda }^{(n1)}(t_1)𝐃_{𝐤_{}+𝐪_{}𝐪_{}^{},\lambda ;𝐤_{}+𝐪_{},l^{}}(q_z^{})\right]𝐄(𝐪^{},t_1).`$
For the diagonal components there is an additional contribution from the relaxation of secondary electrons into the given state out of higher-energy states Rol1 ; Rol2 , i.e., the term with $`𝐪_{}=0`$, $`l^{}=l`$ in Eq. (54). The diagonal components can be written in the implicit form

$`\rho _{\mathrm{irr};𝐤_{}l;𝐤_{}l}^{(n)}(t)={\displaystyle \frac{e}{i\mathrm{}}}{\displaystyle _{\mathrm{}}^t}dt_1{\displaystyle \underset{𝐪^{}\lambda }{}}e^{\mathrm{\Gamma }_{𝐤_{}l;𝐤_{}l}(tt_1)}[𝐃_{𝐤_{}l;𝐤_{}+𝐪_{}^{},\lambda }(q_z^{})\rho _{\mathrm{irr};𝐤_{}+𝐪_{}^{},\lambda ;𝐤_{}l}^{(n1)}(t_1)`$ (56)
$`\rho _{\mathrm{irr};𝐤_{}l;𝐤_{}𝐪_{}^{},\lambda }^{(n1)}(t_1)𝐃_{𝐤_{}𝐪_{}^{},\lambda ;𝐤_{}l}(q_z^{})]𝐄(𝐪^{},t_1)+{\displaystyle }_{\mathrm{}}^tdt_1{\displaystyle }_{𝜿_{}\lambda }^{}\gamma _{𝐤_{}l;𝜿_{}\lambda }\rho ^{(n)}_{\mathrm{irr};𝜿_{}\lambda ;𝜿_{}\lambda }(t_1).`$
If the contribution of secondary electrons is small, Eq. (55) also applies for $`𝐪_{}=0`$, $`l^{}=l`$.

The photoelectron yield for momentum $`𝐤`$ and spin $`\sigma `$ is given by the time integrated photoelectron current. Equivalently, it can be written as the occupation of the appropriate vacuum state,

$$𝒩(𝐤,\sigma )=\rho _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{out}}(t\mathrm{}),$$
(57)
since electrons are assumed not to leave the states $`|𝐤\sigma ,\mathrm{out}`$ again. For the electron states outside of the crystal

$$\frac{d}{dt}\rho _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{out}}=\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}\rho _{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}},$$
(58)
since their occupation can only change due to electrons leaving the crystal. With Eq. (57) we find

$$𝒩(𝐤,\sigma )=\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}_{\mathrm{}}^{\mathrm{}}𝑑t\rho _{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}(t).$$
(59)
The irreducible contribution to order $`n`$ is obtained by inserting the irreducible part of $`\rho ^{(n)}`$, Eq. (56), into this equation. After changing the order of integrals we obtain

$`𝒩_{\mathrm{irr}}^{(n)}(𝐤,\sigma )={\displaystyle \frac{e}{i\mathrm{}}}{\displaystyle \frac{\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}}{\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t{\displaystyle \underset{𝐪\lambda }{}}`$ (60)
$`\times [𝐃_{𝐤\sigma ,\mathrm{in};𝐤_{}+𝐪_{},\lambda }(q_z)\rho _{\mathrm{irr};𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}^{(n1)}(t)`$
$`\rho _{\mathrm{irr};𝐤\sigma ,\mathrm{in};𝐤_{}𝐪_{},\lambda }^{(n1)}(t)𝐃_{𝐤_{}𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}(q_z)]𝐄(𝐪,t),`$
where we have used that there is no relaxation into the states above $`E_{\mathrm{vac}}`$ in the solid. Note, Eq. (60) describes photoemission out of any (possibly non-equilibrium) state.

For ordinary photoemission, $`𝒩^{(2)}`$, we have to calculate the density operator to first order, $`\rho ^{(1)}`$, which is purely irreducible. Consequently, the full ordinary photoelectron yield is obtained by inserting Eq. (55) for $`n=1`$ into Eq. (60),

$$𝒩^{(2)}(𝐤,\sigma )=\underset{𝐪}{}𝑑t_1𝑑t_2\eta _{ij}(𝐪;t_1,t_2;𝐤,\sigma )E_i(𝐪,t_1)E_j(𝐪,t_2)$$
(61)
with the response function

$`\eta _{ij}(𝐪;t_1,t_2;𝐤,\sigma )={\displaystyle \frac{e^2}{\mathrm{}^2}}{\displaystyle \frac{\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}}{\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}}}{\displaystyle \underset{\lambda }{}}`$ (62)
$`\times D_{𝐤\sigma ,\mathrm{in};𝐤_{}+𝐪_{},\lambda }^i(q_z)\mathrm{exp}\left[i{\displaystyle \frac{E_{𝐤_{}+𝐪_{},\lambda }E_{𝐤\sigma ,\mathrm{in}}}{\mathrm{}}}(t_2t_1)\right]`$
$`\times e^{\mathrm{\Gamma }_{𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}|t_2t_1|}f(E_{𝐤_{}+𝐪_{},\lambda })`$
$`\times D_{𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}^j(q_z).`$
It is useful to write our results in the frequency domain. For ordinary photoemission this yields

$$𝒩^{(2)}(𝐤,\sigma )=\underset{𝐪}{}_{\mathrm{}}^{\mathrm{}}𝑑\omega \eta _{ij}(𝐪,\omega ;𝐤,\sigma )E_i(𝐪,\omega )E_j(𝐪,\omega )$$
(63)
with

$`\eta _{ij}(𝐪,\omega ;𝐤,\sigma )={\displaystyle \frac{2\pi ie^2}{\mathrm{}}}{\displaystyle \frac{\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}}{\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}}}`$ (64)
$`\times {\displaystyle \underset{\lambda }{}}D_{𝐤\sigma ,\mathrm{in};𝐤_{}+𝐪_{},\lambda }^i(q_z)`$
$`\times ({\displaystyle \frac{1}{\mathrm{}\omega +E_{𝐤_{}+𝐪_{},\lambda }E_{𝐤\sigma ,\mathrm{in}}+i\mathrm{}\mathrm{\Gamma }_{𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}}}`$
$`{\displaystyle \frac{1}{\mathrm{}\omega +E_{𝐤_{}+𝐪_{},\lambda }E_{𝐤\sigma ,\mathrm{in}}i\mathrm{}\mathrm{\Gamma }_{𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}}})`$
$`\times f(E_{𝐤_{}+𝐪_{},\lambda })D_{𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}^j(q_z).`$
The third-order contribution, $`𝒩^{(3)}`$, is of interest since the irreducible third-order response appears in the reducible contributions to 2PPE. To calculate $`𝒩_{\mathrm{irr}}^{(3)}`$ from Eq. (60), we need the off-diagonal elements of $`\rho _{\mathrm{irr}}^{(2)}`$ only, i.e., the polarization of the electron system, which can be obtained from Eq. (55) alone. The result is

$`𝒩_{\mathrm{irr}}^{(3)}(𝐤,\sigma )={\displaystyle \underset{\mathrm{𝐪𝐪}^{}}{}}{\displaystyle 𝑑t_1𝑑t_2𝑑t_3\eta _{ijk}^{(3)}(𝐪,𝐪^{};t_1,t_2,t_3;𝐤,\sigma )}`$ (65)
$`\times E_i(𝐪,t_1)E_j(𝐪^{},t_2)E_k(𝐪𝐪^{},t_3)`$
with

$`\eta _{ijk}^{(3)}(𝐪,𝐪^{};t_1,t_2,t_3;𝐤,\sigma )={\displaystyle \frac{ie^3}{\mathrm{}^3}}{\displaystyle \frac{\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}}{\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}}}{\displaystyle \underset{\lambda \lambda ^{}}{}}D_{𝐤\sigma ,\mathrm{in};𝐤_{}+𝐪_{},\lambda }^i(q_z)D_{𝐤_{}+𝐪_{},\lambda ;𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}}^j(q_z^{})D_{𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{};𝐤\sigma ,\mathrm{in}}^k(q_zq_z^{})`$ (66)
$`\times (\mathrm{\Theta }(t_1t_2)\mathrm{\Theta }(t_2t_3)e^{i\mathrm{\Omega }_{𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}(t_1t_2)}e^{i\mathrm{\Omega }_{𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{};𝐤\sigma ,\mathrm{in}}(t_2t_3)}[f(E_{𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}})]`$
$`\mathrm{\Theta }(t_1t_3)\mathrm{\Theta }(t_3t_2)e^{i\mathrm{\Omega }_{𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}(t_1t_3)}e^{i\mathrm{\Omega }_{𝐤_{}+𝐪_{},\lambda ;𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}}(t_3t_2)}\left[f(E_{𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}})f(E_{𝐤_{}+𝐪_{},\lambda })\right]`$
$`\mathrm{\Theta }(t_3t_1)\mathrm{\Theta }(t_1t_2)e^{i\mathrm{\Omega }_{𝐤\sigma ,\mathrm{in};𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}}(t_3t_1)}e^{i\mathrm{\Omega }_{𝐤_{}+𝐪_{},\lambda ;𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}}(t_1t_2)}\left[f(E_{𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}})f(E_{𝐤_{}+𝐪_{},\lambda })\right]`$
$`+\mathrm{\Theta }(t_3t_2)\mathrm{\Theta }(t_2t_1)e^{i\mathrm{\Omega }_{𝐤\sigma ,\mathrm{in};𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}}(t_3t_2)}e^{i\mathrm{\Omega }_{𝐤\sigma ,\mathrm{in};𝐤_{}+𝐪_{},\lambda }(t_2t_1)}f(E_{𝐤_{}+𝐪_{},\lambda })).`$
Here, $`\mathrm{\Omega }_{𝐤l;𝐤^{}l^{}}(E_{𝐤l}E_{𝐤^{}l^{}})/\mathrm{}i\mathrm{\Gamma }_{𝐤l;𝐤^{}l^{}}`$ is a complex transition frequency. The four terms in Eq. (66) correspond to different time orders of interactions with the electric field. In itself, $`𝒩^{(3)}`$ is usually negligible compared to $`𝒩^{(2)}`$ for the following reason: The three frequencies of incoming photons have to add up to zero so that one has to be the negative of the sum of the other two. However, then the sum frequency is already present in the exciting laser pulse and ordinary photoemission dominates the signal.

Finally, the irreducible contribution to fourth order has the general form

$`𝒩_{\mathrm{irr}}^{2\mathrm{P}\mathrm{P}\mathrm{E}}(𝐤,\sigma )={\displaystyle \underset{\mathrm{𝐪𝐪}^{}𝐪^{\prime \prime }}{}}{\displaystyle 𝑑t_1𝑑t_2𝑑t_3𝑑t_4}`$ (67)
$`\times \eta _{ijkl}^{2\mathrm{P}\mathrm{P}\mathrm{E}}(𝐪,𝐪^{},𝐪^{\prime \prime };t_1,t_2,t_3,t_4;𝐤,\sigma )E_i(𝐪,t_1)E_j(𝐪^{},t_2)`$
$`\times E_k(𝐪^{\prime \prime },t_3)E_l(𝐪𝐪^{}𝐪^{\prime \prime },t_4).`$
$`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ can be found by inserting Eqs. (55) and (56) into Eq. (60). The photoelectron yield is determined by the off-diagonal components of $`\rho ^{(3)}`$, which in turn depend on *all* components of $`\rho ^{(2)}`$, including the diagonal ones. New physics enters here: The 2PPE current depends on both the polarization and the change of occupation to second order.

If the increase of the occupation due to secondary electrons is small to second order, we can use Eq. (55) to calculate all components of $`\rho ^{(2)}`$. The change of occupation due to dipole transitions and to relaxation out of excited states is included in Eq. (55). Then the response function $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ reads

$`\eta _{ijkl}^{2\mathrm{P}\mathrm{P}\mathrm{E}}(𝐪,𝐪^{},𝐪^{\prime \prime };t_1,t_2,t_3,t_4;𝐤,\sigma )={\displaystyle \frac{e^4}{\mathrm{}^4}}{\displaystyle \frac{\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}}{\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}}}{\displaystyle \underset{\lambda \lambda ^{}\lambda ^{\prime \prime }}{}}D_{𝐤\sigma ,\mathrm{in};𝐤_{}+𝐪_{},\lambda }^i(q_z)D_{𝐤_{}+𝐪_{},\lambda ;𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}}^j(q_z^{})`$ (68)
$`\times D_{𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{};𝐤_{}+𝐪_{}+𝐪_{}^{}+𝐪_{}^{\prime \prime },\lambda ^{\prime \prime }}^k(q_z^{\prime \prime })D_{𝐤_{}+𝐪_{}+𝐪_{}^{}+𝐪_{}^{\prime \prime },\lambda ^{\prime \prime };𝐤\sigma ,\mathrm{in}}^l(q_zq_z^{}q_z^{\prime \prime })`$
$`\times [F(t_1t_2,t_2t_3,t_3t_4;1,0;2,0;3,0)F(t_1t_2,t_2t_4,t_4t_3;1,0;2,0;2,3)`$
$`F(t_1t_4,t_4t_2,t_2t_3;1,0;1,3;2,3)+F(t_1t_4,t_4t_3,t_3t_2;1,0;1,3;1,2)`$
$`F(t_4t_1,t_1t_2,t_2t_3;0,3;1,3;2,3)+F(t_4t_1,t_1t_3,t_3t_2;0,3;1,3;1,2)`$
$`+F(t_4t_3,t_3t_1,t_1t_2;0,3;0,2;1,2)F(t_4t_3,t_3t_2,t_2t_1;0,3;0,2;0,1)]`$
with the auxilliary function

$`F(\mathrm{\Delta }t_1,\mathrm{\Delta }t_2,\mathrm{\Delta }t_3;n_1,n_2;n_3,n_4;n_5,n_6)\mathrm{\Theta }(\mathrm{\Delta }t_1)\mathrm{\Theta }(\mathrm{\Delta }t_2)\mathrm{\Theta }(\mathrm{\Delta }t_3)`$ (69)
$`\times e^{i\mathrm{\Omega }_{n_1,n_2}\mathrm{\Delta }t_1}e^{i\mathrm{\Omega }_{n_3,n_4}\mathrm{\Delta }t_2}e^{i\mathrm{\Omega }_{n_5,n_6}\mathrm{\Delta }t_3}\left[f(E_{n_6})f(E_{n_5})\right],`$
where the states $`|n_i`$ are defined as

$`|0|𝐤\sigma ;\mathrm{in},`$ $`|1|𝐤_{}+𝐪_{},\lambda ,`$
$`|2|𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{},`$ $`|3|𝐤_{}+𝐪_{}+𝐪_{}^{}+𝐪_{}^{\prime \prime },\lambda ^{\prime \prime }.`$ (70)
Thus the first term in the brackets in Eq. (68) reads

$`\mathrm{\Theta }(tt_1)\mathrm{\Theta }(t_1t_2)\mathrm{\Theta }(t_2t_3)e^{i\mathrm{\Omega }_{𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}(tt_1)}e^{i\mathrm{\Omega }_{𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{};𝐤\sigma ,\mathrm{in}}(t_1t_2)}`$ (71)
$`\times e^{i\mathrm{\Omega }_{𝐤_{}+𝐪_{}+𝐪_{}^{}+𝐪_{}^{\prime \prime },\lambda ^{\prime \prime };𝐤\sigma ,\mathrm{in}}(t_2t_3)}\left[\underset{=0}{\underset{}{f(E_{𝐤\sigma ;\mathrm{in}})}}f(E_{𝐤_{}+𝐪_{}+𝐪_{}^{}+𝐪_{}^{\prime \prime },\lambda ^{\prime \prime }})\right]`$
etc. $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ has the same structure as $`\eta ^{(3)}`$, only with more terms due to more possible time orders. One should exclude terms from Eq. (68) that correspond to processes for which the system returns to the equilibrium state after two of the four interactions. These processes only contribute a small correction to the numerical prefactor of the ordinary photoelectron yield.

If, in addition, the typical time scale of the experiment, e.g., the delay time in the pump-probe case, is long compared to the dephasing times, 2PPE can be described by the change of occupation alone. Interference effects are then absent. In this limit, only the terms $`F`$ in Eq. (68) with $`n_3=n_4`$ contribute, where $`n_3`$ is an excited state reachable by a single interaction out of the Fermi sea. Then $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ simplifies to

$`\eta _{ijkl}^{2\mathrm{P}\mathrm{P}\mathrm{E}}(𝐪,𝐪^{},𝐪^{\prime \prime };t_1,t_2,t_3,t_4;𝐤,\sigma )={\displaystyle \frac{e^4}{\mathrm{}^4}}{\displaystyle \frac{\gamma _{𝐤\sigma ,\mathrm{out};𝐤\sigma ,\mathrm{in}}}{\mathrm{\Gamma }_{𝐤\sigma ,\mathrm{in};𝐤\sigma ,\mathrm{in}}}}{\displaystyle \underset{\lambda \lambda ^{}}{}}`$ (72)
$`\times \delta _{𝐪_{}^{}+𝐪_{}^{\prime \prime },0}D_{𝐤\sigma ,\mathrm{in};𝐤_{}+𝐪_{},\lambda }^i(q_z)D_{𝐤_{}+𝐪_{},\lambda ;𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{}}^j(q_z^{})`$
$`\times D_{𝐤_{}+𝐪_{}+𝐪_{}^{},\lambda ^{};𝐤_{}+𝐪_{},\lambda }^k(q_z^{})D_{𝐤_{}+𝐪_{},\lambda ;𝐤\sigma ,\mathrm{in}}^l(q_z)`$
$`\times [F(t_1t_4,t_4t_2,t_2t_3;1,0;1,1;2,1)`$
$`+F(t_1t_4,t_4t_3,t_3t_2;1,0;1,1;1,2)`$
$`F(t_4t_1,t_1t_2,t_2t_3;0,1;1,1;2,1)`$
$`+F(t_4t_1,t_1t_3,t_3t_2;0,1;1,1;1,2)].`$
Note, $`\eta ^{2\mathrm{P}\mathrm{P}\mathrm{E}}`$ is now proportional to $`\mathrm{exp}(\mathrm{\Gamma }_{11}\mathrm{\Delta }t)`$ from the second exponential in Eq. (69), where $`\mathrm{\Gamma }_{11}=\tau _1^1`$ is the energy relaxation rate of state $`|1=|𝐤_{}+𝐪_{},\lambda `$ and $`\mathrm{\Delta }t`$ is the time between the second and third interaction. This is easy to understand: After the second interaction the electron is in the pure state $`|1`$, which decays with the rate $`\mathrm{\Gamma }_{11}`$.

Reducible contributions to the photoelectron current result from nonlinear optical effects in the solid. They are obtained by replacing the effective electric field $`𝐄`$ in Eqs. (63) and (65) by the electric field to second order, see Eq. (13). For 2PPE we obtain the contributions given in Eqs. (20) and (21) and shown diagrammatically in Fig. 7(b) and (c). There are also contributions given by $`\eta `$ times a product of $`𝐄`$ and the third-order polarization $`𝐏^{(3)}`$. These are only significant if the incident light contains a frequency component large enough to allow ordinary photoemission and we neglect them.

Finally, it is also possible to describe (ordinary) photoemission out of a general non-equilibrium state described by the density matrix $`\rho _{\mathrm{neq}}`$. The irreducible contribution is obtained from Eq. (60) for the photoelectron yield $`𝒩^{(n)}`$ by expressing $`\rho _{\mathrm{irr}}^{(n1)}`$ in terms of the lower-order $`\rho _{\mathrm{irr}}^{(n2)}`$ with the help of Eq. (55). This is exact, since only off-diagonal components of $`\rho _{\mathrm{irr}}^{(n1)}`$ are needed. Then $`\rho ^{(n2)}`$ is replaced by $`\rho _{\mathrm{neq}}`$. There is also a reducible part: Two photons can be converted into a single one at the sum frequency, which for non-equilibrium, when a finite polarization exists, can lead to a change of occupation of states above $`E_{\mathrm{vac}}`$ and thus to photoemission.