## 5 Discussion  
In case of the explicit scheme, a good guess for the time step requirement can be obtained from the [PERSON] stability criterion [20], \(\Delta t\leq\frac{(\Delta x)^{2}}{2D_{th}}\), where \(D_{th}\) is thermal diffusion coefficient, which is proportional to thermal conductivity \(k_{e}\) and inversely proportional to the carrier heat capacity \(C_{e-h}\). Under initial (prepulse) conditions, in the absence of free carriers, the latter tends to vanish (see Equation (10)), whereas the former is limited (see Table 1). After the laser irradiation starts, a quick increase of \(T_{e}\) at initially low \(n\) (see also Figure 2 in [25]) leads to an abrupt increase of \(D_{th}\), which influences the von Neumann stability criterion and limits the maximum possible time step for explicit integration methods. Consequently, if one applies an explicit finite-difference scheme for the numerical solution, the stability of Equation (8) limits the maximum possible time step to \(0.5\times 10^{-19}\,\mathrm{s}\) for the set of parameters published in [19]. Unfortunately, a mathematical error in [19] did not allow us to directly compare the results (specifically, in Equations (18) and (19) therein). For the numerical setup presented here, according to our test calculations, the explicit scheme requires time steps as small as \(1\times 10^{-24}\,\mathrm{s}\), making calculations too expensive. In contrast, the proposed semi-implicit numerical integration scheme provides a stable solution for time step as high as \(1\times 10^{-16}\,\mathrm{s}\) with the energy conservation about \(0.16\,\mathrm{\char 37}\) per simulation, Figure 3. At the same time, in case the calculation speed is critical, increasing the time step even higher is possible: \(1\times 10^{-15}\,\mathrm{s}\)  
Figure 2: Electron/lattice energy densities (divided by the melting energy density) and carrier density dynamics, according to the \(n\)TTM model, at the surface of silicon target of \(800\,\mathrm{nm}\) width, followed by the \(130\,\mathrm{fs}\) laser pulse at the incident fluence of \(0.26\,\mathrm{J/cm^{2}}\). The total energy density, averaged through the whole sample, is shown with black solid line. The laser pulse shape, shown with black dotted line, is not in scale. provides the energy conservation within 1.6 %. Thus, the increase in the calculation speed of up to \(10^{9}\) times has been achieved, compared with the explicit finite-difference integration scheme. Figure 4 shows the evolution of root mean square error (RMSE) in the electronic temperature for different time steps. The calculation with time step of \(1\times 10^{-20}\) s was used as a reference. Note that the chosen precision of the predictor-corrector (\(10^{-6}\) K) is relatively low and likely influences RMSE at low time steps, preventing further improvement in RMSE for time steps below \(\Delta t=1\times 10^{-17}\) s.  
The time step is of course limited by the characteristic times of the involved physical processes, such as laser pulse duration, electron-phonon coupling time, and carrier recombination time. As long as it is much smaller than those mentioned above, the presented integration scheme is tested to be unconditionally stable. Figure 4: Root mean square error in electronic temperature depending on time for different time steps. The calculation with time step of \(1\times 10^{-20}\) s was used as a reference. Figure 3: Maximal relative error in energy conservation as a function of simulation time for six different time steps. The quantity is calculated from the difference between total absorbed fluence and the total energy in the system. This approach has been successfully applied earlier in order to investigate and improve the presented \(n\)TTM model [24]. In [25], we used the described scheme for the solution of the continuum part of the atomistic-continuum model MD-\(n\)TTM. The atomistic-continuum model describing the dynamics of gold targets under the ultrashort-pulse lasers also benefited from using the presented approach [23].