## 3 Numerical Solution Scheme
The function \(f_{i}^{k}\) can be found from:  
\[(C_{e-h})_{i}^{k}f_{i}^{k}= (S_{u})_{i}^{k}-\frac{W_{i}^{k}-W_{i-1}^{k}}{\Delta z}-\frac{(C_{ e-h})_{i}^{k}}{(\tau_{e-p})_{i}^{k}}\left[T_{i}^{k}-(T_{a})_{i}^{k}\right] \tag{20}\] \[-\left(\frac{\partial n}{\partial t}\right)_{i}^{k}\left\{(E_{g} )_{i}^{k}+\frac{3}{2}k_{B}T_{i}^{k}\left[H_{\frac{1}{2}}^{\frac{3}{2}}(\eta_{e })+H_{\frac{1}{2}}^{\frac{3}{2}}(\eta_{h})\right]_{i}^{k}\right\}\] \[-n_{i}^{k}\left\{\left(\frac{\partial E_{g}}{\partial n}\right)_{ i}^{k}\left(\frac{\partial n}{\partial t}\right)_{i}^{k}+\left(\frac{ \partial E_{g}}{\partial T_{a}}\right)_{i}^{k}\left(\frac{\partial T_{a}}{ \partial t}\right)_{i}^{k}\right\}\] \[-\frac{3}{2}k_{B}T_{i}^{k}n_{i}^{k}\left(\frac{\partial n}{ \partial t}\right)_{i}^{k}\times\left\{\left(\frac{\partial\eta_{e}}{\partial n }\right)_{i}^{k}\left[1-\left(H_{\frac{1}{2}}^{\frac{3}{2}}(\eta_{e})\right)_{ i}^{k}\left(H_{\frac{1}{2}}^{-\frac{1}{2}}(\eta_{e})\right)_{i}^{k}\right]\right.\] \[+\left.\left(\frac{\partial\eta_{h}}{\partial n}\right)_{i}^{k} \left[1-\left(H_{\frac{1}{2}}^{\frac{3}{2}}(\eta_{h})\right)_{i}^{k}\left(H_{ \frac{1}{2}}^{-\frac{1}{2}}(\eta_{h})\right)_{i}^{k}\right]\right\},\]  
where \(W_{i}^{k}\) is defined according to Figure 1 (between the cells) and Equation (6):  
\[ \begin{split} W_{i}^{k}=\left\{(E_{g})_{i+\frac{1}{2}}^{k}+2k_{B} T_{i+\frac{1}{2}}^{k}[H_{0}^{1}(\eta_{e})+H_{0}^{1}(\eta_{h})]_{i+\frac{1}{2}}^{k} \right\}\times I_{i}^{k}\\ -(k_{e}+k_{h})_{i+\frac{1}{2}}^{k}\frac{T_{i+1}^{k}-T_{i}^{k}}{ \Delta z}\end{split}  \tag{21}\]  
Figure 1: The finite-difference grid mesh for the solution of the \(n\)TTM model equations. Symbol “\(\bullet\)” indicates the grid points for \(n\), \(Te\), \(T_{a}\), and \(E_{g}\) (\(i=1,2,...,N\)); symbol “\(\times\)” indicates the grid points for \(\vec{J}\), \(\vec{W}\), \(D\), and \(\ k_{a}\frac{\partial T_{a}}{\partial z}\) (\(j=1,2,...,N+1\)). Analogously, according to Figure 1 and Equation (4), the carrier current is  
\[ \begin{split} J_{i}^{k}=-D_{i}^{k}\Bigg{[}\frac{n_{i+1}^{k}-n_{i}^{k} }{\Delta z}+\frac{n_{i+\frac{1}{2}}^{k}}{k_{B}T_{i+\frac{1}{2}}^{k}}\left\{ \left[H_{-\frac{1}{2}}^{\frac{1}{2}}(\eta_{e})+H_{-\frac{1}{2}}^{\frac{1}{2}}( \eta_{h})\right]_{i+\frac{1}{2}}^{k}\right\}^{-1}\frac{(E_{g})_{i+1}^{k}-(E_{ g})_{i}^{k}}{\Delta z}\\ +\frac{n_{i+\frac{1}{2}}^{k}}{T_{i+\frac{1}{2}}^{k}}\left\{2\frac {\left[H_{0}^{1}(\eta_{e})+H_{0}^{1}(\eta_{h})\right]_{i+\frac{1}{2}}^{k}}{ \left[H_{-\frac{1}{2}}^{-\frac{1}{2}}(\eta_{e})+H_{\frac{1}{2}}^{-\frac{1}{2}} (\eta_{h})\right]_{i+\frac{1}{2}}^{k}}-\frac{3}{2}\right\}\left\{\frac{T_{i+1} ^{k}-T_{i}^{k}}{\Delta z}\right\}\end{split}  \tag{22}\]  
with  
\[D_{i}^{k}=\frac{k_{B}T_{i+\frac{1}{2}}^{k}}{q_{e}}\frac{\mu_{e}\mu_{h}(H_{ \frac{1}{2}}^{0}\left(\eta_{e}\right))_{i+\frac{1}{2}}^{k}(H_{\frac{1}{2}}^{0 }\left(\eta_{h}\right))_{i+\frac{1}{2}}^{k}}{\mu_{e}(H_{\frac{1}{2}}^{0} \left(\eta_{e}\right))_{i+\frac{1}{2}}^{k}+\mu_{h}(H_{\frac{1}{2}}^{0}\left( \eta_{h}\right))_{i+\frac{1}{2}}^{k}}\left[H_{-\frac{1}{2}}^{\frac{1}{2}} \left(\eta_{e}\right)+H_{-\frac{1}{2}}^{\frac{1}{2}}\left(\eta_{h}\right) \right]_{i+\frac{1}{2}}^{k}. \tag{23}\]  
Any function in between cells can be found by averaging  
\[A_{i+\frac{1}{2}}=\frac{1}{2}(A_{i}+A_{i+1}). \tag{24}\]  
The Fermi-Dirac integrals were calculated using GNU Scientific Library [40] and stored in the tables to speed up the calculations.