A new way for decreasing of amplitude of wave reflected from artificial boundary condition for 1D nonlinear Schrödinger equation

In this report, we focus on computation performance enhancing at computer simulation of a laser pulse interaction with optical periodic structure (photonic crystal). Decreasing the domain before the photonic crystal one can essentially increases a computation performance. With this aim, firstly, we provide a computation in linear medium and storage the complex amplitude at chosen section of coordinate along which a laser light propagates. Then we use this time-dependent value of complex amplitude as left boundary condition for the 1D nonlinear Schrodinger equation in decreased domain containing a nonlinear photonic crystal. Because a part of a laser pulse reflects from faces of photonic crystal layers, we use artificial boundary condition. To decrease amplitude of the wave reflected from artificial boundary we introduce in consideration some additional number of waves related with the equation under consideration.


1.
Introduction. Computation performance enhancing for a problem of laser pulse interaction with optical periodic structure (photonic crystal (PC)) is considered in many papers. As rule, a length of PC element is about 0.1-1 µm or less. In this case, the length of domain before PC, which corresponds to linear propagation of a laser pulse with duration about 20-100 fs is many orders of magnitude larger than the PC length. As a consequence, one has to solve the corresponding Schrödinger equation in this domain and this solution is not of interest for practice. Therefore, decreasing size of the domain before the PC one can essentially increase a computation performance. With this aim, firstly, we provide a computation of femtosecond laser pulse propagation in a linear medium and storage the complex amplitude at chosen sections of coordinate, along which the laser pulse propagates. It should be noted that the corresponding amplitude distributions may be defined also explicitly using analytical formula for certain case of incident laser pulse shape. Then, we use one of these time-dependent distributions of complex amplitude as left boundary condition for nonlinear Schrödinger equation in decreased domain containing a nonlinear PC. Because we reformulate the problem we have to take into account two waves propagating in both directions in dependence of relation between left boundary of the domain and the section, in which the boundary condition is written. Taking into account that a part of laser pulse energy reflects from the faces of PC layers, therefore we use also artificial boundary condition (see, for example [1]- [5]) to increase computer simulation efficiency. To decrease amplitude of the reflected wave we introduce in consideration some additional number of waves related with wave propagating in the PC.
2. Problem statement. Let's consider 1D nonlinear Schrödinger equation, describing femtosecond laser pulse propagation in PC, which consists of two types of layers with identical length for each type of layers, for example. The equation is written below in dimensionless variables with initial condition which corresponds to wave propagating in a positive direction of z-axis if a parameter χ is positive (Fig. 1), and zero-value boundary conditions Above, a function A = A(z, t) is a complex amplitude slowly varying in time; t is a time; z denotes a spatial coordinate along which a laser pulse propagates; L t , L z are maximal values of time and space coordinate, L c is a coordinate of the laser beam center at initial time. D, β, χ are real coefficients, which satisfy the conditions D = 1 4πχ , β = πχ. Functions γ(z), ε(z) depends on the spatial coordinate. The PC begins by a layer with a dielectric permittivity ε 1 , then a layer with the dielectric permittivity ε 2 is placed. In the domain z 0 ≤ z ≤ L 0 before the PC a dielectric permittivity of medium is equal approximately to 1. Parameter L str = L 0 + (d 1 + d 2 )N str defines coordinate of the right boundary of the PC. A substrate with the dielectric permittivity ε 3 is displaced after the PC. It should be stressed if the dielectric permittivity of a medium before the PC is equal to 1 then it means that other dielectric permittivities are measured in units of the dielectric permittivity for this domain. N str is the number of layer pairs with dielectric permittivities ε 1 and ε 2 correspondingly.
As the solution of considered problem is interesting for large time interval, then it lead to necessity of formulating the boundaries of spatial domain far away from the PC if zero-value boundary conditions were formulated. Because a domain occupied by PC is much smaller than whole domain, then for multiple solution of the problem it is necessary to solve the equation in total domain and do a lot of useless mathematic operations.
Another feature of this problem concludes in a reflection of a laser pulse from the boundaries of each layer of the PC, therefore a wave propagating in a negative direction along z coordinate appears. This leads to necessity of domain increasing also. Therefore, it is necessary to use artificial boundary conditions, which allows eliminating the influence of reflected wave with certain accuracy. It should be stressed that the increasing domain size is unacceptable obviously for multidimensional problem.
Taking into account the peculiarities mentioned above, in this paper we propose a new way for increasing of the computer simulation efficiency for problem which described above and it is important that this approach can be easily generalized to the multidimensional case.
It consists of several stages. First of all, we should decrease a computational domain size before the PC. To do this let's consider a linear (γ = 0) problem (1)-(3) without the PC (γ = 0, ε = 1). Its size should be sufficient big to provide the laser pulse passing fully through the chosen sections z 0 and z s1 . In these sections the complex amplitude A is stored and let's denote them by A 0 (z 0 , t), A 0 (z s1 , t) correspondingly. The section z 0 will be used below to write the boundary condition for the equation in decreased domain and section z s1 will be necessary to calculate the wave reflected from the PC.
Thus, the equation (1) is solved in decreased domain. In this case there are two possible locations of the section z 0 . In the first case ( Fig. 2) it is not a left boundary of decreased domain. Then, we solve the equation (1) in a domain z 1 < z < L z with zero-value initial condition A| t=0 = 0 (4) and boundary conditions Here, a parameter Ω is a local wave number may be defined, in particular, by initial condition (2): Ω = 2πχ. It should be stressed that the statement of an artificial boundary condition on the left boundary of considered domain is necessary because the wave which defined in section z 0 propagates in both directions. As it was mentioned above the reflection of laser beam from the PC layers leads to appearance of a wave propagating in negative direction of z axis also. Inaccuracy of the artificial boundary condition leads to a reflection of laser pulse from the boundary and while reflected wave with sufficiently large amplitude will reach the PC, then the solution becomes false.
If the section z 0 is the left boundary of the domain (Fig. 3), the boundary conditions for the equation (1) are defined in the form with zero-value initial condition (4).
Function A − (z, t) describes the propagation of wave in the negative direction of z axis. It should be stressed while the propagating wave is not reflected from the PC the amplitude To calculate the complex amplitude A − (z, t) let's formulate an additional problem. To do this, choose a section z s2 , which is placed closer to the section z = 0 than the section z 0 (Fig. 4). Then, a complex amplitude of this wave is a solution of the following problem with zero-value initial condition and with the following boundary conditions where A 0 (z s1 , t) is the complex amplitude stored at previous stage (see above), A(z s1 , t) is the solution of the equation (1) with boundary conditions (6) and zero-value initial condition (4). It should be stressed that the right equality (9) defines in a section z s1 the wave, which reflected from the PC layers and propagates towards the incident wave. The problem for the complex amplitudes A, A − with corresponding boundary and initial conditions formulated above leads to the elimination of the wave A reflected in section z 0 while a wave A − propagating in the positive direction due to its reflection from the boundary in the section z s2 appears. Its appearance may be caused by inaccuracy of the artificial boundary condition (9). This procedure can be repeated arbitrary number of times to remove the reflected wave A − , which propagates in positive direction. To do this we choose the section z s3 in the interval (z s2 , z s1 ) and will use a value of complex amplitude A − (z s3 , t) as a boundary condition for new problem. As a result, in section z s2 also will be absent reflected wave A − until the wave A (2) − , which propagates in positive direction of spatial coordinate and appear because of inaccuracy of the boundary condition in section z s4 , will return to this section. Then the problem for complex amplitudes A, A − , A with initial condition A| t=0 = 0 (11) and zero-value boundary conditions where A − (z, t) is solution of the following problem with zero-value initial condition and with the following boundary conditions and so on. Then equation for A with zero-value initial condition and boundary conditions 3. Conservative finite-difference scheme. For definiteness let's assume section z s2M is greater than 0. Then in the interval 0 ≤ z ≤ L z and 0 ≤ t ≤ L t we introduce uniform grids Let's define the mesh functions − | 2 and corresponding mesh functions ε = ε(z j ) and γ = γ(z j ).
For the problem (10) -(18) we write the following conservative finite-difference scheme with initial condition A 0,j = 0, j = j 0 + 1, N z − 1 (21) and following boundary conditions where z 0 = j 0 h. Choice of different time layers in (22) will be explained below. The scheme for reflected wave A − (z, t) is with initial condition A −0,j = 0, j = j s2 + 1, j s1 − 1 (24) and boundary conditions where z s1 = j s1 h, z s2 = j s2 h. And so on. The scheme for A The problem (20) is nonlinear one that is why for its solution we use an iteration process. For example, this process can be implemented in following way with initial condition A 0,j = 0, j = j 0 + 1, N z − 1 (30) and following boundary conditions (32) The iteration process is stopped if the following condition, for example, is valid Then we continue to calculate with iterations only main equation for function A(z, t). The scheme for other sought functions described by equations (23) -(28) written above.
As it well known, an effective method for solving the equations written above is the Thomas algorithm.

4.
Computer simulation results. We will consider the results of computer simulation for the parameters: z 0 = 10, z s1 = 14, z s2 = 4, z s3 = 6, z s4 = 0, z s5 = 2, z s6 = −4, L z = 40.0, L c = 5, ε 1 = 2, ε 2 = 3, ε 3 = 1, γ 1 = 1,  In the time moment, when the reflected wave from left boundary reaches the PC boundary, error occurs in the solution inside a PC. The value of this error corresponds to amplitude value of reflected wave. From Table 1 we can see that with increasing a number of additional steps M increasing the time of calculation until the moment when the reflected wave reaches selected section. Accordingly, the program calculation time before the error occurred inside the PC increases. Also the computer simulation time for the problem with additional functions (A − , . . . , A (M ) − ) in decreased domain is smaller than the corresponding time for the problem in the full domain.
Computer simulation results indicate that the results obtained in full domain and the results obtained in decreased domain coincide up to a 10 −6 .
Let's consider the intensity profile propagation of function A in PC. In Fig.6 (a)-(c) we can see the occurrence and propagation of the reflected wave generated at the boundary of the PC layers. Note that the wave is not reflected from the left boundary of the domain due to the artificial left boundary conditions described above. In Fig.6 (d) shows that the wave reflected from the left boundary of the last additional function A

5.
Conclusion. In this paper we have proposed a new method for computation performance enhancing at numerical solution of 1D nonlinear Schrödinger equation. Using a computer simulation results it is shown that the restatement of the problem leads to significantly decreasing computer simulation time. With increasing a number of additional problems with respect to new complex amplitudes M the computer time interval, during which a right evolution of beam propagation takes place, is increased. Also with increasing a number of additional stages the maximum intensity of wave reflected from left artificial boundary condition decreases.
Developed approach can be applied with success for multi-dimensional problem.