Photoionization and transient Wannier-Stark ladder in silicon: First principle simulations versus Keldysh theory

Nonlinear photoionization of dielectrics and semiconductors is widely treated in the frames of the Keldysh theory whose validity is limited to small photon energies compared to the band gap and relatively low laser intensities. The time-dependent density functional theory (TDDFT) simulations, which are free of these limitations, enable to gain insight into non-equilibrium dynamics of the electronic structure. Here we apply the TDDFT to investigate photoionization of silicon crystal by ultrashort laser pulses in a wide range of laser wavelengths and intensities and compare the results with predictions of the Keldysh theory. Photoionization rates derived from the simulations considerably exceed the data obtained with the Keldysh theory within the validity range of the latter. Possible reasons of the discrepancy are discussed and we provide fundamental data on the photoionization rates beyond the limits of the Keldysh theory. By investigating the features of the Stark shift as a function of photon energy and laser field strength, a manifestation of the transient Wannier-Stark ladder states have been revealed which become blurred with increasing laser field strength. Finally, it is shown that the TDDFT simulations can potentially provide reliable data on the electron damping time that is of high importance for large-scale modeling.

The Keldysh photoionization theory for atoms and solids published in 1964 [1] plays a vital role in interpretation of ultrafast laser-induced phenomena, the field which continue to rapidly develop and influence our knowledge about fundamental processes in physics, chemistry, and biology. However, its applicability is limited by the requirements of 'sufficiently small' photon energy compared to ionization potential (band gap in solids) [1,2] and not too strong laser fields [3]. For band-gap solids, the theory [1] enables calculations of the number of electrons excited from the valence to the conduction bands per time unit, using a simplified description of the electronic levels reduced to two bands. This simplification has become commonly employed for qualitative simulations of laser-induced materials damage [4][5][6] while experimental observations indicate that the Keldysh theory can considerably underestimate or overestimate photoionization rates, depending on irradiation regime and kind of solid [2,7,8]. Direct experimental investigations of multiphoton inner excitation of band-gap materials was made possible relatively recently [9][10][11][12]. Supported by comparisons with experimental data, several modifications of the Keldysh model for solids were proposed to address discrepancies between experiment and theory [3,8,13,14].
An important feature of the original Keldysh theory is the dependence of the effective ionization potential for atoms and the effective band gap energy in the case of crystals, U eff , on the laser field intensity [1]. Although the increase of the U eff value with intensity was proven experimentally for atomic gases [15,16], an eventual increase of the band gap energy in crystals [4,17] irradi-ated by linearly polarized light was, to our knowledge, not observed for bulk materials. Recently, a clear shift of the energy levels to higher values was found in monolayer WS 2 [18] under irradiation by circularly polarized femtosecond laser pulses while a replication of the electronic levels by laser dressing was well demonstrated in GaAs [19]. Moreover, recent experimental studies performed for Si and ZnO [20,21] supported a reduction of the band gap energy in the laser field. These findings are consistent with the Keldysh theory where, in the multiphotonic regime, U eff increases with the field amplitude (γ 1 with γ to be the Keldysh adiabadicity parameter) whereas in the tunneling regime the contribution of U eff fully fades out (γ 1). Nevertheless, the interpretations of γ and U eff remain elusive, particularly for the cases of the mixed regime (γ ∼ 1).
Keldysh [1] attributed the increase of U eff to the Stark effect. As a whole, the Stark shift can transiently affect the band gap energy at high intensities while the laser dressing may eventually close the band gap, e.g., in regime when γ 1 [22]. However, there is no clarity of the manifestation and interplay of these effects at different irradiation regimes that is of high importance for predicting the consequences of the action of ultrashort laser pulses on band-gap materials under the real experimental conditions [4][5][6]. There is a need to verify the applicability of the Keldysh theory to bulk materials in a wide range of the irradiation parameters and to find ways of determining the photoionization rates beyond its validity limits and beyond the perturbative regime [23].
In this Letter, we first confront the qualitative insights from the Keldysh theory with a Floquet model supported by the density functional theory (DFT). Then, we perform a quantitative comparison of the photoionization rates obtained in the frames of the Keldysh model and in the simulations based on first-principles time-dependent density functional theory (TDDFT) for a wide range of laser intensities and photon energies for silicon as an example. Finally, based on the TDDFT simulations, we analyse the laser energy absorption and compare it with predictions of the Drude model. Means to improve accuracy of the Drude model for macroscopic description of materials optical response are discussed. At the basis of the Keldysh photoionization model, there is temporal averaging of a Hamiltonian which describes the interaction of the electromagnetic field with a two-level system, representing an atom, a molecule or a solid in a simplified manner [24]. When performing the temporal averaging of electronic energy levels over one laser cycle [25,26], the so-called laser dressing manifests as a replication of the energy levels spaced by the photon energy. The discrete set of possible quasi-states appearing transiently in the energy gap is usually referred as the Wannier-Stark ladder (WSL), the effect which can be described qualitatively using a Floquet model [27]. Figure 1(a) presents a schematics of the WSL quasienergy levels of silicon shifted by the laser field under the combined action of laser dressing and the optical Stark effect. The ground-state band structure shown on the left was calculated using density functional theory (DFT) (see Supplemental Material [28]). Depending on the choice of laser polarization, wavelength and intensity, the distances between the quasi-energy levels can either increase or decrease in presence of the strong optical field as schematically shown on the right part of Fig. 1(a). Figure 1(b) shows the calculated energies of the laser dressed states as a function of laser intensity, obtained using the Floquet model for the m = 8 states at the Γ point of Si. Note that the dipolar matrix elements [27] were computed using the DFT based on the local-density approximation (LDA). With increasing laser intensity, the number of replicas, which are associated with the increasing contribution of the anharmonic motion of the electrons, also increases [26]. With n = 2 replicas per ground electronic state, the Floquet model reveals multiple crossings (as well as anti-crossings, see [29]) of the dressed electronic levels that may be manifested as a transient metallization during the laser pulse [29]. The latter effect gives birth to important applications in ultrafast optoelectronics [22,30,31].
Let us next consider the results of numerical simulations for silicon irradiated by ultrashort laser pulses, which have been obtained by solving the Kohn-Sham (KS) equations using the time-dependent density functional theory (TDDFT) [11,32]. The laser field was introduced in the form of a time-dependent vector potential (velocity gauge). The crystal Si <100> is modeled by the method of ab initio pseudo-potentials [33] (see details in Supplemental Material [28]) involving periodic boundary conditions. The band structure calculated based on the LDA yields the direct band gap for Γ → Γ transition to be equal to E Γ g = 2.56 eV [34]. The approach enables to account directly for the effects of the laser field on the valence and conduction electrons from first principles, in particular, at intensities where a perturbative approach is not applicable [35] and with accounting for a realistic band structure [12,36,37].
Using TDDFT, the number of electrons excited from the valence to the conduction bands, n exc , was calculated using [10] where n and n are band indices, k indexes the electron wave-vectors in 3 dimensions, ψ GS n,k (r) are the KS wavefunctions of the ground state, and ψ * n ,k (r, t) are the KS time-and space-dependent wave-functions, where * indi-cates complex conjugation. V is the volume of the simulation box (constant in this work), N tot is the total number of electrons in the simulation volume, expressed by N tot = n,k ψ GS n,k (r) 2 . The summing is performed over the occupied states in the valence bands. Note that in the ground state, before the laser pulse action, the obtained density of electrons in the conduction bands n exc (t = 0) is zero. The gauge-related problems associated with the virtual excitations were avoided by exploiting the simulation results obtained shortly after the end of the laser irradiation [10]. The KS equation was solved numerically for a number of laser parameters by varying the photon energy ω and the laser intensity I. Details of the numerical approach are given in the Supplemental Material [28]. Figure 2 presents the results of TDDFT simulations for 10, 20 and 30 fs laser pulse durations (a "softened" top-hat temporal shape, see details in the Supplemental Material [28]) at the wavelength λ = 3200 nm ( ω = 0.387 eV). The number of photons per electron for crossing the direct band gap can be estimated as n ∼ E Γ g / ω that is 7 photons at this wavelength. Using the scaling law n exc ∼ I n [38], the slope of the function n exc (I) in Fig.  2(a) can be associated with 5 photons that are needed to cross the band gap at relatively low intensities in the range of 2 × 10 10 − 1.5 × 10 11 W/cm 2 (γ 1) [39]. With increasing the laser intensity, the simulations are consistent with experimental trends of a decrease of the n exc (I) slope [38], which can be conventionally fitted with the I n scaling law at lower n numbers. The corresponding intensity ranges with different n are marked in Fig. 2(a) and shadowed by colors. For each n-range of intensity, an effective (pulse averaged) multiphoton ionization coefficient σ n can be estimated using a simplified formula ∂n exc /∂t ≈ n exc (τ p )/τ p = σ n I n /n ω where n exc (τ p ) is the number of the electrons excited to the conduction bands by the end of the laser pulse of duration τ p . For τ p = 10 fs at λ = 3200 nm, the TDDFT results can be fitted with σ 1 = 2.93 × 10 6 m −1 , σ 2 = 2 × 10 −10 m/W, Note that in a purely multiphoton regime, the Keldysh theory predicts increasing n with intensity (see Refs. [1,4,17]) while the total photoionization rate can be considered as a sum of n-photon processes with different n [14]. Based on ab initio simulations, we interpret the reduction of the photon number n necessary for transition across the band gap as a mutual contribution of laser dressing of electronic states and tunneling that reduces the effective band gap during the laser action. At high intensities when γ < 0.1, the quantity of electrons excited to the conduction bands saturates to the number of electrons available for ionizations (4 valence electrons per atom in Si).
The corresponding excitation rates w TDDFT PI are com-pared with the Keldysh theory in Fig. 2(b). In both approaches, the effective electron mass was taken from the DFT calculations of Ref. [40] (m * = 0.2226m e with m e to be the electron mass in vacuum). It should be stressed that we use here the Keldysh theory for solids corrected by Gruzdev [8] and referred further as the KG model (black solid line in Fig. 2(b)). At laser intensities below 10 13 W/cm 2 , the w TDDFT PI values are considerably larger than the KG photoionization rates. This is in line with statements in [2,8] where it was noticed that the Keldysh theory may underestimate the photoionization rate by order(s) of magnitude. At higher intensities, in the TDDFT simulations the saturation regime is achieved, determined by the number of valence electrons available for ionization. In this regime, the KG rate and its tunneling limit (dashed line in Fig. 2(b) [41]) approach and finally exceed w TDDFT PI .
The results for the original atomic Keldysh theory, both its analytical expressions and numerical integration, are also added in Fig. 2(b) (dash-dotted and dotted lines respectively), obtained for a virtual atom with an effective ionization potential of 2.56 eV [28]. Surprisingly, at intensities below 10 11 W/cm 2 , the atomic theory fits reasonably well the TDDFT results although at higher intensities it deviates from the ab initio data. Important is that replacing the saddle-point approximation in the Keldysh theory by an exact integration has a minor effect, contrary to a conjecture in [2]. As noted in [42], further advances of the Keldysh theory of photoinization for solids can be seen in introducing the realistic band gap structure. However, the rigorous TDDFT simulations presented here represent a solid alternative, both from fundamental and practical points of view. As an illustration, the TDDFT results and their comparisons with the Keldysh theory are given for λ = 1600 nm and 800 nm in the Supplemental Material [28].
Comparing w KG PI (I) and w TDDFT PI (I), one can notice that the KG model yields oscillations of the multiphoton excitation rate with I while the TDDFT-calculated rates demonstrate rather smooth intensity dependences ( Fig. 2(b)). This difference can also be explained by the two-band approximation applied in the KG model, while the TDDFT accounts for multiple bands. We note that the zeros of nth order Bessel function correspond to a regular suppression of the effective transition probability called dynamical localization or destruction of tunneling [43,44]. In a two-band model, this effect manifests itself as a swift decrease of the excitation rate at certain intensities, whereas in the case of a multiband description, the Bourget theorem secures that only a single transition may be simultaneously disabled. As a result, the functions w PI (I) provided by the TDDFT multiband simulations are smoother than those obtained in the frames of the KG model.
An important comment should be made on the excitation rates derived from the TDDFT simulations and cal- , black dashed line) are added for comparison. Also for illustration, the data according to the analytical expression for the excitation rate for atoms [1], w at-anal PI , and a more exact numerical calculations of integrals of the Keldysh theory w at-num PI are shown (see Supplemental Material [28]). culated by the Keldysh theory for the saturation regime. In the Keldysh approach, the ionization rate is calculated by averaging the probability over many laser cycles. In our ab-initio numerical calculations, electron excitation follows the instantaneous laser field and, hence, it is strongly varying during the period of the electromagnetic wave. Furthermore, valence electrons, which are pushed out to the conduction bands at the first part of the laser pulse, can partially return back to the valence bands when the electric field of the wave changes direction [30,45]. Our estimations of the excitation rates given in Fig. 2(b) for three pulse durations correspond averaging over one, two, and three laser periods (note that one period for λ = 3200 nm is ∼ 10.7 fs). Thus, in the saturation regime when all or almost all valence electrons are excited to the conduction bands already during the first cycle, the w TDDFT PI value at 10 fs looks to be more reliable although the excitation rates for all three pulse durations give the same result, transfer of all valence electrons to the conduction bands ( Fig. 2(a)). Note that the KG model does not account for decreasing the number of valence electrons available for excitation into the conduction bands at high intensities and its direct use can lead to unphysically large number of electrons in the conduction bands ( Fig. 2(b)). When a considerable fraction of valence electrons has been transferred to the conduction bands by the front part of the laser pulse, the excitation rate in the trailing edge of the pulse should naturally be decreasing by the dynamic factor (1−n exc V /N tot ) (in our notations; see, e.g., [46]). Further still, the validity range of the KG theory is limited to single ionization per atom while the TDDFT simulations are not subjected to this limitation.
Analysis of the TDDFT-calculated excitation rates as a function of the photon energy provides insights into the WSL and the Stark shift in bulk crystals. Figure 3 presents the comparison of the excitation rates w PI calculated by the KG model ( Fig. 3(a)) and by the TDDFT (Fig. 3(b)) for the Γ → Γ transition in Si. The data are given for the electric field amplitudes E of the electromagnetic wave from 0.5 V/nm to 5 V/nm (curves from bottom to top). The corresponding intensity range is 3.32×10 10 − 3.32 × 10 12 W/cm 2 with the maximum intensity close to the saturation regime for λ = 3200 nm ( Fig. 2(a)).
At low electric fields, both the Keldysh formula and the TDDFT simulations reveal regular drops of the electron excitation at the photon energies close to resonances of the band gap at rest (marked by dashed vertical lines at E Γ g /n with n = {2, 3, 4, 5, 6} in Fig. 3). For each n However, the TDDFT minima sometimes fall beyond the resonant photon energy (e.g. the minima seen between 3PA and 4PA in Fig. 3(b)). Another qualitative difference between w KG PI ( ω) and w TDDFT PI ( ω) can be admitted for IR spectral range at high excitation fields. In the TDDFT simulations the resonance drops are vanishing with increasing laser intensity ( Fig. 3(b)). This is explained by a gradual transition to the tunneling ionization regime. In addition, as noticed above, the KG model does not account for decreasing the number of valence electrons available for excitation into the conduction bands at high intensities that should contribute to vanishing the w PI resonances. We also note that coupling between moving charges in the valence and conduction bands [47,48] can influence the excitation. The latter effect is not included in the Keldysh theory but can be captured by other approaches [49].
In pump-probe and angle-resolved photo-emission experiments [19,50], the possibility to observe the effects of quasi-energy levels (the WSL) in bulk band-gap crystals was recently proven. Here, based on the TDDFT simulations of the laser energy absorbed by electrons (see solid lines in Fig. 4 and Supplemental Material [28]), we anticipate a possibility to directly observe the WSL by measuring the change in time-integrated optical transmission when accurately tuning the pump laser wavelength, preferably in the IR spectral range. Another fundamental aspect, which is of high importance for multiphysics large-scale modeling of the interaction of ultrashort laser pulses with band-gap materials, is an adequate description of dynamically changing reflectivity and absorptivity, which are usually treated in the frames of the Drude model [51]. The TDDFT simulations give an excellent opportunity to verify on how accurate can be this simplest optical theory based on the classical equations of electron motion in an optical electric field.
The TDDFT and the Drude model have been compared based on calculations of the laser energy absorption at a wide range of laser parameters (Fig. 4). For this aim, the TDDFT-calculated density of the laser-excited electrons n TDDFT exc was introduced into the Drude formalism coupled with the energy balance equation for the electrons absorbing laser radiation. The latter was integrated for 20-fs laser pulses of the same shape as in the TDDFT simulations (see details in the Supplemental Material [28]). By employing a Drude damping time of τ D = 6 fs for all tested wavelengths and field strengths, it is possible to obtain a reasonable agreement with the TDDFT results (Fig. 4). A discrepancy between the Drude model and the TDDFT becomes pronounced for longer wavelengths and stronger optical fields. This can be attributed to the fact that τ D is dependent on the density of electrons, their energy, and laser wavelength [51,52]. Thus, the TDDFT simulations can provide a possibility to derive the τ D values for a wide range of laser irradiation conditions and for different materials.
In conclusion, based on the TDDFT simulations we have investigated photoionization of crystalline silicon by ultrashort laser pulses in a wide range of laser intensities and for photon energies from UV to mid-IR spectral range. The excitation rates have been derived and compared with the Keldysh theory [1]. The photoionization rates obtained within the the Keldysh approach in its validity range are smaller by more than order of magnitude as compared with the TDDFT simulation results that is in agreement with previous observations [2,8] that can be attributed to several factors, including simplification of the band structure [42]. We anticipate that the TDDFT predicts well the photoionization process both within and beyond the limits of applicability of the Keldysh formalism. The excitation rates obtained in the first principle simulations represent valuable fundamental and practical information and, being tabulated in a wide range of the irradiation parameters, can be directly used in multiphysics large-scale models, such as laser beam propagation in band-gap materials [53][54][55].
By calculating the non-linear absorption as a function of photon energy for the Si crystal, we observed the quasienergy levels known as the Wannier-Stark ladder. The energy of these laser-dressed states are shifted by the Stark effect at high optical field strengths. With further increasing field strength, these levels become less pronounced and finally disappear in the TDDFT simulation results. Finally, we have verified the Drude model based on laser energy absorption calculated in the frames of the TDDFT and a simplified energy balance equation, thus showing that the TDDFT simulations can potentially provide reliable data on the electron damping time as a function of density and energy of electrons and laser wavelength.
It should be underlined that the results of the present work refer to pure photoionization-related processes in spatially homogeneous optical fields (dipolar approximation) which are not masked by electron-phonon coupling. For longer pulse durations than used here, the role of phonons becomes important, leading to indirect transitions, which are not accounted for in the present simulations. Furthermore, at longer timescales, lattice destabilization, non-thermal melting, Auger recombination, transport of hot quasi-free charge carries and electron energy relaxation are important for silicon and other semiconductors [45,51,52,56,57], which are beyond the scope of this paper.
Here, M12 and M21 are the material-dependent dipolar matrix elements characterizing the dipolar transitions between two electronic levels (indexed by 1 and 2); A is the field strength; ω is the photon frequency. Note that this model can easily be extended to n bands exchanging m photons [25]. [28]  The electrons in the crystal are modeled using the Kohn-Sham (KS) auxiliary system that provides a oneto-one correspondence of the total electron density n e (r, t) of the N non-interacting electrons with the electron density of the many-body problem [1,2]. We solve the time-dependent KS equations for the diamond structure of Si crystal along with applying periodic boundary conditions in all directions. The dynamics of the N electrons in the crystal is described by solving a set of N KS equations, expressed by [3] − i 2m e ∇ r + |e| c A (t) where ∇ r describes the real-space gradient operator, is the reduced Planck constant, m e is the bare electron mass, |e| is the elementary charge, and c is the light velocity in vacuum. Note that the non-local contribution of the pseudo-potential to the external potential has been omitted for simplicity.v ion is the ionic potential,v H is the Hartree potential, n e (r, t) is the total density of electrons (i.e., all electrons in the valence and conduction bands) at the position r and instant t, andv xc is the exchangecorrelation potential. ψ n,k (r, t) is the time-dependent KS wave-function of an electron located in a band n with a wave-vector k. The external vector field A (t) is introduced in the velocity gauge and is taken homogeneous in the simulation volume (dipolar approximation, see Ref. [4]). The ground-state KS wave-function ψ GS n,k (r) is obtained by solving self-consistently the static KS equations [4]. The ground state band structure is given by the energy representation of the eigenvalues ε n,k as a function of k taken along the path of high-symmetry points for Si crystal.

B. Calculations of the excited electron density
The excited electron density n exc (r, t) evolving in the conduction bands can be evaluated by a projection of the occupied time-evolved KS orbitals ψ n ,k (r, t) on the ground state KS orbitals ψ GS n,k (r) expressed by [5,6] n,n ,k V is the volume of the simulation box (constant in this work). N tot is the total number of electrons in the simulation box, expressed by N tot = n,k ψ GS n,k (r) 2 .
Note that initially electrons are absent in the conduction bands.

C. TDDFT calculations of the excitation rates w PI
To avoid the gauge-dependence during the laser pulse, we use the results obtained by the end of the laser pulse of duration τ p , and define a pulse-averaged excitation rate w TDDFT PI based on the density of excited electrons, expressed via n exc (t = 0) = 0 in our case.

D. Numerical details
The presented computation results were obtained for a bulk Si sample. We employed the primitive cell of Si composed of two atoms with an experimental lattice constant, a 0 = 5.431 Å. Non-orthogonal periodic boundary conditions in all directions were used to describe the bulk crystal. Calculations of the ground state Si were performed using a real-space discretization (with grid spacing of 0.227 Å) and converged using the local density approximation (LDA) for the exchange-correlation functional. Several TDDFT calculations were performed using a more accurate but computationally demanding TB09 meta-generalized gradient approximation for the functional [7]. Assessment of the validity of the LDA and TB09 functionals to describe the transient properties of Si is provided in Refs. [8][9][10][11][12][13]. All calculations were carried out using the open source code Octopus [14][15][16]. The modeling of induced fields was disabled. The atomic potential of Si is modeled using a norm-conserving pseudopotential [17]. The k-grid was refined until reaching a constant quantity of excited electrons. The convergence was reached with using a 24 × 24 × 24 grid in the k-space. The integration of the KS equation was performed using the enforced time-reversal symmetry (ETRS) algorithm [18]. The time step was reduced until convergence of the excited electron density below 1% variation achieved with a time step of 6.8 attoseconds for both functionals. Note that, upon using the TB09 functional, the temporal integration was performed based on predictor-corrector method [13]. The TB09-based TDDFT results do not considerably differ from those obtained based on the LDA potential. However, further investigations are needed in this direction.

E. Introducing laser pulses with a top-hat temporal profile
The electric field of the laser pulse is introduced in the KS equation using the dipolar approximation [4]. As a consequence, the time-dependent laser field is homogeneous inside the simulation cell. The vector potential A (t) is linked with the electric field E (t) via If to approximate this integral using the slowly varying envelope approximation, the following relation can be obtained A (t) ∼ icE(t)/ω. We are interested in the excitation rate of the electrons transferred from the valence bands to the conduction bands. In order to provide this value on a similar basis as in the Keldysh theory where the ionization rate is averaged over laser cycles in a constant-amplitude field [19], we introduce a "softened" top-hat (STH) laser pulse of a duration τ p . The main part of STH pulse represents a plateau of the duration τ p with rising and decay phases, both of the duration of τ r , which are short compared with τ p . The pulse is expressed as: The laser frequency is given by ω = 2πc/λ and φ is the carrier envelope phase (CEP), which was set to 0. The rising time τ r for the STH pulses is used here to avoid a temporal discontinuity in the field amplitude when the laser pulse starts and terminates. The τ r value was optimized to obtain the minimum final excited electron density. This optimum was reached for τ r = 4π/ω.

A. Keldysh-Gruzdev excitation model for crystals
To calculate the number of excited electrons, several theories of electron excitation for solids were proposed [20]. Gulley et al. summarized the Keldysh theory without amending modifications [21][22][23][24], whereas corrections were proposed by Gruzdev et al. [25][26][27][28], to account for spin degeneracy. More recently, McDonald et al. [29,30] have used a model based on semiconductor Bloch equations [31] to further study the effect of band dispersion in electron excitation where not only the envelope of the pulse but also the phase of the laser pulse can be accounted for.
In this Section the analytical model is detailed, which was used to calculate the w KG PI values which are shown in Figs. 2-3 of the main manuscript. The Kane band structure, which is applicable for narrow band gap materials [27], is considered in this model. In the general form, the nonlinear photoionization rate w PI (I) in a constantamplitude field is expressed as According to the Keldysh-Gruzdev excitation model, the photoionization rates are expressed as .
Here m * = 0.2226m e is the electron effective mass with m e to be the electron mass in vacuum [32]; n = U eff / ω is the number of photons per electron for overcoming the effective potential barrier. Expression (9) employs the Dawson integral G (z) = z 0 dy e y 2 −z 2 , which is calculated numerically.K (x) andÊ (x) are the complete elliptic integrals of the first and second kind respectively, which have the formsK (x) = π/2 0 1 − x 2 sin 2 θ −1/2 dθ and E (x) = π/2 0 1 − x 2 sin 2 θdθ. The effective ionization potential U eff accounts for the energy shift induced by the Stark effect that can be written as U eff = 2Eg πF1Ê (F 2 ). F 1 and F 2 are the functions of the adiabaticity parameter γ: F 1 = γ/ 1 + γ 2 , F 2 = F 1 /γ. The γ value is given by It characterizes irradiation regimes via the peak intensity I peak or the electric field amplitude E peak . When γ 1, the electron excitation from the valence band to the conduction band occurs via tunneling ionization whereas γ 1 corresponds to multiphoton ionization. For γ 1, the Keldysh model for crystals describes the tunneling mechanism of photoionization in the form where E g is the bare band gap energy [27]. Note that this formula is not applicable at low fields where ionization is governed by the multiphoton mechanism [33].

B. Keldysh ionization model for atomic gases: analytical and numerical integration
In this section, the formulas of the Keldysh theory for ionization of an atomic gas are provided. We have employed the analytical atomic Keldysh model given by Eqs. (16)-(18) from Ref. [19] to compare this theory with both the TDDFT and the Keldysh theory for excitation of band gap solids (see Fig. 2 of the main manuscript): The Dawson integral G (z) is given in Section II A. Numerical integration. We use an estimative adaptation of material photoionization based on the Keldysh theory for atoms [19]. For this aim, we consider a virtual hydrogen-like atom with the energy of its ground state equal to the material band gap, E Γ g = 2.56 eV [12]. The electron of the atom has mass m * , which is taken to be the same as in both the Keldysh theory for band gap solids (Eqs. (7)-(9)) and the TDDFT simulations. This gives formulas for the ionization probability of atoms analogous to Ref. [19]. Note that the ionization rate (Eq. (12)) is multiplied by the atomic density of the solid ρ (silicon in our case) to express the rate in m −3 s −1 .
When solving Eq. (12), we calculated the integrals numerically instead of using the saddle-point method. The integrals have the following form To calculate the integrals, a fine temporal grid t k = k × τ (k = 0, ..., N t ; τ = 2π Ntω ) was used. Within the time segments t k−1 < t < t k+1 , the function F is approximated by the polynomial For η, the linear expansion is used, . Consequently, we have Such calculations of integrals are time-efficient and sufficiently precise, provided that the time step τ is much smaller than the laser cycle.

III. THE TDDFT RESULTS FOR DIFFERENT WAVELENGTHS; FITTING OF THE KG MODEL
According to our TDDFT simulations, the first principles approach yields photoionization rates, which are at least an order of magnitude higher as compared to the Keldysh theory for solids. This has systematically been studied at different laser wavelengths, see Figs. 1-3. We remind that, in our study, the fields induced by the movement of charges in the time-propagation of KS orbitals were disregarded (Eq. (1)). Under this assumption, the laser pulse inside a bulk material is compressed with the intensity multiplied by the material refractive index and it should be the same both in TDDFT simulations and in the Keldysh formulas. However, the TDDFT approach includes several important factors, which are absent in the Keldysh theory and thus can be responsible for the observed discrepancy.
One of the most important factors is the realistic bandgap structure introduced in the TDDFT, which is dynamically varying upon electron excitation into the conduction bands with corresponding distortion of interatomic potential. As a result, the electron wave functions are subjected to the action of the field of the laser wave superimposed with the dynamic interatomic field. One can anticipate that the local field acting on the electronic component of the crystal is enhanced similarly to predictions of Gaier et al. [34] while not necessarily in a fluctuation manner. Another factor, also connected with the dynamic band structure, is the so-called laser dressing of the electronic states or, by other words, appearance of transient quasi-states usually referred as the Wannier-Stark ladder [35]. We can also mention the Abraham-Minkowski problem [36] connected with the momentum of photons inside material, which still calls for investigations. Although this goes beyond the scope of the present study, a particular attention on the role of induced fields in the conservation of the momentum at interfaces is envisioned [6].

Figures 2(b) and 3(b)
show that the KG model underestimates the results of the TDDFT simulations where the band structure is more realistic and laser dressing is naturally addressed. Interesting is to find a factor ζ of "laser field amplification", E → ζE, at which the KG model would fit the TDDFT results. We have performed this procedure for different wavelengths and the results are presented in Figs. 1(b), 2(c), and 3(c) respectively for 3200, 1600, and 800 nm for both the KG model and its tunneling limit. By choosing the ζ value, it was surprisingly found that the accurate enough fits were achieved at ζ = n(λ) where n is the refractive index at the corresponding wavelength λ. It is not clear yet if it is a pure coincidence or it hides a physical cause. Also interesting is that for all studied wavelengths the atomic Keldysh theory applied for our virtual atom agrees reasonably with the TDDFT simulation results at low intensities (in the multiphoton regime, see Figs. 1-3). (b) The same as in (a) but with fitting the KG and tunneling rates to the TDDFT results using normalization of the laser field in the KG formulas using the factor ζ = n(λ) (see text).

IV. LASER DRESSING, EFFECT OF RABI OSCILLATIONS AND THE ROLE OF BESSEL FUNCTIONS IN STRONG FIELD EXCITATION
During laser action, response of a bandgap material to irradiation can be described by shifting the energy levels of the band structure as compared with the ground state structure [37,38]. This shift is known as laser dressing and it originates from the coupling of laser light to the electrons, which lasts the time of the laser illumination. As a result of the light-induced periodic electron motion, the dressed band structure does not represent anymore purely electronic states. Instead, the dressed energy levels are considered as quasi-particles named polaritons. This vision is not new and it was already employed in a number of formalisms, e.g., already in the Keldysh theory [19], in atomic physics [39] and in the non-equilibrium solid state physics [40,41]. Laser dressing may lead to ultrafast transient metallization [42,43], a phenomenon that has enabled the development of novel applications in ultrafast optoelectronics [44].
As stated in [42], a qualitative description of the laser dressing effect is possible from the knowledge of the electronic ground state. We underline that the approach employed here is based on a simplified model Hamiltonian. In particular, the importance of the term scaling in A 2 in the employed Hamiltonian, which is disregarded in our  approach, was investigated in a series of recent publications [41,45]. Also the dipolar matrix elements should be affected by intense laser fields, a phenomenon that is not described in the simplified method we have employed to prepare Fig. 1 of the main manuscript, which therefore should be considered as a contextual illustration. In Fig. 3(a,b) of the main manuscript, the minima of the excitation rates w KG PI ( ω) provided by the KG theory are outlined by grey lines. These reductions of w PI ( ω) are pronounced in the Keldysh theory and are also visible in the TDDFT simulation results, though mildly. In literature, the origin of such reductions of the amplitude probability w PI ( ω) was attributed to the suppression of tunneling [46][47][48] (equivalently, these are described in the real-space as Wannier-Stark localization [35,49,50]). The tunneling drop interpretation can be illustrated by an analytical solution constructed from a simplified Hamiltonian model describing the effect of the laser field on the conduction electronic subsystem [46][47][48]. Using a rigid band viewpoint, a series of the Bessel functions of n-th order appear in the corresponding analytical solution as a multiplicative term for each replication order n [46][47][48]. Since the Bessel function is inside the sum of possible wave-functions, an eventual suppression of interband transitions depends on the number of available electronic levels.
In a two-band description as in the Keldysh theory, suppression of the transition from one band to the other may take place when the Rabi frequency is twice the laser photon frequency [46][47][48] J n 2Ω Rabi ω laser = 0.
However, when several electronic energy levels are available, a simultaneous suppression of all possible transitions (a transition is denoted by its dipolar matrix element d i→j ) may not be feasible. Assuming that the total transition probability can be reasonably described by a series of Bessel functions in the physical reality, the Le Bourget theorem [51] suggests that only one transition could be possibly disabled for a given set of laser parameters (electric field of the wave, wavelength). This mathematical argument provides a physical explanation for rather smooth reduction of total transition probability w PI observed in multiband description such as TDDFT ( Fig. 3(b)) while a two-band description such as the Keldysh model reveals a clear suppression of transitions at given resonant frequencies [ Fig. 3(a)]. Then, since the Rabi frequency is proportional to Ω Rabi ∝ E · d i,j (E is the field amplitude), Eq. (13) can be rewritten as where J n (x) denotes the Bessel function of n-th order. The dipolar transition matrix elements d i,j = ψ i |r|ψ j were obtained from the DFT computation. Figure 4 depicts the value of the dipolar matrix elements for Si in atomic units. The elements indexed from 0 to 3 correspond to the valence band states, the rest correspond to the conduction states. Since the description of the Hamiltonian is hermitian, the matrix elements are symmetric above and below the matrix diagonal, evidencing the equivalent transition probabilities for excitation and recombination in the linear regime. Knowing the nodes of the Bessel functions of n-th order, an illustrative mapping of the transition to be selectively disabled can be provided (Fig. 5).  (t), the formula given in the Appendix H of Ref. [52] was used (it is also given in Chapter 5 of Ref. [4]).
Electron excess energy. The energy density de (r, t) introduced into the sample by the laser at each time interval dt can be expressed using the total current density j (r, t) and the laser electric field E (r, t) as Usually the total current is separated to contributions from the polarization (originating from interband transitions) and the free carrier current j free (t) related to intraband processes. The dynamics of the absorbed energy density e (r, t) can be described as interband absorption where w PI (r, t) is the instantaneous photoionization rate and n = Egap ω is the number of photons that are required to overcome the bare band gap energy E gap using the electric field of the laser wave at frequency ω. In order to compare with the total energy, calculated by solving the KS equation (Eq. 1), we integrate Eq. (16) over time in a volume V as The density n exc is calculated by Eq. (2). Using the dipolar approximation for the external field, we obtain the energy ξ el received by electrons from the laser field, which is defined as Time-dependent description of the free-carrier absorption The energy transferred from the laser light to the electrons via intraband absorption can be computed from the excited electron density n exc (t) and the conduction band electron current J free (t) by using of the generalized Ohm law, expressed by J free (t) = dt σ free (t, t ) E (t ) .
The contribution of the free electron currents can be described using a time-dependent conductivity σ free (ω; t) = −iωε 0 [ε Drude (ω; t) − 1], where the dielectric permittivity ε Drude (t, ω) is given by the Drude model for solids [53] ε Drude [n exc (t)] (ω) = 1 − e 2 m e m eff ε 0 ω 2 × n exc (t) 1 + i ν ω . (19) The electron collision (or damping) time τ D = ν −1 is usually considered to be from one to several femtoseconds. In this work, the effective electron mass is taken as given by the DFT [32], m eff = 0.2226. The collision time τ D = 6 fs was adjusted to obtain the best fit to the TDDFT simulation results (see Fig. 4 of the main manuscript).
(20) In the case of a quasi-continuous wave of frequency centered at ω 0 , one has E (ω) → δ (ω ± ω 0 ) E (ω) that yields To calculate the absorbed electron energy at a time moment t in a volume V , we integrate the above expression over time from t = 0 to t that gives The later expression depends on the excited electron density n exc (t) and on three free parameters, namely the electron collision frequency ν, the electron effective mass m eff , and the band gap energy E gap . We used Eq. (21) to compare it with the electron excess energy ξ TDDFT el obtained in our TDDFT simulations.

VI. EXCITATION PROBABILITIES AS A FUNCTION OF LASER PARAMETERS
From the calculated dependences n TDDFT exc (I peak ), the effective multiphoton rates σ n associated with n−photon transitions can be extrapolated using the multiphoton approximation expressed by ∂n TDDFT exc ∂t = σ n I n n ω .
We remind that the LDA functional can underestimate the bandgap energy of crystals, which is usually smaller as compared to more sophisticated modeling approaches for band structure calculations [12]. Therefore, we have to note that here the effective multiphoton excitation rates σ n for silicon were calculated for the direct bandgap energy, which is somewhat smaller than the experimentally measured value. Using Eq. (22), one can fit the TDDFT results from Figs. 1(a), 2(b), and 3(b), thus deriving the σ n values. The data are summarised in Table  I for the intensities below the saturation regimes. Where possible, we provide comparison of our data with the multiphoton excitation rates available in the literature. As one can see, a reasonable agreement is achieved although in experimental studies pure multiphoton excitation can be masked by other processes involved at longer pulse durations such as collisional ionization and phonon-assisted indirect transitions.
We notice that for 800 nm wavelength, the value of σ 2 obtained using first-principle simulations somewhat decreases when increasing the laser pulse duration and the same is observed for other wavelengths and other σ n . This may originate from the dynamic behavior of the excitation. Once states at the bottom of the conduction bands are populated, the excitation probability decreases rapidly, an effect known as the Burstein-Moss effect [60]. Although the multiphoton ionization rates σ 2 obtained in the first-principles calculations show a dependence on pulse duration τ p , the σ n values mostly remain within the same order of magnitude. Also we have calculated the single-and two-photon absorption rates σ 1 and σ 2 using the LDA and TB09 functionals at wavelengths of 484 nm and 407.58 nm respectively, the latter in order to adjust the laser wavelength for the resonant excitation of single-photon transition, see Table I. In the tunneling regime where the photoionization rate obtained in the TDDFT simulations scales linearly with the peak intensity, an effective tunneling ionization rate σ 1 (using Eq. (22), n = 1) can be estimated as a function of pulse duration. For λ = 800 nm, σ 1 = 1.447 × 10 6 m −1 at τ p = 10 fs; σ 1 = 7.985 × 10 6 m −1 at τ p = 20 fs; σ 1 = 5.754 × 10 6 m −1 at τ p = 30 fs (not reported in the Table).
Note that the direct comparison of the TDDFT photoionization rates with the experimental values is not straightforward since the most experimental measurements involve also indirect transitions (Γ → L). Therefore, the comparison with experimental measurements is only qualitative here. Note also that indirect band gap transition rates can be computed using the TDDFT by employing a localized electric field, as was very recently shown by Noda et al [61].
The full database of pulse-averaged photoionization rates obtained by our TDDFT simulations is available at http://www.quantumlap.eu/ photo-ionization-database/.