1,706
Views
4
CrossRef citations to date
0
Altmetric
Original Articles

Conservative finite-difference scheme for the 2D problem of femtosecond laser pulse interaction with kink structure of high absorption in semiconductor

, &
Pages 207-244 | Received 19 Sep 2017, Accepted 01 Jun 2018, Published online: 07 Jul 2018

ABSTRACT

The problem of high-intense laser pulse interaction with a semiconductor under the condition of a light energy nonlinear absorption, which results in high absorption domains formation, is considered. Such interaction allows reaching a construction of an element for all-optical data treatment. For its adequate description we propose new mathematical model taking into account the longitudinal and transverse diffraction effects. The longitudinal diffraction induced a reflection of laser radiation from boundaries of the high absorption domains that results in changing of their spatial structure. The conservative finite-difference scheme (FDS) is developed for numerical computation of the complicated nonlinear processes. The property of the FDS conservatism is proved. For the proposed FDS realization the two-stage iteration method is proposed. Computer simulation results are presented. We show the uniform boundedness of the mesh functions at the iterations and the convergence of the iteration process. We show also positiveness and boundedness of the mesh function corresponding to free-charged particles. We discuss some properties of differential problem including the problem invariants, positiveness of the free-charged concentrations.

2010 MSC CLASSIFICATIONS:

1. Introduction

Laser radiation interaction with a semiconductor is a very modern problem for the past decade. One of the important application of such investigation is a creation of an optical bistable element using various nonlinear responses of a semiconductor at a laser radiation exposure. As it is well-known, the optical bistability (OB) is a very promising phenomenon for creation and developing of all-optical data processing and many authors pay their attention to OB phenomenon research since pioneer works [Citation10,Citation11] and continued, for example, in [Citation8,Citation9,Citation13,Citation15,Citation19,Citation23,Citation28,Citation31,Citation32,Citation33,Citation38,Citation39]. In the case of OB existence, the hysteresis loop of semiconductor characteristics in dependence on a laser pulse intensity takes place. As a result, the domains of charged particles high concentration as well as kink structures of a high absorption could be formed, which results in the appearance of contrast spatial structure of free-charged particle concentration as well as explosive absorption of laser energy. In the present paper, we consider the 2D resonatorless model of the absorption OB.

The femtosecond (1015c) laser pulse interaction with a semiconductor under the condition of absorption OB occurrence is described by the set of nonlinear PDEs [Citation9]. For numerical solution of such problems, the split-step method is widely used [Citation2,Citation4,Citation6,Citation17,Citation29,Citation30,Citation36,Citation37,Citation40]. For example, in [Citation40] this method is used for a numerical solution of the nonlinear Schrödinger equation. It is shown, how a conservation law and instability are reflected in the numerical finite-difference schemes (FDS). In [Citation37] an efficient time-splitting compact finite-difference method for Gross–Pitaevskii equation is proposed. In [Citation29] a simulation of the nonlinear Schrödinger equation in 1D, 2D and 3D cases is studied on the base of a time-splitting method. The FDS conservatism is investigated on the base of some numerical experiments. In [Citation4] a compact finite-difference method is proposed for the nonlinear Schrödinger equations with constant and variable coefficients.

Much attention is paid to a question of the split-step method properties, such as the conservatism, accuracy and stability. In [Citation2] it is shown, that using of energy-conserving methods for the nonlinear Schrödinger equation, able to conserve a discrete counterpart of the Hamiltonian functional, confers more robustness on the numerical solution of a problem. In [Citation17] the problem of numerical instability of the split-step method applied to the generalized nonlinear Schrödinger equation is widely studied. It should be stressed, that the Fourier method or its modifications applied to such problems [Citation6,Citation40], but this approach requires the certain boundary conditions (BCs) for the problem under consideration. Our aim is to develop the finite-difference method for the problems with arbitrary BCs.

Earlier, at theoretical analysing of the resonatorless OB, in particular in our papers [Citation32,Citation33,Citation34,Citation35] also, the charged particles concentration, a potential of the electric field induced by a laser radiation and a laser pulse intensity changing due to a nonlinear absorption at its propagation in a semiconductor have been considered as a rule. However, it is well-known that the free-charged particle high concentration domains formation is observed under the condition of the hysteresis loop realization. Therefore, a laser beam interaction with these domains may result in partial light reflection from their boundaries. To describe this process, we propose a new mathematical model of the problem. Its main feature consists in a replacement of the equation concerning the laser pulse intensity by the Schrödinger equation with respect to the envelope (complex amplitude, slowly varying in time only) of a wave packet. In our opinion, the diffraction effects could play a significant role for changing the OB element characteristics, such as a stability of the bistable states and time of switching from one state to another one, for example. In this regard, taking into account the diffraction effects at the analysis of switching waves is the urgent problem.

For computer simulation of the problem, we develop a conservative FDS similar to that, which we proposed in [Citation34] for the 2D problem without considering diffraction effects. It is a nonlinear implicit one, so to realize it we proposed an original two-stage iteration process. Computer simulation results confirmed that the FDS also possesses an asymptotic stability property that is very important because we should provide computation during the long time interval.

2. Problem statement

In the present paper, we consider the 2D problem of femtosecond laser pulse interaction with a semiconductor taking into account a longitudinal and transverse diffraction of a laser beam. The semiconductor has a rectangular form and could be placed in the external electric field (). Characteristic time of changing the problem parameters is about several tens femtoseconds. A mathematical model for the problem consists from a set of four 2D dimensionless PDEs concerning semiconductor characteristics – free-electron concentration n(x,y,t), ionized donors concentration N(x,y,t), electric field potential ϕ(x,y,t) [Citation9,Citation26] and a laser pulse complex amplitude A(x,y,t), which is slowly varying in time only: (1) Nt=G(n,N,|A|2)R(n,N),0<x<Lx,0<y<Ly,t>0,(1) (2) nt=κxxnxμxnϕx+κyynyμynϕy+G(n,N,|A|2)R(n,N),(2) (3) 2ϕx2+2ϕy2=γ(nN),(3) (4) At+iDAx2Ax2+iDAy2Ay2+βAδ02δ(n,N)A=0,0<x<Lx,Ly<y<Ly,t>0.(4)

Figure 1. Setup for computer simulation of the femtosecond pulse interaction with a semiconductor. Ex, Ey – external electric field strength (x, y are dimensionless spatial coordinates).

Figure 1. Setup for computer simulation of the femtosecond pulse interaction with a semiconductor. Ex, Ey – external electric field strength (x, y are dimensionless spatial coordinates).

The following notations are used: variables x, y are the dimensionless spatial coordinates, parameters Lx, Ly denote their maximal values, correspondingly. Variable t denotes a dimensionless time. The electron diffusion coefficients κx, κy, the diffraction coefficients DAx,DAy and the electron mobility coefficients μx, μy are non-negative parameters. Parameter γ depends, in particular, on the maximal concentration of free-charged particles. δ0 denotes a maximal value of the semiconductor laser energy absorption. χ is a wave number, βA=πχ, DAy=1/4πχ. According to physical sense of the problem, the concentrations n(x,y,t), N(x,y,t) and the absorption coefficient δ(n,N) must be non-negative. Therefore, the concentration N(x,y,t) must be less or equal to unity, as it is shown in Section 2.1.

The homogenous BCs concerning the free-elector flow correspond to the electric current absence through a semiconductor surface: (5) nxμxnϕxx=0,Lx=nyμynϕyy=0,Ly=0.(5) The non-zero electric field strength Ex, Ey corresponds to the external electric field presence: (6) ϕxx=0,Lx=Ex,ϕyy=0,Ly=Ey.(6) We stress that below we consider a case of the external electric field is absent: (7) Ex=Ey=0.(7) The zero-value BCs with respect to a complex amplitude correspond to its finite distribution: (8) A|x=0,Lx=A|y=Ly,Ly=0,(8) and we suppose, that the function A and its nth derivatives on x-coordinate and y-coordinate decrease exponentially when the corresponding coordinate tends to the domain boundaries. Some remarks concerning Equation (4) are discussed in Section 2.4.

Gaussian beam falls on a semiconductor. We stress that the spatial distribution of a laser pulse complex amplitude along a propagation coordinate coincides with its time distribution because a laser pulse propagates in a linear medium before its interaction with a semiconductor: (9) ϕ|t=0=0, n|t=0=N|t=0=n0,A|t=0=A0=e(xLxc)2/a2(yLyc)102iπχ(yLyc),(9) where Lxc,Lyc are coordinates of the incident laser beam centre position, a is a beam radius on x-coordinate. As the initial laser beam profile, which is shown in Figure , is a symmetric one, and then the problem solution should also keep symmetry if an external electric field is absent. It is an important feature for numerical method efficiency assessment.

Figure 2. Square root from the intensity distribution of incident pulse (a) and a distribution of square root from the intensity along y-coordinate at beam centre on x-coordinate (x=Lxc) (b).

Figure 2. Square root from the intensity distribution of incident pulse (a) and a distribution of square root from the intensity along y-coordinate at beam centre on x-coordinate (x=Lxc) (b).

Functions G and R describe charged particle’s generation and recombination [Citation9,Citation26] in the area of semiconductor, correspondingly: (10) G(n,N,|A|2)=q0|A|2δ(n,N),(10) (11) R(n,N)=nNn02τR,0xLx, 0yLy, t>0(11) otherwise, these functions are equal to zero: (12) G(n,N,|A|2)=R(n,N)=0,0xLx, Lyy<0, t>0.(12) We will keep in mind this definition of the functions below. Above q0 is a dimensionless maximum value of the laser pulse intensity, parameter n0 is the equilibrium value of free-electron concentration and ionized donor concentration before the laser pulse action, τR characterizes a recombination time (in present work we provide computations for τR=1). The absorption coefficient δ(n,N) could be approximated by different ways. Below we consider its following approximation [Citation9]: (13) δ(n,N)=(1N)eψ(1ξn),(13) where ξ, ψ are non-negative constants. Boundedness of the absorption coefficient is discussed in Section 2.2. We also write in Section 2.7 some assessments for the electric field potential, which is a solution of Equation (3) and there are some remarks concerning Equation (2) solution existence in Section 2.8.

2.1. Analysis of changing for the function N

Let us prove that the concentration of the ionized donors is less than unity: N1. First, we stress that the derivatives from the generation function on the concentrations n,N are: (14) Gn=ξψG0,if G0,(14) (15) GN=eψ(1ξn)|A|20.(15) Consequently, the function G increases with increasing of the concentration n and decreases with increasing of the concentration N.

The corresponding derivatives from the recombination function on the concentrations functions are: (16) Rn=N,RN=n(16) and this function increases with the concentrations of growth.

At the initial time moment, the functions G and R are equal to zero because the laser intensity is equal to zero and the concentrations are equal to their equilibrium values. If a laser pulse penetrates in a semiconductor then the function G becomes positive and the function N (as well as n) becomes grow. Let us suppose that in certain time moment the function G-R possesses the negative value: the right part of Equation (1) will be negative. Consequently, there is a value N(and it is less than unity), at which the right part of Equation (1) is equal to zero. It means that further increase of the function N stops. Thus, the function N cannot be greater than unity.

2.2. Boundedness of the absorption coefficient

Taking into account that the value of the concentration N does not exceed unity, the concentration n is greater than zero, and then the following assessment for the function δ(n,N) can be written: (17) δ(n,N)δmax=eψ+ξψnmax,nmax=maxx,y,tn(x,y,t)(17) for enough small time interval t[0,tε]. More weak assessment can be obtained taking into account the equation [Citation18,Citation22]: (18) t(n+(1N))=DxJxx+DyJyy,(18) where Jx=nxμxnϕx, Jy=nyμynϕy are the electric current flows.

We see that changing of the concentration sum is due to free-electron flow. Obviously, in small time interval t[0,tε], depending on diffusion coefficients, this flow is close to zero approximately and the derivative on time coordinate in (18) also is equal to zero. In this case, we can believe that the free-electron concentration n is equal to the concentration of ionized donors N (this statement is confirmed by numerical simulation because the diffusion velocity is bounded by the diffusion coefficients and spatial distribution of the free-electron concentration). Therefore, at first stage of laser pulse interaction with a semiconductor (during small time interval) and for qualitative analysis of the absorption coefficient from the free-electron concentration n, we analyse the following function: (19) δ(n,N)=δ(n)=(1n)eψ(1ξn).(19) To find the maximum of this function we take a first derivative, which is (20) δn=(ψξ(1n)1)eψ(1ξn).(20) At a value of the concentration n=11/(ψξ), the first derivative (20) is equal to zero. The second derivative from the function δ(n) is negative at this value of the free-electron concentration (21) 2δ(n)n2=ψξeψ+ψξ1<0.(21) Thus, in a dependence of value for the parameter ψξ there are three various dependences of the absorption coefficient from the free-electron concentration, which depicted in . We see in Figure  if the parameter ψξ is greater than unity that there is the same value of the function δ(n) for two values of free-electron concentration. It is easy to see that a maximal value of the absorption coefficient is defined by the following formula: (22) δmax=1ψξeψ+ψξ,ψξ>1,eψ+ψξ,ψξ1.(22)

Figure 3. Qualitative dependence of the absorption coefficient from the free-electron concentration at various values of the parameters multiplication ψξ.

Figure 3. Qualitative dependence of the absorption coefficient from the free-electron concentration at various values of the parameters multiplication ψξ.

2.3. Invariants of the problem and some estimation for the problem solution (Equations (1) and (2))

As it is known, the charge conservation law takes place for this problem (23) Q(t)=0Ly0Lx(n(x,y,t)N(x,y,t))dxdy=const,(23) if the charge flow through the semiconductor boundaries is absent. It is easy to see from Equations (1) and (2). To obtain formula (23) it is necessary to subtract Equation (1) from Equation (2) and then to take an integral over the semiconductor domain. Therefore, the FDS conservatism consists of a validity of this invariant difference analogue.

In the same way, one can write one more integral from the concentrations of free-electron and ionized donors: (24) Q¯(t)=0Ly0Lx(n(x,y,t)+(1N(x,y,t)))dxdy=0Ly0Lx(n0+(1n0))dxdy=LxLy.(24) If we define a norm L1 for non-negative function as ||f||L1=0Ly0Lxf(x,y)dxdy, and suppose non-negativeness of the functions n(x,y,t),1N(x,y,t) (for N this assumption was considered above) then one can re-write quality (24) in the following way: (25) ||n||L1+||1N||L1=||n0||L1+||1n0||L1=LxLy.(25) Let notice that taking into account the initial conditions (9) with respect to charged particles concentrations one can re-write the integral (23) as: (26) ||n||L1=||N||L1.(26) Thus, if the norm for one of the functions n(x,y,t),N(x,y,t) is bounded then the norm for other function is bounded also.

2.4. Estimation for the Schrödinger equation solution (Equation (4))

For the set of Equations (1)–(4) we can write an integral relation with respect to the solution of the Schrödinger equation (4). For this purpose, we multiply Schrödinger equation concerning complex amplitude A by the conjugated amplitude A, and the conjugated Schrödinger equation – by the A and then summarize these two equations. After that we integrate obtained equation: (27) tLyLy0Lx|A|2dxdy+iDAxLyLy0LxA2Ax2A2Ax2dxdy+iDAyLyLy0LxA2Ay2A2Ay2dxdy+βAδ0LyLy0Lxδ(n,N)|A|2dxdy=0.(27) Integrating the second and the third terms of Equation (10) by parts and considering the zero-value BCs for the function A we obtain the equalities: (28) LyLy0LxA2Ax2A2Ax2dxdy=0,LyLy0LxA2Ay2A2Ay2dxdy=0.(28) Thus, we obtain the following integral equality (energy invariant): (29) LyLy0Lxt|A|2+βAδ0δ(n,N)|A|2dxdy=0.(29) We can write an assessment for the norm L2 from the function A using (29). Indeed, let us suppose the existence of a solution for the problem (1)–(4) and a validity of the following inequality for absorption function δ(n,N) (see also (17)): (30) δminδ(n,N)δmax,δmin,δmax>0.(30) In this case, from Equation (13), one can write two inequalities: (31) tLyLy0Lx|A|2dxdy+βAδ0δminLyLy0Lx|A|2dxdy0,(31) (32) tLyLy0Lx|A|2dxdy+βAδ0δmaxLyLy0|A|2dxdy0.(32) If we introduce the norm L2for the complex amplitude A as (33) ||A||L22=(A,A)=LyLy0LxAAdxdy=LyLy0Lx|A|2dxdy,(33) which characterizes the pulse energy, then we obtain the following assessment of its evolution in time (34) ||A||L22||A0||L22eβAδ0δmint,(34) (35) ||A||L22||A0||L22eβAδ0δmaxt.(35) Therefore, the laser pulse energy decreases with time, but its value is always greater than zero.

2.5. Integral relations for the pulse impulse

It is possible to write another two integral relations from Equation (4), which characterize the pulse energy motion in an x-coordinate or in y-coordinate. With this aim, we multiply the Schrödinger equation concerning laser pulse complex amplitude A by a derivative of the conjugated amplitude on x-coordinate (A/x), and the conjugated Schrödinger equation by the (A/x) and then subtract these two equations. After integrating the obtained equation we can get: (36) LyLy0LxAtAxAtAxdxdy+iDAxLyLy0Lx2Ax2Ax+2Ax2Axdxdy+iDAyLyLy0Lx2Ay2Ax+2Ay2Axdxdy+0.5βAδ0LyLy0Lxδ(n,N)AAxAAxdxdy=0.(36) Using the partial integration for the second term and the third term in (36) and taking into account the zero-value of the function and its derivatives if the corresponding coordinate tends to the domain boundaries we obtain the following relation: (37) tLyLy0LxImAAxdxdy+βAδ0LyLy0Lxδ(n,N)ImAAxdxdy=0,(37) or (38) LyLy0LxImAAxdxdy+βAδ00tLyLy0Lxδ(n,N)ImAAxdxdydη=LyLy0LxImA0A0xdxdy,(38) where a symbol Im denotes an imaginary part of complex function.

Let us represent the function A using Euler's formula: (39) A=A¯eiS,(39) here introduced functions A¯,S are real ones. Then, we obtain the following integral relation: (40) LyLy0LxA¯2Sxdxdy+βAδ00tLyLy0Lxδ(n,N)A¯2Sxdxdydη=LyLy0LxA¯02S0xdxdy.(40) Thus, if initial distribution of the pulse phase S0 does not contain the odd powers from xLxc and initial function distribution A¯0 is symmetrical one on x-coordinate with respect to the point Lxc, then the integral in the right part of equality (40) is equal to zero and the first integral in (40) also will be equal to zero: (41) LyLy0LxA¯2Sxdxdy=0.(41) because the functions δ(n,N) and A¯ are symmetrical ones on x-coordinate with respect to this point Lxc in this case. This last statement is easy to prove by changing x on –x and taking into account the symmetry of the initial distribution of the complex amplitude in this coordinate. The integral (41) is used to control the computer simulation results.

In the same way one can obtain the following integral relation: (42) LyLy0LxImAAydxdy+βAδ00tLyLy0Lxδ(n,N)ImAAydxdydη=LyLy0LxImA0A0ydxdy(42) or in the representation (39) for the complex amplitude: (43) LyLy0LxA¯2Sydxdy+βAδ00tLyLy0Lxδ(n,N)A¯2Sydxdydη=LyLy0LxA¯02S0ydxdy.(43) For initial condition (9) this equality is re-written as: (44) LyLy0LxA¯2Sydxdy+βAδ00tLyLy0Lxδ(n,N)A¯2Sydxdydη=2πχLyLy0LxA¯02dxdy=const.(44) Let us remind that the first derivative from wavefront distribution on spatial coordinate characterizes the laser beam motion along this coordinate. Therefore, if we take into account the optical reflection from the semiconductor face and from boundaries of the high absorption domains, which leads to changing of the sign of the derivative for part of the pulse, then a modulus of the similar derivative for the other part of the laser beam must increase. This will lead to acceleration of laser pulse motion. This phenomenon is confirmed by the computer simulation results (see below).

2.6. Integral relation for the third invariant of Equation (4)

Another integral relation, which is similar to the Hamiltonian of Equation (4), can be obtained by multiplying the Schrödinger equation concerning laser pulse complex amplitude A by a derivative from the conjugated amplitude on t coordinate (A/t), and the conjugated Schrödinger equation is multiplied by the (A/t) and then subtract these two equations. After integrating of the obtained equation we can get: (45) iDAxLyLy0Lx2Ax2At+2Ax2Atdxdy+iDAyLyLy0Lx2Ay2At+2Ay2Atdxdy+0.5βAδ0LyLy0Lxδ(n,N)AAtAAtdxdy=0.(45) Using the partial integration for two terms in (36) and taking into account the zero-value of the function and its derivatives if the corresponding coordinate tends to the domain boundaries we obtain the following relation: (46) itLyLy0LxDAxAx2+DAyAy2dxdy+0.5βAδ0×LyLy0Lxδ(n,N)AAtAAtdxdy=0(46) or (47) LyLy0LxDAxAx2+DAyAy2dxdy=I30+βAδ00tLyLy0Lxδ(n,N)ImAAηdxdydη,(47) where I30=LyLy0LxDAxA0x2+DAyA0y2dxdy.

Using the representation (39), we write the equality (47) in the form: (48) LyLy0LxDAxA¯x2+DAyA¯y2+A¯2DAxSx2+DAySy2dxdy=I30βAδ00tLyLy0Lxδ(n,N)A¯2Sηdxdydη.(48) To simplify this equality, we need to write the equations with respect to A¯,S from Equation (4). Finally, they are written in the form: (49) A¯St+DAx2A¯x2A¯Sx2+DAy2A¯y2A¯Sy2=0,(49) (50) A¯2t2DAxxA¯2Sx2DAyyA¯2Sy+βAδ0δ(n,N)A¯2=0.(50) Using (49), one can transform equality (48) to the form: (51) LyLy0LxDAxA¯x2+A¯2A¯x2+DAyA¯y2+A¯2A¯y2+A¯2Stdxdy=I30βAδ00tLyLy0Lxδ(n,N)A¯2Sηdxdydη.(51) Thus, using the partial integration in (51) and taking into account the zero-value of the function A¯ if the corresponding coordinate tends to the domain boundaries we obtain the following relation: (52) LyLy0LxA2Stdxdy=I30βAδ00tLyLy0Lxδ(n,N)A¯2Sηdxdydη(52) or (53) tLyLy0LxA2Stdxdy+βAδ0LyLy0Lxδ(n,N)A¯2Stdxdy=0.(53) This integral relation can be used also to control computer simulation results. If we introduce the function (54) F=0tδ(n,N)A¯2Sηdη,(54) then one can rewrite Equation (52) in the form (55) LyLy0Lx1δ(n,N)t(FeβAδ0t)dxdy=I30eβAδ0t.(55) Because the right part of Equation (55) is positive then the integral in the left part of (55) must be positive also. The function δ(n,N) is positive as well. So, using the inequality (30) for the absorption coefficient and after some algebra, we obtain the following assessment: (56) LyLy0Lxδ(n,N)A¯2StdxdyI30δmaxeβAδ0t.(56) Moreover, from (56) it is as follows: (57) LyLy0LxA¯2StdxdyδmaxδminI30eβAδ0t.(57) Taking into account the inequality (34)–(35) and positiveness of the function A¯2 one can claim that the change of the pulse phase S will be bounded as in the norm L2 as well as in the norm C.

2.7. Equation (3)

Let us write some assessments for the potential of the electric field, which is a solution of the Poisson equation. As this equation is a linear one with respect to the electric field potential and the domain is a rectangle, then an existence of its smooth solution is investigated (see, e.g. [Citation7,Citation16,Citation20,Citation24,Citation27]). We suppose that the functions in the right part of the equation are known.

Let us multiply Equation (3) on the function ϕ and then take an integral over the domain occupied by a semiconductor. As a result, we obtain: (58) 0Ly0Lx(ϕΔϕ)dxdy=γ0Ly0Lx((Nn)ϕ)dxdy(58) or, taking into account the BCs with respect to the electric field potential: (59) 0Ly0Lx(ϕ)2dxdy=γ0Ly0Lx((Nn)ϕ)dxdy,ϕ=ϕx2+ϕy2.(59) Consequently, using the Cauchy–Bunkyakovsky inequality and ε inequality we write an assessment: (60) ||ϕ||L22γ|((Nn),ϕ)|γε||nN||L22+14ε||ϕ||L22,ε>0,(60) where the L2 norm from the electric field potential, or charged particles concentrations is defined in the domain, occupied by a semiconductor, as (61) ||ϕ||L22=(ϕ,ϕ)=0Ly0Lxϕ2dxdy.(61) Next, multiplying Equation (3) on the function (N–n) and then take an integral over the domain occupied by a semiconductor, we obtain: (62) ((Nn),Δϕ)=γ||Nn||L22(62) or (63) ((Nn),ϕ)=γ||Nn||L22.(63) Further, take modulus from the left part of Equation (63) and using the Cauchy–Bunkyakovsky inequality and ε inequality, we write an inequality: (64) γ||Nn||L22ε1||(Nn)||L22+14ε1||ϕ||L22,ε1>0.(64) From (60) and (64) one can write two assessments: (65) ||ϕ||L221ε4ε1ε1ε||(Nn)||L22+γ4ε||ϕ||L22,(65) (66) ||Nn||L22(1ε4ε1)ε1γ||(Nn)||L22+14εε1||ϕ||L22.(66) Consequently, a norm L2 of the electrical field potential gradient is bounded if the similar norm from the gradient of the concentrations difference is bounded. A norm L2 of the electrical field potential is bounded because a solution of the linear Poisson equation exists. These assessments can be used for proving the problem solution in small time interval as well as for proving the solution of the corresponding difference problem. Obviously, in the last case, we have to use the difference analogues of the inequalities written above.

2.8. Remarks about Equation (2) solution existence

The most difficult problem consists in proving the existence and uniqueness of the solution for Equation (2). Let us stress that, of course, it is necessary to consider the solution of the problem (1)–(4) with corresponding BCs and initial conditions. As our aim in this paper is not the proving of this solution we discuss only a set of possible steps for this proving by using already well-known approaches or theorems containing infamous books because PDEs involving in the problem under consideration is a subject of various investigations.

Analysing Equation (2) we can see that this equation has a structure which is similar to the hydrodynamics equations with taking into account the thermo-diffusion flows [Citation5]. The similar operators can meet for description of the gas mixture dynamics and they are similar to the Fokker–Planck equation [Citation12]. In particular, in [Citation5] the existence and uniqueness of the solution is demonstrated for the set of equations which possess the similar differential operator to Equation (2). Therefore, one can assume to use the results obtained in this paper for proving the solution of the set of equations (1)–(4).

There is another way for proving of the existence and uniqueness of the solution for the problem (1)–(4) during small time interval t[0,tε]. To reach this it is necessary to use multi-stages iteration process similar to the one used by us below for solution of the nonlinear difference equations. Let note that Equation (2) can transform to more convenient form: (67) nt=Dxxeμxϕx(neμxϕ)+Dyyeμyϕy(neμyϕ)+GR.(67) In Equation (67) we see the nonlinear dependence of the diffusion coefficients from the electric field potential. These coefficients are positive and do not equal to zero. In literature, the heat conduction equation with the nonlinear dependence of heat diffusion was investigated in [Citation7,Citation25]. But this equation was considered in the 1D case. That is why we propose to use two-stage iteration process.

It should be emphasized that in two particular cases it is easy to use the results of the book [Citation7]. First of them corresponds to the validity of the equality between the mobility coefficients: μx=μy. In this case Equation (67) reduces to the heat conduction equation with the nonlinear dependence of the diffusion coefficient and thermal capacity. The second one corresponds to the case, at which there is a relation between diffusion coefficients and electron mobility coefficients: Dxμx=Dyμy. Under this condition Equation (67) is similar to the Fokker–Planck equation.

Let us notice that when we use the iteration process for proving of the existence and uniqueness of the problem solution for Equation (2) (or of the whole set of equations (1)–(4)) then at each of the iterations, the equations with respect to the functions n, N, ϕ will be linear ones and the solutions of these equations (except Schrödinger equation) will exist due to theorems proved in [Citation7]. On the other hand, the existence and uniqueness of the solution for nonlinear heat conduction equations, including the nonlinear BCs, were proved in [Citation7].

Thus, summing the discussion mentioned above one can state that the solution of the problem (1)–(4) exists and this solution unique. We stress that for proving of this, one can use theorems containing in [Citation5,Citation7,Citation25] and applying the iteration process.

Under consideration of the FDS, we will assume that the functions G and R possess the Lipschitz continuity property.

3. Finite-difference scheme

To solve the differential initial-boundary problem (1)–(6) numerically we approximate it by the set of finite-difference equations. For this aim, we introduce in the domain G¯={0xLx}×{0yLy}×{0tLt} the uniform time and space grids: (68) ωx={xk=khx,i=0,Px¯,hx=Lx/Px},ωy={yj=jhy,j=0,Py¯,hy=Ly/Py},(68) ωy={yj=jhy,j=Py,Py¯,hy=Ly/Py},ωt={tm=mτ,τ=0,Pt¯,τ=Lt/Pt}, where hx, hy, τ are the grid steps on spatial coordinates and on time, correspondingly, Px, Py, Pt denote the number of grids nodes. Using the 1D grids defined in (68) we define temporal–spatial grids (69) Ω1=ωx×ωy×ωt,Ω2=ωx×ωy×ωt.(69) Let us define the mesh functions nh, Nh, ϕh, Gh, Rh, δh on the grid Ω1 in the following way: nkjm=n(xk,yj,tm), Nkjm=N(xk,yj,tm), ϕkjm=ϕ(xk,yj,tm), Gkjm=G(xk,yj,tm), Rkjm=R(xk,yj,tm), δkjm=δ(xk,yj,tm) and the mesh function Ah is defined on a grid Ω2, which is constructed using the grid ωy with the extended area along y coordinate: Akjm=A(xk,yj,tm) ().

Figure 4. Template of the mesh function definition.

Figure 4. Template of the mesh function definition.

Grid parameter δ0h (maximal value of the semiconductor laser energy absorption) for the parameter δ0 is defined on the grid Ω2 by the rule: δ0h=0,j=Py,1¯0.5δ0,j=0δ0,j=1,Py¯. For brevity, below we use the following index-free notations: (70) f=fh=fkjm=f(xk,yj,tm),fk±1=fk±1jm,fj±1=fkj±1m,fˆ=fkjm+1,(70) where f is one of the defined mesh functions and omit index h. We also define some useful notations using (14): (71) fk±0.5=0.5(fk+fk+1),fj±0.5=0.5(fj+fj+1),f0.5=0.5(fm+fm+1),(71) G0.5=q0|A|20.5δ0.5,R0.5=0.5(R(nˆ,Nˆ)+R(n,N)),δ0.5=0.5(δ(nˆ,Nˆ)+δ(n,N)),|A|20.5=0.5(|Aˆ|2+|A|2) The first and the second difference derivatives are defined in a standard way and they are notated as follows: fx, fx¯, fx¯x, fy, fy¯, fy¯y, ft.

For notation brevity, we introduce the finite-difference operators: Zx¯x(n)ϕ=μxhx(nk+0.5ϕxnk0.5ϕx¯), Zy¯y(n)ϕ=μyhy(nj+0.5ϕynj0.5ϕy¯), and the 2D Laplace difference operator: Λf=fx¯x+fy¯y.

It should be stressed one more, that below we provide all computations for zero-value of the external electric field Ex=Ey=0. However, the developed FDS is also an effective one for arbitrary BCs.

At FDS construction we use the Crank–Nicolson FDS. As a result, we propose the following FDS for the set of equations (1)–(4): (72) NˆNτ=G0.5R0.5,k=0,Px¯, j=0,Py¯,(72) (73) nˆnτ=κxn0.5x¯x0.5(Zx¯x(nˆ)ϕˆ+Zx¯x(n)ϕ)+κyn0.5y¯y0.5(Zy¯y(nˆ)ϕˆ+Zy¯y(n)ϕ)+G0.5R0.5,k=1,Px1¯,j=1,Py1¯,(73) (74) Λϕˆ=γ(nˆNˆ),k=1,Px1¯, j=1,Py1¯,(74) (75) AˆAτ+iDAxA0.5x¯x+iDAyA0.5y¯y+βAδ0h2δ0.5A0.5=0,k=1,Px1¯, j=Py+1,Py1¯.(75)

Let us note, that δ0.5kj=0 if j=Py+1,1¯.

The BCs (5) for the charged particles concentrations and electric field potential are approximated with the first order: (76) ϕˆx|k=0j=ϕˆx¯|k=Pxj=0, nˆx|k=0j=nˆx¯|k=Pxj=0, j=0,Py¯,(76) ϕˆy|kj=0=ϕˆy¯|kj=Py=0, nˆy|kj=0=nˆy¯|kj=Py=0, k=0,Px¯, Aˆ|x=0j=Aˆ|x=Pxj=0, j=Py,Py¯, Aˆ|kj=Py=Aˆ|kj=Py=0, k=0,Px¯.

The corresponding initial conditions for electron and donor concentrations are written as: (77) nkj0=Nkj0=n0,ϕkj0=0, k=0,Px¯, j=0,Py¯.(77) The complex amplitude initial condition is defined in the following way: (78) Akj0=e(xkLxc)2a2(yjLyc)102iπχ(yjLyc), k=0,Px¯, j=Py,Py¯.(78) At writing the initial condition (78), we assume that the corresponding mesh nodes coincide with coordinates Lxc, Lyc. This assumption does not restrict our consideration.

Theorem 3.1:

FDS (72)–(75) possesses the second order of approximation on spatial coordinates and on time coordinate in inner grid nodes concerning the point (xk,yj,(tm+τ/2)) on sufficient smooth solution of the problem (1)–(4). BCs (76) possess the first order of approximation.

For brevity, we omit the proof of Theorem 3.1 because this fact is easy to prove using the standard differential derivatives expansion in a Taylor series on sufficient smooth solution of the problem under consideration.

3.1. Conservativeness of the FDS with respect to the charge

Let us notice, that for the problem under consideration the FDS conservatism means a validity of the difference analogue of the conservation law (11). For numerical computation, we use the following its approximation: (79) Q(tm)=j=1Py1k=1Px1hyhx(nkjmNkjm).(79) At BCs approximation, we should take into account the requirement of the FDS conservatism. For this purpose, we have formulated and proved

Theorem 3.2:

The FDS (72)–(75) is a conservative one at the BCs approximation with the first order on spatial coordinates (76). It means, that a charge, defined in (79), is constant: Q(tm)=const.

Proof:

We can write below the sum, following from the equations concerning the free-electron concentration and ionized donor concentration: (80) j=1Py1k=1Px1(ntNt)hyhx=j=1Py1k=1Px1((nˆNˆ)(nN))hyhxτ=j=1Py1k=1Px1hyhx(κxn0.5x¯x+κyn0.5y¯y)0.5j=1Py1k=1Px1hyhx(κx(Zx¯x(nˆ)ϕˆ+Zx¯x(n)ϕ)κy(Zy¯y(nˆ)ϕˆ+Zy¯y(n)ϕ)).(80) Let us calculate the sums entering the right part of relation (80). For example, we see that two equalities are valid: k=1Px1hxnˆx¯x=1hxk=1Px1(nˆk12nˆ+nˆk+1)=(nˆx¯,Pxjnˆx,0j), k=1Px1hxZx¯x(nˆ)ϕˆ=μxi=1Px1nˆk+0.5ϕˆk+1ϕˆkhxnˆk05ϕˆkϕˆk1hx=μx(nˆPx0.5jϕˆx¯,Pxjnˆ0.5jϕˆx,0j). Obviously, that at the BCs approximation (76), the following equalities are valid: (81) j=1Py1k=1Px1hyhxnˆx¯x=j=1Py1hy(nˆx¯,Pxjnˆx,0j)=0,(81) j=1Py1k=1Px1hyhxZx¯x(nˆ)ϕˆ=μxj=1Py1hy(nˆPx0.5jϕˆx¯,Pxjnˆ0.5jϕˆx,0j)=0.

Similar equalities are valid also for other sums entering to the right part of expression (80): (82) j=1Py1k=1Px1hyhxnˆy¯y=k=1Px1hx(nˆy¯,kPynˆx,k0)=0,(82) j=1Py1k=1Px1hyhxZy¯y(nˆ)ϕˆ=μyk=1Px1hx(nˆPyk0.5ϕˆy¯,kPynˆk0.5ϕˆy,k0)=0. Therefore, the right part of equality (80) is equal to zero. Thus, the charge conservation law (79) takes place and FDS (72)–(75) is a conservative one.

Corollary 3.1:

Important consequence follows from Theorem 3.2 which consists in constancy of the sum (83) j=1Py1k=1Px1hyhx(nˆNˆ)=j=1Py1k=1Px1hyhx(nN)=0.(83) If we take into account the initial conditions (77) then the equality (83) can also be written in the form (84) ||nˆ||L1=||Nˆ||L1,(84) where mesh norm L1for the positive mesh function is defined as (85) ||u||L1=j=1Py1k=1Px1hyhxukj.(85) Applying similar way as we provide for equalities (80) one can obtain the following relation: (86) ||nˆ||L1+||1Nˆ||L1=LxLy=const.(86) At writing (86), we suppose that mesh function Nˆ does not exceed the unity. We see that sum of the norms for two mesh functions is bounded by constant. Thus, the mesh solution of the equations concerning free-electron concentration and ionized donor concentration will be stable in the norm L1.

3.2. Proving the positivity property for Nˆ

From difference equation (72) with respect to the ionized donor concentration on upper layer on time coordinate, we obtain the following expression with respect to the mesh functionNˆ: (87) Nˆ1+0.5τq0|A|20.5eψ+ξψnˆ+nˆτR=N10.5τq0|A|20.5eψ+ξψn+nτR+τ0.5q0|A|20.5eψ(eξψnˆ+eξψn)+n02τR.(87) Let us suppose that the mesh functions n and N are positive at previous time layer. In this case, the right part of Equation (87) will be positive if the following condition (88) τ2q0|A|max20.5eψ+ξ||n||C+||n||CτR,|A|max20.5=maxk,j|A|20.5k=1,Px1¯, j=1,Py1¯(88) is valid.

With respect to the left part of Equation (87), it is necessary to emphasize that the expression in big parentheses is positive at any sign of the mesh function nˆ due to choosing the grid step on time coordinate. Indeed, if the mesh function nˆ is positive then this term is positive also. If a value of nˆ is negative, then there are two situations. The first one corresponds to zero-value of a term in small parentheses is equal to zero then the mesh function Nˆ will be positive. For the second one to achieve a positive value of the mesh function Nˆ, it is sufficient choosing of the mesh step on time coordinate as: (89) τ2τR|minnˆ|,minnˆ=minj,knˆ,if nˆ<0.(89) Therefore, independent from a sign of the mesh function nˆ the value of Nˆwill be positive at a validity of the conditions (88), (89). We stress that these conditions are sufficient.

Theorem 3.3:

Function Nˆ is positive if the conditions (88)–(89) are valid in the hypothesis that N is positive.

3.3. Proving the inequality Nˆ<1

Let us re-write Equation (87) as (90) (1Nˆ)1+0.5τq0|A|20.5eψ+ξψnˆ=0.5τnˆNˆτR+(1N)10.5τq0|A|20.5eψ+ξψn+nτR+0.5τn2n02τR.(90) Because of choosing of the function N normalization, a value n0 is less than unity and at the initial time moment n=n0, the inequality n>n02 is valid at previous time layer. We suppose that the mesh function describing the concentration of the ionized donors at the low time layer is less than unity: N<1. Then, discussing in a similar way as we discuss a validity of positiveness for the mesh function Nˆ in (87), we conclude about a validity of the mesh function property Nˆ<1.

Theorem 3.4:

Function Nˆ is bounded by unit.

3.4. Assessment for the mesh function A in the mesh norm L2. Equation (75)

Let us introduce a scalar product for a mesh complex function with zero-value in boundary's nodes as (91) (u,v)=k=1Px1j=Py+1Py1hxhyuv=j=Py+1Py1k=1Px1hyhxukjvkj.(91) Therefore, the mesh norm L2is defined as (92) ||u||L22=(u,u)=k=1Px1j=Py+1Py1hxhyuu=k=1Px1j=Py+1Py1hxhy|u|2.(92) Multiplying Equation (75) on Akj0.5 and the equation, conjugated to Equation (75), – on Akj0.5, summing the result and multiplying it on the grid steps on spatial coordinates and then taking a sum over grid nodes, we obtain the following equation: (93) j=Py+1Py1k=1Px1hyhx|Aˆ|2|A|2τ+βAδ0hδ0.5|A|20.5=0(93) from which one can write (94) ||Aˆ||L22=||A||L22τβAδohj=Py+1Py1k=1Px1hyhxδ0.5|A|20.5.(94) We have to stress that for writing expression (93) we use the difference Green formulas for the mesh functions (see, e.g. [Citation24]): (95) (ux¯x,v)=(ux¯,vx¯)+uvx¯|Pxvux|0.(95) As we assumed above the mesh function δ(n,N) is bounded from below by the constant δmin=minj,kδ(n,N,nˆ,Nˆ). Therefore, the next inequality follows from (94): (96) ||Aˆ||L221+τ2βAδohδmin||A||L221τ2βAδohδmin.(96) Thus, we obtain that the norm L2 from the mesh complex function A at current layer on the time coordinate mesh is bounded by its value for initial distribution: (97) ||Aˆ||L22qm||A0||L22||A0||L22,q=1τ2βAδ0hδmin1+τ2βAδ0hδmin1,(97) if the following condition (98) τ<τ0=2βAδ0hδmin(98) is valid. It means Equation (75) solution stability.

Theorem 3.5:

The norm L2 of Equation (75) difference solution is bounded if the mesh step on time coordinate satisfies inequality (98).

Thus, the solution of the difference equation (75) is stable in the norm L2.

4. Iteration process for developed FDS

To solve the obtained set of 2D nonlinear difference equations we use two-stage iteration process. We construct it in such a way, that the FDS becomes a conservative one on each of iterations. It is an important feature of our approach because it allows avoiding disadvantages which arise from using the split-step method. It should be emphasized that in [Citation6] for the 2D problem without taking into account diffraction effects we had showed, that the split-step method using accumulates computing mistakes at calculating on a large time interval and, therefore, asymptotic stability property violation takes place.

Below we write these two stages separately with corresponding BCs. For iteration process, we use the following notations: upper indices s and s+1 mean the iteration number.

The first stage of the iteration process is: (99) Nˆs+1Nτ=G0.5s,s+1R0.5s,s+1,k=0,Px¯, j=0,Py¯,(99) (100) nˆs+1nτ=κxn0.5s+1x¯x+κyn0.5sy¯y+G0.5s,s+1R0.5s,s+10.5μxκxZx¯x(nˆs)ϕˆs+Zx¯x(n)ϕ+μyκyZy¯y(nˆs)ϕˆs+Zy¯y(n)ϕ,k=1,Px1¯, j=1,Py1¯,(100) (101) G0.5s,s+1=q0|A0.5s|2δ0.5s,s+1,R0.5s,s+1=0.5Rnˆs,Nˆs+1+R(n,N),(101) δ0.5s,s+1=0.5δnˆs,Nˆs+1+δ(n,N),|A|20.5s=0.5|Aˆ|2s+|A|2 (102) Λϕˆs+1=γnˆs+1Nˆs+1,k=1,Px1¯, j=1,Py1¯,(102) (103) Aˆs+1Aτ+iDAxAx¯x0.5s+1+iDAyAy¯y0.5s+βAδoh2δ0.5s,s+1A0.5s=0,k=1,Px1¯, j=Py+1,Py1¯,(103) (104) ϕˆs+1xk=0=ϕˆs+1x¯k=Px=0, j=0,Py¯, ϕˆs+1yj=0=ϕˆs+1y¯j=Py=0, k=0,Px¯,(104) nˆs+1xk=0=nˆs+1x¯k=Px=0, j=0,Py¯,Aˆ0js+1=AˆPxjs+1=0, j=Py,Py¯.

The second stage is: (105) Nˆs+2Nτ=G0.5s+1,s+2R0.5s+1,s+2,k=0,Nx¯, j=0,Ny¯,(105) (106) nˆs+2nτ=κxn0.5x¯xs+1+κyn0.5y¯ys+2+G0.5s+1,s+2R0.5s+1,s+20.5μxκxZx¯x(nˆs+1)ϕˆs+1+Zx¯x(n)ϕ+μyκyZy¯y(nˆs+1)ϕˆs+1+Zy¯y(n)ϕ,k=1,Px1¯, j=1,Py1¯,(106) (107) G0.5s+1,s+2=q0|A0.5s+1|2δ0.5s+1,s+2,R0.5s+1,s+2=0.5Rnˆs+1,Nˆs+2+R(n,N),(107) δ0.5s+1,s+2=0.5δnˆs+1,Nˆs+2+δ(n,N),|A|20.5s+1=0.5|Aˆ|2s+1+|A|2. (108) Λϕˆs+2=γnˆs+2Nˆs+2,k=1,Px1¯, j=1,Py1¯,(108) (109) Aˆs+2Aτ+iDAxAx¯x0.5s+1+iDAyAy¯y0.5s+2+βAδoh2δ0.5s+1,s+2A0.5s+1=0,k=1,Px1¯, j=Py+1,Py1¯,(109) (110) ϕˆs+2xk=0=ϕˆs+2x¯k=Px=0, j=0,Py¯, ϕˆs+2yj=0=ϕˆs+2y¯j=Py=0, k=0,Px¯,(110) nˆs+2yj=0=nˆs+2y¯j=Py=0, AˆkPxs+2=AˆkPys+2=0, k=0,Px¯. The iterations are stopped if the following inequalities (111) nˆs+2nˆsε1nˆs+ε2, Nˆs+2Nˆsε1Nˆs+ε2, ϕˆs+2ϕˆsε1ϕˆs+ε2,(111) Aˆs+2Aˆsε1Aˆs+ε2, ε1>0,ε2>0 are valid in all grid nodes. We check criterion (111) only after both iteration stages are made.

As initial approach for the iteration process we use the function values, obtained on the previous time layer: (112) nˆs=0=n, Nˆs=0=N, ϕˆs=0=ϕ, Aˆs=0=A.(112) It should be stressed, that the iteration process (99)–(110) possesses the conservatism property on each of the iterations.

Theorem 4.1:

The iteration process (99)–(110) possesses a conservativeness property on each of the iterations. It means, that a charge, defined in (79), is constant: Q(tm)=const.

4.1. Remark to the Poisson equation solution

Let us discuss the problem of 2D Poisson equation solution. For zero-value Dirichlet BCs, one can apply the split-step method in combination with FFT method. However, for the arbitrary BCs, the FFT cannot be used, obviously. In the case of the electric field potential evolution describing by linear 2D Poisson equation, the split-step method with a combination of another algorithm, which is caused by BCs, is an effective tool for solving the linear problem. Therefore, below we use this method also for solving the Poisson equation (102) with BCs (104) and describe a method application only at the first stage (99)–(104) of the main iteration process.

For the second stage of the iteration process, the same algorithm is used.

Because there are many approaches for multidimensional split-step method construction and they possess the same properties such as an accuracy order or the stability conditions [Citation21,Citation27], then in the present paper we use the stabilization method, which is written in the following manner for our problem: (113) Fp+1Fpτ¯=Fx¯xp+1+Fy¯ypγnˆs+1Nˆs+1,k=1,Px1¯, j=1,Py1¯,(113) (114) Fxp+1|k=0=Fx¯p+1|k=Px=0, j=0,Py¯,(114) (115) Fp+2Fp+1τ¯=Fx¯xp+1+Fy¯yp+2γnˆs+1Nˆs+1,k=1,Px1¯, j=1,Py1¯(115) (116) Fyp+2j=0=Fy¯p+2j=Py=0, k=0,Px¯,(116) where F is an auxiliary mesh function defined on the spatial grid Ω¯=ωx×ωy, τ¯ is an iteration parameter, which is not equal to time grid step τ, p is an iteration number. The initial approach for this iteration process is (117) F0=ϕˆs.(117) This method has the second order of an accuracy concerning the iteration parameter τ¯.

We use the Thomas algorithm to solve each of the difference equations (41) and (43). The convergence assessment for the iteration process (41)–(44) is given by the inequality (118) Fx¯xp+1+Fy¯yp+2γnˆs+1Nˆs+1<ε3, ε3>0,k=1,Px1¯, j=1,Py1¯(118) and it is a discrepancy for the difference Poisson equation (31). Using this criterion allows us to calculate the electric field potential with high accuracy and with good high-speed performance. If the criterion (46) is satisfied, then the iteration process (41)–(44) is stopped and we compute the electric field potential at s+1 iteration: (119) ϕˆs+1=Fp+2.(119)

Theorem 4.2:

The iteration method (113)–(117) converges if the following inequality τ¯τ0 is valid, where τ0=min{hx2,hy2}.

4.2 Proving the positivity property for 0<Nˆs+1<1

From Equation (99), one can write the following expression for the mesh function Nˆ at s+1iteration (Nˆs+1): (120) Nˆs+11+0.5τq0|A|20.5seψ+ξψnˆs+nˆsτR=N10.5τq0|A|20.5seψ+ξψn+nτR+0.5τq0|A|20.5seψeξψnˆs+eξψn+n02ττR.(120) Supposing that the mesh function N at previous time layer is positive and bounded: (121) 0<NQN,QN=const>0,(121) the mesh functions n, nˆs are positive and bounded also: (122) 0<nQn,Qn=const>0,0  nˆsQn+εn,εn>0,(122) and the mesh function A on previous time layer as well as on upper layer at s iteration is bounded also (123) ||A||C2QA||Aˆ||C2sQA+εA,(123) where (124) ||A||C=maxj,k|A|, j=Py+1,Py1¯k=1,Px1¯,(124) then from Equation (120) the positiveness of the mesh function Nˆs+1 on upper time layer occurs if the following inequality is valid (125) τ<2eψq0ImeξψQn(1+eξψnεn),(125) where (126) Im=0.5(2QA+εA).(126) For proving the inequality Nˆs+1<1 we re-write Equation (120) in the form: (127) (1Nˆs+1)1+0.5τq0|A|20.5seψ+ξψnˆs+nˆsτR=(1N)10.5τq0|A|20.5seψ+ξψn+nτR+0.5τnˆs+n2n02τR.(127) We see from (127) that at inequality (125) validity of the mesh function for ionized donors at s+1 iteration is less than unity (Nˆs+1<1). As above (see Section 3.3) a value n0 is less than unity and at the initial time moment n=n0, then the inequality n>n02 is valid at previous time layer. We suppose that the mesh function describing the concentration of ionized donors at the low time layer is less than unity: N<1. It should be stressed if even the last term in (127) becomes negative (let suppose this) then choosing of the grid step τ in appropriate way one can increase a value of the first term in the right part of Equation (127) to achieve a positiveness of the right part for Equation (127). Thus, we conclude about a validity of the mesh function property Nˆs+1<1.

Theorem 4.3:

Value of the mesh function Nˆ on upper time layer at s+1 iteration belongs to interval (0,1) if the ττ0 (see inequality (125)).

4.3 Uniform boundedness of Nˆs+1

From the previous paragraph, it follows that the mesh function Nˆs+1 is bounded by unity. Below we show the uniform boundedness of this function with respect to the iteration process. For this aim let us suppose that Equation (121) is valid. Let us suppose that the mesh function describing the concentration of the ionized donors on upper time layer at s-th iteration is bounded by (128) NˆsQN+εN.(128) Let us show that the next inequality (129) Nˆs+1QN+εN(129) is valid.

With this aim, Equation (120) is re-written in the form (130) Nˆs+1=N10.5τq0|A|20.5seψ+ξψn+nτR+0.5τq0|A|20.5seψ(eξψnˆs+eξψn)+n02ττR1+0.5τq0|A|20.5seψ+ξψnˆs+nˆsτR.(130) Thus, there is an inequality (131) Nˆs+1QN+τn02τR+0.5q0Im+0.5εAeψ(1ξQn)(1+eψξεn)QN+εN,(131) where instead (122) we introduce a new notation (132) ||n||CQn,||nˆ||CsQn+εn,(132) if does not suppose the positiveness of the mesh function for free-electron concentration at previous time layer and its value at the previous iteration. Therefore, if the grid step on time coordinate satisfies the inequality (133) τεNn02τR+0.5q0Im+0.5εAeψ(1ξQn)(1+eψξεn)2εNeψ(1ξQn)q0Im(1+eψξεn)(133) at a validity of inequality (125), the mesh function Nˆs+1is uniformly bounded.

This means a validity of the

Theorem 4.4:

The mesh function Nˆ is uniformly bounded on each of the iterations.

4.4. Difference of values at the iterations s+1, s and s+2, s+1 for the mesh function Nˆ

From Equation (99), one can write the following expression: (134) Nˆs+1Nˆs=0.5ττRNˆs+1nˆsNˆsnˆs10.5τq0eψNˆs+1|A|20.5seψξnˆsNˆs|A|20.5s1eψξnˆs1|A|20.5seψξnˆs|A|20.5s1eψξnˆs1.(134) Let us transform the terms in the right part of Equation (134). First of all, we transform the first two terms in (134). So, one can write (135) Nˆs+1nˆs±NˆsnˆsNˆsnˆs1=Nˆs+1Nˆsnˆs+Nˆsnˆs+1nˆs.(135) In a similar way, we transform the second term: (136) Nˆs+1|A|0.5seψξnˆsNˆs|A|20.5s1eψξnˆs1=Nˆs+1Nˆs|A|20.5s1eψξnˆs+Nˆseψξnˆs|A|20.5s|A|20.5s1eψξnˆs+Nˆs|A|20.5s1eψξnˆseψξnˆs1.(136) The last two terms can be transformed as: (137) |A|20.5seψξnˆs±|A|20.5s1eψξnˆs1|A|20.5s1eψξnˆs1=eψξnˆs|A|20.5s|A|20.5s1+|A|20.5s1eψξnˆseψξnˆs1.(137) Thus, Equation (134) could be written as (138) Nˆs+1Nˆs1+0.5τeψ+1τR=0.5τq0eψ((1Nˆs))|A|20.5s1eξψnˆseξψnˆs1+eξψnˆs|A|20.5s|A|20.5s10.5ττRnˆsnˆs1.(138) If we use Taylor's series for the exponential term (or a property of Lipchitz continuity) then (139) eξψnˆseξψnˆs1=Mn(nˆsnˆs1),(139) where Mn is a constant of the series. Finally, we obtain the inequality (140) |Nˆs+1Nˆs|0.5τq0eψ1+0.5τeψ+1τRnˆsnˆs1ImMn+eξψ(Qn+εn)|A|20.5s|A|20.5s10.5τq0eψnˆsnˆs1ImMn+eξψ(Qn+εn)|A|20.5s|A|20.5s1.(140) For the next neighbour iterations s + 2 and s + 1, we obtain the inequality which is similar to inequality (140): (141) Nˆs+2Nˆs+10.5τq0eψ|nˆs+1nˆs|ImMn+eξψ(Qn+εn)||A|20.5s+1|A|20.5s|.(141) Inequalities (140) and (141) will be used below for writing the conditions for the iteration process convergence.

4.5. Uniform boundedness of ||Aˆ||L22s+1

Let us remind that in Section 3.4 we showed the validity of the inequality for the norm L2 from the mesh function corresponding to the complex amplitude (142) ||Aˆ||L22CA||A0||L22,(142) if the grid step on time coordinate satisfies inequality (97). Because all norms in finite-dimensional spaces are equivalent [Citation1] then we proposed the validity of boundedness for the mesh function A on low layer in time coordinate and for the mesh function Aˆs on upper layer in time coordinate at s iteration (see (123)). Below we show that the assessment (143) ||Aˆ||Cs+1QA+εA(143) is valid for the mesh function Aˆs+1. It means the uniform boundedness of this function on upper grid layer on time coordinate. For this aim, let us write Equation (103) in the following form: (144) Aˆs+1=Bx1BxAiτDAyBx1Λy¯yA0.5sτβAδ04Bx1(1Nˆs)eψ(1ξnˆs)+(1N)eψ(1ξn)A0.5s,(144) where we introduce the following operators: (145) Bx=E+iτ2DAxΛx¯x,Bx=Eiτ2DAxΛx¯x.(145) Multiplying Equation (144) by the operator Bx1 from the left and taking into account the assessments for the operators [Citation14] (146) ||Bx1Bx||C1,||Bx1||C<1,||Λy¯y||Cc1hy2,c1=const>0(146) and boundedness of the mesh functions n, nˆs, N, Nˆs in this norm as well as ||Nˆs||c<1 we obtain: (147) ||Aˆ||Cs+1||A||C+τ(2QA+εA)c12hy2DAy+βAδ04eψ(1ξQn)(1+eψξεn).(147) Therefore, in Equation (143) will be valid if the mesh step on time coordinate satisfies to the inequality: (148) τ4hy2εA(2QA+εA)(c1DAy+hy2βAδ0eψ(1ξQn)(1+eψξεn)).(148) Carrying out the similar algebra for the mesh function Aˆs+2 we obtain from (109) the following equation: (149) Aˆs+2=By1ByAiτDAxBy1Λx¯xA0.5s+1τβAδ04By1(1Nˆs+1)eψ(1ξnˆs+1)+(1N)eψ(1ξn)A0.5s+1,(149) where we introduce the following operators: (150) By=E+iτ2DAyΛy¯y,By=Eiτ2DAyΛy¯y.(150) Multiplying Equation (149) by the operator By1 from the left and taking into account the assessments for the operators [Citation14] (151) ||By1By||C1,||By1||C<1,||Λx¯x||Cc2hx2,c2=const>0(151) and supposing that this assessment (152) ||nˆ||Cs+1Qn+εn(152) is valid (we discuss this statement below), we obtain (153) ||Aˆ||Cs+2||A||C+τ(2QA+εA)c22hx2DAx+βAδ04eψ(1ξQn)(1+eψξεn).(153) Therefore, the inequality (154) ||Aˆ||Cs+2QA+εA(154) will be valid if the mesh step on time coordinate satisfies to the inequality: (155) τ4hx2εA(2QA+εA)(c2DAx+hx2βAδ0eψ(1ξQn)(1+eψξεn)).(155) Thus, the following theorem is valid:

Theorem 4.5:

The mesh function Aˆ is uniformly bounded on iterations if inequalities (148) and (155) are valid.

4.6. Difference of values at the iterations s+1, s and s+2, s+1 for the mesh function Aˆ

Using the difference equation (103) for the mesh function Aˆs+1 and the difference equation, which is similar to (109), for the mesh function Aˆs we obtain the following difference equation with respect to difference of these mesh functions (156) Aˆs+1Aˆs+i0.5τDAx(Aˆx¯xs+1Aˆx¯xs+Aˆx¯xsAˆx¯xs1)+τβAδ02(δ0.5sA0.5sδ0.5s1A0.5s1)=0(156) or (157) Bx(Aˆs+1Aˆs)+i0.5τDAx(Aˆx¯xsAˆx¯xs1)+0.5τβAδ0(δ0.5sA0.5sδ0.5s1A0.5s1)=0.(157) To write the final difference equation we need to transform the last difference in Equation (157), which can be written as (158) (AˆsAˆs1)eψξnˆs(1Nˆs1)+(1N)eψξn+Aˆs1(1Nˆs1)(eψξnˆseψξnˆs1)(NˆsNˆs1)eψξnˆsAˆs+A(158) or using Equation (139), expression (158) transforms to (159) (AˆsAˆs1)eψξnˆs(1Nˆs1)+(1N)eψξn+(A+Aˆs1)(1Nˆs1)Mn(nˆsnˆs1)(NˆsNˆs1)eψξnˆs(Aˆ+A)(159) Thus, Equation (157) is written in the form: (160) Bx(Aˆs+1Aˆs)+0.5τiDAxΛx¯x+βAδ04eψξnˆs(1Nˆs1)+(1N)eψξn(AˆsAˆs1)+τβAδ08(A+Aˆs1)(1Nˆs1)Mn(nˆsnˆs1)(Aˆs+A)eψξnˆs(NˆsNˆs1)=0.(160) In a similar way, one can transform the second difference between the mesh function Aˆ at s+2 and s+1 iterations: (161) By(Aˆs+2Aˆs+1)+i0.5τDAyΛy¯y(Aˆy¯ys+1Aˆy¯ys)+τβAδ02(δ0.5s+1A0.5s+1δ0.5sA0.5s)=0.(161) After applying some algebra, we write the equation, which is similar to Equation (160), (162) By(Aˆs+2Aˆs+1)+0.5τiDAyΛy¯y+βAδ04eψξnˆs+1(1Nˆs)+(1N)eψξn(Aˆs+1Aˆs)+τβAδ08(A+Aˆs)(1Nˆs)Mn(nˆs+1nˆs)(Aˆs+A)eψξnˆs+1(Nˆs+1Nˆs)=0.(162) The next step consists of writing the inequalities for a difference of the mesh function Aˆ values at s+1 and s iterations, and s+2 and s+1 iterations. With this aim we re-writing Equations (160) and (162) by multiplying them by Bx1,By1 from the left, correspondingly, we obtain (163) (Aˆs+1Aˆs)+0.5τBx1iDAxΛx¯x+βAδ04eψξnˆs(1Nˆs1)+(1N)eψξn(AˆsAˆs1)+τβAδ08Bx1(A+Aˆs1)(1Nˆs1)Mn(nˆsnˆs1)(Aˆs+A)eψξnˆs(NˆsNˆs1)=0,(163) (164) (Aˆs+2Aˆs+1)+0.5τBy1iDAyΛy¯y+βAδ04eψξnˆs+1(1Nˆs)+(1N)eψξn(Aˆs+1Aˆs)+τβAδ08By1(A+Aˆs)(1Nˆs)Mn(nˆs+1nˆs)(Aˆs+A)eψξnˆs+1(Nˆs+1Nˆs)=0.(164) Then, taking into account the assessments for the corresponding operators, we write: (165) (Aˆs+1Aˆs)C0.5τDAxc2hx2+βAδ04eψξQn(1+eψξεn)(AˆsAˆs1)C+τβAδ08(2QA+εA)Mn(nˆsnˆs1)C+eψξ(Qn+εn)(NˆsNˆs1)C,(165) (166) (Aˆs+2Aˆs+1)C0.5τDAyc1hy2+βAδ04eψξQn(1+eψξεn)||(Aˆs+1Aˆs)||C+τβAδ08(2QA+εA)Mn(nˆs+1nˆs)C+eψξ(Qn+εn)(Nˆs+1Nˆs)C.(166) The last two inequalities together with inequalities (140) and (141) are used for a convergence condition of the differences between the mesh functions values at neighbour iterations.

4.7. Uniform boundedness of ||nˆs||C

Let us take a sum from both parts of Equation (100) multiplying them by hx, hy. Taking into account the equations corresponding to the boundary nodes, we obtain: (167) j=1Py1k=1Px1hxhynˆs+1nτ=0.5j=1Py1k=1Px1hxhy((1Nˆs+1)eψ+ξψnˆs+(1N)eψ+ξψn)|A|20.5s0.5Nˆs+1nˆsn02τR+Nnn02τR.(167) Let us emphasize that a value of the mesh function Nˆs+1 can be expressed from Equation (99) via values of the mesh function at the previous iteration. Moreover, in Section 4.3 we showed the uniform boundedness of this mesh function. We take into account below that the mesh function Nˆs+1 is positive and bounded. In this case, from (167) the inequality for nˆs+1 follows: (168) ||nˆ||L1s+1||n||L1+0.5τeψeξψnˆsL1+eξψnL1Im+ττRn02LxLy+0.5ττR||n||L1+||nˆ||sL1.(168) At writing the ||nˆ||L1s+1 we suppose that nˆs+1 is positive. Some discussion about this statement will be present in the end of Section 4. In general case, it is necessary to use modulus from this mesh function.

Because in expression (122) we suppose the boundedness of the mesh functions n,nˆs in the norm C then we suppose that these mesh functions will be bounded in the norm L1 also in accordance with the equivalence property of all finite-dimensional spaces [Citation1]. Therefore, we state: (169) ||n||L1Q¯n,Q¯n=const>0,||nˆs||L1Q¯n+ε¯n,ε¯n>0.(169) Let notice that other reason according to which the norm ||n||L1is bounded arises from a validity of the invariant Q(tm)=const (see (79)) and boundedness of the norm ||N||L1. With accordance to Theorem 4.1 (see Section 4) and boundedness of the norm ||Nˆs||L1, the ||nˆ||L1s is bounded also.

Thus, the norm ||nˆ||L1s+1 is bounded by the same constant Q¯n+ε¯n, if the mesh step on time coordinate satisfies the inequality: (170) τ2eψε¯neξψnˆsL1+eξψnL1(2QA+εA)+2τRn02LxLy+1τR2Q¯n+ε¯n.(170) The similar assessment can be obtained with respect to the mesh function nˆs+2 (106): (171) ||nˆ||L1s+2||n||L1+0.5τeψeξψnˆs+1L1+eξψnL1(2QA+εA)+ττRn02LxLy+0.5ττR||n||L1+nˆs+1L1.(171) From this inequality, one can conclude that (172) ||nˆ||L1s+2Q¯n+ε¯n(172) if the grid step on time coordinate satisfies the inequality: (173) τ2eψε¯neξψnˆs+1L1+eξψnL1(2QA+εA)+2τRn02LxLy+1τR(2Q¯n+ε¯n).(173) Finally, the following theorem is valid:

Theorem 4.6:

The mesh function nˆ at iterations will uniformly bounded in the norm L1, if the mesh step on time coordinate is satisfied both inequalities (170) and (173).

4.8. Difference of values at the iterations s+1, s and s+2, s+1 for the mesh function nˆ

From the iteration process (99)–(104), the following equalities follows (174) j=1Py1k=1Px1hxhy(nˆs+1nˆs)=j=1Py1k=1Px1hxhy(Nˆs+1Nˆs)=τj=1Py1k=1Px1hxhy(G0.5s,s+1R0.5s,s+1),(174) (175) j=1Py1k=1Px1hxhy(nˆs+2nˆs+1)=j=1Py1k=1Px1hxhy(Nˆs+2Nˆs+1)=τj=1Py1k=1Px1hxhy(G0.5s+1,s+2R0.5s+1,s+2),(175) in accordance to Theorem 4.1 (see Section 4). We use these relations for assessment of the geometric ratio for the iteration process (99)–(110). Using these equalities, one can write the following equations: (176) (10.5τκxΛx¯x)(nˆs+1nˆs)=0.5τκxΛx¯x(nˆsnˆs1)+((Nˆs+1Nˆs)(NˆsNˆs1))0.5τμxκxZx¯x(nˆs)ϕˆsZx¯xnˆs1ϕˆs1+μyκyZy¯ynˆsϕˆsZy¯ynˆs1ϕˆs1,(176) and (177) (10.5τκyΛy¯y)(nˆs+2nˆs+1)=0.5τκyΛy¯y(nˆs+1nˆs)+((Nˆs+2Nˆs+1)(Nˆs+1Nˆs))0.5τμxκxZx¯xnˆs+1ϕˆs+1Zx¯xnˆsϕˆs+μyκyZy¯ynˆs+1ϕˆs+1Zy¯ynˆsϕˆs.(177) If we introduce the operators (178) Bnx=10.5τκxΛx¯x,Bny=10.5τκyΛy¯y,(178) Equations (176), (177) can be re-written as (179) (nˆs+1nˆs)=0.5τκxBnx1Λx¯x(nˆsnˆs1)+Bnx1((Nˆs+1Nˆs)(NˆsNˆs1))0.5τBnx1μxκxZx¯xnˆsϕˆsZx¯xnˆs1ϕˆs1+μyκyZy¯ynˆsϕˆsZy¯ynˆs1ϕˆs1,(179) (180) (nˆs+2nˆs+1)=0.5τκyBny1Λy¯y(nˆs+1nˆs)+Bny1((Nˆs+2Nˆs+1)(Nˆs+1Nˆs))0.5τBny1μxκxZx¯xnˆs+1ϕˆs+1Zx¯xnˆsϕˆs+μyκyZy¯ynˆs+1ϕˆs+1Zy¯ynˆsϕˆs.(180) Thus, in the norm C, we obtain (181) (nˆs+1nˆs)C=0.5τκxc2hx2(nˆsnˆs1)C+Nˆs+1NˆsC+NˆsNˆs1C+0.5τ(Qn+εn)×μxκxhx(Q¯n+ε¯n+Q¯N+ε¯N)+μyκyhy(Q¯n+ε¯n+Q¯N+ε¯N),(181) (182) (nˆs+2nˆs+1)C=0.5τκyc1hy2(nˆs+1nˆs)C+Nˆs+2Nˆs+1C+Nˆs+1NˆsC+0.5τ(Qn+εn)×μxκxhx(Q¯n+ε¯n+Q¯N+ε¯N)+μyκyhy(Q¯n+ε¯n+Q¯N+ε¯N).(182) Summing all conditions for uniform boundedness of the mesh functions at the iterations and for the differences of the mesh functions we obtain the relation between the mesh steps (183) τminC1hy2,C2hx2,C1,C2>0,(183) where the parameters C1,C2 depend on the problem parameters.

Theorem 4.7:

The iteration process (99)–(110) converges with the geometric ratio q=minτ{(C1/hy2),(C2/hx2)} if inequality (183) is satisfied.

We showed that the mesh functions are uniformly bounded at each of iterations and the iteration process converges. It means validity of the following theorem:

Theorem 4.8:

The solution of the difference problem (1)–(4) existed on upper time layer.

Remark:

A positiveness of the mesh function nˆ can be stated by using the approach developed in [Citation3]. This requires additional transform of the mesh equation with respect to mesh concentration nˆ and, as we believe, it is a subject of other paper.

5. Computer simulation results

As an illustration of the FDS efficiency we show in Figures the laser pulse intensity distribution evolution (Figures , , ) and the free-electron concentration evolution (Figure ) computed for a set of the problem parameters la=0.1,δ0=0.25,γ=1000,τp=1,q0=4,n0=0.01,ξ=3,ψ=2.553,κx=κy=105,μx=μy=1,DAx=105,χ=5,Lxc=0.5,Lyc=2,Lx=1,Ly=30. Grid steps on spatial coordinates and on time are chosen as hx=hy=102, τ=103. The iteration parameters are chosen as ε1=104,ε2=106. It should be stressed, that we provide our computations for extended area on y-coordinate, but in figures we use a decreased area, to optimize the computer simulation results visualization.

Figure 5. Distribution of square root from laser pulse intensity |A|2 at time moment t = 0 (a),1.4 (b), 3.6(c), 5(d), 10(e), 20(f).

Figure 5. Distribution of square root from laser pulse intensity |A|2 at time moment t = 0 (a),1.4 (b), 3.6(c), 5(d), 10(e), 20(f).

Figure 6. Square root from beam profile |A|2 along y-coordinate at the beam axis on x-coordinate (x = 0.5) at time moment t = 0 (a), 1.4 (b), 3.6(c), 5(d), 10(e), 20(f).

Figure 6. Square root from beam profile |A|2 along y-coordinate at the beam axis on x-coordinate (x = 0.5) at time moment t = 0 (a), 1.4 (b), 3.6(c), 5(d), 10(e), 20(f).

Figure 7. Distribution of free-electron concentration n at time moment t = 3.2 (a), 3.4 (b), 3.6 (c), 4 (d), 5 (e), 10 (f).

Figure 7. Distribution of free-electron concentration n at time moment t = 3.2 (a), 3.4 (b), 3.6 (c), 4 (d), 5 (e), 10 (f).

Figure 8. Square root from beam profile |A|2 along y-coordinate at the beam axis on x-coordinate (x = 0.5) at time moment t = 1.4 (a), 3 (b), 3.6(c), 5(d), 10(e), 20(f).

Figure 8. Square root from beam profile |A|2 along y-coordinate at the beam axis on x-coordinate (x = 0.5) at time moment t = 1.4 (a), 3 (b), 3.6(c), 5(d), 10(e), 20(f).

We see in Figure  that a part of laser beam reflects from a semiconductor face in time moment t=3.6. The laser beam transmitted through a boundary of a semiconductor with air acquires very complicated profile. The travelling part of the beam penetrated in a semiconductor becomes with intensity minimum at the x-axis near the face of a semiconductor (Figure (b,c)). With time increasing, the beam profile becomes much more disturbed one in comparison with the incident Gaussian profile. The reflected beam possesses also complicated profile in the x-direction as well as in y-direction. (Figure (d)). A part of the reflected beam moves faster than other ones (Figure ). For the beam propagating in a semiconductor, the pulse front moves faster than a basic part of the beam. For illustration of the complicated spatial–temporal intensity distribution, in Figure  the beam profile along y-coordinate at beam axis on x-coordinate is depicted at various time moments. We see a formation of the laser sub-pulses propagating in both directions along y-coordinate.

Analysis of the free-electron spatial–temporal distribution (Figure ) shows that the laser beam diffraction action results in multi-domains formation for high absorption (Figure (b,c)). We see an appearance of the free-electron high concentration domains on both coordinates (Figure (b)). For example, such domains form in the section y = 0.12 which cannot appear without taking into account the beam diffraction.

Figure  illustrates the influence of the reflected optical pulse from the high concentration domain. As a result of this, we see that the pulse intensity increases near the face of the semiconductor with increasing time.

In the values of invariant (79) and the difference analogue of integral equality (29) (184) Q¯¯(tm)=j=Py+1Py1k=1Px1hyhx|Aˆkjm|2|Akjm|2τ+βAδ0hδ0.5|Akjm|20.5(184) are shown. It should be noticed, that invariant (29) approximated by formula (184) with the second order on spatial and time coordinates. As we can see, both characteristics Q(tm) and Q¯¯(tm) preserve their values with high accuracy. Thus, a conservatism property which was proved analytically in Theorem 3.1 is confirmed also by computer simulation. It means, that our numerical method allows providing computations on long time interval without accumulating the rounding errors and this method possesses the asymptotic stability property.

Table 1. Invariants Q(tm) and Q¯¯(tm) values computed using formulas (79) and (184), correspondingly.

6. Conclusions

A new mathematical model for the resonatorless OB scheme based on a semiconductor is proposed. The main feature of this model consists in taking into account the longitudinal diffraction of the optical beam. Its action results in a laser beam reflection from the boundaries of high absorption domains. The transverse diffraction induces the laser beam reshaping.

To solve the 2D nonlinear Schrödinger equation together with the set of the equations describing an electron–hole plasma evolution in a semiconductor, the conservative FDS on the base of Crank–Nicolson scheme is proposed and its conservatism is proved. For its realization, the two-stage iteration process is used. Such iteration method allows to achieve the conservatism property on each of the iterations.

We showed the uniform boundedness of the mesh functions at the iterations, the convergence of the iteration process. We showed also positiveness and boundedness of the mesh function corresponding to free-charged particle concentrations. We discuss some properties of the differential problem including the problem invariants, positiveness of the free-charge concentrations.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This paper was made using the support of the Russian Science Foundation [grant number 14-21-00081].

References

  • G. Bachman and L. Narici, Functional Analysis, Dover Publications, Inc., Mineola, NY, 1998.
  • L. Barletti, L. Brugnano, G. Frasca Caccia, and F. Iavernaro, Energy-conserving methods for the nonlinear Schrödinger equation, Appl. Math. Comput. 318 (2018), pp. 3–18.
  • O.S. Bondarenko, Mathematical modeling of some problems of optical bistability on the basis of semiconductors, PhD thesis (1997) (in Russian).
  • M. Dehghan and A. Taleei, A compact split-step finite difference method for solving the nonlinear Schrödinger equations with constant and variable coefficients, Comput. Phys. Commun. 181 (2010), pp. 43–51. doi: 10.1016/j.cpc.2009.08.015
  • J.I. Dıaz and G. Galiano, Existence and uniqueness of solutions of the Boussinesq system with nonlinear thermal diffusion, topological methods in nonlinear analysis, J. Juliusz Schauder Center 11 (1998), pp. 59–82.
  • D. Dockery and J.R. Kuttler, An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation, IEEE Trans. Antennas Propag. 44 (1996), pp. 1592–1599. doi: 10.1109/8.546245
  • A. Friedman, Partial Differential Equations of Parabolic Type, Dover Publications, Inc., Mineola, NY, 2008.
  • A. Ghadi and S. Mirzanejhad, Two-photon absorption effect on semiconductor microring resonators, Optik (Stuttg) 126 (2015), pp. 1645–1649. doi: 10.1016/j.ijleo.2015.05.058
  • H.M. Gibbs, Optical Bistability: Controlling Light with Light, Academic Press, New York, 1985.
  • H.M. Gibbs, S.L. McCall, and T.N.C. Venkatesan, Optical bistable devices – the basic components of all-optical systems, Opt. Eng. 19(4) (1980), pp. 463–468. doi: 10.1117/12.7972544
  • H.M. Gibbs, S.L. McCall, T.N.C. Venkatesan, A.C. Gossard, A. Passner, and W. Wiegmann, Optical bistability in semiconductors, Appl. Phys. Lett. 35(6) (1979), pp. 451–453. doi: 10.1063/1.91157
  • W. Huang, M. Ji, Z. Liu, and Y. Yi, Steady states of Fokker–Planck equations: I. Existence, J. Dynam. Differential Equations 27 (3–4) (2015), pp. 721–742.doi: 10.1007/s10884-015-9454-x
  • A. Hurtado, M. Nami, I. D. Henning, M.G. Adams, and L.F. Lester, Bistability patterns and nonlinear switching with very high contrast ratio in a 1550 nm quantum dash semiconductor laser, Appl. Phys. Lett. 101 (2012), pp. 161117. doi: 10.1063/1.4761473
  • Yu N. Karamzin, A.P. Sukhorukov, and V.A. Trofimov, Mathematical Modeling in Nonlinear Optics, Moscow University Publishing, Moscow, 1989, p. 156 (in Russian).
  • P.K. Kwan and Y.Y. Lu, Computing optical bistability in one-dimensional nonlinear structures, Optics Commun. 238 (2004), pp. 169–175. doi: 10.1016/j.optcom.2004.04.025
  • O. Ladyzhenskaya and N.N. Uraltseva, Linear and Quasilinear Elliptic Equations, Nauka, Moscow, 1973 (in Russian).
  • T.I. Lakoba, Instability of the finite-difference split-step method applied to the generalized nonlinear Schrödinger equation. III. External potential and oscillating pulse solutions, Numer. Methods Partial Differential Equations, 33 (2017) pp. 633–650. doi: 10.1002/num.22071
  • A.W. Leung, Nonlinear Systems of Partial Differential Equations: Applications to Life and Physical Sciences, Word Scientific Publishing Co. Pvt. Ltd, Hackensack, NJ, 2009.
  • L. Li, Optical bistability in semiconductor lasers under intermodal light injection, J.f Quantum Electronics 32 (1996), pp. 248–256. doi: 10.1109/3.481872
  • J. Loustau, Numerical Differential Equations: Theory and Technique,ode Methods, Finite Differences, Finite Elements and Collocation, WordScientific Publishing Co. Pte. Ltd., New Jersey, 2016.
  • G.I. Marchuk, Methods of Numerical Mathematics, New York, NY, Springer, 1975.
  • J.D. Murray, Mathematical Biology II: Spatial Models and Biomedical Application, Springer-Verlag, New York, 2003.
  • R. Nasehi, Ultra-low threshold optical bistability and multi-stability in dielectric slab doped with semiconductor quantum well, Commun. Theor. Phys. 66 (2016), pp. 129–132. doi: 10.1088/0253-6102/66/1/129
  • A.A. Samarskii, Difference Schemes Theory, Nauka, Moscow, 1977 (in Russian).
  • A.A. Samarskii, V.A. Galaktionov, S.P. Kurdumov, and A.P. Mikhailov, The Modes with Aggravation in Problems for the Quasilinear and Parabolic Equations, Nauka, Moscow, 1987 (in Russian).
  • R.A. Smith, Semiconductors, Cambridge Univ. Press, Cambridge, 1978.
  • G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (3) (1968), pp. 506–517. doi: 10.1137/0705041
  • A. Suryanto, E. van Groesen, and M. Hammer, Finite element analysis of optical bistability in one-dimensional nonlinear photonic band gap structures with a defect, J. Nonlinear Optic. Phys. Mater. 12 (2003), pp. 187. doi: 10.1142/S0218863503001328
  • A. Taleei and M. Dehghan, Time-splitting pseudo-spectral domain decomposition method for the soliton solutions of the one- and multi-dimensional nonlinear Schrödinger equations, Comput. Phys. Commun. 185 (2014), pp. 1515–1528. doi: 10.1016/j.cpc.2014.01.013
  • Z. Toroker and M. Horowit, Optimized split-step method for modeling nonlinear pulse propagation in fiber Bragg gratings, JOSA B 25 (2008), pp. 448–457. doi: 10.1364/JOSAB.25.000448
  • S. K. Tripathy and S. Swain, Optical bistable switching in semiconductor heterostructure containing a quantum dot layer: The effect of phonons, Optik (Stuttg) 124 (2013), pp. 2723–2726. doi: 10.1016/j.ijleo.2012.08.088
  • V.A. Trofimov, M.M. Loginova, and V.A. Egorenkov, 2D wave structure induced by femtosecond laser pulse in semiconductor, Adv. Optoelectronics Lasers, Kharkov (2013), pp. 296–298.
  • V.A. Trofimov, M.M. Loginova, and V.A. Egorenkov, Developing of 2D helical waves in semiconductor under the action of femtosecond laser pulse and external electric field, Proc. SPIE, 9586 (2015), 95860 K. doi: 10.1117/12.2189308
  • V.A. Trofimov, M.M. Loginova, and V.A. Egorenkov, New two-step iteration process for solution of semiconductor plasma generation problem with arbitrary boundary conditions in 2D case, WIT trans. Model. Simul. 59 (2015), pp. 85–96.
  • V.A. Trofimov, M.M. Loginova, and V.A. Egorenkov, Conservative finite-difference scheme for computer simulation of field optical bistability, Lecture Notes Comput. Sci. 10187 (2017), pp. 682–689. doi: 10.1007/978-3-319-57099-0_78
  • H. Wang, Numerical studies on the split-step finite difference method for nonlinear Schrödinger equations, Appl. Math. Comput. 170 (2005), pp. 17–35.
  • H. Wang, X. Ma, J. Lu, and W. Gao, An efficient time-splitting compact finite difference method for gross–Pitaevskii equation, Appl. Math. Comput. 297 (2017), pp. 131–144.
  • H. Wang, H.-T. Zhang, and Z.-P. Wang, Optical bistability via incoherent pumping fields in semiconductor quantum wells, Mod Phys Lett B 25 (2011), pp. 97.
  • D. Weaire and J. P. Kermode, Dispersive optical bistability – numerical methods and definitive results, JOSA B 3 (1986), pp. 1706–1711. doi: 10.1364/JOSAB.3.001706
  • J.A.C. Weideman and B.M. Herbst, Split-step methods for the solution of the nonlinear Schrödinger equation, SIAM J. Numer. Anal. 23 (1986), pp. 485–507. doi: 10.1137/0723033