## 2 Methods  
We describe all steps of the photoemission process from excitation to detection by the TDSE in one dimension and in a single-active-electron picture,  
\[i\hbar\frac{\partial\psi(z,t)}{\partial t}=\left(-\frac{\hbar^{2}}{2m_{\rm e}} \frac{\partial^{2}}{\partial z^{2}}+V_{0}(z)+V_{\rm IR}(z,t)+V_{\rm EUV}(p,t) \right)\psi(z,t), \tag{3}\]  
where \(V_{0}(z)\) is the time-independent effective one-electron potential, and \(V_{\rm IR}(z,t)\) and \(V_{\rm EUV}(p,t)\) are the IR and EUV interaction potentials in length and velocity gauge, respectively. All calculations shown here are based on an approach introduced by [PERSON] and [PERSON] to describe atomic photoionization [33], as well as by [PERSON] and [PERSON], which was used to describe photoemission dynamics from a metal surface [16]. Note that the present implementation of the 1D model does not account for angular momentum effects and thus cannot be used to describe atomic-photoemission processes. The central idea of this approach is to express the wave function as a sum of initial state and final state to separate fast phase oscillations. In the following, it is assumed that the IR field vanishes inside the solid and that EUV interaction is such a small perturbation that depopulation of the initial and final states by EUV excitation can be neglected. Then, the initial state wave function can be treated as time-independent, and the full wave function can be expressed as:  
\[\psi(z,t)=\Phi_{0}(z)+e^{-i(\frac{\Phi_{0}}{\hbar}+\omega)t}\Phi(z,t), \tag{4}\]  
where \(\Phi_{0}(z)\) is the initial state wave function that fulfills  
\[\left(-\frac{\hbar^{2}}{2m_{\rm e}}\frac{\partial^{2}}{\partial z^{2}}+V_{0}( z)\right)\Phi_{0}(z)=\varepsilon_{0}\Phi_{0}(z). \tag{5}\]  
Photoelectron wave function \(\Phi(z,t)\) is zero before the EUV pulse arrives. The rapidly oscillating exponential is separated for computational efficiency. Here, \(\omega\) is the center frequency of EUV light. Inserting Equation (4) into Equation (3), and using Equation (5) as well as the rotating-wave approximation, results in a differential equation for photoelectron wave function [16; 10]:  
\[i\hbar\frac{\partial\Phi(z,t)}{\partial t}=\left(-\frac{\hbar^{2}}{2m_{\rm e} }\frac{\partial^{2}}{\partial z^{2}}+V_{0}(z)+V_{\rm IR}(z,t)-(\hbar\omega+ \varepsilon_{0})\right)\Phi(z,t)+i\frac{\hbar e}{2m_{\rm e}}A_{\rm EUV}^{0}( t)\frac{\partial\Phi_{0}(z)}{\partial z}, \tag{6}\]  
where \(A_{\rm EUV}^{0}(t)\) is the envelope of the EUV pulse vector potential. Note that this approach is nonperturbative with respect to solid potential and IR potential, i.e., it does not apply strong field approximation, but describes the laser field and the solid's effective potential on equal footing. The unperturbed solid is described by the sum of jellium-type potential and electron-hole interaction. Jellium potential is given by:  
\[V_{\rm JP}(z)=\frac{V_{\rm JP}}{2}\left(\tan\hbar\left(\frac{z-z_{\rm s}}{w_{ \rm s}}\right)-1\right), \tag{7}\]  
where \(z_{\rm s}\) is the surface position, \(w_{\rm s}\) is the width of the transition region, and \(V_{\rm JP}\) is the depth of the potential well. The screened electron-hole interaction is described by a one-dimensional Yukawa potential:  
\[V_{\rm h}(z)=-e\frac{\exp\left(-\frac{|z-z_{0}|}{\mu}\right)}{\sqrt{\left(z-z_ {0}\right)^{2}+a^{2}}}. \tag{8}\]Here, \(z_{0}\) is the effective initial position, \(\upmu\) is the screening length, and \(a\) is a regularization parameter [16]. The spatial distribution of IR potential is assumed as a smooth but narrow step at the surface because modeling of IR field penetration is beyond the scope of this article. Damping of the photoelectron wave function due to inelastic scattering processes is neglected here. Instead, damping enters the effective initial position via the inelastic mean free path. This approach can be used to calculate streaking spectra as well as RABBITT spectra.