Iteration of semiconductor Bloch equations for ultrashort laser pulse propagation

The numerical propagation of intense laser pulses through bulk material requires the recurrent calculation of the nonlinear material response. To describe the optical Kerr effect and the current in the conduction band for macroscopic propagation distances, very simplified models are typically used. Recent studies of the response of dielectrics to intense few-cycle pulses have revealed that ionization does not accumulate monotonically, but conduction bands are populated both irreversibly and reversibly during a laser cycle. The reversible (or transient or virtual) population of the conduction bands is not captured by simple response models. Here, an efficient iteration based on the semiconductor Bloch equations for three bands is developed, which consistently captures the laser cycle resolved interband polarization and intraband current. The full calculation of the nonlinear material response at each propagation step is avoided, instead only the incremental modification of the previous propagation step is calculated. The iteration is particularly well-suited for very short pulses and can be applied for intensities above the critical value at which perturbation theory does not converge. Furthermore, it is shown that virtual currents and dynamic Bloch oscillations are mechanisms which are missing in the Drude model, but these two mechanisms prevail for short intense pulses. Therefore, a generalized Drude model is derived from the SBEs, which is capable to account for arbitrary band shapes and both real and virtual ionization.


Introduction
The numerical methods employed in ultrafast nonlinear optics, especially in laser pulse propagation, have reached a degree of maturity at which standard tools are available [1]. Pulse propagation equations at various levels of approximations of Maxwell's equations are available, which can be used in tandem with material response models. The material response is typically divided into the nonlinear polarization (also called Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Kerr-related component) and the current (also called plasmarelated component). Usually, the two components of the model are independently parametrized, thus not reflecting the fact that these effects have a common origin and are therefore intimately related [2].
The ionization is typically modeled by a rate, and as such it neglects the history of the material response in the time dependent laser field and results in a monotonically growing density of charge carriers. This does not fit well to recent experiments of attosecond spectroscopy in dielectrics, which have revealed that the electron density grows in an oscillatory fashion, consisting of both an irreversible and a reversible component [3][4][5]. The irreversible component associated with a monotonically growing electron density prevails in the case of electron tunneling in atoms, which evolves in steps at the crests of the laser electric field strength [6]. By contrast, the reversible component associated with a transient population exchange between valence and conduction bands prevails in dielectrics [7]. The reversible (or transient or virtual) population of the conduction bands oscillates at a frequency twice that of the laser frequency. The oscillations of the virtual conduction band population can reach peak values well above the critical value at which optical damage occurs for real populations [8].
How to efficiently incorporate this type of material response into models for laser pulse propagation is the topic of this article. Simple propagation models, where the Drude model is used for the current and the ionization-related energy losses are modeled by a purely phenomenological current [2], seem inherently inappropriate, because only bound and free states are considered, whereas correlated electronic states (coherences) are neglected. The coherences, however, are indispensable for virtual populations. Their dynamics are driven by interband and intraband mechanisms. Interband dynamics result from coupling between the electronic energy bands of a crystal. Intraband dynamics, in contrast, refer to the motion of electrons and holes within the bands, which is described by the evolution of the electron density along the crystal quasi-momentum k. Realistic numerical modeling can be achieved on the basis of semiconductor Bloch equations (SBEs) [9]. This allows the calculation of the microscopic material response to a given electric field of a laser pulse, but it remains rather challenging to calculate the reshaping of the laser pulse during macroscopic propagation, where the material response needs to be calculated recurrently along the laser propagation direction. Coupling SBEs with Maxwell's equations has been demonstrated for micrometer propagation distances [10], but the thickness of bulk dielectrics used in experiments often exceeds 100 μm. The propagation effects over such distances are rather pronounced for strong-field pulses; recently the comparison of high-order harmonic generation (HHG) in transmission and reflection geometry revealed that propagation effects overwhelm the phase variation of the microscopic response [11].
Here, a method for the efficient iteration of SBEs is developed. The aim is not to describe the material as realistically as possible, but rather to create a model that is as simple as possible and still captures the mechanisms relevant to laser pulse propagation. The model contains three electronic bands, a valence band and a conduction band, which are coupled to a third band. It is shown that the mechanisms missing in the Drude model are virtual currents and dynamic Bloch oscillations, which describe phenomena that occur when the electrons are reflected at an edge of the Brillouin zone [12]. A generalization of the Drude model is derived from the SBEs, which is capable to account for arbitrary band shapes and both real and virtual ionization. It is shown that both real and virtual ionization seed the intraband current, but the influence of the virtual ionization prevails for very short pulses. The equations of motion of the model are solved by a mixed time-domain/frequency-domain iteration similar to perturbation theory, which is shown to be particularly well-suited for very short pulse durations. The full calculation of the nonlinear response at each propagation step is avoided, instead the nonlinear response at the previous propagation step is used as the initial value of the iteration. This allows to apply the iteration for intensities above the critical value at which perturbation theory does not converge. The determination of the free parameters of the model is discussed using fused silica as an example.

Fourier iteration of semiconductor Bloch equations
The SBEs consistently describe the coupled interband and intraband dynamics induced by a laser pulse [9]. Intraband dynamics refer to the changes that a light field induces in the crystal quasi-momentum k. In conventional optics, these changes are very small and can often be neglected. In the regime of extreme nonlinear optics, the field-induced shifts of the wave vector are very large and can exceed the width of the Brillouin zone [13]. On the other hand, the Coulomb interaction among the carriers, which is in the order of the exciton binding energy and important in conventional optics, can be regarded as a small perturbation in the regime of extreme nonlinear optics [14]. Omitting the Coulomb interaction, the SBEs in the length gauge read (atomic units are used) The diagonal elements ρ k nn of the density matrix are the k-dependent populations of the electronic bands with energies ω n , the off-diagonal elements ρ k nm (n = m) are the coherences between the states with transition energies ω nm = ω m − ω n . While the terms proportional to the dipole matrix elements d k nm describe interband transitions driven by the electric field E, the terms proportional to ∇ k lead to changes in k which represent the intraband transitions. The scattering terms ∂ t ρ k nm scatt summarize all couplings to higher-order correlations and will be treated in the Markov approximation as phenomenological damping terms.
The coordinate transform k → k − A, where A = −∂ t E is the vector potential, realizes the transition into a time dependent k-space, where dipole matrix elements and band energies become time dependent. Applying this transform to (1) yields which is the velocity gauge representation of the SBEs [10,15]. The time-domain integration of the SBEs is readily performed for the calculation of microscopic responses, but the repeated integration that is necessary for macroscopic propagation is mostly avoided. The aim of the following is to find an approximation that is as simple as possible but still reproduces the essence of the SBEs important for strong-field laser pulse propagation, such as transient conduction band populations and dynamic Bloch oscillations. The restriction to two bands, a valence and a conduction band, is often employed, but the consistent treatment of linear absorption and nonlinear refraction may require coupling to a third band, see section 4. Therefore, another band is coupled to the valence band (band 1) and to the conduction band (band 2) as depicted in figure 1. To avoid resonance effects that depend on the details of the band structure, ω k 3 = ω k 1 will be assumed later in the calculation. Alternatively, ω k 3 ω k 1 may be assumed. Only the conduction band energy is considered to be k-dependent (ω k 32 = ω 32 ; ω k 13 = ω 13 ; d k 32 = d 32 ; d k 13 = d 13 ); only off-diagonal dipole matrix elements are considered, which are assumed to be real ( d nn = 0 and d nm = d mn ). For systems with inversion symmetry, d 32 = 0 will be assumed.
For dielectrics, the population transfer into the conduction band is only a small fraction of the valence band population when irreversible material changes are avoided. For 50 fs pulses at 800 nm, a free-carrier density < 10 20 cm −3 was observed to cause damage in fused silica [16], which is small in proportion to the density of valence band electrons > 10 23 cm −23 . This justifies the assumptions ρ k which effectively decouples the diagonal and off-diagonal elements of the SBEs. Applying the approximations to (2), the equations for the off-diagonal elements read where T 2 is the interband dephasing time. Here it is assumed that T 2 is identical for all coherences and independent on k, but a more realistic treatment may be achieved using individual and k-dependent dephasing times [7,10,17,18]. Mutual coupling between all bands is considered in (3), because it can be used to model the generation of even harmonics for systems without inversion symmetry [19,20]. In most cases (for systems with inversion symmetry or when complex transition dipoles are used for even harmonic generation [18]) d 32 = 0 can be assumed, which eliminates the direct coupling of ρ 12 and ρ 13 , see section 4.
For an efficient calculation, where the elements D nm are defined by here, the time-dependent transition energy ω k− A nm is separated into a static part ω nm = ω k=0 nm and the dynamic part The transition energy and dipole matrix element between bands 1 and 3 is assumed to be constant as depicted in figure 1. When the elements D nm are set to zero, then (4) is the solution in the limit of linear optics. The elements D nm are the corrections due to the nonlinearity. A recursive method similar to perturbation theory is used for their calculation: In the nth step of iteration, the elements ρ nm in (4) are calculated using the elements D nm of the (n − 1)th step. Each step requires (inverse) Fourier transforms and time-domain multiplications. Three options for the initial value are considered here: In option 1, D nm = 0 is assumed in step 0. Using this option, the iteration is very similar to perturbation theory. Accordingly there is an upper limit for the electric field above which the iteration does not converge. In the limit of a monochromatic field with frequency ω 0 , the Keldysh parameter γ = ω 0 √ 2ω 12 | E| [21] appears in the convergence criterion: In (6), n denotes the step of iteration and the following simplifications have been applied in addition to the assumption of a monochromatic field: a two-band model ( d 13 = 0 and d 32 = 0); parabolic band evaluated at In this approximation, the Fourier iteration of (4) and (5) will not converge for γ 1. This is not surprising, because the Fourier iteration is a power-law expansion in orders of the electric field, and the Keldysh parameter describes the separation of the multiphoton regime (γ > 1) and the tunneling regime (γ 1). However, the convergence criterion is relaxed for short pulses, because the ratio decreases as the laser pulse is shortened. In the numerical implementation in section 4, it was found that the application of a filter in time domain and frequency domain (in the simplest case setting all ρ nm and D nm to zero above and below threshold values of time and frequency) is very important for fast convergence. The shorter the laser pulse, the better is the performance of the Fourier iteration. Option 2 of the iteration can be applied for intensities above the convergence limit of option 1. Option 2 is identical to option 1, but the electric field in step 0 is assumed to be lower than the convergence limit of option 1. In subsequent steps the electric field is slowly increased to the actualvalue.
Option 3 of the iteration can be applied for recurrent calculations of laser pulse propagation. If the propagation of an electric field in direction z is calculated with step-size dz, then the initial values of the elements D nm at position z 0 are taken from the prior calculation at z 0 − dz. At the beginning of the laser pulse propagation, either a time-domain integration of (3) or option 2 of the iteration can be performed. Option 3 can be applied for intensities above the convergence limit of option 1 if the change of the electric field per propagation step is sufficiently small. To fulfill this condition, the laser pulse propagation should be performed in the local pulse frame [1].
After the calculation of the off-diagonal density matrix elements, the diagonal elements can easily be calculated. Since the intraband current is zero for bands without curvature (bands 1 and 3), not all diagonal elements, but only the conduction band population needs to be calculated: The calculation can be performed as a simple time-domain integration. Damping is included via the lifetime T 1 and the collision time T c . The lifetime in conduction bands of dielectrics usually exceeds 100 fs, justifying the assumption T 1 = ∞. A damping mechanism that causes a decay of the currents, even though neglected in many recent studies, cannot be neglected, because Drude collision times are in the few-femtosecond range. In (7), the terms proportional to 1/T c preserve the total conduction band population but damp its asymmetry. This approach follows reference [20] and is in accordance with the qualitative picture of excited electrons that first undergo rapid momentum relaxation and thereafter energy relaxation on a longer timescale [22]. An alternative for faster calculations is to replace 22 , where indicates the average over all k-values.
From the density matrix elements, the polarization P and the current J, which are the source terms in Maxwell's equations, as well as the total conduction band population ρ can be calculated by integration over the Brillouin zone (BZ):

The limiting case of the Drude model
The connection between the SBEs and the Drude model is revealed in the length-gauge representation of (7): With the assumption T 1 = ∞, the application of BZ ∇ k ω k 2 to both sides of the equation yields (10) Periodic boundary conditions of the BZ have been assumed in order to carry out a partial integration of the first integral expression. The first integral describes the Drude current, which is seeded by real conduction band populations. The second integral describes the virtual current, seeded by the coherences, which is not captured by the Drude model. The link to the Drude model is most clearly revealed by considering a free electron gas with parabolic bands. In this case, ∇ 2 k ω k 2 = 1, and the first integral is simply the total conduction band population ρ. If the approximation of parabolic bands with an effective mass is applied to realistic bands, then one will find that the effective mass averaged over all charge carriers strongly depends on the amplitude of the laser pulse, as reported in reference [23].
The second integral vanishes when the coherences of the density matrix are symmetric with respect to k, specifically Im(ρ k nm ) ) = Im(ρ −k nm ) = for n = m, because ∇ 2 k ω k 2 = −∇ 2 k ω − k 2 and d k nm = d − k nm . The k-asymmetry is induced by the vector potential of the laser pulse and oscillates during the course of a laser cycle. Therefore, the second integral can be neglected when only cycle-averaged currents are considered, but not when subcycle currents due to virtual population are considered. In summary, for parabolic bands without population at the Brillouin zone edges, equation (10) for cycle-averaged currents reads This is the response model that is frequently used in laser pulse propagation and for the calculation of Brunel harmonics. It is often simply called the Drude model [1]. Strictly speaking, the formulation with a time-dependent density ρ goes beyond the traditional Drude model and was derived in two different ways. Brunel derived it from semiclassical electron trajectories in his original paper [24]. In the field of laser pulse propagation, it was derived from fluid equations [25,26]. The derivation from the SBEs seems to be missing so far. Comparison of (10) and (11) reveals the mechanisms that go beyond the Drude model: (i) the effective mass approximation breaks down; (ii) Bloch oscillations are not captured; (iii) virtual currents are not are not captured.

Numerical application for fused silica
To demonstrate a simple implementation, the SBE iteration is performed in one dimension of k and adapted to fused silica. Even though fused silica is not a crystal, its short range order is believed to prevail on the microscopic length scale of semiclassical electron trajectories [3]. A substantial simplification arises from the inversion symmetry, expressed as d 32 = 0. As a result, there is no direct coupling of ρ 12 and ρ 13 in (4) and (5), because D 12 does not explicitly depend on ρ 13 . Therefore it is a valid approximation to neglect the term F d k−A 12 Eρ k−A 13 , reducing equations (4) and (5) to and It suffices to calculate ρ k−A 13 after the iteration and not for each step of the iteration, because ρ k−A 13 does not appear on the righthand side of (12) or (13).
The exact relation between the dephasing time T 2 and the Drude collision time T c is not known. Experimental studies were reported where very short dephasing times < 10 fs yield the best agreement of models and data [20,[27][28][29], which is increasingly discussed in recent theoretical studies [7,10,30]. A Drude collision time as low as 0.4 fs was reported for multiphoton ionization of fused silica [16], but values used for pulse propagation vary broadly in the region 1-100 fs [31]. A consistent theoretical model should also reproduce the direct measurement of the electron drift mobility, where a temperature dependent mean free path < 4 nm was reported for fused silica [32]. The phenomenological picture suggests that if scattering occurs to an electron at position k, its interband-and intrabandcoherences are likewise destroyed, which is expressed by the relation T 2 = 2T c . With the additional assumption T 1 = ∞, equation (7) becomes For the conduction band, a tight-binding band structure is assumed with a bandwidth Δ = 4 eV and a lattice constant a = 0.49 nm. These values are consistent with band structure  calculations of SiO 2 [33]. Although there is no long range order in amorphous silicon dioxide, its electronic structure is similar to the electronic structure of crystalline quartz [34]. The dipole matrix elements are calculated as [9] d k 12 = d k=0 For laser pulse propagation it is essential that the calculation of the macroscopic polarization and current is correct on the quantitative level, rather than only qualitatively. Since the three-band approximation does not reflect all contributing electrons, the macroscopic observables are scaled to an effective electron number density N e by performing the integration in (8) as BZ → N e k δk, where δk is the spacing in the kgrid. For the following calculations, a k-grid with 29 points and a t-grid with 50 001 points in the interval [−150, 150] fs are used. After these considerations, six parameters remain to be determined: N e , T 2 , ω k=0 12 , d k=0 12 , ω 3 and d 13 . While the former four parameters affect both the linear and the nonlinear response, the latter two affect only the nonlinearity. To first adapt the linear response, the linear susceptibility is calculated from (12) by setting D 12 = 0, which in turn determines the complex refractive index n = 1 + 4πχ (1) .
The position of the absorption edge determines ω k=0 12 , which is set to 10.5 eV by comparison with numerical refractive index data (see figure 1). The absorption in the region [10,15] eV is affected by both N e and d k=0 12 , but only marginally by T 2 . A rather good agreement with the reference data is achieved for N e d k=0 12 2 = 0.035. In the following d k=0 12 = 0.1 and accordingly N e = 3.5 is used. The numerical results discussed below depend only weakly on d k=0 12 or N e as long as the product N e d k=0 12 2 is fixed. The dephasing time T 2 determines the broadening of the absorption below the bandgap. The calculation and the reference data cannot be matched perfectly, which is reasonable considering that the Sellmeier equations for dielectrics typically include two resonances in the ultraviolet and one resonance in the infrared. T 2 = 10 fs reproduces the refractive index reasonably well and will be used in the following to model pronounced virtual populations. A second set of parameters (d k=0 12 = 0.1; N e = 0.326; T 2 = 1 fs) is used for contrasting simulations with weak virtual populations. Both parameter sets, however, strongly overestimate the absorption in the optical regime. The transmission trough 1 mm at 800 nm is only 0.024% for the first and 0.003% for the second parameter set. Calculations for realistic laser pulse propagation should reproduce the linear response as precisely as possible. Therefore, it is advisable to separate linear and nonlinear polarization P (NL) in the calculation, and use numerical refractive index data for the linear polarization and P (NL) from the three-band calculation. In practice, the linear polarization is incorporated through the refractive index in pulse propagation equations [1].
The most common measure for the nonlinear response is the nonlinear refractive index n 2 . The assessment of more than 30 literature values of n 2 [35], mostly obtained using the zscan technique [36], suggests n 2 = 2.73 × 10 −20 m 2 W −1 at 800 nm. With a linear refractive index of 1.45, this corresponds to the nonlinear susceptibility χ (3) = 4.27. Since the literature values have been obtained using long laser pulses where 18% of χ (3) are due to the delayed Raman response [37], the effective value for the nonlinear electronic response is χ (3) = 3.50 (that is 1.66 × 10 −22 m 2 V −2 in SI units). The nonlinear electronic polarization, responsible for the optical Kerr effect (OKE), is usually modeled as an instantaneous response, P (NL) = χ (3) E 3 , which is used here to determine the interband coupling parameters. The restriction to two bands by setting d 13 = 0 matches the instantaneous response model with χ (3) = 2.14, which severely underestimates the OKE (not shown). Therefore, coupling to a third band needs to be included.  pronounced for a pulse at 1600 nm, where P (NL) of the threeband model is slightly diminished in comparison with the instantaneous response model. This wavelength-scaling is consistent with the literature values [35].
The performance of the three-band model with regard to multiphoton absorption is tested by a scan of the central laser frequency ω 0 , see figure 3. As expected, the final conduction band population after the laser pulse ρ(t = ∞) peaks when multiples of ω 0 equal the transition energy ω 12 . The base level below the multiphoton peaks is the tail of the one-photon absorption process. The shorter the pulses are, the broader the multiphoton peaks are and the lower their height. For high orders of absorption, the multiphoton peaks are buried in the one-photon absorption tail. The shorter the pulses are, the less pronounced are the multiphoton peaks. This behavior is consistent with the fact that the critical field strength for convergence of perturbation theory scales proportional to the square root of the pulse duration, see section 2. For the 5 fs pulse, the Fourier iteration shows almost perfect agreement with the time-domain integration after 3 iterations.
The conduction band density ρ k 22 transiently peaks during the cycles of the laser electric field and is shifted in k-space because of the interband coupling ( figure 4). The transient peaks are distinct for longer dephasing times T 2 and washedout for shorter dephasing times. The interband coupling causes a current, which can be divided into the Drude current and the virtual current. For the calculation of the individual contributions to the total current, equation (10) is used. As displayed in figure 4, the Drude current and the virtual current are of similar magnitude and largely cancel out for long dephasing times, but the Dude current becomes dominant for short dephasing times. The Fourier iteration with only 1 iteration yields already good agreement with the time-domain integration in terms of the conduction band population, but at least two iterations are required to reproduce the current ( figure 6).
Laser pulse propagation over macroscopic distances requires the recurrent calculation of the polarization and current. The calculation using the Fourier iteration method is much faster compared to the time domain integration and can readily be applied to macroscopic distances of multiple hundreds of micrometers. This is demonstrated here in one spatial dimension (the propagation direction z) by integrating the equation where c is the speed of light. This equation results from any carrier-resolving unidirectional propagation equation [1] when the lateral spatial dimensions are neglected. The propagation through 200 μm of fused silica (step-size 0.25 μm) is performed using the split-step method. For the laser pulse propagation, a coarser t-grid with 6001 points in the interval [−250, 250] fs is used. The larger interval is needed because of pulse broadening, but the grid spacing may be much wider for the Fourier iteration compared to the RK4 integration. At the first spatial step, option 2 of the iteration is used to calculate P (NL) and J, because the laser intensity (10 TW cm −2 ) is chosen above the convergence limit. At subsequent steps, option 3 with 30 iterations is used. Numerical refractive index data is used for the linear polarization. The calculation is done with a time-frame that travels at the group velocity at 800 nm. The two parameter sets lead to different laser pulse shapes ( figure 7). The parameter set for distinct transient populations (T 2 = 10 fs) shows more spectral broadening but less absorption. This is understandable because the dephasing time was linked to the Drude collision time.

Conclusion
Based on the SBEs for three bands, equations of motion were derived that can be regarded as the minimal-effort equations for a consistent calculation of the interband polarization and the interband current. The resulting interband current includes virtual currents and the influence of dynamical Bloch oscillations in addition to what is expected from the Drude model. To avoid the expansive time-domain integration, equations for a Fourier iteration were presented, which can be applied for intensities above the critical value at which perturbation theory does not converge. The iteration is particularly well-suited for laser pulse propagation, because the full calculation of the nonlinear response at each propagation step is avoided, instead only the incremental modification of the previous propagation step is calculated. This might open new opportunities, because strong-field propagation effects may be harnessed for novel phase-matching schemes [38]. Furthermore it will enable the interpretation of spectroscopic data obtained with methods where propagation effects cannot be neglected, such as noncollinear spectroscopy [27,39].