Violation of Neumann Problem’s Solvability Condition for Poisson Equation, Involved in the Non-Linear PDEs, Containing the Reaction-Di ﬀ usion-Convection-Type Equation, at Numerical Solution by Direct Method

: In this paper, we consider the 3D problem of laser-induced semiconductor plasma generation under the action of the optical pulse, which is governed by the set of coupled time-dependent non-linear PDEs involving the Poisson equation with Neumann boundary conditions. The main feature of this problem is the non-linear feedback between the Poisson equation with respect to induced electric ﬁ eld potential and the reaction-di ﬀ usion-convection-type equation with respect to free electron concentration and accounting for electron mobility (convection’s term). Herein, we focus on the choice of the numerical method for the Poisson equation solution with inhomogeneous Neumann boundary conditions. Despite the ubiquitous application of such a direct method as the Fast Fourier Transform for solving an elliptic problem in simple spatial domains, we demonstrate that applying a direct method for solving the problem under consideration results in a solution distortion. The reason for the Neumann problem’s solvability condition violation is the computational error’s accumulation. In contrast, applying an iterative method allows us to provide ﬁ nite-di ﬀ erence scheme conservativeness, asymptotic stability, and high computation accuracy. For the iteration technique, we apply both an implicit alternating direction method and a new three-stage iteration process. The presented computer simulation results con ﬁ rm the advantages of using iterative methods.


Introduction
In the current study, we consider a computer simulation of optical pulse propagation in a semiconductor, which is described by a set of coupled non-linear 3D PDEs. This set involves time-dependent equations describing the evolution of electron and ionized donor concentrations, the equation describing the optical pulse evolution, and the Poisson equation with Neumann boundary conditions (BCs) with respect to the laser-induced electric field potential [1][2][3]. As it is well-known, when solving the Poisson equation, it is necessary to account for the Neumann problem's solvability condition [4,5], which is a consequence of the Gauss theorem [6]. In the problem under consideration, there is a mutual non-linear influence of the charged particles' concentrations and the electric field potential because accounting for the electron mobility and describing the evolution of the electron concentration via the reaction-diffusion-convection-type equation. Because of the problem's nonlinearity, the accuracy of the Poisson equation's solution is a key issue, and our research aims to address the choices of numerical methods for solving it by comparing the direct method (DM) and the iteration method (IM).
The most popular DM for solving linear elliptic problems stated in the domain of simple (rectangular) geometry is the Fast Fourier Transform (FFT) algorithm. After its development in the second half of the 20th century [7,8], it has been applied ubiquitously in almost every scientific discipline because of its accuracy and cost-effectiveness [9][10][11][12][13][14][15][16][17][18]. Moreover, there are many commercially available packages that implement the FFT algorithm. In the current study, we use Intel MKL Fast Poisson Solver Routines for solving the Poisson equation with Neumann-type inhomogeneous BCs in a 3D case.
However, to our best knowledge, no study has investigated the effectiveness and applicability of DM if the elliptic equation with Neumann BCs is involved in a set of nonlinear time-dependent PDEs, which implies the mutual feedback between the Poisson equation solution and a solution of the reaction-diffusion-convection-type equation. As demonstrated below, due to the violation of the Neumann problem's solvability condition, the DMs (FFT) become ineffective. It is important that the Neumann problem's solvability condition coincides with the conservation law for the problem under consideration. Therefore, the problem's invariant is not preserved with high accuracy at using DM, and its value grows over time. Consequently, the Neumann problem's solvability condition is also violated, and the numerical solution does not converge with those of the differential problem.
To overcome this, we propose to solve the Neumann problem for the Poisson equation using IM. It is known that most IMs are unstable in 3D cases. The stable method for 3D problems, named the implicit alternating direction method (ADM), was proposed almost simultaneously by Peaceman, Rachford, and Douglas [19,20]. It is applied as a finitedifference scheme (FDS) to the parabolic equation and as an iteration technique for solving stationary elliptic problems. Further, in [21], the implicit ADM for 3D problems, which is unconditionally stable and has the second-order approximation both on time coordinates and spatial coordinates, was developed. Let us note that nowadays, ADM and its modifications are extensively used for solving multidimensional non-linear problems [22][23][24][25][26][27][28]. We also propose another IM both for solving the 3D Poisson equation and for the implementation of the conservative non-linear FDS developed for the numerical solution of other equations of the problem based on Crank-Nicolson-type FDS-the original threestage iteration process (TSIP) [29]. The main advantages of TSIP consist of providing conservativeness property on the iterations and asymptotic stability of FDS. These are crucial features for solving non-linear time-dependent problems.
Below, we compare the effectiveness of using the DM (FFT) and IM (ADM or TSIP) for solving the Neumann problem for the Poisson equation, involved in the set of PDEs, which contains a convection-type term. It should be emphasized that the obtained results and conclusions are regardless of the particular method's implementation and have a general character because they arise from the solvability condition's violation caused by the round-off errors accumulation at using DM. It should be emphasized that in [30], the crucial influence of a convection-type term related to electron mobility on the computation's effectiveness was demonstrated, and at zero-value electron mobility, both methods (DM and IM) possessed the same high accuracy. However, in [30], the problem with homogeneous BCs or a particular case of the external electric field (EEF) action along the direction that is transverse to the optical pulse propagation direction was only considered. On the contrary, below, we study the problem with inhomogeneous BCs corresponding to the action of EEF in one or two directions, including the optical pulse propagation direction. This paper is organized as follows. In Section 2, we present a mathematical model of the optical pulse propagation in a semiconductor as well as the problem of funding the initial function's distribution under the condition of the inhomogeneous BCs. Section 3 is devoted to the numerical methods of the problem's solution. Here, we discuss the DM and IM for a solution to the Neumann problem for the Poisson equation. In Section 4, extensive computer simulation results are given to illustrate the theoretical propositions regarding the effectiveness of the methods for the Poisson equation solution and the influence of EEF on the applicability of these methods. Some concluding remarks are presented in Section 5.

Semiconductor Plasma Evolution Equations
The 3D problem of semiconductor optical-induced plasma generation is governed by a set of well-known time-dependent, dimensionless, non-linear PDEs [1]: Above, the following notations are used: variables x, y, and z are the dimensionless spatial coordinates, and the parameters x L , y L , and z L denote their maximal values, respectively. Variable t represents a dimensionless time unit. The optical pulse propagates along the z-coordinate, and it is named a longitudinal coordinate; the x-coordinate and ycoordinate are transverse ones. Function ( , , , ) x y z t ϕ is a potential of the optical-induced electric field, function ( , , , ) N x y z t is an ionized donor concentration, function ( , , , ) n x y z t is a free electron concentration in the conduction band of a semiconductor, A semiconductor is assumed to be a rectangular parallelepiped: , , The BCs for the Equation (2) are the following: Constants x E , y E , and z E are responsible for the absence (in the case of their zerovalue) or presence of EEF along the corresponding spatial coordinate. The conditions (6) correspond to the absence of an electric current through the semiconductor faces. The intensity distribution of the incident optical radiation has a Gaussian profile of the beam and time-dependent pulse shape: In the initial time moment, the optical pulse intensity is equal to zero: The initial distributions for other functions are defined by the BCs, and they are discussed below in Section 2.2.
The absorption coefficient fundamentally influences the regimes occurring in a semiconductor. Below, we consider its non-linear dependence on free-charged particles concentrations [31][32][33][34]: Such an absorption coefficient can result in occurring OB, called a concentration OB. The following functions describe the generation and the recombination of free-charged particles: where a parameter 0 q characterizes the incident optical pulse maximal intensity, the parameter eq n denotes an equilibrium value of the free-charged carrier concentration, and the parameter R τ characterizes a free-electron recombination time.

Remark 1.
As it is known, the Neumann problem is solved only up to constant, and the multiplicity of its solution occurs. So, changing the electric field potential on a constant does not influence the problem's (1)- (8) solution if the absorption coefficient is defined by (9). However, if the absorption coefficient depends on the electric field potential, for example, as follows: which takes place if a semiconductor is undergoing the action of a high intensive laser pulse [35][36][37][38], it is necessary to provide re-normalization of the electric field potential computed by solving the Neumann problem for the Poisson equation [39].
We want to stress that due to the presence of electron mobility (coefficients x μ , y μ , z μ ) in the term describing convective transfer in Equation (3), the non-linear feedback between the charged particle concentrations and the electric field potential occurs: the electric field potential influences the free electron concentration evolution and vice versa. Evidently, this non-linear feedback is "switched off" at zero-value electron mobility. In [30], we demonstrated that the presence of such a term plays a crucial role in the Neumann problem's solvability condition violation using DM. Let us notice that the existence and uniqueness of the problem (1)- (11) solution was discussed in [40], where the laser pulse propagation was described by the non-linear Schrödinger equation, which is more complicated than Equation (4). We described how to prove the solution's existence. However, it should be taken into account that the uniqueness of the problem solution can be absent in general cases because the laser beam interaction with the semiconductor under certain conditions possesses the property of optical bistability, which means that for one value of the beam intensity, there are two stable states of electron concentrations and ionized donors' concentrations.

The Problem of Finding the Semiconductor Characteristics' Initial Distributions
For the problem's numerical solution, it is necessary to compute the initial distributions of the free-electron and ionized-donor concentrations as well as a spatial distribution of the electric field potential. For this purpose, we solve the PDEs (1)-(3) in a steady-state at zero-value optical pulse intensity. In this case, the free-charged carriers' generation is absent, and the stationary problem is written in the following way: with BCs (5)- (6). We denote the solution of the problem (12)-(14), (5)-(6) as: Let us note that if the EEF is absent ( )

The Charge Conservation Law and Its Correlation with the Neumann Problem's Solvability
For the problem under consideration, the charge conservation law occurs if the charge flows through the semiconductor surfaces (conditions (6)) are absent.
Because a semiconductor is electrically neutral before the optical pulse action, the problem invariant is equal to zero: ( ) (0) 0. Q t Q = = We follow this invariant when providing the computation because it is crucial for a computer simulation of the process: changing the invariant ( ) Q t leads to no correspondence of the numerical solution of the problem to its physical one. In turn, as it is wellknown [4,5], the necessary and sufficient condition of the solvability for the Neumann problem for the Poisson equation is the integral equality, which follows from the Gauss theorem [6]: where ( , , ) f x y z is the right part of the Poisson equation defined in the domain Ω , and ( , , ) g x y z is the right part of Neumann's BCs. For the problem (1), (5), condition (18) takes the form: where So, the right part of (19) is equal to zero if the BCs (5) are stated, and therefore, the Neumann problem's solvability condition (18) coincides with the charge conservation law (17). Thus, the FDS's conservativeness is the necessary condition for providing the Neumann problem (1), (5) solvability condition.

Conservative Finite-Difference Scheme and Three-Stage Iteration Process for Its Implementation
To solve the set of time-dependent non-linear PDEs in the rectangular 3D domain, we develop the FDS. For this aim, we introduce uniform time and spatial meshes: The mesh function h I is defined on the mesh ˆ′ ′ Ω with shifted nodes both on spatial coordinate z and on time coordinate: Approximation of the beam intensity on a shifted mesh is caused by the FDS secondorder accuracy reaching. For brevity, we use below the index-free notations: where f is one of the mesh functions , , h h h n N ϕ , and the following index-free notations for the mesh function h I : We also introduce other notations in the following way: The first and the second difference derivatives are defined in a standard way, and they are notated as follows: , , , , , , , , , ence operator is also stated in a standard way: For brevity, we introduce the finite-difference operators: and similar operators on other spatial coordinates.
Based on the Crank-Nicolson scheme and using introduced notations, we write the FDS for the problem (1)-(8): The BCs are approximated as follows: The initial intensity distribution is approximated in the following manner: Let us note that the FDS (20)-(24) has the second order of approximation in the inner mesh nodes, while the BCs are approximated with the first order. This inconsistency is caused by reaching the FDS conservativeness [29]. The difference analog of the invariant (17) is computed by using the formula: which possesses the second order of approximation with respect to the mesh node ( ) , assuming the solution of the differential problem is smooth enough.

Three-Stage IM for FDS Implementation
As it is well-known, the splitting methods are widely used for the numerical solution of multidimensional linear and non-linear problems. Most of them are effective for 1D and 2D problems, but their generalization for the 3D case is not obvious and often leads to the method's instability. In [29], we compared the efficiency of the three-stage iterative process and the splitting method and demonstrated the drawbacks of the split-step method for the problem under consideration. The main drawback is a violation of the FDS's conservativeness. Another one is an absence of asymptotic stability (long-time stability) due to the accumulation of computational error. This manifests in the numerical solution's symmetry violation contrary to the stated conditions of the process. To overcome this issue, it is necessary to decrease the mesh steps. As the non-stationary process is under consideration, the machine time could become enormous, especially for the 3D problem solution. Contrary to the split-step method, TSIP eliminates these disadvantages, as it provides conservativeness and asymptotic (long-term) stability [29]. In turn, both methods are economic ones.
Below, at the writing iteration stages, we use notations: where s denotes an iteration number. Similar notations are used for the operators yy zz Z n Z n ϕ ϕ as well as for the free-electron concentration and the pulse intensity at the corresponding iteration stages. The stages of the TSIP are the following: The first stage: The second stage: The third stage: We use values of the mesh functions at the previous time layer as the initial approach for those at the next time layer: Non-linear difference Equations (28), (35), and (40), with respect to the free-electron concentration with the corresponding BCs, are solved using the Thomas algorithm. The key issue is the choice of a numerical method for solving the Poisson equation concerning the electric field potential, and it is discussed below in Section 3.3. In turn, the numerical approach for computing the mesh function's initial distribution is presented in Section 3.4.
The criterion for the iteration process convergence is based on the relative error estimation: The computation is continued until all inequalities (46) are satisfied in all mesh nodes.

Numerical Methods for the Neumann Problem Solution
According to the research's purpose, we solve the Poisson equation at each of the stages by using the DM (FFT) and the IM (ADM). We use the ADM below because the Poisson equation is linear, and, in this case, the three-stage IM may not be efficient in comparison with ADM if the mesh number is not big enough. Then, we compare the computer simulation results obtained using both methods. As mentioned above, the Neumann problem solvability condition coincides with the charge conservation law for the problem under consideration. If we use DM for the Poisson equation solution, then the difference analogs of the ratio (19) should also be valid at the problem's numerical solution. Evidently, when providing numerical experiments, the conservation law is preserved in the best case within the accuracy of the problem approximation and for the problem under consideration, this is enough for the growth of the invariant's value because of non-linear feedback between the equations concerning the electron concentration and the electric field potential due to the electron mobility. Therefore, the Neumann problem solvability condition could not be fulfilled. As a consequence, the numerical solution will not approximate a solution to the differential problem.
Consequently, the numerical solution will not approximate a solution of the differential problem.
If we solve the Poisson equation on each of the stages by applying an IM, then it transforms into the difference "heat-conduction equation" by introducing the artificial time step τ as an iteration parameter: At each of the time layers, one can treat this equation as the difference analog of the differential Helmholtz equation The Equation (48) with the Neumann BCs is well known to be uniquely solvable if the following inequality 0 η > is valid [5]. Because the parameter τ is always positive, the Equation (47)  To apply the ADM [21], we introduce the auxiliary mesh functions F and Φ , defined on the spatial mesh x y z ω ω ω × × . Function F corresponds to the electric field potential for the iteration process stages (27) each of the iteration stages. As a result, the iteration process for the Poisson equations on each of the stages is written in the following manner: 1, 1, 1, 1, As it is well-known, such kinds of FDSs are unconditionally stable [21]. Each of the Equations (49)-(51) is solved by using the Thomas algorithm. The process is stopped if the estimation for the residual of the Poisson equation is satisfied: for all inner mesh nodes. As a result, the field potential at the corresponding iteration of the multi-stage process is equal to the function

Remark 2. Alternatively, it is possible to solve the Poisson equation on each of the stages by applying the TSIP:
( ) with the convergence criterion: Section 4 presents some relevant computer simulation results obtained using TSIP (54)-(57). These computations require more computational time than modified ADM (49)-(51); they have the same number of iterations but different operation numbers.

Numerical Method for Finding the Initial Distributions of the Mesh Functions
Accounting for the solvability condition of the Neumann problem for the difference Poisson equation, we replace the set of Equations (12)- (14) with equivalent ones by introducing a pseudo-time in each of the equations: where t is a positive parameter, which characterizes the artificial time, and constant r is a regularization parameter. The t -dependent BCs are stated for the electric field potential: is introduced to provide smooth reaching a steady-state of the electric field strength, and ϕ τ is a positive parameter characterizing a rate of exponent tending to zero. The BCs for the free-electron concentration are stated as in (6). The initial conditions for the Equations (58)-(60) are: So, to obtain initial distributions of the semiconductor characteristics at the arbitrary BCs, we approximate the system (58)-(61), (6) by FDS and realize it by using the iteration process: The iteration process convergence criterion is given by the formulas:   x in x y z in eq N x y z n e μ ϕ − =

Computer Simulation Results
Below, we compare the computer simulation results obtained using DM (FFT) and IM for solving the Neumann problem for the Poisson equation involved in the set (1)-(4) of the equations and accounting for the action of EEF on a semiconductor at statement of the BCs. We evaluate the method's efficiency in dependence on the EEF strength and directions of the electric field action. We also follow the invariant (26) to assess the method's accuracy.
For the DM implementation, the FFT code belonging to the Intel Math Kernel Library (Poisson Library Routines) is applied to solve the Poisson equation. The solver uses the forward and backward FFT. The description of the library application stresses that "For Poisson equation with Neumann BCs the right-hand side of the problem cannot be arbitrary. Poisson Library verifies that the classical solution exists (up to rounding errors). In any case, Poisson Library computes the normal solution, that is, the solution that has the minimal Euclidean norm of residual" [41] ADM (49)-(51) as well as TSIPs (27)- (45) are realized as the computer codes developed using C++ language.
The algorithm for solving the problem (20)-(25) is the following. Firstly, we compute the functions at the initial time moment based on the Formulas (63)-(65). Further, we use them to start the iteration process (27)- (44). On each of the iterations, we compute the 3D Poisson equation using one of the discussed methods (ADM or FFT). We provide computations for the following cases: EEF acts only along the longitudinal coordinate or along one of the transverse coordinates (x-axis, for example), or acts simultaneously along longitudinal and transverse ones. The problem (1)-(4) parameters have the following values: The mesh steps and the iteration process parameters are taken as follows:  In Figure 1, we show 2D projections of the electric field potential and free-electron concentration distributions in a moment of time t = 100. It is easy to see the strong influence of the EEF on them. If the EEF is absent and the incident beam profile is defined in the form (7), then the computed functions are symmetrical with respect to the center of the semiconductor's illuminated area. Preserving the solution symmetry is an important feature that allows us to estimate the numerical method's correctness. However, if the EEF acts on a semiconductor along the transverse coordinate, then the function distributions become asymmetric, complicating the evaluation of the obtained numerical solution. The distribution of the free-electron concentration represents a complicated spatial-temporal oscillating structure (Figure 1, right column). It should be stressed that these distributions depend significantly on the choice of numerical method. 10,
As mentioned above, conservativeness is one of the most important properties of the FDS on par with numerical method stability. Let us consider how a Poisson equation's solution method affects invariant preservation accuracy. As follows from Figure 2, at using IM, the invariant is preserved with very high accuracy: ( ) Q t order is about of magnitude if EEF is absent or it acts in one direction (along the optical pulse propagation or in the opposite direction). The invariant value increases by three orders up to a magnitude of 8 10 − if the EEF acts in two directions. Thus, the 2D EEF influences the preservation of the invariant more strongly than an action of the 1D EEF. Such ( ) Q t preservation accuracy is enough for providing the conservativeness property of the FDS if accounting for the order of the FDS approximation. Moreover, as follows from Table 1, the invariant value practically does not grow over time after reaching a certain value.  In contrast, at using DM for the Neumann problem solution, the conservativeness of the FDS is violated: the invariant's value grows over time even at the zero-value EEF (see Figure 3 and Table 1). There is also a dependence of the invariant's growth on the direction of EEF action. For example, if . In this case, as follows from Figure 3, the invariant's growth occurs approximately linearly, and evidently, it continues in time. The reason for such invariant evolution at using DM for the Poisson equation solution is the Neumann problem solvability condition violation. Due to the nonlinear feedback between the equations concerning the charged particle concentrations and the electric field potential, the round-off errors accumulate and, consequently, induce errors in computing the right part of the Poisson equation.
Computer experiments show that the EEF presence (a case of the inhomogeneous Neumann BCs) essentially influences the gradient and magnitude of the invariant's growth. This could be avoided by using a mesh with decreased steps. This can be seen in Figure 3b by comparing the black solid and the red dotted lines. However, decreasing the mesh steps leads to an unacceptable computational cost when solving multidimensional time-dependent problems.   (Table 1). At providing computation using IM, the invariant's ( ) Q t value increases up to several orders of magnitude (nevertheless maintaining a high accuracy) in the very short time interval at the beginning of computation in comparison with its value at the initial time moment. Then, it practically does not change. The situation becomes the opposite at using DM: regardless of the invariant's computation accuracy at the initial time moment, its value increases in time.
We stress the absence of the problem solution's stability if the Poisson equation is solved by using the DM ( Table 2). Increasing the computation accuracy for computation of the initial functions' distribution by increasing the severity of the criterion (66) (  Table 2, white lines). In contrast, when applying IM, the increase in the accuracy of the initial functions' distribution computation leads to the corresponding increase in accuracy of the preservation Q(t) ( Table 2, gray lines).  We see in Table 3 that using DM for the Poisson equation solution with the EEF acting in one direction the decreasing of mesh step z h by one order of magnitude improves the invariant conservation accuracy by two orders of magnitude. However, it is still six orders of magnitude greater than the invariant preservation accuracy reached when using IM. It should be noted that with the same mesh steps, DM computes faster than the IM: computation of up to 100 dimensionless time units takes about 63 h. In contrast, the iteration process requires 190 h. However, decreasing the spatial step in one of the coordinates by an order increases the computation time extremally at using DM: up to 570 h.
 Similar norms are used for an assessment of a difference between the free-electron concentrations computed using IM and DM for the numerical solution of the Poisson equation.
As one can see from the graphs depicted in Figure 6, the EEF plays a fundamental role in the DM's efficiency. If the EEF action is absent, then the C-norm for the difference between IM ϕ and DM ϕ has the order of magnitude 2 10 − , and the 2 L -norm is about 3 10 − (Figure 6a,b). The norm values dramatically grow (up to two orders) if the EEF is non-zero: C-norm value exceeds the unit for some cases. Analyzing Figure 6c,d, we can state that the difference between numerical solutions is maximal in the case of 2D EEF or in the case of EEF acting along longitudinal coordinate. Thus, the Poisson equation solution obtained using DM differs strongly in these cases. In Figure 7, the free-electron concentration distribution on spatial coordinates is depicted. We see a pronounced difference in both spatial distributions at using IM or DM. They differ both in maximal values of the free-electron concentration and in spatial coordinates of their appearance. We also depict the norms of the difference between functions which is defined by summing both the numbers of iteration processes that are necessary to fulfill the conditions (53), (57) and the number of iterations performed to satisfy the criterion (46). As one can see in Figure 9, at the beginning of the computation, the ADM requires significantly fewer iterations than the TSIP. With time increasing, it S becomes almost the same for both methods. It should be emphasized that the iteration process convergence is strongly determined by the regimes occurring in a semiconductor. Because the physical process under consideration is accompanied by forming complicated contrast dynamical structures, a number of the iterations, which are required for achieving sufficient solution accuracy, can vary significantly according to the problem parameters. From Figure 10, it follows that both methods give numerical solutions that coincide with high accuracy. To assess the difference between the solutions, we introduce additional notations:

Conclusions
We demonstrate the inefficiency of DM application for solving the Poisson problem with inhomogeneous Neumann BCs if it is involved to the set of non-stationary non-linear PDEs with realized non-linear feedback. This is a consequence of the violation of the solvability condition of the Neumann problem. We stress that the conservation law of the problem under consideration coincides with the corresponding solvability condition. Therefore, a violation of this law for the FDS leads to a violation of the corresponding solvability condition. Using the IM for the Poisson equation solution allows us to avoid a violation of the Neumann problem solvability condition.
The difference in the results obtained using DM and IM becomes more pronounced at the action of the EEF on a semiconductor. The direction of the EEF action significantly influences the difference between numerical solutions obtained by using IM or DM. When solving the Neumann problem using IM, the problem's invariant is preserved with high accuracy: it does not grow in time. By decreasing the spatial mesh steps at using DM for a solution of the Poisson equation, it is possible to reach a coincidence of the spatial distribution of the electric field strength obtained using IM with essentially big values of these mesh steps. However, this coincidence occurs during a limited time interval.
We demonstrated that the TSIP (54)-(56) is a promising alternative to the implicit ADM (49)-(51) for a numerical solution of 3D PDEs contained the terms describing reaction-diffusion-convection (in the case under consideration, the generation and recombination of the free charged particles, and electron diffusion and mobility). It possesses high accuracy and asymptotic stability. TSIP is easy to realize and can be used for parallelization algorithms. We highlighted that it is an economical one, as well as the ADM.
We stress that the approach presented in the study could be directly applied to many other scientific problems described by similar PDE sets, for example, in modeling cell chemotaxis or organic molecular systems [42][43][44][45][46].