## 3 Numerical Solution Scheme  
To solve the system of Equations (7) to (9), we utilize the finite difference grid mesh presented in Figure 1. Sample is divided into cells as indicated, and the local thermodynamic parameters are calculated in each cell. The spatial derivatives of \(n\), \(T_{e}\), \(T_{a}\), \(J,W,k_{a}\frac{\partial T_{a}}{\partial z}\), and \(E_{g}\) at the interior points are approximated with the central differences, and those at the boundaries are evaluated with the first-order approximation. Equations (7) and (9) are initially solved explicitly (\(T\equiv T_{e}\)):  
\[\frac{n_{i}^{k+1}-n_{i}^{k}}{\Delta t}+\frac{l_{i}^{k}-l_{i-1}^{k}}{\Delta x} =(S_{n})_{i}^{k}-\gamma(n_{i}^{k})^{3}+\delta_{i}^{k}n_{i}^{k}, \tag{17}\] \begin{table}
\begin{tabular}{l l r} \hline \hline
**Parameter Name** & **Value** & **Citation** \\ \hline Initial carrier density & \(n_{0}=1\times 10^{16}\,\text{m}^{-3}\) & [29] \\ \hline Initial lattice and carrier temperature & \(T_{0}=300\,\text{K}\) & \\ \hline Lattice specific heat & \(C_{a}\ =\ 1.978\times 10^{6}\) + \(3.54\times 10^{2}T_{a}-3.68\times 10^{6}/T_{a}^{2}\), \(\text{J}/(\text{m}^{3}\text{K})\) & [30] \\  & (\(T_{a}\) in K) & \\ \hline Lattice thermal conductivity & \(k_{a}=1.585\times 10^{5}\times T_{a}^{-1.23}\), \(\text{W}/(\text{m}\cdot\text{K})\) (\(T_{a}\) in K) & [30] \\ \hline Carrier thermal conductivity & \(k_{e}=k_{h}=-3.47\times 10^{18}\)+ \(4.45\times 10^{16}T_{e}\), \(\text{eV}/(\text{s m K})\) & [31] \\ \hline Indirect band gap & \(E_{g}=1.170-4.73\times 10^{-4}T_{a}^{2}/(T_{a}+636)-1.5\times 10^{-10}n^{1/3}\) & \\  & if \(1.170-4.73\times 10^{-4}T_{a}^{2}/(T_{a}+636)-1.5\times 10^{-10}n^{1/3}\geq 0\) & [32] \\  & and 0 otherwise, eV (\(T_{a}\) in K, \(n\) in m\({}^{-3}\)) & [33] \\ \hline Interband absorption & \(\alpha=1.34\times 10^{5}\exp{(T_{a}/427)}\), \(\text{m}^{-1}\) & [34] \\ (taken from 694 nm laser) & \(\alpha=1.34\times 10^{5}\exp{(T_{a}/427)}\), \(\text{m}^{-1}\) & [34] \\ \hline Two-photon absorption & \(\beta=15\,\text{cm}/\text{GW}\) & [25] \\ \hline Reflectivity & \(R=0.329+5\times 10^{-5}(T_{a}-300)\) (\(T_{a}\) in K) & [35] \\ \hline Auger recombination coefficient & \(\gamma=3.8\times 10^{-43}\), \(\text{m}^{6}/\text{s}\) & [36] \\ \hline Impact ionization coefficient & \(\delta=3.6\times 10^{10}\exp{\left(-1.5E_{g}/k_{B}T_{e}\right)}\), \(\text{s}^{-1}\) & [37] \\ \hline Free-carrier absorption cross section & \(\Theta=2.91\times 10^{-22}T_{a}/300\), \(\text{m}^{2}\) (\(T_{a}\) in K) & [38] \\ \hline Electron-phonon relaxation time & \(\tau_{e-p}=0.5\times 10^{-12}\left[1+n/(2\times 10^{27})\right]\), \(\text{s}\) (\(n\) in \(\text{m}^{-3}\)) & [31] \\ \hline Electron effective mass & \(m_{e}^{*}=0.36m_{e}\) & [39] \\ \hline Hole effective mass & \(m_{h}^{*}=0.81m_{e}\) & [39] \\ \hline Mobility of electrons (taken at 1000 K) & \(\mu_{e}=0.0085\,\text{m}^{2}/\text{V}\cdot\text{s}\) & [38] \\ \hline Mobility of holes (taken at 1000 K) & \(\mu_{h}=0.0019\,\text{m}^{2}/\text{V}\cdot\text{s}\) & [38] \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Model parameters. \[(C_{a})_{i}^{k}\frac{(T_{a})_{i}^{k+1}-(T_{a})_{i}^{k}}{\Delta t}= \frac{1}{(\Delta z)^{2}}\left[(k_{a})_{i+\frac{1}{2}}^{k}(T_{i+1}^ {k}-T_{i}^{k})-(k_{a})_{i-\frac{1}{2}}^{k}(T_{i}^{k}-T_{i-1}^{k})\right] \tag{18}\] \[+\frac{(C_{e-h})_{i}^{k}}{(\tau_{e-p})_{i}^{k}}\left[T_{i}^{k}-(T _{a})_{i}^{k}\right],\]  
where index \(i\) is connected to the cell number (see Figure 1) and \(k\) to the moment of time. Therefore, before solving Equation (8) we already have the predictions for \(n^{k+1}\) and \((T_{a})^{k+1}\). The approach is based on the Crank-Nicolson semi-implicit scheme [21, 22]. Equation (8) can be rewritten in the following finite-difference form:  
\[\frac{T_{i}^{k+1}-T_{i}^{k}}{\Delta t}=(1-\psi)f_{i}^{k}+\psi f_{i}^{k+1}. \tag{19}\]  
The right-hand side contains parameter \(\psi\), which can be 0 for explicit scheme, 1 for implicit, and \(\frac{1}{2}\) for semi-implicit.