Simulation of attosecond transient soft x-ray absorption in solids using generalized Kohn–Sham real-time time-dependent density functional theory

Time-dependent density functional theory (TDDFT) simulations of transient core-level spectroscopies require a balanced treatment of both valence- and core-electron excitations. To this end, tuned range-separated hybrid exchange–correlation functionals within the generalized Kohn–Sham scheme offer a computationally efficient means of simultaneously improving the accuracy of valence and core excitation energies in TDDFT by mitigating delocalization errors across multiple length-scales. In this work range-separated hybrid functionals are employed in conjunction with the velocity-gauge formulation of real-time TDDFT to simulate static as well as transient soft x-ray near-edge absorption spectra in a prototypical solid-state system, monolayer hexagonal boron nitride, where excitonic effects are important. In the static case, computed soft x-ray absorption edge energies and line shapes are seen to be in good agreement with experiment. Following laser excitation by a pump pulse, soft x-ray probe spectra are shown to exhibit characteristic features of population induced bleaching and transient energy shifts of exciton peaks. The methods outlined in this work therefore illustrate a practical means for simulating attosecond time-resolved core-level spectra in solids within a TDDFT framework.


Introduction
Over the past twenty years, rapid advances have been realized in the ability to generate extremely short X-ray laser pulses in a variety of settings ranging from table-top highharmonic generation (HHG) [1,2] setups to large-scale synchrotron and X-ray freeelectron laser (XFEL) [3,4] facilities. This in turn has led to the development of a wide variety of new time-domain X-ray spectroscopic techniques capable of investigating electron and ion dynamics in matter with picosecond to attosecond time-resolution [5,6,7]. Due to their characteristic element-specificity and unprecedented temporal resolution, ultrafast extreme-UV (XUV) and X-ray spectroscopies are emerging as indispensable tools in chemistry and materials science [8,7]. In particular, X-ray transient absorption spectroscopy (XTAS) which deploys short X-ray pulses to probe ultrafast modulations in the X-ray absorption near-edge structure (XANES) is now increasingly employed to investigate electron dynamics in molecular and solid-state systems [6,7,9]. In the case of static X-ray absorption spectroscopy (XAS) [10] which is widely used for characterizing the ground state electronic structure of materials, the interpretation of experimentally measured spectra has been greatly facilitated over the years by the development of theoretical tools for atomistic simulation of XANES using quantum mechanical methods [11]. A wide variety of first-principles approaches now exist for static XANES simulations and have been reviewed in the literature [10,11,12,13,14,15,16,17,18]. The advent of XTAS similarly requires the concomitant development of theoretical techniques and tools to aid the interpretation of time-resolved X-ray spectra that probe the dynamics of quantum mechanical excited states of matter [9]. In this context, the framework of time-dependent density functional theory (TDDFT) [19,20,21] which has long been a workhorse method for excited-state simulations in quantum chemistry represents an attractive option for its potential to balance computational feasibility and accuracy. For simulating the response of materials to intense ultra-short laser pulses relevant to pump-probe experiments in particular, real-time TDDFT (RT-TDDFT) [20,21,22] offers a convenient formalism. Hence in recent years, RT-TDDFT has been deployed for modeling femto-to attosecond timescale transient absorption in both molecules and solids for energies in the infrared (IR) to XUV range [23,24,25,26,27]. To date however, RT-TDDFT simulations of XTAS involving core-excitations in the soft X-ray energy range have not been reported especially in extended systems. In this work therefore, the velocity-gauge formulation of RT-TDDFT (VG-RT-TDDFT) [28,29,30,31] is employed to simulate XTAS in monolayer hexagonal boron nitride (h-BN) which is a prototypical 2D-periodic insulator characterized by strong excitonic effects [32,31].
XTAS pump-probe experiments [9,7] typically involve excitation by photon energies on the order of ∼1 eV from intense IR-UV pump laser pulses to drive valence electron dynamics in a sample while XUV/X-ray probe pulses investigate the pump-induced excited-state dynamics of valence electrons by recording time-dependent modulations of core-electron excitation features in the ∼10-1000 eV range. Therefore, a balanced theoretical description of XTAS needs among other things, to provide a reliable treatment of both valence and core electronic excitations simultaneously. This necessarily translates in a TDDFT context, to a need for choosing appropriate exchange-correlation (XC) potentials [20,21] which primarily dictate the accuracy of the simulations. As noted by Zhang et al [33], similar considerations apply in emerging nonlinear X-ray spectroscopy protocols where nonlinear excitation of coreelectrons by intense coherent X-ray pulses can be used to in turn drive valence electron dynamics. While the choice of appropriate XC functionals is generally an open problem in TDDFT, several useful guidelines have been proposed in the literature over the past two decades [34,13,35,14,36,37,38,39,40]. Firstly, in the case of valence excitations, local, semi-local and short-range nonlocal hybrid functionals even where reliable for ground state total energy calculations, are prone to significant errors in their TDDFT description of charge-transfer excitations in molecular systems [34,35,36,37] and excitonic effects in semiconducting or insulating solids [38,39]. To mitigate these shortcomings, long-range corrected (LRC) XC functionals [34,35,36,38,39] within the generalized Kohn-Sham (GKS) approach [41] have been proposed and applied successfully in many molecular and solid-state systems. However, LRC-XC functionals that provide a reliable description of valence electronic excitations by improving the long-range asymptotic behavior of the XC potential [38,39] do not necessarily improve the description of localized core-electrons and their excitations relevant to XUV and Xray spectroscopies [14,37,31]. In particular, as shown recently [31], LRC-XC functionals can significantly underestimate core-excitonic effects in extended systems. In the context of molecular systems, a number of groups have proposed short-range corrected (SRC) functional forms that are either tailored specifically to improve the description of coreexcitations [14] or more generally to provide a balanced treatment of both core and valence electronic excitations simultaneously [35,37]. In this work a simple procedure is proposed to incorporate short-range corrections for core-electrons in an extended periodic system such as 2D h-BN. It is shown that by augmenting a screened LRC-XC functional with a SRC component, the description of core-level excitonic effects is significantly improved providing a balanced treatment of valence and core-excitations in VG-RT-TDDFT.

Velocity-gauge 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] realtime TDDFT framework has been described in detail elsewhere [30,31]. Within this approach the time-dependent generalized Kohn-Sham (TDGKS) equations for electron dynamics with φ i being single-particle orbitals andĤ GKS the VG GKS Hamiltonian, are propagated in real-time. The VG GKS Hamiltonian is given bŷ and features a kinetic term that includes coupling to the time-dependent vector potential − → A (t) representing the applied external fields, the VG electron-nuclear interaction termV ion , the Hartree potential and a non-multiplicative XC operatorV XC as a functional of the instantaneous single-particle GKS density matrix ρ( − → r , − → r , t) [43].
Time propagation of equation 1 yields at every time-step, relevant quantities such as the GKS 1-body density matrix which in turn is obtained from the time-dependent current density − → j ( − → r , t) given as and includes the generalized momentum 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].

Tuned range-separated hybrid functionals
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) exchangecorrelation (XC) functional form that provides improved accuracy compared to semilocal 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 ∆ self-consistent field (∆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: 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 where the E SR HFX , E LR 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 LDAC is set to be that of the LDA [47]. Thus the functional form adopted features fractions (α + β SR ), (α + β 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 α, β LR , ω LR are chosen such that which ensures the correct long-range 1 ∞r behavior of the XC potential in the r → ∞ asymptotic limit while tuning α, ω 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 β SR , ω SR on the other hand predominantly determines the core-level GKS  Figure 1: Charged supercell calculation of the absolute binding energy in eV (black line) of the N 1s core-level in 2D h-BN with extrapolation to the infinite supercell size limit under the scaling L → α L of linear cell dimensions. α = 1 corresponds to a 18.927×16.392×18.927 a.u 3 orthorhombic supercell. α −1 → 0 corresponds to the relevant infinite cell size limit and the extrapolated binding energy is indicated by the red circle. (inset) Isosurface depicting the charge density difference between a N 1s coreionized configuration and the neutral ground state of 2D h-BN. Green and blue spheres represent B and N sites. Structure visualization is performed using VESTA 3 [49] eigenenergies. In principle both β SR and ω SR can be tuned separately so that the following condition is satisfied for a chosen set of core-levels: where ε 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 coreorbital 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 ∆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 ∆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 Table 1: Core-electron binding energies and the band-gap at Γ (∆ Γ G ) inferred from KS LDA or GKS SLRC functional eigenvalues are compared against target energies inferred from ∆SCF or GW [32] calculations. Within the GKS scheme [41], the SLRC functional eigenvalues reproduce the target single-particle excitation energies for both core-and valence-electrons to within 2% for the tuned parameter set α = 0, β LR = 1, ω LR = 0.129 Bohr −1 , β LR = 1, ω SR = 1.40 Bohr −1 . The standard KS system is not required to reproduce single-particle excitation energies except for the highest occupied orbital [55].
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 || ∞ → 1 [56] dictates a choice of α + β LR = 1 from eqn. 8 above. Physical considerations similarly suggest unscreened Fock exchange in the very short-range [57] which requires α + β SR = 1. While α is in principle a free parameter, in this work we set α = 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 β SR and the range-separation parameters ω SR , ω LR . The choice of α = 0 fixes both β LR = 1, β SR = 1 leaving ω SR , ω LR as the two parameters to be tuned in eqn. 6. The two coreelectron energy eigenvalues ε B 1s , ε N 1s and the Γ → Γ GKS band-gap ∆ Γ G 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 ε B 1s , ε N 1s , ∆ Γ G and their comparison to target data are shown in Table 1 for the optimized parameter set α = 0, 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 details
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 Γ-centered 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.

Linear response
To evaluate the tuned SLRC functional, firstly the response of 2D h-BN to a weak Dirac-δ 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.  Table 2: Computational parameters used for valence and core excitation simulations of 2D h-BN. The LCAO basis set is specified using nl-ζ notation where n, l are principal and azimuthal quantum numbers respectively and ζ represents the number of functions of each nl type. Additional ghost basis functions of s character are used in the vacuum region surrounding the 2D h-BN layer within the LCAO approach. Hartree and semilocal exchange-correlation potentials are evaluated by representing the charge density over a real-space mesh [58] for which an equivalent planewave cutoff is listed.  Table 3: Energies of key valence and core-excitations in 2D h-BN involving as obtained from RT-SLRC TDDFT are compared with reference values. The % error in absolute RT-SLRC energies is the range of 2%.
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, realtime 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 1 ∞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 The current fluctuations include low and high frequency components that correspond to valence and core-level dipole excitations respectively. (b) Valence optical absorption spectrum in 2D h-BN from the RT-SLRC approach (red) compared to BSE spectra (blue) digitized from reference [32]. (c) Boron and (d) Nitrogen K-edge X-ray absorption spectra in 2D h-BN from the RT-SLRC approach (red) compared with reference experimental data (blue) digitized from references [61] and [62] respectively. Adiabatic LDA (RT-LDA) spectra are also shown (gray) for both the valence and core-excitation cases.
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 Focklike 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 → π * 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 X-ray absorption
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 − → I (t) (Eqn. 3) within a material that is under the influence of external laser fields described by the vector potential obtained from an electric field − → E (t). In the frequency domain, the macroscopic current − → I (ω) associated with a weak probe pulse considered to be in the linear response regime and described by a given external perturbation − → E (ω) can be written as where σ is the conductivity tensor. In a pump-probe scheme for a certain choice of pump and probe fields described by a field profile − → E pump-probe (t) one can calculate through propagating the GKS equation 1 the associated macroscopic current − → I pump-probe (t) using equation. 3. Similarly, for the pump field − → E pump (t) alone, once can calculate the corresponding current − → I pump (t). Then the probe current in the pump-probe scheme is defined as This separation of the total current into pump and probe parts is sensible for probes of any intensity but − → 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 P as well as the pump-probe delay ∆T . Additionally, for weak probe pulses, equation 11 can be generalized to in terms of the transient frequency dependent conductivity tensor σ P (ω, ∆T ) which now depends on the pump parameters P and the delay ∆T . For comparison to XTAS experiments, it is useful to calculate the transient dielectric tensor 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 pumpprobe 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 ∼ 10 12 W/cm 2 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 → π * and 1s → σ * excitations. With this set up, the transient dielectric functions for different pump intensities are calculated and its imaginary part (ε 2 ) is plotted in Fig. 5 for both the B and N K-edges. By comparing the transient dielectric function (tr-ε 2 ) after the pump with its ground state counterpart (gs-ε 2 ) calculated with no pump, it is apparent that tr-ε 2 is modulated in several ways at the B and N K-edges. Firstly, it is seen that relative to gsε 2 the 1s → π * peak in tr-ε 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-2p orbital derived π and B-2p orbital derived π * characters respectively (see Fig 6). Valence excitations in the energy range of the pump pulse therefore exhibit π → π * character with the electronic part of the electron-hole wavefunction being predominantly π * like [32]. This is equivalent within a projected single-particle picture to a transient filling of the π * CBM states [29] following laser excitation by the pump. As the probe pulse that immediately follows the pump drives 1s → π * transitions, a certain fraction of these are blocked in proportion to the transient occupation of the π * CBM states. Thus the intensity of the 1s → π * peak in tr-ε 2 is reduced. Note that the higher lying 1s → σ * peaks occurring near ∼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 σ * states. Meanwhile, the pump also transiently induces a hole population near the VBM which is of π character. Hence from the perspective of the probe pulse certain 1s → π 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 ∼ 397.5 eV peak observed in the N K-edge tr-ε 2 spectrum which is missing in gs-ε 2 . A similar feature is not observed at the B K-edge however because the B 1s → π oscillator strength is much weaker relative to the B 1s → π * strength given that the VBM has very little B-2p character. This difference between In addition to the intensity modulations discussed in the previous paragraph, trε 2 peaks also exhibit certain energy shifts on the order of 0.1-0.4 eV. In particular we note that the B 1s → π * peak at 188.1 eV is transiently red-shifted (see Fig 5) while the N 1s → π * peak at 403.3 eV is blue-shifted. Higher lying 1s → σ * 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-2p and B-2p 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-ε 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.

Conclusions
In summary, this work investigated the use of range-separated hybrid functionals employing both short-and long-range corrected (SLRC) nonlocal screened exchange contributions in the context of real-time TDDFT simulations of solid state systems. It was shown that in a prototypical strongly excitonic system like 2D h-BN, the use of SLRC potentials within a generalized Kohn Sham TDDFT framework can enable a balanced treatment of both valence and core excitations simultaneously improving absolute spectral energy positions as well as overall lineshapes relative to semi-local XC approximations. The use of such SLRC potentials to directly simulate attosecond transient absorption spectra of core-excitations was illustrated. Following laser driven electron-hole excitations at energies near the first valence excitonic peak of 2D h-BN, nonequilibrium absorption spectra at the B and N K-edges were shown to exhibit characteristic features consistent with population induced state blocking and transient core-level energy shifts. The findings in this work suggest that the range-separated hybrid functional form in spite of certain limitations has sufficient flexibility to provide a practical and computationally efficient methodology for improving the accuracy of TDDFT based time-domain spectroscopy simulations in condensed matter.