{'fraction_non_alphanumeric': 0.03615488518685277, 'fraction_numerical': 0.012877082395317424, 'mean_word_length': 4.8098352079518705, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-Markovian pure dephasing in a dielectric excited by a few-cycle laser pulse

We develop the theory of pure dephasing in a solid exposed to an ultrashort laser pulse beyond the commonly used Markov approximation. This approach takes into account the finite cutoff energy of the bath and can be applied to both many-particle and phonon environments. With numerical simulations performed with the time-dependent Hartree-Fock equations, we investigate how the excitation probability and high-harmonic generation are described by different models of decoherence. It is shown that the time-dependent rates allow for temporally high dephasing to successfully reproduce the main features of high-harmonics spectrum and avoid an overestimation of the charge carrier population after the pulse, which is a common problem of the Markov approximation.

The majority of modern theoretical treatments of these systems are employing the Markov approximation describing the relaxation phenomena by constant rates 1/T 1 and 1/T 2 . This approach was successfully applied in the recent studies of high harmonic generation (HHG) in solids [5,[11][12][13][14] provided experimental methods allowing for a distinction of the inter-and intraband components of polarization via analysis of the spectrum, waveform and group delay of the emitted radiation. Notably, the HHG measurements in thin films of fused silica have shown that the group delay and scaling of the cutoff in high-frequency plateaus with the field amplitude demonstrate the features of an intraband current [13,15], whereas the measurements in a semiconductor (GaSe) demonstrate the leading role of interband transitions and their interference [16,17].
The state-of-the-art quantum-mechanical models significantly overestimate the interband polarization and require very short dephasing time T 2 ≈ 0.4 − 4 fs for reproducing the experimental data in SiO 2 [13,15] and ZnO [11,12]. This assumption is incompatible with simulations of carrier-envelope-phase (CEP) control of the current in dielectrics and semiconductors [4,18,19] and the dynamic Franz-Keldysh effect [20][21][22], where ultrafast dephasing was not required for reaching an agreement with experiment.
The recent multiscale ab initio simulations of HHG in a diamond [23], considered additional averaging by intensity distribution in the laser beam spot and emphasized the role of propagation effects in the build-up of a smooth harmonic spectrum. Nevertheless, these additional considerations still cannot reproduce the absence of a group delay dispersion in the emitted radiation [13,15], which is peculiar to the intraband current, and simulations with insufficiently high decoherence rates still predict the dominant contribution of interband polarization in the high-frequency plateaus of the spectrum [24]. On one hand, assumption of ultrafast scattering times ∼ 0.1 − 4 fs previously reported in the semiclassical simulations of a high-field transport in SiO 2 [25,26] and fully microscopic simulations of nonlinear spectroscopy in semiconductors [27,28] increases the intraband component of the current density to the reproduces both the high-energy plateaus of experimental HHG spectra and group delay [5,13,15,29]. On the other hand, as will be shown below, an ultrafast dephasing rate results in the overestimated spectral broadening and wrong scaling of the carrier population with the field intensity due to the opening of an artificial single-photon excitation channel.
The recently reported experimental results on the optically-controlled current suggest that the total charge induced by the laser pulse in a circuit scales close to the perturbative result at low intensities scales when the system is still in the multiphoton regime [4,19]. Thus, a more appropriate model of dephasing in dielectrics needs to be developed.
In this work, we theoretically investigate pure dephasing in a dielectric beyond the Markov approximation [30][31][32][33]. We show that it is possible to introduce the dephasing rate as a timedependent function with a slowly-varying envelope accounting for a non-trivial bath spectral function. Starting from the model Hamiltonian, we derive the semiconductor Bloch equations applicable for a description of pure dephasing beyond the Markov approximation. Both environments are analyzed in the framework of the harmonic oscillator model, where finite spectral cutoff energy is taken into account. The strong dependence of dephasing rate on the cutoff energy explains why the phonon bath is approximately Markovian even on a fewfemtosecond time scale and why the many-particle environment features a temporally high dephasing rate, which has not been observed in the Markovian limit. These results are illustrated by numerical simulations in both independent-particle and time-dependent Hartree-Fock approximations. It is shown that the time-dependent rate given by the excitationinduced dephasing with the envelope of harmonic oscillator improves the field amplitude scaling of excitation probability and produces high-harmonic spectra in a good agreement with experimental data.

II. APPROXIMATIONS AND APPLICABILITY LIMITS
We start from a discussion of material parameters in dielectrics and semiconductors to validate the standard approximations applied in simulations of ultrafast spectroscopy. Fig. 1 shows the state-of-the-art ab initio simulation of α-SiO 2 band structure in the self-consistent quasiparticle GW approximation. The circles show relative contributions of the Bloch orbitals to the excitonic 1S-state determined from the solution of the Bethe-Salpeter equations. The direct quasiparticle bandgap at the Γ point is E QP g = 10.1 eV, and the optical bandgap E g = 8.9 eV is shifted from it by the exciton binding energy E ex = 1.2 eV obtained by solving of the Bethe-Salpeter equations. Remarkably, the exciton binding energy in α-SiO 2 is by 2-3 orders of magnitude larger than that in the commonly studied semiconductors, e.g., E GaAs ex = 4 meV [34], E Si ex = 15 meV [35]. The value of E ex = 1.2 eV is larger than the Rabi frequency ω max R = F 0 max k |ξ nm (k)| ≈ 0.74 eV even at the field amplitude of F 0 = 2 V/Å close to the damage threshold, where max k |ξ nm,k | ≈ 0.37Å. [36] Therefore, electron-electron interaction plays a significant role in electron dynamics of the wide bandgap dielectrics, even though the excitonic peaks are broadened and not visible in the absorption spectrum at high intensities. This statement is also true for some other wide bandgap materials, e.g., CaF 2 [37,38], where the exciton binding energy is also ∼ 1 eV.  d dt E nK(t) and τ AS ≡ 2π/ max u nK(t) | d dt |u mK(t) = 2π/ max k |F 0 · ξ nm,K(t) |, respectively, K(t) = k +A(t) is the time-dependent crystal momentum given by the acceleration theorem [42], is the vector potential, and ξ nm,k = u nk |i∂ k |u mk is the matrix element of a coordinate operator in the crystal momentum representation.
Our ab initio simulations show that in α-SiO 2 and other dielectrics with large effective masses of carriers, the matrix elements ξ nm,k are slowly-varying functions of k, and thus τ AE τ AS . In the materials with small effective masses of electrons in the conduction band, e.g., GaAs, [43] the optical matrix element change more rapidly with k than the band energies, so the opposite situation (τ AE τ AS ) might be realized as well.

III. NON-MARKOVIAN MASTER EQUATIONS
If the system's evolution time is much longer than the correlation decay time, the scattering events can be viewed as instantaneous in comparison to evolution time, which leads to the Markov approximation with constant dephasing rates ∝ 1/T 2 . As shown in the previous section (see Table II), the many-particle environment does not satisfy this condition, since the correlation decay time is comparable to the evolution time τ corr. t max . The fat band plot in Fig. 1 shows that the 1S excitonic state is formed primarily by the Bloch orbitals at k = 0 and suggests that it is possible to adiabatically separate the Hilbert space. The fast-evolving single-particle states participate in the intraband motion and dynamic Bloch oscillations, while the slow carriers stay in the middle of BZ and form the many-particle states.
The full Hamiltonian is written as where is the Hamiltonian of the single-particle states occupying the accelerated Bloch (Houston) state at the band n and initial crystal momentum k, a nk and a † nk are their ladder operators, H L (t) is the non-adiabatic part of interaction with an external field responsible for interband transitions, is the bath Hamiltonian, where the ladder operators b † pq and b pq satisfy the Bose commutation rules in the cases of phonons and many-particle states including an even number of carriers (excitons) and the Fermi commutation rules for odd number of carriers (defects with trapped charge carriers, trions). The bath states are described by the quantum number p and quasimomentum q. For simplicity of notation, we describe them by a single index α ≡ {p, q}.
In the case of the fast system and slow environment, the model Hamiltonian of systembath interaction can be written as follows are the operators acting only on the system and bath states, g α is the amplitude of systembath interaction determined by intrinsic bath properties, and f nk,α (t) is the time-dependent part of the interaction amplitude depending on the rapidly changing parameters, e.g., charge carrier density or kinetic energy.
To obtain the equations of motion beyond the Markov approximation, we employ the time-convolutionless (TCL) projection operator technique [30,31,44], which yields the following equation for the reduced density matrix in the interaction representation Substituting (2) into (5) and replacing the t 1 with τ = t − t 1 in the integrand, one obtains where ∆t = t−t 0 is the time elapsed since the initial time is the bath correlation function, and B α (τ ) = g α b † α e iωατ . In the Markov approximation, one assumes that the bath correlation function C αβ decays much faster than ∆t, which allows extending the integration limit to infinity (∆t → ∞). We set t 0 = 0 and keep the finite limit of integration over τ to consider the non-Markovian case.
Eq. (5) is time-local only in terms of the density matrix ρ S (t), but it is not in the Lindblad form because the convolution between C αβ (τ ) and S mk,β (t − τ ) is still required. As shown below, in the case of pure dephasing, it can be reduced to the Lindblad form.
Following [44], we rewrite S mk (t − τ ) as where the evolution operator is separated in two parts.
The part describing the evolution from t − τ to t can be approximated as Assuming the completeness of the instantaneous eigenbasis, one can decompose the operator S β (t) into summation over all instantaneous Bohr frequencies ω(t) in the system: This leads to the following approximation for the operator Here, the convolution between C αβ (τ ) and S mk, is the spectral correlation tensor connected with the correlation function via the finite Fourier transform.
Pure dephasing is dominated by elastic collisions described by energy-conserving terms with zero Bohr frequencies ω(t) = 0. Then from (8) it follows that S nk,α (t − τ ) ≈ S nk,α (0, t), the relaxation superoperator (9) takes the Lindblad form, and the spectral correlation tensor is connected with the correlation function simply via time integration We neglect the non-diagonal elements of the correlation function and use a single index to enumerate the diagonal ones: C α ≡ C αα . Transforming Eq. (5) back from the interaction to the Schrödinger picture and using the definitions (3), we obtain the master equation in the Redfield form: where (11) is the time-dependent pure dephasing rate.
Thus, the pure dephasing rate can be written as a product of the slowly varying envelope dτ depending on the bath correlation function, and the rapidly varying function F Evolution of a quasiparticle interacting with an environment according to master equation (10) can be described by the effective non-Hermitian Hamiltonian where the electron energies E nk are replaced by the complex-valued time-dependent quasiparticle energies E nk (t).
Considering the field-matter interaction in the length gauge, one obtains the system of partial differential equations to the semiconductor Bloch equations [16,45]. Furthermore, applying the method of characteristics [46] to the partial differential equations, one derives the Bloch acceleration theorem ∂ t K(t) = −F (t) and the following system of ordinary differential equations similar to those in Refs. [47][48][49]: where is the matrix element of the field-matter interaction, is the change of a total quantum phase between the Houston states in the bands n and m, is the quasiparticle energy describing the electron interacting with an environment, E nK (t) = E nK + F (t) · ξ nn,K are the modified band energies accounting for the geometric phase contribution.
Including the Coulomb interaction between electrons and keeping only the first-order terms, one obtains the well-known semiconductor Bloch equations [50] Note that the diagonal matrix elements of the self-energy Σ TDHF nn (t) are real-valued quantities. Imaginary part of quasiparticle energy describing the damping of single-particle states due to interaction with the many-particle environment appears only in the higher-order approximations, e.g., the GW or coupled-clusters, which are very computationally expensive and currently applicable for real-time simulations of simple atomic systems [52]. Therefore, a reasonable non-Markovian model for γ nK (t) is still required to describe pure dephasing.

IV. BATH OF HARMONIC OSCILLATORS
As was originally demonstrated by Feynman [53], Caldeira and Leggett [54], the interaction with any structured environment can be rigorously mapped onto a bath of harmonic oscillators, if the interaction is sufficiently weak and perturbation theory is applicable. In this model, the bath is characterized by the spectral function describing the distribution of oscillator's energy levels and their coupling to the system.
To simplify further analysis, we assume the ohmic spectral density with an exponential cutoff [32,55] J(ω) = J 0 ωe −ω/ωc (20) allowing for analytical expressions of both the correlation function and relaxation rate. Here, J 0 = γ 0 /ω c is a dimensionless constant, γ 0 is the dephasing rate amplitude, and ω c is the cutoff energy. For the many-particle environment, the cutoff is defined by the exciton binding energy E ex , and for the phonon environment, it is given by the highest phonon energy LO .
If the bath in thermal equilibrium before interaction with the system and approaches it afterward, its correlation function is given by (Ref. [55], p. 181): where the factor coth ω 2k B T = 1 + 2N (ω, T ), N (ω, T ) = 1 e ω/k B T − 1 appears due to the Bose-Einstein population distribution N (ω, T ). Here, the spectral function is extended to negative frequencies as J(−ω) = −J(ω).
To make the integration in (21) analytical, we assume that the thermal energy k B T α of the environment is much smaller than the cutoff frequency ω c . For α-SiO 2 at the room temperature (T = 300 K, k B T ≈ 25 meV), this condition is fully satisfied for both phonon (E LO ≈ 150 meV) and many-particle environments (E ex ≈ 1.2 eV). Thus the dephasing rate envelope in the harmonic oscillator model is given by the following expression In general, the dephasing rate envelope in the harmonic oscillator model is determined For the phonon bath, it is comparable to the peak value due to initial quantum fluctuations, but for the many-particle bath, it is smaller by more than an order of magnitude, as shown in Fig. 2d. This property explains why the phonon bath is approximately Markovian even on a few-femtosecond time scale. On the other hand, to reproduce the behavior of manyparticle bath with the Markov approximation, one has to assume very short dephasing times T 2 , which was done in previous simulations of high-harmonics spectroscopy [5,12,13]. In the previous works on semiconductors exposed to terahertz laser pulses [56][57][58], the experimentally measured dephasing rate was well described within the excitation-induced dephasing (EID) model, where the scattering rate is inversely proportional to the mean inter-particle distance r(t) = ρ −1/3 (t) Here, ρ(t) is the density of excited charge carriers, and γ 0 is the dephasing rate accounting for level broadening due to intrinsic lattice defects and frozen phonons, and γ 1 is the dephasing rate due to electron-electron scattering.
The EID model can be obtained as a particular case of Eq. (11), where F   In Fig. 4, we compare the Markovian and excitation-induced dephasing models with a constant (EID) and harmonic oscillator envelopes (HOEID). The switch-on time for the many-particle environment is synchronized with the main optical cycle of the laser pulse, where the majority of charge carriers is excited, and response to the electric field becomes significantly nonlinear. As one can see from Eqs. (13), the band populations are determined by the integral of dephasing rate. The HOEID model has a smallest area under the curve, so it should have a lower excitation probability than the other models, while allowing for a temporally high dephasing rate. This hypothesis will be numerically tested in the next section.

V. NUMERICAL RESULTS AND DISCUSSION
In this section, we present the numerical simulations for a bulk α-SiO 2 interacting with the few-cycle IR laser pulse (Fig. 5a)  with four valence and four conduction bands. Both models give qualitatively similar results, but the TDHF model reduces the interband current and transient populations due to the coupling of density matrix elements at different k points and renormalization of interband interaction energy (18). On the other hand, the residual population has increased by nearly 12% in the TDHF simulation (see Fig. 5b). After the laser pulse, the current density predicted by IPA simulation demonstrates the rephasing effect (see Figs. 5c and 5d). The total interband coherence is partially restored and oscillates with a period of ∼ 2.5 fs. This effect is not observed in the TDHF simulation partially including electron-electron interaction.
In Fig. 6, we compare numerical simulations of the charge carrier population on the amplitude of the NIR pulse with the cos 4 envelope, where pure dephasing is described by the HOEID model (solid curve) with the other approximations: fully coherent TDHF equations, the EID model, and the Markovian constant decoherence rate. As expected, the simulation without dephasing follows the perturbative scaling law ∝ F 10 0 at low field amplitudes (F 0 < 0.8 V/Å). At higher fields, the scaling law changes due to closing of the lowest multiphoton channel predicted by the Keldysh theory and its modern generalization [19,[61][62][63]. By contrast, the numerical simulation with a constant dephasing time T 2 = 4 fs, which is required for reproduction of the experimental HHG spectrum, shows a quadratic scaling of the population in entire range of the field amplitudes. This corresponds to an artificial single-photon excitation channel due to spectral broadening, which is not typical for solids, As shown in Fig. 6, the unphysical scaling of excitation probability can be corrected by using the time-dependent dephasing rates. The outcome of HOEID model approaches closer to the coherent one at high field amplitudes F 0 0.8 V/Å, where the ponderomotive energy becomes sufficiently large (U p > ω 0 ) to overcome the spectral broadening introduced by pure dephasing. In the recent CEP current control measurements [65], scaling of the transferred charge Q ∝ F N 0 with powers smaller than the perturbative result N = 11 were observed. This observation can also be explained by a non-trivial dependence of pure dephasing rate on laser field and material parameters. A rigorous analysis of similar measurements and the high-harmonic generation spectroscopy with simulations based on Eq. (13) can be used for determining of the material-specific time-dependent dephasing rates. Finally, we compare the simulations of high harmonic spectra with different models of dephasing. Fig. 7 shows that simulation obtained with the constant dephasing time T 2 = 4 fs (dashed curve) still gives the best agreement with experimental results, where the cutoff is extended beyond 35 eV [5,13]. The HOEID model approaches closer to the Markovian result than the EID. Note that both the shape and intensity of high harmonics are sensitive to the time dependence of dephasing rate, which suggests that the high-harmonic spectroscopy can be used for reconstruction of time-dependent dephasing rates containing information on interaction with phonon and many-particle environments.

VI. CONCLUSIONS
To summarize, we developed the non-Markovian theory of pure dephasing in a dielectric excited by an ultrashort IR/visible laser pulse. It is shown that in the case of fast singleparticle states and slow environment the adiabatic separation of system and bath results in the time-dependent dephasing rate with a slowly-varying envelope defined by the bath spectral function and rapidly-varying part determined by system's interaction with an external field.
We studied both phonon and many-particle baths within the harmonic oscillator model and showed that the spectral function cutoff significantly changes time-dependent envelope of the dephasing rate as well as its peak and thermalized values. This explains why the phonon bath is approximately Markovian even on a few-femtosecond time scale and why the many-particle bath features unusually high values of the dephasing rate, which were not observed in the experiments with much longer laser pulses.
Numerical simulations show that the time-dependent dephasing rate with the envelope derived from the harmonic oscillator model significantly improves the problem of overestimated excitation probability at high intensities and allows for a temporally high dephasing rate, which is necessary for reproducing the experimental HHG spectrum of α-SiO 2 .