Simulation of attosecond transient soft X-ray absorption in solids using generalized Kohn-Sham real-time TDDFT

The velocity-gauge (VG) formulation of real-time TDDFT (VG-RT-TDDFT) which is well suited for the treatment of laser-driven electron dynamics in Bloch-periodic systems has been previously discussed by several authors [22, 29, 42, 30]. In particular, the numerical approach employed in this work, which is based on a linear-combination of atomic-orbitals (LCAO) implementation of the generalized Kohn-Sham (GKS) [43] real-time TDDFT framework has been described in detail elsewhere [30, 31]. Within this approach the time-dependent generalized Kohn-Sham (TDGKS) equations for electron dynamics

$\displaystyle\imath\hbar\frac{\partial}{\partial t}{\phi}_{i}(\overrightarrow{r},t)=\hat{{H}}_{GKS}{\phi}_{i}(\overrightarrow{r},t)$

(1)

with ${\phi}_{i}$ being single-particle orbitals and $\hat{{H}}_{GKS}$ the VG GKS Hamiltonian, are propagated in real-time. The VG GKS Hamiltonian is given by

	$\displaystyle\hat{{H}}_{GKS}=$	$\displaystyle\frac{1}{2m}\left[\overrightarrow{\hat{p}}+\frac{e}{c}\overrightarrow{A}(t)\right]^{2}+\hat{\tilde{V}}_{ion}+$			(2)
		$\displaystyle\int d\overrightarrow{r}^{\prime}\frac{e^{2}}{|\overrightarrow{r}-\overrightarrow{r}^{\prime}|}n(\overrightarrow{r}^{\prime},t)+\hat{V}_{XC}[\rho(\overrightarrow{r},\overrightarrow{r}^{\prime},t)]$

and features a kinetic term that includes coupling to the time-dependent vector potential $\overrightarrow{A}(t)$ representing the applied external fields, the VG electron-nuclear interaction term $\hat{\tilde{V}}_{ion}$, the Hartree potential and a non-multiplicative XC operator $\hat{V}_{XC}$ as a functional of the instantaneous single-particle GKS density matrix $\rho(\overrightarrow{r},\overrightarrow{r}^{\prime},t)$ [43]. Time propagation of equation 1 yields at every time-step, relevant quantities such as the GKS 1-body density matrix

$$\rho(\overrightarrow{r},\overrightarrow{r}^{\prime},t)=\sum_{i}^{N_{occ}}{\phi}_{i}(\overrightarrow{r},t){\phi}_{i}^{*}(\overrightarrow{r}^{\prime},t)$$

electron density,

$$n(\overrightarrow{r},t)=\rho(\overrightarrow{r},\overrightarrow{r},t)=\sum_{i}^{N_{occ}}|\tilde{\psi}_{i}(\overrightarrow{r},t)|^{2}$$

and macroscopic current

$$\overrightarrow{I}(t)=-\frac{e}{\Omega}\int_{\Omega}d\overrightarrow{r}\overrightarrow{j}(\overrightarrow{r},t).$$

(3)

which in turn is obtained from the time-dependent current density $\overrightarrow{j}(\overrightarrow{r},t)$ given as

$$\overrightarrow{j}(\overrightarrow{r},t)=\sum_{i}\frac{e}{2m}\left\{{\phi}^{*}_{i}(\overrightarrow{r},t)\overrightarrow{\pi}{\phi}_{i}(\overrightarrow{r},t)+c.c\right\}$$

(4)

and includes the generalized momentum

$$\overrightarrow{\pi}=\frac{m}{\imath\hbar}[\overrightarrow{r},\hat{{H}}_{GKS}].$$

(5)

Following numerical time propagation of the TD GKS equation 1, frequency domain observables related to the time-dependent electron density or current are readily calculated via Fourier transformation[20, 28, 29, 30].

As outlined in the introduction, one of the aims of this work is to investigate in the context of extended systems, the use of a range-separated hybrid (RSH) exchange-correlation (XC) functional form that provides improved accuracy compared to semi-local approximations with regards to the treatment of both core and valence electrons simultaneously. To this end an RSH functional which combines both long-range (LR) and short-range (SR) corrections is employed according to the following strategy: For an adequate treatment of valence excitations including excitonic effects, LR corrections along the lines proposed by Refaely-Abramson et al [38] and previously demonstrated in the case of monolayer h-BN [31] are employed. Additional SR corrections which primarily affect core-electron energies are then applied based on satisfying an approximate Koopman condition [44] where by core-level GKS eigenenergies are required to match core-electron removal energies from $\Delta$ self-consistent field ($\Delta$SCF) [45] calculations. Such a procedure is analogous in the case of finite systems to enforcing the linearity condition for orbital energies (LCOE) as outlined previously by a number of authors [37, 44, 46].

To incorporate SR and LR corrections, the following partitioning of the Coulomb operator is employed featuring the error- and complementary-error functions:

	$\displaystyle\frac{1}{|\overrightarrow{r}-\overrightarrow{r}^{\prime}|}=$	$\displaystyle\frac{\alpha+\beta_{LR}~{}\mathrm{erf}(\omega_{LR}|\overrightarrow{r}-\overrightarrow{r}^{\prime}|)+\beta_{SR}~{}\mathrm{erfc}(\omega_{SR}|\overrightarrow{r}-\overrightarrow{r}^{\prime}|)}{|\overrightarrow{r}-\overrightarrow{r}^{\prime}|}+$			(6)
		$\displaystyle\frac{1-\{\alpha+\beta_{LR}~{}\mathrm{erf}(\omega_{LR}|\overrightarrow{r}-\overrightarrow{r}^{\prime}|)\}+\beta_{SR}~{}\mathrm{erfc}(\omega_{SR}|\overrightarrow{r}-\overrightarrow{r}^{\prime}|)}{|\overrightarrow{r}-\overrightarrow{r}^{\prime}|}$

The first and second terms on the right hand side respectively are then used to evaluate the nonlocal Hartree-Fock (HFX) like and local density approximation (LDAX) like contributions to the exchange energy. This leads to an effective partitioning of the total XC energy as

	$\displaystyle E_{\mathrm{XC}}=$	$\displaystyle(\alpha+\beta_{SR})E^{\mathrm{SR}}_{\mathrm{HFX}}+(1-\alpha-\beta_{SR})E^{\mathrm{SR}}_{\mathrm{LDAX}}+$			(7)
		$\displaystyle(\alpha+\beta_{LR})E^{\mathrm{LR}}_{\mathrm{HFX}}+(1-\alpha-\beta_{LR})E^{\mathrm{LR}}_{\mathrm{LDAX}}+E_{\mathrm{LDAC}}$

where the $E^{\mathrm{SR}}_{\mathrm{HFX}},E^{\mathrm{LR}}_{\mathrm{HFX}}$ respectively stand for the short-range (SR) and long-range (LR) Hartree-Fock exchange (HFX) contributions with similar notation for the LDA counterparts. The correlation contribution to the total energy $E_{\mathrm{LDAC}}$ is set to be that of the LDA [47]. Thus the functional form adopted features fractions $(\alpha+\beta_{SR}),(\alpha+\beta_{LR})$ of HFX in the SR and LR respectively. For brevity, the functional form implied by eqns. 6, 7 is henceforth refereed to as the short- and long-range corrected (SLRC) form. The SR and LR contributions to LDAX are calculated using the Coulomb attenuated form due to Toulouse et al [48]. The parameters $\alpha,\beta_{LR},\omega_{LR}$ are chosen such that

$\displaystyle\alpha+\beta_{LR}=\frac{1}{\epsilon_{\infty}}$

(8)

which ensures the correct long-range $\frac{1}{\epsilon_{\infty}r}$ behavior of the XC potential in the $r\rightarrow\infty$ asymptotic limit while tuning $\alpha,\omega_{LR}$ so that the GKS band-gap of the extended system matches the quasiparticle band-gap obtained from benchmark GW calculations [32]

The choice of $\beta_{SR},\omega_{SR}$ on the other hand predominantly determines the core-level GKS eigenenergies. In principle both $\beta_{SR}$ and $\omega_{SR}$ can be tuned separately so that the following condition is satisfied for a chosen set of core-levels:

$\displaystyle\varepsilon^{A}_{c}=-(E^{A+}_{c}-E^{GS})$

(9)

where $\varepsilon^{A}_{c}$ is the GKS eigenvalue of the core-level $c$ on atom $A$, $E^{GS}$ is the total energy of the neutral ground state of the system of interest and $E^{A+}_{c}$ is the total energy of the positively ionized system resulting from the removal of one electron form the core-orbital $c$ localized on atom A. Thus $(E^{A+}_{c}-E^{GS})$ represents a core-electron removal energy including final-state electronic screening effects which can be accessed from a constrained-DFT $\Delta$SCF [50] or a quasiparticle GW calculation [51]. In this work, $(E^{A+}_{c}-E^{GS})$ is calculated as a LDA total-energy difference for 2D periodic h-BN using a supercell $\Delta$SCF approach including 2D image charge corrections to $E^{A+}_{c}$ from finite size scaling [52, 53](see Fig. 1). Because of the localized nature of the core hole and associated short-range screening response (see Fig.  1) supercell total energy methods for charged defects can be used to approximate $E^{A+}_{c}$ [53, 54]. Since in h-BN, the B $1s$ and N $1s$ orbitals represent the core-states of interest, the removal energy $(E^{A+}_{c}-E^{GS})$ is calculated separately for each core hole case .i.e., $c=1s$ for $A=$B, N (see Fig. 1 and Table  1).

The parameters on the right hand side of eqn. 6 are then determined as follows: In monolayer h-BN, the long-range asymptotic limit of the macroscopic dielectric constant $\epsilon_{\infty}^{||}\rightarrow 1$ [56] dictates a choice of $\alpha+\beta_{LR}=1$ from eqn. 8 above. Physical considerations similarly suggest unscreened Fock exchange in the very short-range [57] which requires $\alpha+\beta_{SR}=1$. While $\alpha$ is in principle a free parameter, in this work we set $\alpha=0$ to be consistent with a previously reported [31] valence-electron LRC functional parameterization for h-BN without explicit treatment of core-electrons. Thus the effect of including core-electrons in the description is intended to be captured by $\beta_{SR}$ and the range-separation parameters $\omega_{SR},\omega_{LR}$. The choice of $\alpha=0$ fixes both $\beta_{LR}=1$, $\beta_{SR}=1$ leaving $\omega_{SR},\omega_{LR}$ as the two parameters to be tuned in eqn. 6. The two core-electron energy eigenvalues $\varepsilon^{B}_{1s}$, $\varepsilon^{N}_{1s}$ and the $\Gamma\rightarrow\Gamma$ GKS band-gap $\Delta_{G}^{\Gamma}$ in 2D h-BN represent three quantities to be fit to the corresponding quasiparticle energy differences as outlined above (see eqn. 9) by minimizing the percentage error in a least squares sense. The resulting values for $\varepsilon^{B}_{1s}$, $\varepsilon^{N}_{1s}$, $\Delta_{G}^{\Gamma}$ and their comparison to target data are shown in Table 1 for the optimized parameter set $\alpha=0$, $\beta_{LR}=1$, $\omega_{LR}=0.129~{}\mathrm{Bohr^{-1}}$, $\beta_{LR}=1$, $\omega_{SR}=1.40~{}\mathrm{Bohr^{-1}}$.

The variation in the effective fraction of Fock exchange as a function of inter-electron separation for the tuned SLRC functional is plotted in Figure 2. The corresponding GKS band structure and independent-particle approximation (IPA) optical response of 2D h-BN are also shown. Thus the SLRC functional yields a description in which both core and valence electron GKS eigenvalues approximate single-particle addition/removal energies of the system. In the next section RT-TDDFT results for both equilibrium and non-equilibrium optical response employing this SLRC functional are discussed.

Numerical simulations in this work are carried out using a linear-combination of atomic-orbital (LCAO) framework for Velocity-gauge real-time TDDFT (VG-RT-TDDFT) [30, 31] based on a modification of the SIESTA [58] density functional platform. Accordingly, for 2D h-BN, electron-ion interactions are modeled through norm-conserving nonrelativistic LDA pseudopotentials generated using the Troullier-Martins [59] scheme while the GKS Hamiltonian, density matrix and wavefunctions are represented in an LCAO basis. Non-local Fock exchange integrals are evaluated using an auxiliary Gaussian basis set representation as described previously [31]. Pseudopotential and basis set parameters are listed in Table 2. Brillouin zone sampling is carried out using $\Gamma$-centered $\mathbf{k}$-point grids also as listed in Table 2. Time propagation of the GKS [43, 31] equations is carried out using a Crank-Nicholson scheme [60] with a predictor-corrector step.

To evaluate the tuned SLRC functional, firstly the response of 2D h-BN to a weak Dirac-$\delta$ like time-localized electric field perturbation is simulated and fluctuations of the induced macroscopic current are calculated in real-time as plotted in Figure 3(a). The corresponding frequency domain dielectric response is obtained through Fourier transformation of the time-dependent macroscopic current [29]. The low and high frequency fluctuations apparent in the time-evolution of the current respectively correspond to the valence- and core-electronic optical response. The imaginary part $\epsilon_{2}$ of the dielectric function is plotted in figure 3 over energy ranges relevant to valence as well as B, N $1s$ (K-edge) core-excitations in h-BN. For the sake of comparison, real-time adiabatic LDA (ALDA) results are also shown. As has been extensively reviewed previously [40, 63], Kohn-Sham ALDA does not capture excitonic features in the optical response either in the valence or core-excitation case [31]. In ALDA the energy onset of UV absorption in 2D h-BN is underestimated by over 1 eV and the strong excitonic enhancement of absorption at 5.58 eV predicted by benchmark BSE simulations [32] is not observed. In the absence of explicit self-interaction corrections [47], the ALDA XC potential asymptotically decays exponentially which translates to an over-screening of the Coulomb interactions necessary to stabilize excitons. In contrast because of the incorporation of the correct $\frac{1}{\epsilon_{\infty}r}$ long-range asymptotic behavior of the XC potential via the screened exchange component of the nonlocal SLRC functional, the RT-SLRC time dynamics reproduces the strong valence excitonic features expected in 2D h-BN to within 0.2 eV of BSE results. In the case of X-ray excitations at the B and N K-edges in the hundreds of eV range, ALDA once again underestimates the onset of the absorption edge by over $9\%$ and $6\%$ respectively while also missing prominent core-excitonic features leading to an overall spectral line shape that bears little resemblance to experimental spectra. The short-range component of the nonlocal Fock-like exchange incorporated into the SLRC functional overcomes these shortcomings. Firstly, it improves absolute core-excitation energies by mitigating delocalization errors on the length-scale of localized core electrons towards restoring the piece-wise linearity condition on core-orbital energies [37]. Concomitantly, because of reduced screening of short-range electron-hole interactions relative to ALDA, the density evolution under the SLRC XC potential encodes core-excitonic effects. As apparent from Fig. 3(b), 3(c) and table 3, with the RT-SLRC approach absolute energies of the B, N K-edge positions are reproduced to within 2$\%$ of experiment. In particular the excitonic $1s\rightarrow\pi^{*}$ transitions at the B and N K-edges are correctly described while overall spectral features in a 20 eV energy range above the onset are also in reasonable agreement with experiment considering that the frequency-dependence of lifetime broadening effects is not included within the adiabatic XC approximation used here. Thus with time-propagation of the GKS equations for a delta perturbation, the tuned SLRC XC potential employed in this study is able to provide a balanced description of electron-hole excitations across a wide energy range in strongly excitonic 2D h-BN.

Attosecond Transient Absorption Spectroscopy (ATAS)  [9, 7, 64, 65, 23, 27] typically involves a pump-probe protocol wherein the absorption of an ultrashort probe pulse in a material excited by a pump pulse is measured. The intensity, duration and shape of the pump pulse as well as the pump-probe delay provide external control parameters that can be tuned. The simulation, within a TDDFT framework, of ATAS employing low energy IR-XUV probe pulses in molecules [23, 25] and solids [66, 27] has been discussed in the literature. Soft X-ray transient absorption (XTAS) [9, 7, 65] is similar in principle except that higher frequency (100 - 1000 eV) wide band-width X-ray probe pulses tuned to core-excitation edges are typically employed yielding an element specific perspective of transient electron dynamics with attosecond time-resolution Within the real-time velocity gauge formulation of TDDFT [29], one has access to the time-dependent macroscopic current $\overrightarrow{I}(t)$ (Eqn. 3) within a material that is under the influence of external laser fields described by the vector potential

$$\overrightarrow{A}(t)=-c\int^{t}\overrightarrow{E}(t^{\prime})dt^{\prime}$$

(10)

obtained from an electric field $\overrightarrow{E}(t)$. In the frequency domain, the macroscopic current $\overrightarrow{I}(\omega)$ associated with a weak probe pulse considered to be in the linear response regime and described by a given external perturbation $\overrightarrow{E}(\omega)$ can be written as

$$\overrightarrow{I}(\omega)=\overline{\overline{\sigma}}(\omega).\overrightarrow{E}(\omega)$$

(11)

where $\overline{\overline{\sigma}}$ is the conductivity tensor.

In a pump-probe scheme for a certain choice of pump and probe fields described by a field profile $\overrightarrow{E}_{pump\mh probe}(t)$ one can calculate through propagating the GKS equation 1 the associated macroscopic current $\overrightarrow{I}_{pump\mh probe}(t)$ using equation. 3. Similarly, for the pump field $\overrightarrow{E}_{pump}(t)$ alone, once can calculate the corresponding current $\overrightarrow{I}_{pump}(t)$. Then the probe current in the pump-probe scheme is defined as

$$\overrightarrow{I}_{probe}(t)\equiv\overrightarrow{I}_{pump\mh probe}(t)-\overrightarrow{I}_{pump}(t)$$

(12)

This separation of the total current into pump and probe parts is sensible for probes of any intensity but $\overrightarrow{I}_{probe}(t)$ has a dependence on the parameters describing the pump pulse such as its shape, intensity, duration, polarization etc indicated collectively by the set $\mathcal{P}$ as well as the pump-probe delay $\Delta T$. Additionally, for weak probe pulses, equation 11 can be generalized to

$$\overrightarrow{I}_{probe}(\omega)=\overline{\overline{\sigma}}_{\mathcal{P}}(\omega,\Delta T).\overrightarrow{E}_{probe}(\omega)$$

(13)

in terms of the transient frequency dependent conductivity tensor $\overline{\overline{\sigma}}_{\mathcal{P}}(\omega,\Delta T)$ which now depends on the pump parameters $\mathcal{P}$ and the delay $\Delta T$. For comparison to XTAS experiments, it is useful to calculate the transient dielectric tensor

$$\overline{\overline{\epsilon}}_{\mathcal{P}}(\omega,\Delta T)=1+\frac{4\pi\imath\overline{\overline{\sigma}}_{\mathcal{P}}(\omega,\Delta T)}{\omega}$$

(14)

whose imaginary part is related to X-ray absorption.

To illustrate the simulation of XTAS using the RT-SLRC hybrid functional approach within GKS-TDDFT, the pump-probe set up shown in Fig. 4 is considered in this study. Accordingly a 5 femtosecond (fs) pump pulse with a a $Sin^{2}$ envelope function and a central frequency of 5.67 eV tuned to the first valence exciton in 2D h-BN is employed. Because of the finite time duration of the pulse, its frequency band-width is roughly 2 eV. Different intensities of this pulse impinging on 2D h-BN with its electric field polarized in-plane are simulated. The chosen intensity range of $\sim 10^{12}$ W/cm² corresponds to that typically employed [64] in attosecond pump-probe experiments of electron dynamics. The pump pulse acts on the sample to excite electron-hole pairs at a density which is roughly 2% of the valence electron density in 2D h-BN. Concomitantly, it also deposits on the order of 1 eV per unitcell of energy during the course of its interaction with the material. Immediately following the passage of the pump pulse, the material is probed using a weak delta perturbation. While a weak probe of any shape could be employed, a delta perturbation corresponds to the infinite band-width limit which is convenient in this context as it provides a probe-shape independent view of the transient dielectric function. The electric field of the probe is polarized along the diagonal of the plane perpendicular to the in-plane pump electric field. The probe therefore couples to both $1s\rightarrow\pi^{*}$ and $1s\rightarrow\sigma^{*}$ excitations. With this set up, the transient dielectric functions for different pump intensities are calculated and its imaginary part ($\mathrm{\varepsilon_{2}}$) is plotted in Fig. 5 for both the B and N K-edges.

By comparing the transient dielectric function (tr-$\mathrm{\varepsilon_{2}}$) after the pump with its ground state counterpart (gs-$\mathrm{\varepsilon_{2}}$) calculated with no pump, it is apparent that tr-$\mathrm{\varepsilon_{2}}$ is modulated in several ways at the B and N K-edges. Firstly, it is seen that relative to gs-$\mathrm{\varepsilon_{2}}$ the $1s\rightarrow\pi^{*}$ peak in tr-$\mathrm{\varepsilon_{2}}$ at both the B (188.1 eV) and N (403.3 eV) K-edge exhibits a reduction in intensity or bleaching that increases with increasing pump fluence. This effect can be understood as an excitation induced Pauli-blocking [67, 68, 64, 9]. In 2D h-BN the valence band maximum (VBM) and conduction band minimum (CBM) are primarily of N-2$p$ orbital derived $\pi$ and B-2$p$ orbital derived $\pi^{*}$ characters respectively (see Fig 6). Valence excitations in the energy range of the pump pulse therefore exhibit $\pi\rightarrow\pi^{*}$ character with the electronic part of the electron-hole wavefunction being predominantly $\pi^{*}$ like [32]. This is equivalent within a projected single-particle picture to a transient filling of the $\pi^{*}$ CBM states [29] following laser excitation by the pump. As the probe pulse that immediately follows the pump drives $1s\rightarrow\pi^{*}$ transitions, a certain fraction of these are blocked in proportion to the transient occupation of the $\pi^{*}$ CBM states. Thus the intensity of the $1s\rightarrow\pi^{*}$ peak in tr-$\mathrm{\varepsilon_{2}}$ is reduced. Note that the higher lying $1s\rightarrow\sigma^{*}$ peaks occurring near $\sim$194.8 eV and 410.1 eV at the B and N K-edges respectively, do not show such a reduction in intensity as the pump excitation does not lead to significant occupation of the conduction band $\sigma^{*}$ states. Meanwhile, the pump also transiently induces a hole population near the VBM which is of $\pi$ character. Hence from the perspective of the probe pulse certain $1s\rightarrow\pi$ transitions of core electrons to VBM hole states become allowed. These would otherwise be blocked in the ground state when the valence band is fully filled. This kind of Pauli-unblocking of core-excitations is readily apparent as a new $\sim$ 397.5 eV peak observed in the N K-edge tr-$\mathrm{\varepsilon_{2}}$ spectrum which is missing in gs-$\mathrm{\varepsilon_{2}}$. A similar feature is not observed at the B K-edge however because the B $1s\rightarrow\pi$ oscillator strength is much weaker relative to the B $1s\rightarrow\pi^{*}$ strength given that the VBM has very little B-2$p$ character. This difference between the B and N K-edges shows the orbital selectivity of features in tr-$\mathrm{\varepsilon_{2}}$.

In addition to the intensity modulations discussed in the previous paragraph, tr-$\mathrm{\varepsilon_{2}}$ peaks also exhibit certain energy shifts on the order of 0.1-0.4 eV. In particular we note that the B $1s\rightarrow\pi^{*}$ peak at 188.1 eV is transiently red-shifted (see Fig 5) while the N $1s\rightarrow\pi^{*}$ peak at 403.3 eV is blue-shifted. Higher lying $1s\rightarrow\sigma^{*}$ peaks show similar trends being transiently red- and blue-shifted at the B and N K-edges respectively. These can be understood in terms of valence excitation induced core-level energy shifts. Note that in h-BN, because of the predominant N-2$p$ and B-2$p$ orbital characters of the VBM and CBM respectively, an across the band-gap excitation has a charge-transfer like character whereby N sites are incrementally oxidized and B sites are reduced. The amount of fractional oxidation or reduction scales with the density of electron-hole pairs created by the pump. Localized core-electronic levels on atoms are sensitive to the local oxidation state because of its effect on electronic screening of the attractive nuclear potential. Thus core-electron energy-levels on transiently oxidized species (N) see their binding energies increase while those on transiently reduced species (B) find themselves at lower binding energies relative to the vacuum level (see Fig. 6). This behavior is exploited in time-resolved X-ray photoemission spectroscopy [69] that directly tracks transient changes in core-electron binding energies. In XTAS simulated with RT-TDDFT we do not directly calculate the transient core-electron single particle energy shifts but nevertheless the signature of these binding energy shifts is observed in the probe induced core to valence dipole transition energies. Accordingly, pump intensity dependent peak shifts that modulate tr-$\mathrm{\varepsilon_{2}}$ are seen in Fig 5.

While the XTAS spectra in Fig 5 are analyzed qualitatively in terms of pump excitation induced state blocking and core-level shifts, in general depending upon the excitation regime, transient particle-hole response in semiconductors and insulators incorporates additional effects such as band gap renormalization [70, 71] and nonequilibrium screening induced exciton energy shifts [64, 72]. For the specific mode of excitation and the low 1-2% particle-hole density considered here, these effects, to the extent that they are modeled within adiabatic GKS theory, do not seem prominent. The decomposition of the full electron-hole transient response available from XTAS into separate one-particle and two-particle contributions [64] can in principle be carried out within a GKS RT-TDDFT framework by also simulating separately, time-resolved core and valence [73] photoemission spectra. Such simulations could enable a comparison of GKS RT-TDDFT to nonequilibrium quasi-particle theories [74, 71] and will be a subject of future studies. Nevertheless, the lack of explicit memory dependence in the XC potentials studied here also imposes certain limitations with regards to the description of excited state energy relaxation and decoherence [75, 21, 20]. Additionally, it should be noted that for time-overlapped pump and probe pulses transient spectra exhibit light-field induced modulations such as the dynamical Franz-Keldysh [66, 26] and Stark [76] effects among others. Although not considered here as the pump and probe pulses used do not overlap in time, these types of nonlinear optical phenomena are in principle directly accessible from both KS and GKS RT-TDDFT subject as always to the accuracy limits of specific XC approximations employed.