# Space Charge Accumulation and Decay in Dielectric Materials with Dual Discrete Traps
## 2 Dynamic Responses of Charge Trapping and De-trapping
### Instantaneous Charge Injection  
When a certain fixed number of free charges \(q_{total}\) is assumed instantaneously, it is difficult to obtain the analytical solutions for the charge trapping/de-trapping Equations (3) and (4). Nevertheless, when trapped charge densities, \(q_{trap1}(t)\), are always much smaller than the densities of trap centers in a dielectric material, \(q_{trap1}(t)\)/\(Q_{Ti}\) can be approximated by zero, and the ordinary differential equations can be solved analytically. The first order charge trapping/de-trapping equations in the dielectric materials containing two discrete trap levels can be simplified as the following equations. \[\frac{\mathrm{d}q_{trap1}}{\mathrm{d}t}\ =\ -(P_{T1}+P_{D1})q_{trap1}-P_{T1}q_{trap2}+P_{T1}q_{ total} \tag{8}\]  
\[\frac{\mathrm{d}q_{trap2}}{\mathrm{d}t}\ =\ -P_{T2}q_{trap1}-(P_{T2}+P_{D2})q_{trap2}+P_{T2}q_{ total} \tag{9}\]  
Given \(P_{Ti}\), \(P_{Di}\), and \(q_{total}\), Equations (8) and (9) are linear system of differential equations with constant coefficients. The undetermined function \(q_{trap2}(t)\) and the first-order derivative \(\mathrm{d}q_{tmp2}(t)\)/\(\mathrm{d}t\) can be expressed as functions of \(q_{trap1}(t)\), \(\mathrm{d}q_{trap1}(t)\)/\(\mathrm{d}t\), and \(\mathrm{d}^{2}q_{trap1}(t)\)/\(\mathrm{d}t^{2}\) by Equation (8). We can eliminate \(q_{trap2}(t)\) and \(\mathrm{d}q_{trap2}(t)\)/\(\mathrm{d}t\) in Equation (9) to obtain the linear differential equation of second order with constant coefficients of \(q_{trap1}(t)\). Similarly, we can find the second order linear differential equation with constant coefficients of \(q_{trap2}(t)\). \[\frac{\mathrm{d}^{2}q_{trap1}}{\mathrm{d}t^{2}}+(P_{T1}+P_{D1}+P_{T2}+P_{D2}) \frac{\mathrm{d}q_{trap1}}{\mathrm{d}t}+(P_{T1}P_{D2}+P_{T2}P_{D1}+P_{D1}P_{D 2})q_{trap1}-P_{T1}P_{D2}q_{total}\ =\ 0 \tag{10}\]  
\[\frac{\mathrm{d}^{2}q_{trap2}}{\mathrm{d}t^{2}}+(P_{T1}+P_{D1}+P_{T2}+P_{D2}) \frac{\mathrm{d}q_{trap2}}{\mathrm{d}t}+(P_{T1}P_{D2}+P_{T2}P_{D1}+P_{D1}P_{D 2})q_{trap2}-P_{T2}P_{D1}q_{total}\ =\ 0 \tag{11}\]  
The second order linear differential Equations (10) and (11) have the same characteristic equation. \[r^{2}+(P_{T1}+P_{D1}+P_{T2}+P_{D2})r+P_{T1}P_{D2}+P_{T2}P_{D1}+P_{D1}P_{D2}\ =\ 0 \tag{12}\]  
Solving the characteristic Equation (12), we can obtain two unequal real roots, namely \(\lambda_{1}\) and \(\lambda_{2}\). Then, we can solve the second order linear inhomogeneous differential Equations (10) and (11). The time dependent trapped charge densities \(q_{trap1}(t)\) and \(q_{trap2}(t)\) are given by,  
\[q_{trap1}(t) =\ C_{11}\exp(\lambda_{1}t)+C_{12}\exp(\lambda_{2}t)+P_{T1}P_{D2}q _{total}/(P_{T1}P_{D2}+P_{T2}P_{D1}+P_{D1}P_{D2}) \tag{13}\] \[q_{trap2}(t) =\ C_{21}\exp(\lambda_{1}t)+C_{22}\exp(\lambda_{2}t)+P_{T2}P_{D1}q _{total}/(P_{T1}P_{D2}+P_{T2}P_{D1}+P_{D1}P_{D2}) \tag{14}\]where \(C_{11}\), \(C_{12}\), \(C_{21}\), and \(C_{22}\) are coefficients dependent on the initial conditions. The analytical solutions of coefficients \(C_{11}\), \(C_{12}\), \(C_{21}\), and \(C_{22}\) are shown in Appendix A. Substituting Equations (13) and (14) into \(q_{free}(t)=q_{total}-q_{trap1}(t)-q_{trap2}(t)\), we can calculate the time dependent free charge density. \[q_{free}(t)\ =\ P_{D1}P_{D2}q_{total}/(P_{T1}P_{D2}+P_{T2}P_{D1}+P_{D1}P_{D2})-(C_{ 11}+C_{21})\exp(\lambda_{1}t)-(C_{12}+C_{22})\exp(\lambda_{2}t) \tag{15}\]  
In the case of constant total charge density, the time dependent effective carrier mobility, \(\mu_{eff}(t)\), has a similar function as \(q_{free}(t)\). \[ \begin{array}{ll}\mu_{eff}(t)=&\mu_{0}P_{D1}P_{D2}/\left(P_{T1}P_{D2}+P_{T2} P_{D1}+P_{D1}P_{D2}\right)\\ &-\mu_{0}[(C_{11}+C_{21})exp(\lambda_{1}t)-(C_{12}+C_{22})\exp(\lambda_{2}t)] /q_{total}\end{array}  \tag{16}\]  
When the approximation, \(q_{trapi}(t)/Q_{Ti}\approx 0\), is not valid for \(i=1\) or \(2\), alternatively, we can solve the charge trapping/de-trapping equations numerically by an unconditionally stable nonstandard finite differential (NSFD) method [29].