## 5 Numerical Analysis  
The numerical analysis will be based on the system of ordinary differential Equation (27) rewritten in terms of the corresponding dimensionless values. Figure 3: The density of the polarization current in a constant field (57) with parameters of Figure 2.  
We will investigate the response of the system to four electric field models: the constant electric field (57), the Eckart - Sauter field model  
\[E(t)=E_{a}\cosh^{-2}(t/T),\;\;A(t)=-E_{a}\tanh(t/T), \tag{58}\]  
the harmonic function with a constant amplitude  
\[E(t)=E_{a}\sin(\omega t),\;\;A(t)=\frac{E_{a}}{\omega}-\frac{E_{a}}{\omega} \cos(\omega t), \tag{59}\]  
where \(\omega\) is the angular frequency, and the \"laser field\" model [28]  
\[E(t) = E_{a}\cos(\omega t)\exp(-t^{2}/2\tau^{2}), \tag{60}\] \[A(t) = -\sqrt{\frac{\pi}{8}}E_{a}\tau\exp{(-\sigma^{2}/2)}\;\text{erf}( \frac{t}{\sqrt{2}\tau}-i\frac{\sigma}{\sqrt{2}})+c.c.,\]  
where \(\sigma=\omega\tau\). In all the cases in this section we assume that the electric field is directed along the first coordinate axis. We start with the most convenient model of the field (58). From Equation (27) follows that the speed of the filling process of the conduction band is determined by the amplitude of the transitions (10). In the denominator of (10) the quasienergy (8) takes zero values on the Dirac line (26). This feature of the amplitude of transitions should be reflected in the behavior and properties of the distribution function. From the form of the evolution of the vector potential (58) it follows that the Dirac line in this case should be represented in the momentum space by a segment with the endpoints determined by \(A(t\rightarrow-\infty)\) and \(A(t\rightarrow\infty)\) in accordance with the conditions (26). In Figure 4 we demonstrate the presence of such characteristic features of the distribution function. On the left panel the Dirac line has the end point coordinates \(p_{1}\mp 0.1\), \(p_{2}=0.0\) while on the right panel the pulse duration is five times larger so that the coordinates of the end points are \(p_{1}\mp 0.5\), \(p_{2}=0.0\). The Dirac line itself cuts in the distribution function a very thin canyon that is not visible in this figure owing to the selected scale of the numerical calculations. In the next step, we consider the constant field (57) at \(t\geq 0\). The distribution function at the time \(t=10.0\)\(t_{0}\) of the field action is presented on the left panel of Figure 5. Results of the field action with five times longer duration are represented on the right panel of Figure 5.  
Another frequently used model is the harmonic electric field (59). The procedure of switching on at \(t=0\) and off at \(t_{m}=2\pi m/\omega\) can be realized with sufficient accuracy in the numerical calculations. The shape of the distribution function and its change in time (\(m=1,2,4\) and \(10\)) for the field (59) are presented in Figure 6.  
Figure 4: The distribution function in the planar momentum space after the action of the Eckart-Sauter pulse (58). **Left panel**: \(E_{a}=0.01\)\(E_{0}\) (\(1.088\times 10^{8}\) V/m) and \(T=10\)\(t_{0}\) (\(2.46\times 10^{-15}\) s), **Right panel**: \(E_{a}=0.01\)\(E_{0}\) (\(1.088\times 10^{8}\) V/m) and \(T=50\)\(t_{0}\) (\(1.23\times 10^{-14}\) s). Finally, let us consider the more realistic model (60) of the \"laser pulse\". In this case the vector potential and the field strength are changed smoothly at any moment of time and do not bear any problems for the numerical calculations. The shape of the distribution function and its dependence on the pulse width determined by parameter \(\sigma\) are presented in Figure 2.0  
The above results correspond to very short time intervals from \(T=10\ t_{0}\ (2.46\times 10^{-15}\ \mathrm{s})\) to \(T=50\ t_{0}\ (1.13\times 10^{-14}\ \mathrm{s})\) of the electric field action for the models (57) and (58) and for a very high frequency of oscillations \(t_{0}/T=0.1(\approx\)400 THz) for the models (59) and (60).