Conservative finite-difference scheme for the 2D

International Journal of Computer Mathematics

ISSN: 0020-7160 (Print) 1029-0265 (Online) Journal homepage: http://www.tandfonline.com/loi/gcom20

Conservative finite-difference scheme for the 2D problem of femtosecond laser pulse interaction with kink structure of high absorption in semiconductor Vyacheslav A. Trofimov, Maria M. Loginova & Vladimir A. Egorenkov To cite this article: Vyacheslav A. Trofimov, Maria M. Loginova & Vladimir A. Egorenkov (2018): Conservative finite-difference scheme for the 2D problem of femtosecond laser pulse interaction with kink structure of high absorption in semiconductor, International Journal of Computer Mathematics, DOI: 10.1080/00207160.2018.1492117 To link to this article: https://doi.org/10.1080/00207160.2018.1492117

© 2018 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group. Accepted author version posted online: 02 Jul 2018. Published online: 07 Jul 2018. Submit your article to this journal

Article views: 53

View Crossmark data

Full Terms & Conditions of access and use can be found at http://www.tandfonline.com/action/journalInformation?journalCode=gcom20

INTERNATIONAL JOURNAL OF COMPUTER MATHEMATICS https://doi.org/10.1080/00207160.2018.1492117

Conservative finite-difference scheme for the 2D problem of femtosecond laser pulse interaction with kink structure of high absorption in semiconductor Vyacheslav A. Trofimov, Maria M. Loginova and Vladimir A. Egorenkov Faculty of Computational Mathematics & Cybernetics, Lomonosov Moscow State University, Russian Federation ABSTRACT

ARTICLE HISTORY

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.

Received 19 September 2017 Revised 4 May 2018 Accepted 1 June 2018 KEYWORDS

Conservative finite-difference scheme; femtosecond laser pulse; nonlinear absorption; diffraction of laser beam; kink structure 2010 MSC CLASSIFICATIONS

65N06; 65N12; 35Q60

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 [10,11] and continued, for example, in [8,9,13,15,19,23,28,31,32,33,38,39]. 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 (10−15 c) laser pulse interaction with a semiconductor under the condition of absorption OB occurrence is described by the set of nonlinear PDEs [9]. For numerical solution of such problems, the split-step method is widely used [2,4,6,17,29,30,36,37,40]. For example, in [40] this method is used for a numerical solution of the nonlinear Schrödinger equation. It is shown, CONTACT Vyacheslav A. Trofimov

[email protected]

© 2018 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group. This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivatives License (http:// creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited, and is not altered, transformed, or built upon in any way.

2

V. A. TROFIMOV ET AL.

how a conservation law and instability are reflected in the numerical finite-difference schemes (FDS). In [37] an efficient time-splitting compact finite-difference method for Gross–Pitaevskii equation is proposed. In [29] 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 [4] 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 [2] 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 [17] 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 [6,40], 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 [32,33,34,35] 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 [34] 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 (Figure 1). 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) [9,26] and a laser pulse complex amplitude A(x, y, t), which is slowly varying in time only: ∂N = G(n, N, |A|2 ) − R(n, N), 0 < x < Lx , 0 < y < Ly , t > 0, ∂t ∂n ∂ ∂n ∂ ∂n ∂ϕ ∂ϕ = κx − μx n + κy − μy n + G(n, N, |A|2 ) − R(n, N), ∂t ∂x ∂x ∂x ∂y ∂y ∂y ∂ 2ϕ ∂ 2ϕ + 2 = γ (n − N), 2 ∂x ∂y

(1) (2) (3)

INTERNATIONAL JOURNAL OF COMPUTER MATHEMATICS

3

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).

∂A ∂ 2A ∂ 2 A β A δ0 + iDAx 2 + iDAy 2 + δ(n, N)A = 0, ∂t ∂x ∂y 2

0 < x < Lx ,

−Ly < y < Ly ,

t > 0. (4)

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:

∂n ∂ϕ ∂n ∂ϕ = n = 0. − μx n − μ y ∂x ∂x x=0,Lx ∂y ∂y y=0,Ly

(5)

The non-zero electric field strength Ex , Ey corresponds to the external electric field presence: ∂ϕ = −Ex , ∂x x=0,Lx

∂ϕ = −Ey . ∂y y=0,Ly

(6)

We stress that below we consider a case of the external electric field is absent: Ex = Ey = 0.

(7)

The zero-value BCs with respect to a complex amplitude correspond to its finite distribution: 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

4


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).

laser pulse propagates in a linear medium before its interaction with a semiconductor: ϕ|t=0 = 0, n|t=0 = N|t=0 = n0 ,

A|t=0 = A0 = e−(x−Lxc )

2 /a2 −(y−L

yc )

10 −2iπ χ (y−L ) yc

,

(9)

where Lxc , Lyc are coordinates of the incident laser beam centre position, a is a beam radius on xcoordinate. As the initial laser beam profile, which is shown in Figure 2, 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. Functions G and R describe charged particle’s generation and recombination [9,26] in the area of semiconductor, correspondingly: G(n, N, |A|2 ) = q0 |A|2 δ(n, N), nN − n20 , 0 ≤ x ≤ L x , 0 ≤ y ≤ Ly , t > 0 τR otherwise, these functions are equal to zero: R(n, N) =

G(n, N, |A|2 ) = R(n, N) = 0,

0 ≤ x ≤ Lx , −Ly ≤ y < 0, t > 0.

(10) (11)

(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 [9]: δ(n, N) = (1 − N)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


5

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: N ≤ 1. First, we stress that the derivatives from the generation function on the concentrations n, N are: ∂G = ξ ψG ≥ 0, ∂n

if G ≥ 0,

(14)

∂G (15) = −e−ψ(1−ξ n) |A|2 ≤ 0. ∂N 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: ∂R ∂R = N, =n (16) ∂n ∂N 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: δ(n, N) ≤ δmax = e−ψ+ξ ψnmax ,

nmax = max n(x, y, t) x,y,t

(17)

for enough small time interval t ∈ [0, tε ]. More weak assessment can be obtained taking into account the equation [18,22]: ∂Jy ∂ ∂Jx (n + (1 − N)) = Dx + Dy , ∂t ∂x ∂y

(18)

where Jx =

∂n ∂ϕ − μx n , ∂x ∂x

Jy =

∂n ∂ϕ − μy n ∂y ∂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

6


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

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: δ(n, N) = δ(n) = (1 − n)e−ψ(1−ξ n) .

(19)

To find the maximum of this function we take a first derivative, which is ∂δ = (ψξ(1 − n) − 1)e−ψ(1−ξ n) . ∂n

(20)

At a value of the concentration n∗ = 1 − 1/(ψξ ), 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 ∂ 2 δ(n∗ ) = −ψξ e−ψ+ψξ −1 < 0. ∂n2

(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 Figure 3. We see in Figure 3 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: ⎧ ⎨ 1 e−ψ+ψξ , ψξ > 1, ψξ δmax = (22) ⎩e−ψ+ψξ , ψξ ≤ 1.

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 Ly Lx Q(t) = (n(x, y, t) − N(x, y, t))dx dy = const, 0

(23)

0

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


7

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: Ly Lx Ly Lx ¯ Q(t) = (n(x, y, t) + (1 − N(x, y, t)))dxdy = (n0 + (1 − n0 ))dxdy = Lx Ly . 0

0

0

0

(24)

If we define a norm L1 for non-negative function as ||f ||L1 =

Ly

0

Lx

0

f (x, y)dxdy,

and suppose non-negativeness of the functions n(x, y, t), 1 − N(x, y, t) (for N this assumption was considered above) then one can re-write quality (24) in the following way: ||n||L1 + ||1 − N||L1 = ||n0 ||L1 + ||1 − n0 ||L1 = Lx Ly .

(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: ||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: ∂ ∂t

∂ 2 A∗ A |A| dx dy + iDAx − A 2 dx dy ∂x2 ∂x −Ly 0 −Ly 0 Ly Lx Ly Lx ∂ 2A ∂ 2 A∗ + iDAy A∗ 2 − A 2 dx dy + βA δ0 δ(n, N)|A|2 dx dy = 0. ∂y ∂y −Ly 0 −Ly 0

Ly

Lx

Ly

2

Lx

∗∂

2A

(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:

Ly −Ly

Lx 0

A∗

∂ 2A ∂ 2 A∗ −A 2 2 ∂x ∂x

dx dy = 0,

Ly −Ly

Lx 0

A∗

∂ 2A ∂ 2 A∗ −A 2 2 ∂y ∂y

dx dy = 0.

(28)

Thus, we obtain the following integral equality (energy invariant):

Ly −Ly

Lx 0

∂ |A|2 + βA δ0 δ(n, N)|A|2 dx dy = 0. ∂t

(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

8


function δ(n, N) (see also (17)): δmin ≤ δ(n, N) ≤ δmax ,

δmin , δmax > 0.

(30)

In this case, from Equation (13), one can write two inequalities: ∂ ∂t

∂ ∂t

Ly

Lx

−Ly

0

Ly

−Ly

|A| dx dy + βA δ0 δmin

Lx

0

Ly

2

Lx

−Ly

0

|A|2 dx dy + βA δ0 δmax

Ly −Ly

0

|A|2 dx dy ≤ 0,

(31)

|A|2 dx dy ≥ 0.

(32)

If we introduce the norm L2 for the complex amplitude A as ||A||2L2 = (A, A) =

Ly

−Ly

Lx 0

AA∗ dx dy =

Ly

−Ly

Lx 0

|A|2 dx dy,

(33)

which characterizes the pulse energy, then we obtain the following assessment of its evolution in time ||A||2L2 ≤ ||A0 ||2L2 e−βA δ0 δmin t ,

(34)

||A||2L2 ≥ ||A0 ||2L2 e−βA δ0 δmax t .

(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: Ly Lx 2 ∂A ∂A∗ ∂ A ∂A∗ ∂A∗ ∂A ∂ 2 A∗ ∂A − dx dy + iDAx + dx dy ∂t ∂x ∂t ∂x ∂x2 ∂x ∂x2 ∂x −Ly 0 −Ly 0 Ly Lx 2 ∂ A ∂A∗ ∂ 2 A∗ ∂A + iDAy + dx dy ∂y2 ∂x ∂y2 ∂x −Ly 0 Ly Lx ∂A∗ ∗ ∂A δ(n, N) A −A dx dy = 0. + 0.5βA δ0 ∂x ∂x −Ly 0

Ly

Lx

(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: ∂ ∂t

Ly −Ly

Lx 0

∂A∗ Im A ∂x

dxdy + βA δ0

Ly −Ly

0

Lx

∂A∗ δ(n, N)Im A ∂x

dxdy = 0,

(37)


or

t Ly Lx ∂A∗ ∂A∗ Im A δ(n, N)Im A dxdy + βA δ0 dxdydη ∂x ∂x −Ly 0 0 −Ly 0 Ly Lx ∂A∗0 Im A0 dxdy, = ∂x −Ly 0

Ly

9

Lx

(38)

where a symbol Im denotes an imaginary part of complex function. Let us represent the function A using Euler’s formula: A = A¯ · eiS , ¯ S are real ones. Then, we obtain the following integral relation: here introduced functions A, Ly Lx t Ly Lx Ly Lx ∂S0 2 ∂S ¯A2 ∂S dxdy + βA δ0 ¯ A¯ 20 δ(n, N)A dxdydη = dxdy. ∂x ∂x ∂x −Ly 0 −Ly 0 0 −Ly 0

(39)

(40)

Thus, if initial distribution of the pulse phase S0 does not contain the odd powers from x − Lxc 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: Ly Lx ∂S A¯ 2 dxdy = 0. (41) ∂x −Ly 0 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: Ly Lx t Ly Lx ∂A∗ ∂A∗ Im A δ(n, N)Im A dxdy + βA δ0 dxdydη ∂y ∂y −Ly 0 0 −Ly 0 Ly Lx ∂A∗ Im A0 0 dxdy (42) = ∂y −Ly 0 or in the representation (39) for the complex amplitude: Ly Lx t Ly Lx Ly Lx ∂S ∂S ∂S0 A¯ 2 dxdy + βA δ0 A¯ 20 δ(n, N)A¯ 2 dxdydη = dxdy. ∂y ∂y ∂y −Ly 0 −Ly 0 0 −Ly 0 For initial condition (9) this equality is re-written as: t Ly Lx Ly Lx ∂S ∂S A¯ 2 dxdy + βA δ0 δ(n, N)A¯ 2 dxdydη ∂y ∂y −Ly 0 0 −Ly 0 Ly Lx A¯ 20 dxdy = const. = −2π χ −Ly

(43)

(44)

0

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).

10


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: Ly Lx 2 Ly Lx 2 ∂ A ∂A∗ ∂ A ∂A∗ ∂ 2 A∗ ∂A ∂ 2 A∗ ∂A iDAx + dx dy + iD + dx dy A y ∂x2 ∂t ∂x2 ∂t ∂y2 ∂t ∂y2 ∂t −Ly 0 −Ly 0 Ly Lx ∂A ∂A∗ + 0.5βA δ0 δ(n, N) A − A∗ dx dy = 0. (45) ∂t ∂t −Ly 0 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: ∂ −i ∂t

Ly Lx −Ly 0

Ly Lx × −Ly 0

or

Ly −Ly

Lx 0

2 2

∂A ∂A DAx + DAy dxdy + 0.5βA δ0 ∂x ∂y

∂A∗ ∗ ∂A δ(n, N) A −A dxdy = 0 ∂t ∂t

(46)

2 2

t Ly Lx ∂A ∂A ∂A∗ δ(n, N)Im A dxdydη, DAx + DAy dxdy = I30 +βA δ0 ∂x ∂y ∂η 0 −Ly 0 (47)

where

I30 =

Ly

−Ly

Lx 0

∂A0 2 ∂A0 2 + DA dxdy. DAx y ∂x ∂y

Using the representation (39), we write the equality (47) in the form: ¯ 2 ¯ 2 2 2

Ly Lx ∂A ∂A ∂S ∂S + DAy + A¯ 2 DAx + DA y DAx dxdy ∂x ∂y ∂x ∂y −Ly 0 t Ly Lx ∂S δ(n, N)A¯ 2 dxdydη. = I30 − βA δ0 ∂η 0 −Ly 0

(48)

¯ S from Equation (4). Finally, To simplify this equality, we need to write the equations with respect to A, they are written in the form: 2

2

2A 2A ¯ ¯ ∂S ∂S ∂S ∂ ∂ A¯ − A¯ − A¯ + D Ax + DA y = 0, (49) ∂t ∂x2 ∂x ∂y2 ∂y ∂ ¯ 2 ∂S ∂ ¯ 2 ∂S ∂ A¯ 2 − 2DAx A − 2DAy A + βA δ0 δ(n, N)A¯ 2 = 0. ∂t ∂x ∂x ∂y ∂y

(50)


11

Using (49), one can transform equality (48) to the form:

Ly

−Ly

Lx

DAx

0

= I30 − βA δ0

t

2

∂ A¯ ∂x Ly

−Ly

0

∂ 2 A¯ + A¯ 2 ∂x

Lx 0

+ DA y

δ(n, N)A¯ 2

∂ A¯ ∂y

2

∂ 2 A¯ + A¯ 2 ∂y

∂S dx dy + A¯ 2 ∂t

∂S dx dy dη. ∂η

(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:

Ly

−Ly

or ∂ ∂t

Lx

A2

0

Ly −Ly

Lx

t Ly Lx ∂S ∂S δ(n, N)A¯ 2 dx dy dη dx dy = I30 − βA δ0 ∂t ∂η 0 −Ly 0 A2

0

Ly Lx ∂S ∂S δ(n, N)A¯ 2 dx dy = 0. dx dy + βA δ0 ∂t ∂t −Ly 0

(52)

(53)

This integral relation can be used also to control computer simulation results. If we introduce the function t ∂S δ(n, N)A¯ 2 dη, F= (54) ∂η 0 then one can rewrite Equation (52) in the form

Ly

−Ly

Lx 0

1 ∂ (F · eβA δ0 t )dx dy = I30 eβA δ0 t . δ(n, N) ∂t

(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: Ly Lx ∂S 2 δ(n, N)A¯ dx dy ≤ I30 δmax e−βA δ0 t . −Ly 0 ∂t Moreover, from (56) it is as follows: Ly Lx ∂S δ max 2 A¯ I e−βA δ0 t . dx dy ≤ −Ly 0 δmin 30 ∂t

(56)

(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. [7,16,20,24,27]). We suppose that the functions in the right part of the equation are known.

12


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:

Ly 0

0

Lx

(ϕϕ)dx dy = −γ

Ly

0

Lx 0

((N − n)ϕ)dx dy

(58)

or, taking into account the BCs with respect to the electric field potential:

Ly 0

0

Lx

Ly

2

(∇ϕ) dx dy = γ

0

0

Lx

((N − n)ϕ)dx dy,

∇ϕ =

∂ϕ ∂x

2 +

∂ϕ ∂y

2 .

(59)

Consequently, using the Cauchy–Bunkyakovsky inequality and ε inequality we write an assessment: ||∇ϕ||2L2

≤ γ |((N − n), ϕ)| ≤ γ

ε||n − N||2L2

1 2 + ||ϕ||L2 , ε > 0, 4ε

(60)

where the L2 norm from the electric field potential, or charged particles concentrations is defined in the domain, occupied by a semiconductor, as ||ϕ||2L2 = (ϕ, ϕ) =

0

Ly

0

Lx

ϕ 2 dx dy.

(61)

Next, multiplying Equation (3) on the function (N–n) and then take an integral over the domain occupied by a semiconductor, we obtain: ((N − n), ϕ) = −γ ||N − n||2L2

(62)

(∇(N − n), ∇ϕ) = γ ||N − n||2L2 .

(63)

or

Further, take modulus from the left part of Equation (63) and using the Cauchy–Bunkyakovsky inequality and ε inequality, we write an inequality: γ ||N − n||2L2 ≤ ε1 ||∇(N − n)||2L2 +

1 ||∇ϕ||2L2 , ε1 > 0. 4ε1

(64)

From (60) and (64) one can write two assessments: ε γ ≤ ε1 ε||∇(N − n)||2L2 + ||ϕ||2L2 , ||∇ϕ||2L2 1 − 4ε1 4ε ||N − n||2L2 (1 −

ε ε1 1 ) ≤ ||∇(N − n)||2L2 + ||∇ϕ||2L2 . 4ε1 γ 4εε1

(65)

(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.


13

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 [5]. The similar operators can meet for description of the gas mixture dynamics and they are similar to the Fokker–Planck equation [12]. In particular, in [5] 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: ∂n ∂ ∂ μx ϕ ∂ −μx ϕ μy ϕ ∂ −μy ϕ ) + Dy ) + G − R. = Dx e (ne e (ne ∂t ∂x ∂x ∂y ∂y

(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 [7,25]. 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 [7]. 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 [7]. On the other hand, the existence and uniqueness of the solution for nonlinear heat conduction equations, including the nonlinear BCs, were proved in [7]. 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 [5,7,25] 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 ¯ = {0 ≤ x ≤ Lx } × {0 ≤ y ≤ Ly } × {0 ≤ t ≤ Lt } G

14


Figure 4. Template of the mesh function definition.

the uniform time and space grids: ωx = {xk = khx , i = 0, Px , hx = Lx /Px }, ωy = {yj = jhy , j = −Py , Py , hy = Ly /Py },

ωy = {yj = jhy , j = 0, Py , hy = Ly /Py },

(68)

ω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 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: m m nm kj = n(xk , yj , tm ), Nkj = N(xk , yj , tm ), ϕkj = ϕ(xk , yj , tm ), m m Gm kj = G(xk , yj , tm ), Rkj = R(xk , yj , tm ), δkj = δ(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: Am kj = A(xk , yj , tm ) (Figure 4). Grid parameter δ0h (maximal value of the semiconductor laser energy absorption) for the parameter δ0 is defined on the grid 2 by the rule: ⎧ j = −Py , −1 0, ⎪ ⎪ ⎨ . δ0h = 0.5δ0 , j = 0 ⎪ ⎪ ⎩ δ0 , j = 1, Py For brevity, below we use the following index-free notations: f = fh = fkjm = f (xk , yj , tm ),

m fk±1 = fk±1j ,

m fj±1 = fkj±1 ,

fˆ = fkjm+1 ,

(70)


15

where f is one of the defined mesh functions and omit index h. We also define some useful notations using (14): fk±0.5 = 0.5(fk + fk+1 ),

0.5 0.5

0.5

G = q0 |A|2 δ , 0.5

0.5

f = 0.5(f m + f m+1 ),

fj±0.5 = 0.5(fj + fj+1 ),

0.5

ˆ + R(n, N)), R = 0.5(R(ˆn, N)

(71)

0.5

ˆ + δ(n, N)), δ = 0.5(δ(ˆn, N)

2

ˆ + |A|2 ) |A|2 = 0.5(|A| 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)ϕ =

μy μx (nk+0.5 ϕx − nk−0.5 ϕx¯ ), Zy¯ y (n)ϕ = (nj+0.5 ϕy − nj−0.5 ϕy¯ ), hx hy

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): 0.5 0.5 Nˆ − N = G − R, τ

nˆ − n = κx τ

0.5 n x¯ x

0.5

k = 0, Px , j = 0, Py ,

(72)

0.5 − 0.5(Zx¯ x (ˆn)ϕˆ + Zx¯ x (n)ϕ) + κy n y¯ y − 0.5(Zy¯ y (ˆn)ϕˆ + Zy¯ y (n)ϕ)

0.5

+ G − R,

k = 1, Px − 1,

j = 1, Py − 1,

ˆ k = 1, Px − 1, j = 1, Py − 1, ϕˆ = γ (ˆn − N), 0.5 0.5 βA δ0h 0.5 0.5 Aˆ − A + iDAx A x¯ x + iDAy A y¯ y + δ A = 0, τ 2

k = 1, Px − 1, j = −Py + 1, Py − 1.

(73)

(74)

(75)

0.5

Let us note, that δ kj = 0 if j = −Py + 1, −1. The BCs (5) for the charged particles concentrations and electric field potential are approximated with the first order: ϕˆx |k=0j = ϕˆx¯ |k=Px j = 0, nˆ x |k=0j = nˆ x¯ |k=Px j = 0, j = 0, Py ,

(76)

16


ϕˆy |kj=0 = ϕˆy¯ |kj=Py = 0, nˆ y |kj=0 = nˆ y¯ |kj=Py = 0, k = 0, Px , ˆ x=Px j = 0, j = −Py , Py , A| ˆ kj=−Py = A| ˆ kj=Py = 0, k = 0, Px . ˆ x=0j = A| A| The corresponding initial conditions for electron and donor concentrations are written as: n0kj = Nkj0 = n0, ϕkj0 = 0, k = 0, Px , j = 0, Py .

(77)

The complex amplitude initial condition is defined in the following way: A0kj = e

−

(xk −Lxc )2 −(yj −Lyc )10 −2iπχ(yj −Lyc ) a2

, 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: Py −1 Px −1

Q(tm ) =

m hy hx (nm kj − Nkj ).

(79)

j=1 k=1

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: Py −1 Px −1

Py −1 Px −1

(nt − Nt )hy hx =

j=1 k=1

ˆ − (n − N)) ((ˆn − N)

j=1 k=1 Py −1 Px −1

=

0.5

hy hx τ

0.5

hy hx (κx n x¯ x + κy n y¯ y )

j=1 k=1 Py −1 Px −1

− 0.5

hy hx (κx (Zx¯ x (ˆn)ϕˆ + Zx¯ x (n)ϕ) − κy (Zy¯ y (ˆn)ϕˆ + Zy¯ y (n)ϕ)).

j=1 k=1

(80)


17

Let us calculate the sums entering the right part of relation (80). For example, we see that two equalities are valid: P x −1

hx nˆ x¯ x

k=1 P x −1 k=1

hx Zx¯ x (ˆn)ϕˆ = μx

Px −1 1 = (ˆnk−1 − 2ˆn + nˆ k+1 ) =(ˆnx¯ ,Px j − nˆ x,0j ), hx k=1

P x −1 i=1

ϕˆk+1 − ϕˆk ϕˆk − ϕˆk−1 nˆ k+0.5 − nˆ k−05 hx hx

= μx (ˆnPx −0.5j ϕˆx¯ ,Px j − nˆ 0.5j ϕˆx,0j ).

Obviously, that at the BCs approximation (76), the following equalities are valid: Py −1 Px −1

Py −1

hy hx nˆ x¯ x =

j=1 k=1

hy (ˆnx¯ ,Px j − nˆ x,0j ) = 0,

(81)

j=1

Py −1 Px −1

Py −1

hy hx Zx¯ x (ˆn)ϕˆ = μx

j=1 k=1

hy (ˆnPx −0.5j ϕˆx¯ ,Px j − nˆ 0.5j ϕˆx,0j ) = 0.

j=1

Similar equalities are valid also for other sums entering to the right part of expression (80): Py −1 Px −1

hy hx nˆ y¯ y =

j=1 k=1

hy hx Zy¯ y (ˆn)ϕˆ = μy

j=1 k=1

hx (ˆny¯ ,kPy − nˆ x,k0 ) = 0,

(82)

k=1

Py −1 Px −1

P x −1

P x −1

hx (ˆnPy −k0.5 ϕˆy¯ ,kPy − nˆ k0.5 ϕˆy,k0 ) = 0.

k=1

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 Py −1 Px −1

ˆ = hy hx (ˆn − N)

j=1 k=1

Py −1 Px −1

hy hx (n − N) = 0.

(83)

j=1 k=1

If we take into account the initial conditions (77) then the equality (83) can also be written in the form ˆ L1 , ||ˆn||L1 = ||N||

(84)

where mesh norm L1 for the positive mesh function is defined as Py −1 Px −1

||u||L1 =

hy hx ukj .

(85)

j=1 k=1

Applying similar way as we provide for equalities (80) one can obtain the following relation: ˆ L1 = Lx Ly = const. ||ˆn||L1 + ||1 − N||

(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 .

18


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:

0.5 0.5 nˆ n ˆ 2 −ψ+ξ ψ n 2 −ψ+ξ ψn + + = N 1 − 0.5τ q0 |A| e Nˆ 1 + 0.5τ q0 |A| e τR τR

0.5 n2 + τ 0.5q0 |A|2 e−ψ (eξ ψ nˆ + eξ ψn ) + 0 . (87) τR 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 τ≤

0.5

2 0.5 q0 |A|2max e−ψ+ξ ||n||C

, +

||n||C τR

0.5

|A|2max = max |A|2 k = 1, Px − 1, j = 1, Py − 1 k,j

(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: 2τR ˆ if nˆ < 0. (89) τ≤ , minˆn = min n, ˆ |min n| j,k ˆ Therefore, independent from a sign of the mesh function nˆ the value of Nwill 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

0.5 2 −ψ+ξ ψ nˆ

ˆ 1 + 0.5τ q0 |A| e (1 − N) + 0.5τ

0.5 n nˆ Nˆ 2 −ψ+ξ ψn + (1 − N) 1 − 0.5τ q0 |A| e + = 0.5τ τR τR

n − 2n20 . τ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 > n20 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.


19

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 (u, v) =

P x −1

Py −1

Py −1

hx hy uv =

k=1 j=−Py +1

P x −1

∗

j=−Py +1 k=1

∗ hy hx ukj vkj .

(91)

Therefore, the mesh norm L2 is defined as ||u||2L2 = (u, u) =

Py −1

P x −1

hx hy uu∗ =

k=1 j=−Py +1

P x −1

Py −1

hx hy |u|2 .

(92)

k=1 j=−Py +1

0.5

0.5

Multiplying Equation (75) on A∗kj and the equation, conjugated to Equation (75), – on Akj , 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: Py −1

P x −1

hy hx

j=−Py +1 k=1

2

ˆ − |A|2 0.5 0.5 |A| + βA δ0h δ |A|2 τ

=0

(93)

from which one can write Py −1

ˆ 2L ||A|| 2

=

||A||2L2

− τ βA δoh

P x −1

0.5 0.5

hy hx δ |A|2 .

(94)

j=−Py +1 k=1

We have to stress that for writing expression (93) we use the difference Green formulas for the mesh functions (see, e.g. [24]): (ux¯ x , v) = −(ux¯ , vx¯ ) + uvx¯ |Px − vux |0 .

(95)

As we assumed above the mesh function δ(n, N) is bounded from below by the constant δmin = ˆ Therefore, the next inequality follows from (94): ˆ N). min δ(n, N, n, j,k

ˆ 2L 1 + τ βA δoh δmin ≤ ||A||2L 1 − τ βA δoh δmin . ||A|| 2 2 2 2

(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: ˆ 2L ≤ qm ||A0 ||2L ≤ ||A0 ||2L , ||A|| 2 2 2

−1

τ τ q = 1 − βA δ0h δmin 1 + βA δ0h δmin , 2 2

(97)

if the following condition τ < τ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 .

20


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 [6] 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: s+1

s,s+1

s,s+1

0.5 0.5 Nˆ −N = G − R , τ

k = 0, Px , j = 0, Py ,

(99)

s,s+1 s,s+1 s+1 s s s 0.5 0.5 nˆ −n 0.5 0.5 ˆ ϕˆ +Zx¯ x (n)ϕ = κx n x¯ x + κy n y¯ y + G − R −0.5 μx κx Zx¯ x (n) τ s s ˆ ϕˆ +Zy¯ y (n)ϕ , k = 1, Px − 1, j = 1, Py − 1, +μy κy Zy¯ y (n)

s+1

s,s+1 0.5

s

0.5

G = q0 | A |

s,s+1 0.5

s,s+1 2 0.5

ˆ Nˆ δ = 0.5 δ n,

s+1

s+1

+ δ(n, N) ,

nˆ − Nˆ

s+1

s

k=0

k=Px

+ R(n, N) ,

s

0.5 2

s

(101)

ˆ 2 +|A|2 |A| = 0.5 |A|

s+1

s+1

ϕˆ = γ

s s+1

ˆ Nˆ R = 0.5 R n,

δ ,

s s+1

s,s+1 0.5

(100)

k = 1, Px − 1, j = 1, Py − 1,

,

(102)

s,s+1 s

0.5 0.5 βA δoh 0.5 0.5 Aˆ −A + iDAx Ax¯ x +iDAy Ay¯ y + δ A = 0, k = 1, Px − 1, j = −Py + 1, Py − 1, (103) τ 2 s+1 s+1 s+1 s+1 ϕˆ x = ϕˆ x¯ = 0, j = 0, Py , ϕˆ y = ϕˆ y¯ = 0, k = 0, Px , (104)

nˆ x

s+1

k=0

j=0

s+1 = nˆ x¯

= 0, j = 0, Py ,

k=Px

j=Py

s+1

s+1

0j

Px j

Aˆ = Aˆ = 0, j = −Py , Py .

The second stage is: s+2

Nˆ −N = τ

s+2

s+1

s+2

s+1,s+2 0.5

nˆ −n 0.5 0.5 = κx n + κy n + x¯ x y¯ y τ

G

−

s+1,s+2 0.5

G

s+1,s+2 0.5

−

R

,

s+1,s+2 0.5

R

k = 0, Nx , j = 0, Ny , s+1 s+1 −0.5 μx κx Zx¯ x ( nˆ ) ϕˆ +Zx¯ x (n)ϕ

(105)


s+1 s+1

+μy κy Zy¯ y ( nˆ ) ϕˆ +Zy¯ y (n)ϕ

s+1,s+2 0.5

0.5

= q0 | A |

G

s+1,s+2 0.5

δ

s+1,s+2 0.5

s+1

2

δ

s+1 s+2

s+1

k=0

s+2

k=Px

nˆ y

s+2

j=0

s+1 s+2

nˆ , Nˆ

= 0.5 R

⎛

s+1

0.5 2

(106)

+ R(n, N) ,

s+1 2

(107)

⎞

ˆ +|A|2 ⎠ . |A| = 0.5 ⎝|A|

nˆ − Nˆ

s+2 = ϕˆ x¯

+ δ(n, N) ,

0.5 0.5 β A δ oh Aˆ −A + iDAx Ax¯ x +iDAy Ay¯ y + τ 2

s+2 ϕˆ x

k = 1, Px − 1, j = 1, Py − 1,

,

s+2

s+2

ϕˆ = γ

s+2

R

nˆ , Nˆ

s+2

s+1,s+2 0.5

,

= 0.5 δ

21

,

k = 1, Px − 1, j = 1, Py − 1,

s+1,s+2 s+1 0.5 0.5

δ

A = 0,

k = 1, Px − 1, j = −Py + 1, Py − 1, (109)

s+2 = 0, j = 0, Py , ϕˆ y

j=0

s+2 = nˆ y¯

(108)

s+2 = ϕˆ y¯

= 0, k = 0, Px ,

(110)

j=Py

s+2

s+2

k−Px

kPy

= 0, Aˆ = Aˆ = 0, k = 0, Px .

j=Py

The iterations are stopped if the following inequalities s+2 s s nˆ − nˆ ≤ ε1 nˆ + ε2 ,

s s+2 s ˆ N − Nˆ ≤ ε1 Nˆ + ε2 ,

s+2 s s ϕˆ − ϕˆ ≤ ε1 ϕˆ + ε2 ,

(111)

s+2 s s ˆ ˆ A − A ≤ ε1 Aˆ + ε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: s=0

s=0

s=0

s=0

nˆ = n, Nˆ = N, ϕˆ = ϕ, Aˆ = 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.

22


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 [21,27], then in the present paper we use the stabilization method, which is written in the following manner for our problem: F p+1 − F p p+1 p = Fx¯ x + Fy¯ y − γ τ¯ p+1

Fx

nˆ − Nˆ

p+1

p+2 j=0

s+1

s+1

|k=0 = Fx¯

F p+2 − F p+1 p+1 p+2 = Fx¯ x + Fy¯ y − γ τ¯ Fy

,

|k=Px = 0, j = 0, Py ,

(113)

(114)

s+1

s+1

nˆ − Nˆ

p+2

= Fy¯

k = 1, Px − 1, j = 1, Py − 1,

j=Py

,

k = 1, Px − 1, j = 1, Py − 1

= 0, k = 0, Px ,

(115)

(116)

¯ = ωx × ωy , τ¯ is an iteration where F is an auxiliary mesh function defined on the spatial grid parameter, which is not equal to time grid step τ , p is an iteration number. The initial approach for this iteration process is s

F 0 = ϕˆ .

(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

s+1 s+1 p+2 p+1 Fx¯ x + Fy¯ y − γ nˆ − Nˆ < ε3 , ε3 > 0,

k = 1, Px − 1, j = 1, Py − 1

(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: s+1

ϕˆ = F p+2 .

(119)

Theorem 4.2: The iteration method (113)–(117) converges if the following inequality τ¯ ≤ τ0 is valid, where τ0 = min{h2x , h2y }.


23

s+1

4.2. Proving the positivity property for 0 < Nˆ < 1 From Equation (99), one can write the following expression for the mesh function Nˆ at s + 1iteration s+1

( Nˆ ): ⎛

⎛

s+1

s s 0.5 2 −ψ+ξ ψ nˆ

⎜ ⎜ Nˆ ⎝1 + 0.5τ ⎝q0 |A| e

s

+

⎞⎞

⎛

⎛

s 0.5 2 −ψ+ξ ψn

nˆ ⎟⎟ ⎜ ⎜ ⎠⎠ = N ⎝1 − 0.5τ ⎝q0 |A| e τR s 0.5 2 −ψ

+ 0.5τ q0 |A| e

s

e

ξ ψ nˆ

⎞⎞ +

ξ ψn

n ⎟⎟ ⎠⎠ τR

+e

+

n20 τ . (120) τR

Supposing that the mesh function N at previous time layer is positive and bounded: 0 < N ≤ QN , QN = const > 0,

(121)

s

the mesh functions n, nˆ are positive and bounded also: 0 < n ≤ Qn , Qn = const > 0,

s

0 < nˆ ≤ Qn + ε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 s

2

ˆ C ≤ QA + εA , ||A||C 2 ≤ QA , ||A||

(123)

||A||C = max |A|, j = −Py + 1, Py − 1, k = 1, Px − 1,

(124)

where j,k

s+1

then from Equation (120) the positiveness of the mesh function Nˆ on upper time layer occurs if the following inequality is valid τ
n20 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, s+1

we conclude about a validity of the mesh function property Nˆ < 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)). s+1

4.3. Uniform boundedness of Nˆ

s+1

From the previous paragraph, it follows that the mesh function Nˆ 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 s

Nˆ ≤ QN + εN .

(128)

Let us show that the next inequality s+1

Nˆ ≤ QN + εN

is valid. With this aim, Equation (120) is re-written in the form ⎛ ⎛ ⎞⎞ s ⎜ ⎜ N ⎝1 − 0.5τ ⎝q0 s+1

Nˆ =

0.5 |A|2 e−ψ+ξ ψn

+

n ⎟⎟ τR ⎠⎠ + 0.5τ q0

⎛ ⎜ 1 + 0.5τ ⎝q0

s 0.5 s |A|2 e−ψ+ξ ψ nˆ

(129)

s s 0.5 |A|2 e−ψ (eξ ψ nˆ

+ eξ ψn ) +

n20 τ τR

⎞

. (130)

s

+

nˆ ⎟ τR ⎠

Thus, there is an inequality 2 s+1 n0 + 0.5q0 Im + 0.5εA e−ψ(1−ξ Qn ) (1 + eψξ εn ) ≤ QN + εN , Nˆ ≤ QN + τ τR

(131)

where instead (122) we introduce a new notation ||n||C ≤ Qn ,

s

||ˆn|| ≤ Qn + εn ,

(132)

C

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 τ≤

n20 τR

εN 2εN e−ψ(1−ξ Qn ) ≤ q0 Im (1 + eψξ εn ) + 0.5q0 Im + 0.5εA e−ψ(1−ξ Qn ) (1 + eψξ εn ) s+1

at a validity of inequality (125), the mesh function Nˆ is uniformly bounded.

(133)


25

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: s+1

s

0.5τ Nˆ − Nˆ = − τR ⎛

s+1 s

s s−1

Nˆ nˆ − Nˆ nˆ

s s 0.5 2 ψξ nˆ

⎜ − ⎝|A| e

⎛

−ψ

− 0.5τ q0 e

s s+1 0.5

s

s

s−1 0.5

⎜ ˆ 2 ψξ nˆ − Nˆ |A|2 eψξ ⎝ N |A| e

s−1

nˆ

⎞⎞

s−1 s−1 0.5 2 ψξ nˆ ⎟⎟

− |A| e

⎠⎠ .

(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 s+1 s

s s

s s−1

s+1

Nˆ nˆ ± Nˆ nˆ − Nˆ nˆ =

s

s

s

Nˆ − Nˆ nˆ + Nˆ

s+1 s nˆ − nˆ .

(135)

In a similar way, we transform the second term: s s+1 0.5

s

s

s−1 s−1 0.5 2 ψξ nˆ

Nˆ |A| eψξ nˆ − Nˆ |A| e

s+1

s

s−1 s 0.5 2 ψξ nˆ

Nˆ − Nˆ |A| e

=

s

s−1 s 0.5 2 ψξ nˆ

+ Nˆ |A|

e

⎛ s

s 0.5 2

s

⎞

s−1 s 0.5 2 ψξ nˆ ⎟

⎜ + Nˆ eψξ nˆ ⎝|A| − |A| e

s−1

⎠

ψξ nˆ

−e

.

(136)

The last two terms can be transformed as: s s 0.5 2 ψξ nˆ

|A| e

s−1 s−1 0.5 2 ψξ nˆ

± |A| e

s−1 s−1 0.5 2 ψξ nˆ

− |A| e

⎛ s

=e

ψξ nˆ ⎜

s 0.5 2

⎞

s−1 s 0.5 2 ψξ nˆ

s−1 0.5 2⎟

⎝|A| − |A| ⎠ + |A|

e

s−1

ψξ nˆ

−e

.

(137)

Thus, Equation (134) could be written as

s+1

s

Nˆ − Nˆ

1 + 0.5τ e ⎛

s

s 0.5 2

−ψ

1 + τR

⎞⎞

s−1 0.5 2 ⎟⎟

⎛ s−1

0.5τ ⎜ +eξ ψ nˆ ⎝|A| − |A| ⎠⎠ − τR

= 0.5τ q0 e s s−1 nˆ − nˆ .

−ψ

s

0.5

s

2 ξ ψ nˆ ˆ ⎜ ((1 − N)) − eξ ψ ⎝|A| e

s−1

nˆ

(138)

26


If we use Taylor’s series for the exponential term (or a property of Lipchitz continuity) then s

eξ ψ nˆ − eξ ψ

s−1

nˆ

s

s−1

= Mn (nˆ − nˆ ),

(139)

where Mn is a constant of the series. Finally, we obtain the inequality s ⎞ s−1 s s−1 0.5 0.5 0.5τ q0 ⎟ ⎜ ˆ ≤

⎝nˆ − nˆ Im Mn + eξ ψ(Qn +εn ) |A|2 − |A|2 ⎠ | Nˆ − N| 1 1 + 0.5τ e−ψ + τR s ⎞ ⎛ s−1 s s−1 0.5 0.5 ⎟ ⎜ ≤ 0.5τ q0 e−ψ ⎝nˆ − nˆ Im Mn + eξ ψ(Qn +εn ) |A|2 − |A|2 ⎠ . ⎛

s+1

e−ψ

s

(140)

For the next neighbour iterations s + 2 and s + 1, we obtain the inequality which is similar to inequality (140): ⎞ ⎛ s+1 s s+2 s+1 0.5 0.5 s+1 s ˆ ⎟ ⎜ ˆ m Mn + eξ ψ(Qn +εn ) ||A|2 − |A|2 |⎠ . N − Nˆ ≤ 0.5τ q0 e−ψ ⎝| nˆ − n|I

(141)

Inequalities (140) and (141) will be used below for writing the conditions for the iteration process convergence. s+1

ˆ 2 4.5. Uniform boundedness of ||A|| L2 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 ˆ 2L ≤ CA ||A0 ||2L , ||A|| 2 2

(142)

if the grid step on time coordinate satisfies inequality (97). Because all norms in finite-dimensional spaces are equivalent [1] then we proposed the validity of boundedness for the mesh function A on s

low layer in time coordinate and for the mesh function Aˆ on upper layer in time coordinate at s iteration (see (123)). Below we show that the assessment s+1

ˆ C ≤ QA + εA ||A||

(143)

s+1

is valid for the mesh function Aˆ . 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:

s s s 0.5 τ β δ A 0 ˆ ∗ −1 ˆ −ψ(1−ξ n) Aˆ = B−1 + (1 − N)e−ψ(1−ξ n) A , B−1 (1 − N)e x Bx A − iτ DAy Bx y¯ y A − x 4 (144) where we introduce the following operators:

τ τ Bx = E + i DAx x¯ x , B∗x = E − i DAx x¯ x . (145) 2 2 s+1

s 0.5


27

Multiplying Equation (144) by the operator B−1 x from the left and taking into account the assessments for the operators [14] c1 ∗ ||B−1 ||y¯ y ||C ≤ 2 , c1 = const > 0 (146) ||B−1 x Bx ||C ≤ 1, x ||C < 1, hy s

s

s

ˆ N, Nˆ in this norm as well as || Nˆ ||c < 1 we obtain: and boundedness of the mesh functions n, n,

s+1 c1 βA δ0 −ψ(1−ξ Qn ) ψξ εn ˆ ||A||C ≤ ||A||C + τ (2QA + εA ) DA + (1 + e ) . (147) e 2h2y y 4 Therefore, in Equation (143) will be valid if the mesh step on time coordinate satisfies to the inequality: τ≤

4h2y εA (2QA + εA )(c1 DAy + h2y βA δ0 e−ψ(1−ξ Qn ) (1 + eψξ εn ))

.

(148)

s+2

Carrying out the similar algebra for the mesh function Aˆ we obtain from (109) the following equation:

s+1 s+1 s+2 s+1 s+1 0.5 0.5 τ β δ A 0 ˆ −1 ∗ −1 −1 −ψ(1−ξ −ψ(1−ξ n) n ) Aˆ = By By A − iτ DAx By x¯ x A − + (1 − N)e By (1 − Nˆ )e A, 4 (149) where we introduce the following operators:

τ τ (150) By = E + i DAy y¯ y , B∗y = E − i DAy y¯ y . 2 2 Multiplying Equation (149) by the operator B−1 y from the left and taking into account the assessments for the operators [14] c2 ∗ ||B−1 ||x¯ x ||C ≤ 2 , c2 = const > 0 (151) ||B−1 y By ||C ≤ 1, y ||C < 1, hx and supposing that this assessment s+1

||ˆn|| ≤ Qn + εn

(152)

C

is valid (we discuss this statement below), we obtain s+2 ˆ C ≤ ||A||C + τ (2QA + εA ) c2 DAx + βA δ0 e−ψ(1−ξ Qn ) (1 + eψξ εn ) . ||A|| 2h2x 4

(153)

Therefore, the inequality s+2

ˆ C ≤ QA + εA ||A||

(154)

will be valid if the mesh step on time coordinate satisfies to the inequality: τ≤

4h2x εA . 2 (2QA + εA )(c2 DAx + hx βA δ0 e−ψ(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.

28


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

Using the difference equation (103) for the mesh function Aˆ and the difference equation, which is s

similar to (109), for the mesh function Aˆ we obtain the following difference equation with respect to difference of these mesh functions s+1

s

s+1

s

s

s−1

s

s

s−1 s−1

βA δ0 0.5 0.5 0.5 0.5 Aˆ − Aˆ +i0.5τ DAx (Aˆ x¯ x − Aˆ x¯ x + Aˆ x¯ x − Aˆ x¯ x ) + τ (δ A − δ A)=0 2

(156)

or s+1

s

s

s s 0.5 0.5

s−1

s−1 s−1 0.5 0.5

ˆ + i0.5τ DAx (Aˆ x¯ x − Aˆ x¯ x ) + 0.5τ βA δ0 ( δ A − δ A ) = 0. Bx ( Aˆ − A)

(157)

To write the final difference equation we need to transform the last difference in Equation (157), which can be written as

s s−1 s−1 s s s−1 s−1 s−1 ˆ ψξ ψξ n n (1 − Nˆ ) + (1 − N)e (Aˆ − Aˆ ) e + Aˆ (1 − Nˆ )(eψξ nˆ − eψξ nˆ ) s

s−1

s

ψξ nˆ

− (Nˆ − Nˆ )e

s

Aˆ +A

(158)

or using Equation (139), expression (158) transforms to

s s−1 s−1 s s−1 s−1 s s−1 (Aˆ − Aˆ ) eψξ nˆ (1 − Nˆ ) + (1 − N)eψξ n + (A + Aˆ )(1 − Nˆ )Mn (nˆ − nˆ ) s

s−1

s

− (Nˆ − Nˆ )eψξ nˆ (Aˆ + A)

(159)

Thus, Equation (157) is written in the form:

s+1 s s s−1 s s−1 δ β A 0 ˆ ψξ ψξ n n ˆ + 0.5τ iDAx x¯ x + Bx ( Aˆ − A) (1 − Nˆ ) + (1 − N)e e (Aˆ − Aˆ ) 4

s−1 s s s s−1 s−1 s s−1 τ β A δ0 ˆ ψξ n + (Nˆ − Nˆ ) = 0. (A + Aˆ )(1 − Nˆ )Mn (nˆ − nˆ ) − (Aˆ +A)e 8

(160)

In a similar way, one can transform the second difference between the mesh function Aˆ at s + 2 and s + 1 iterations: s+2

s+1

s+1

s

s+1 s+1

s

s

βA δ0 0.5 0.5 0.5 0.5 By ( Aˆ − Aˆ ) + i0.5τ DAy y¯ y (Aˆ y¯ y − Aˆ y¯ y ) + τ ( δ A − δ A ) = 0. 2 After applying some algebra, we write the equation, which is similar to Equation (160),

s+2 s+1 s+1 s s+1 s δ β A 0 ˆ ψξ ψξ n n ˆ ˆ + (1 − N)e By ( Aˆ − Aˆ ) + 0.5τ iDAy y¯ y + (1 − N) e ( Aˆ − A) 4

s s s+1 s+1 s s s+1 s τ β A δ0 ˆ ψξ n ˆ ˆ n ( nˆ − nˆ ) − (Aˆ +A)e ˆ + − N)M ( Nˆ − N) (A + A)(1 = 0. 8

(161)

(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)


29

−1 and (162) by multiplying them by B−1 x , By from the left, correspondingly, we obtain s+1

( Aˆ

s

ˆ + 0.5τ B−1 − A) x

β A δ0 iDAx x¯ x + 4

s

e

ψξ nˆ

s−1

ψξ n

(1 − Nˆ ) + (1 − N)e

s

s−1

(Aˆ − Aˆ )

s−1 s s s s−1 s−1 s s−1 τ βA δ0 −1 ˆ ψξ n + (Nˆ − Nˆ ) = 0, Bx (A + Aˆ )(1 − Nˆ )Mn (nˆ − nˆ ) − (Aˆ +A)e 8 s+2

s+1

( Aˆ − Aˆ

) + 0.5τ B−1 y

β A δ0 iDAy y¯ y + 4

s+1

ψξ nˆ

e

s

ψξ n

ˆ + (1 − N)e (1 − N)

(163)

s+1

s

ˆ ( Aˆ − A)

s s s+1 s+1 s s s+1 s τ βA δ0 −1 ˆ ψξ n ˆ ˆ n ( nˆ − nˆ ) − (Aˆ +A)e ˆ + − N)M ( Nˆ − N) By = 0. (A + A)(1 8

(164)

Then, taking into account the assessments for the corresponding operators, we write: s s−1 s+1 s c δ D β Ax 2 A 0 ˆ ψξ Q ψξ ε n n ˆ ≤ 0.5τ + (1 + e ) (Aˆ − Aˆ ) e ( A − A) h2x 4 C C

s s−1 s s−1 τ β A δ0 ψξ(Q +ε ) n n + (2QA + εA ) Mn (Nˆ − Nˆ ) , (nˆ − nˆ ) + e 8 C

(165)

C

s+1 s s+2 s+1 DAy c1 βA δ0 ψξ Qn ˆ ψξ εn ˆ ˆ − A)|| ˆ C A + (1 + e ) ||( e ( A − A ) ≤ 0.5τ h2y 4 C

s+1 s s+1 s τ β A δ0 ψξ(Qn +εn ) ˆ ˆ ˆ +e + (2QA + εA ) Mn ( nˆ − n) ( N − N) . 8 C

(166)

C

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. s

ˆ C 4.7. Uniform boundedness of ||n|| 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: Py −1 Px −1

j=1 k=1

⎛

s+1

⎞

⎜ nˆ −n ⎟ hx hy ⎝ ⎠ = 0.5 τ

Py −1 Px −1

⎡ s+1

s

s 0.5

⎢ hx hy ⎣((1 − Nˆ )e−ψ+ξ ψ nˆ + (1 − N)e−ψ+ξ ψn ) |A|2

j=1 k=1

⎛ s+1 ⎞⎤ s 2 2 Nn − n0 ⎟⎥ ⎜ Nˆ nˆ −n0 −0.5 ⎝ + ⎠⎦ . τR τR s+1

(167)

Let us emphasize that a value of the mesh function Nˆ 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

30


s+1

boundedness of this mesh function. We take into account below that the mesh function Nˆ is positive s+1

and bounded. In this case, from (167) the inequality for nˆ follows:

s s+1 s τ 2 τ −ψ ξ ψ nˆ ξ ψn ||n||L1 + ||ˆn||L1 . ||ˆn|| ≤ ||n||L1 + 0.5τ e Im + n0 Lx Ly + 0.5 + e e L1 τR τR L1 L1 (168) s+1

s+1

At writing the ||ˆn|| we suppose that nˆ is positive. Some discussion about this statement will be L1

present in the end of Section 4. In general case, it is necessary to use modulus from this mesh function. s

Because in expression (122) we suppose the boundedness of the mesh functions n, nˆ 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 [1]. Therefore, we state: ¯ n, Q ¯ n = const > 0, ||n||L1 ≤ Q

s

¯ n + ε¯ n , ε¯ n > 0. || nˆ ||L1 ≤ Q

(169)

Let notice that other reason according to which the norm ||n||L1 is bounded arises from a validity of the invariant Q(tm ) = const (see (79)) and boundedness of the norm ||N||L1 . With accordance to s

s

Theorem 4.1 (see Section 4) and boundedness of the norm || Nˆ ||L1 , the ||ˆn|| is bounded also. L1

s+1

¯ n + ε¯ n , if the mesh step on time coordinate Thus, the norm ||ˆn|| is bounded by the same constant Q L1

satisfies the inequality: τ ≤ s eξ ψ nˆ

L1

+ eξ ψn L

2eψ ε¯ n

(2QA + εA ) +

1

2 2 τR n0 Lx Ly

+

1 τR

¯ n + ε¯ n 2Q

.

(170)

s+2

The similar assessment can be obtained with respect to the mesh function nˆ (106):

s+1 s+2 ξ ψn τ −ψ ξ ψ nˆ ||ˆn|| ≤ ||n||L1 + 0.5τ e e + e L1 (2QA + εA ) + n20 Lx Ly τ L1 R L1

s+1 τ + 0.5 ||n||L1 + nˆ L . 1 τR

(171)

From this inequality, one can conclude that s+2

¯ n + ε¯ n ||ˆn|| ≤ Q

(172)

L1

if the grid step on time coordinate satisfies the inequality: τ ≤ s+1 eξ ψ nˆ

+ eξ ψn

L1

L1

2eψ ε¯ n (2QA + εA ) +

. 2 2 τR n0 Lx Ly

+

1 ¯ τR (2Qn

(173)

+ ε¯ n )

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).


31

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 Py −1 Px −1

s+1

Py −1 Px −1

s

ˆ = hx hy ( nˆ − n)

j=1 k=1 Py −1 Px −1

s+1

s

ˆ =τ hx hy ( Nˆ − N)

j=1 k=1

s+2

Py −1 Px −1

s+1

hx hy ( nˆ − nˆ ) =

j=1 k=1

s,s+1 0.5

Py −1 Px −1

s,s+1 0.5

hx hy ( G − R ),

(174)

j=1 k=1

s+2

s+1

hx hy ( Nˆ − Nˆ ) = τ

j=1 k=1

Py −1 Px −1

s+1,s+2 0.5

hx hy ( G

−

s+1,s+2 0.5

R

), (175)

j=1 k=1

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: s+1

s

s

s+1

s−1

s

s

s−1

ˆ − (Nˆ − Nˆ )) ˆ = 0.5τ κx x¯ x (nˆ − nˆ ) + (( Nˆ − N) (1 − 0.5τ κx x¯ x )( nˆ − n) s s s−1 s−1 ˆ ϕˆ −Zx¯ x nˆ ϕˆ − 0.5τ μx κx Zx¯ x (n)

(176)

s s s−1 s−1 ϕˆ , +μy κy Zy¯ y nˆ ϕˆ −Zy¯ y nˆ and s+2

s+1

s+1

s

s+2

s+1

s+1

s

ˆ (1 − 0.5τ κy y¯ y )( nˆ − nˆ ) = 0.5τ κy y¯ y ( nˆ − nˆ ) + (( Nˆ − Nˆ ) − ( Nˆ − N)) s+1 s+1 s s ϕˆ −Zx¯ x nˆ ϕˆ − 0.5τ μx κx Zx¯ x nˆ

(177)

s+1 s+1 s s ϕˆ −Zy¯ y nˆ ϕˆ . Zy¯ y nˆ

+μy κy If we introduce the operators

Bnx = 1 − 0.5τ κx x¯ x , Bny = 1 − 0.5τ κy y¯ y ,

(178)

Equations (176), (177) can be re-written as s+1

s

s

s−1

s+1

s

s

s−1

ˆ ˆ ˆ ˆ ˆ − nˆ ) + B−1 ( nˆ − nˆ ) = 0.5τ κx B−1 nx x¯ x (n nx (( N − N) − (N − N )) s s s−1 s−1 ˆ ϕˆ −Zx¯ x nˆ ϕˆ − 0.5τ B−1 nx μx κx Zx¯ x n

(179)

s s s−1 s−1 Zy¯ y nˆ ϕˆ −Zy¯ y nˆ ϕˆ ,

+μy κy s+2

s+1

s+1

s

s+2

s+1

s+1

s

ˆ ˆ ˆ ˆ ˆ − n) ˆ + B−1 ( nˆ − nˆ ) = 0.5τ κy B−1 ny y¯ y ( n ny (( N − N ) − ( N − N)) s+1 s+1 s s nˆ ϕˆ −Zx¯ x nˆ ϕˆ − 0.5τ B−1 ny μx κx Zx¯ x s+1 s+1 s s ϕˆ −Zy¯ y nˆ ϕˆ . +μy κy Zy¯ y nˆ

(180)

32


Thus, in the norm C, we obtain s+1 s s s−1 s+1 s s s−1 c 0.5τ κ x 2 ( nˆ − nˆ ) = (nˆ − nˆ ) + Nˆ − Nˆ + Nˆ − Nˆ + 0.5τ (Qn + εn ) 2 h x C C C C μy κy μx κx ¯ ¯ N + ε¯ N ) + ¯ n + ε¯ n + Q ¯ N + ε¯ N ) , × (Qn + ε¯ n + Q (Q hx hy s+2 s+1 s+2 s+1 s+1 s 0.5τ κy c1 s+1 s ˆ ( nˆ − n) ( nˆ − nˆ ) = ˆ + N − Nˆ + Nˆ − Nˆ + 0.5τ (Qn + εn ) h2y C C C C κ μ μx κx ¯ ¯ N + ε¯ N ) + y y (Q ¯ n + ε¯ n + Q ¯ N + ε¯ N ) . × (Qn + ε¯ n + Q hx hy

(181)

(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 C 1 C2 , (183) τ ≤ min , C1 , C2 > 0, h2y h2x 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 /h2y ), (C2 /h2x )} 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 4.9: A positiveness of the mesh function nˆ can be stated by using the approach developed in [3]. 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 5–8 the laser pulse intensity distribution evolution (Figures 5, 6, 8) and the free-electron concentration evolution (Figure 7) computed for a set of the problem parameters a = 0.1, δ0 = 0.25, γ = 1000, τp = 1, q0 = 4, n0 = 0.01, ξ = 3, ψ = 2.553, κx = κy = 10−5 , μx = μy = 1, DAx = 10−5 , χ = 5, Lxc = 0.5, Lyc = −2, Lx = 1, Ly = 30. Grid steps on spatial coordinates and on time are chosen as hx = hy = 10−2 , τ = 10−3 . The iteration parameters are chosen as ε1 = 10−4 , ε2 = 10−6 . 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. We see in Figure 5 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 5(b,c)). With time increasing, the beam profile becomes much more disturbed one in comparison with the incident Gaussian


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).

33

34


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).


35

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).

profile. The reflected beam possesses also complicated profile in the x-direction as well as in ydirection. (Figure 5(d)). A part of the reflected beam moves faster than other ones (Figure 6). 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 6 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 7) shows that the laser beam diffraction action results in multi-domains formation for high absorption (Figure 7(b,c)). We see an appearance of the free-electron high concentration domains on both coordinates (Figure 7(b)). For example, such domains form in the section y = 0.12 which cannot appear without taking into account the beam diffraction. Figure 8 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.

36


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).

In Table 1 the values of invariant (79) and the difference analogue of integral equality (29)

¯¯ Q(t m) =

Py −1

P x −1

j=−Py +1 k=1

⎛ hy hx ⎝

⎞

2

m 2 |Aˆ m kj | − |Akj |

τ

0.5

+ βA δ0h δ

0.5 2⎠ |Am kj |

(184)

are shown. It should be noticed, that invariant (29) approximated by formula (184) with the sec¯¯ ond order on spatial and time coordinates. As we can see, both characteristics Q(tm ) and Q(t m) 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.


37

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

Invariant Q(tm )

¯¯ ) Invariant Q(t m

1.4 3.6 4 5 10

2.81732e−06 1.56882e−05 1.95336e−05 1.54824e−05 8.78868e−06

1.11856e−07 1.5281e−06 1.09616e−06 6.44219e−07 9.8865e−08

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.

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

References [1] G. Bachman and L. Narici, Functional Analysis, Dover Publications, Inc., Mineola, NY, 1998. [2] 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. [3] O.S. Bondarenko, Mathematical modeling of some problems of optical bistability on the basis of semiconductors, PhD thesis (1997) (in Russian). [4] 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. [5] 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. [6] 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. [7] A. Friedman, Partial Differential Equations of Parabolic Type, Dover Publications, Inc., Mineola, NY, 2008. [8] A. Ghadi and S. Mirzanejhad, Two-photon absorption effect on semiconductor microring resonators, Optik (Stuttg) 126 (2015), pp. 1645–1649. [9] H.M. Gibbs, Optical Bistability: Controlling Light with Light, Academic Press, New York, 1985. [10] 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. [11] 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.

38


[12] 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. [13] 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. [14] Yu N. Karamzin, A.P. Sukhorukov, and V.A. Trofimov, Mathematical Modeling in Nonlinear Optics, Moscow University Publishing, Moscow, 1989, p. 156 (in Russian). [15] P.K. Kwan and Y.Y. Lu, Computing optical bistability in one-dimensional nonlinear structures, Optics Commun. 238 (2004), pp. 169–175. [16] O. Ladyzhenskaya and N.N. Uraltseva, Linear and Quasilinear Elliptic Equations, Nauka, Moscow, 1973 (in Russian). [17] 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. [18] A.W. Leung, Nonlinear Systems of Partial Differential Equations: Applications to Life and Physical Sciences, Word Scientific Publishing Co. Pvt. Ltd, Hackensack, NJ, 2009. [19] L. Li, Optical bistability in semiconductor lasers under intermodal light injection, J.f Quantum Electronics 32 (1996), pp. 248–256. [20] J. Loustau, Numerical Differential Equations: Theory and Technique,ode Methods, Finite Differences, Finite Elements and Collocation, WordScientific Publishing Co. Pte. Ltd., New Jersey, 2016. [21] G.I. Marchuk, Methods of Numerical Mathematics, New York, NY, Springer, 1975. [22] J.D. Murray, Mathematical Biology II: Spatial Models and Biomedical Application, Springer-Verlag, New York, 2003. [23] 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. [24] A.A. Samarskii, Difference Schemes Theory, Nauka, Moscow, 1977 (in Russian). [25] 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). [26] R.A. Smith, Semiconductors, Cambridge Univ. Press, Cambridge, 1978. [27] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (3) (1968), pp. 506–517. [28] 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. [29] 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. [30] 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. [31] 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. [32] 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. [33] 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. [34] 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. [35] 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. [36] H. Wang, Numerical studies on the split-step finite difference method for nonlinear Schrödinger equations, Appl. Math. Comput. 170 (2005), pp. 17–35. [37] 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. [38] 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. [39] D. Weaire and J. P. Kermode, Dispersive optical bistability – numerical methods and definitive results, JOSA B 3 (1986), pp. 1706–1711. [40] 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.