Solving a System of Differential Equations Containing a Diffusion Equation with Nonlinear Terms on the Example of Laser Heating in Silicon

: We present a ﬁnite-difference integration algorithm for solution of a system of differential equations containing a diffusion equation with nonlinear terms. The approach is based on Crank–Nicolson method with predictor–corrector algorithm and provides high stability and precision. Using a speciﬁc example of short-pulse laser interaction with semiconductors, we give a detailed description of the method and apply it to the solution of the corresponding system of differential equations, one of which is a nonlinear diffusion equation. The calculated dynamics of the energy density and the number density of photoexcited free carriers upon the absorption of laser energy are presented for the irradiated thin silicon ﬁlm. The energy conservation within 0.2 % has been achieved for the time step 10 8 times larger than that in case of the explicit scheme, for the chosen numerical setup. The implemented Fortran source code is available in the Supplementary Materials. We also present a few examples of successful application of the method demonstrating its beneﬁts for the theoretical studies of laser–matter interaction problems. Finally, possible extension to 2 and 3 dimensions is discussed.


Introduction
Many phenomena occurring in nature for their investigation can be described via mathematical models based on time-dependent nonlinear diffusion equations [1]. Examples include genetics [2,3], image processing [4], quantum mechanics [5], and laser-material interactions [6]. Although during the last decades a big effort has been undertaken to find efficient numerical schemes for solution of the corresponding mathematical problem, some of the applications are still a challenging task. Specifically, the efforts in model implementation as well as their demands on the computational power during processing can substantially hinder the theoretical interpretation of the investigated problem. In this work, we consider an application of the nonlinear parabolic diffusion equation to describe the response of solids to an ultrashort laser pulse irradiation. Apart from the insights into the material structure, this topic is important for the description of laser machining [7][8][9] and nanostructuring experiments [10][11][12] with applications in biotechnology [13] and IT [14,15]. For metals, the problem may be mathematically formulated in the form of frequently used Two-Temperature Model (TTM) [16], whereas for semiconductors, a similar TTM-like approach has been proposed [17]. The latter is based on the system of partial differential equations, reflecting the conservation laws in the atomic subsystem of a solid and its electronic subsystem. Although it is relatively simple to apply an explicit finite-difference numerical scheme to solve such systems in metals [18] or semiconductors [19], the corresponding stability criteria demand the integration time steps to be small, causing high computational costs as a result. The main restriction on the time step often comes from the nonlinear diffusion equation describing the carrier heat conduction process [20]. One of the possibilities to increase the time step of diffusion equation is to use implicit or semi-implicit integration schemes. For instance, the Crank-Nicolson semi-implicit scheme [21,22] provides unconditionally stable solution when applied to linear diffusion equations. However, this approach is not directly applicable when nonlinear terms play an important role.
In this work, we present a semi-implicit finite-difference method for the solution of a system of differential equations-one of which is a diffusion equation with nonlinear terms-and apply it to model short laser pulse interaction with silicon. The presented approach is based on Crank-Nicolson method with predictor-corrector algorithm and provides high stability and precision. It has been already successfully applied for the investigation of ultrashort laser interaction with metals [23] and semiconductors [24][25][26]. Section 2 describes the theoretical model and presents the system of equations where a nonlinear diffusion equation results in strong restriction on the time step for the explicit integration algorithm. In Section 3, we give a detailed description of semi-implicit numerical solution scheme, which was modified with predictor-corrector algorithm to account for the nonlinearity in the diffusion equation. Further, in Section 4, the calculation results for a particular set of parameters are presented and the energy conservation versus the applied iteration time step is investigated. Section 5 mentions the existing works, in which this approach has been successfully utilized, and suggests possible improvements for the application of the presented approach in two-dimensional (2D) and three-dimensional (3D) cases. Finally, in Section 6, we give a summary of our results.

Model Description
To demonstrate the application of the finite-difference scheme, we use an example of laser irradiation on silicon. Below, we present the full set of nonlinear differential equations for the continuum description of electron density and electron/phonon energy density dynamics in silicon under ultrashort laser irradiation. For the derivation of the following expressions we refer to the works in [17] and [25]. Due to laser pulse irradiation (in this example Ti:Sapphire laser at 800 nm wavelength), free carriers are generated in the material, electrons in the conduction band, and holes in the valence band. Both types of carriers are assumed to quickly equilibrate in the corresponding bands and move together due to the Dember field preventing charge separation [27]. To each of them, we apply separate Fermi-Dirac distributions with different chemical potentials, φ e and φ h , for the electrons and holes, respectively, but with shared carrier density n and temperature T e (two-chemical potentials model). The reduced chemical potentials are defined as follows, where E C and E V are the conduction and valence band energy levels, respectively, so the energy gap is E g = E C − E V . The integration of the carrier distribution functions over the energy leads to the expressions for the carrier density (parabolic bands are assumed): where subscript c stands as e for electrons and h for holes. The Fermi-Dirac integral is defined as The carrier current is the sum of contributions from the electrons and the holes: with q e the elementary charge. Ambipolar energy flow is the sum of diffusion and thermal energy currents inside the carrier subsystem and can be written as The dynamics of semiconductors under the irradiation of ultrashort laser pulses can be modeled with the system of three continuum equations [17,28]: continuity equation for free carrier density and two coupled energy balance equations, one for the carriers and one for atoms: The meanings of symbols in Equations (7) to (9) are the following, S n is the source of new carriers (excitation rate of new carriers by the laser), S u describes the energy source (rate of laser energy absorption), and T a is atomic temperature. The terms on the right hand side of Equation (7) are responsible for the carrier generation, Auger recombination, and impact ionization, consequently. Equation (8) describes the energy balance in the photoexcited electron-hole pairs and is a nonlinear diffusion equation. The last Equation (9) describes the energy balance in the atomic subsystem. C e−h is specific heat capacity of the electron-hole pairs: The total energy of electron-hole pairs consists of the energy gap and the kinetic energy of electrons and holes (taking into account the Fermi statistics), The parameters used in the calculations as well as the meanings of other symbols are presented in Table 1.
To present an example of the model application, we use the following source terms. The rate of free carrier density growth, S n , and the corresponding rate of their energy increase, S u , are given by S u = αI abs ( r, t) + βI 2 abs ( r, t) + ΘnI abs ( r, t) .
In the last two equations, the first and second terms on the right hand side represent the influence of one-and two-photon absorption, respectively, and the third term in the second equation represents the laser energy absorption by the excited free carriers.
One-dimensional (1D) laser heating problem is analyzed in this work. The laser is focused on the material surface. The corresponding form of laser intensity at the surface (z = 0) in this case is where Φ inc is the incident fluence, ς = 4 ln 2, and R(T a ) is the reflectivity function (see Table 1). In the present calculations, the center of Gaussian pulse is shifted in time from the initial time, t = 0, to 3 pulse duration times, 3 t p , which in turn is defined as the pulse width at the half of maximum. The dependence of the laser pulse intensity, I abs , on depth can be found upon the solution of differential equation of the attenuation process: where z is the depth into sample; the terms on the right side are responsible for one-and two-photon absorption, and for the free-carrier absorption processes. Thus, from the system of Equations (7) to (9), we can fully determine the dynamics of n, T e , and T a in 1D using the following initial and boundary conditions, suitable for a free standing film: n (z, 0) = n eq = 1 × 10 16 m −3 , ref. [29], where L is the thickness of the sample. Owing to its similarity with an ordinary well-known TTM model [16], but with an additional equation for free carriers density n, we refer to the described approach as nTTM model, as it was suggested in [24].

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 ∂T a ∂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 ≡ T e ): 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: The right-hand side contains parameter ψ, which can be 0 for explicit scheme, 1 for implicit, and 1 2 for semi-implicit. The function f k i can be found from: where W k i is defined according to Figure 1 (between the cells) and Equation (6): Analogously, according to Figure 1 and Equation (4), the carrier current is with Any function in between cells can be found by averaging The Fermi-Dirac integrals were calculated using GNU Scientific Library [40] and stored in the tables to speed up the calculations. ∂η c ∂n can be found by derivative of Equation (2) by the carrier density: ∂η c ∂T e can be found by taking the derivative of equation Equation (2) by the electronic temperature: The boundary conditions can be rewritten in the finite-difference form as follows, T k+1 and on the other edge n k+1 All the other equations and connections between variables at the boundaries stay the same and can be straightforwardly obtained by substituting i = 1 and i = N into Equations (17), (21), (25) and (26).
At the current time step, k, we do not have any information about the following parameters from the future time step k + 1, Initially we set them all except the first three to be equal to the corresponding old values (at time step k): Here, (0) means the 0 th step of the corrector. With this assumption, Equation (19) becomes which can be represented as a tridiagonal system of equations, where for i = 2, ..., N − 1, and the boundary conditions are 1 ∂T a ∂t k+1 1 (45) and (48) Such a system can be resolved with respect to T k+1 i N i=1 by using the well-known tridiagonal matrix algorithm [41]. We denote the electronic temperature calculated with assumption (36) as (T k+1 i ) (1) , showing with "(1)" the first correction step. This result allows to calculate the corrected new values of functions in list (35): For an improved precision, the corrected new values of n and T a can be calculated using the semi-implicit approach (instead of the explicit scheme, Equations (17) and (18), used for the predictor): In turn, the corrected values allow to calculate (T k+1 i ) (2) from Equation (38) and so on. Owing to its similarity with predictor-corrector methods, we call it a "predictor-corrector" algorithm. With this approach, Equation (19) can be rewritten in the following form, where index (l) shows the current step of correction and (l) = (0) means the value is old, i.e., taken at time step k. This procedure continues until the difference between last two corrected values of electronic temperature is less than the demanded precision: It takes up to 200 corrections to reach the chosen precision of ε = 10 −6 K during the laser pulse action, whereas when the laser is ended, 5 corrections are usually enough. The chosen numerical setup is discussed in the next section.

Calculation Example
As an example of application of our algorithm to the described system of equations Equations (7) to (9), we perform the simulations of 800 nm thick silicon target's response to ultrashort laser pulse irradiation. The parameters of the irradiation are 130 fs duration, 800 nm wavelength, and 0.26 J/cm 2 incident fluence. For these conditions, the experimental melting threshold fluence is 0.27 J/cm 2 [42], which is in agreement with the result of the nTTM model [43]. The value of fluence is chosen to be just below the melting threshold, providing the applicability of this simple model in the absence of phase transition processes. The sample was divided into 160 cells according to Figure 1. In Figure 2, we show the dynamics of electron-hole carrier density and electronic and atomic energy densities at the silicon surface. The shown energy density is scaled to the melting energy density, which is found to be 3.86 × 10 9 J/m 3 , according to the simulations. Though Equations (8) and (9) are written in terms of temperature of carriers and atoms, we plotted the corresponding energy densities instead, because their dynamics represents the energy flow between the subsystems and allows plotting the same scale for electrons and atoms, whereas electronic temperature is much higher than the atomic one (see also Figure 2 in [25]). In addition, this choice provides a possibility to show the energy conservation with the total average energy density of the sample (shown with black solid line).
The initial increase in the carrier density followed by the laser pulse is connected to the excitation of new carriers by one-and two-photon absorption processes. With time, the increase changes to the decay due to strong Auger recombination and diffusion processes. The strong peak in the electronic energy density is mostly connected to the free-carrier absorption. Finally, the thermal energy from electron-hole carriers is transferred to the atomic subsystem of the sample leading to gradual increase in the lattice energy density upon the electron-phonon equilibration.

Discussion
In case of the explicit scheme, a good guess for the time step requirement can be obtained from the von Neumann stability criterion [20], ∆t ≤ (∆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 × 10 −19 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 × 10 −24 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 × 10 −16 s with the energy conservation about 0.16 % per simulation, Figure 3. At the same time, in case the calculation speed is critical, increasing the time step even higher is possible: 1 × 10 −15 s 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.     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.
This approach has been successfully applied earlier in order to investigate and improve the presented nTTM model [24]. In [25], we used the described scheme for the solution of the continuum part of the atomistic-continuum model MD-nTTM. The atomistic-continuum model describing the dynamics of gold targets under the ultrashort-pulse lasers also benefited from using the presented approach [23]. The high speed and precision of the scheme allowed to significantly decrease the computational costs of the corresponding simulations. Recently, Rathore et al. utilized the presented finite-difference scheme to simulate temporal evolution of photoinduced thermal strain in InSb [26].
In the mentioned applications, the corresponding system was solved in 1D, based on the assumption of wide laser spot in comparison with the lateral sizes of the computational setup [25]. Whenever it is not the case, one needs to solve the corresponding problem (the vector system of Equations (7) to (9)) in 2D or 3D case. According to the work in [44], the diffusion equation in 3D case can be solved in 3 subsequent steps, each of which involves implicit solution in only one direction (X, Y, or Z) and explicit scheme in the other two directions. In other words, one can use 1D implicit scheme three times: for X, Y, and Z directions separately and consequently. This approach is called alternating direction implicit method (ADI) and is also widely applied for the corresponding 2D problems [45]. Therefore, with appropriate modifications, the presented scheme should be applicable for the considered nonlinear problem in 2D and 3D cases as well.

Conclusions
We proposed the semi-implicit integration scheme for the solution of nonlinear diffusion equations or systems of equations containing nonlinear diffusion equations. The scheme is based on the Crank-Nicolson finite-difference integration method, modified with a predictor-corrector algorithm, according to Equation (52). The modification resulted in a possibility to solve nonlinear diffusion equations with high stability and precision. The implemented Fortran source code is available in the Supplementary Materials under the terms of the GNU General Public License version 3 or later.
In the presented example of the scheme application, we reached the speed up of the calculations (by the increase of the integration time step) by up to 10 8 times compared with the explicit scheme, keeping the error in energy conservation below 0.2 %. This error increases linearly with the time step. The algorithm is applicable in case the time step is much smaller than all the characteristic times of the involved physical processes. The existing applications that use the proposed scheme are mentioned and the possible extension for 2D and 3D cases is suggested.