EFFICIENT FINITE DIFFERENCE METHODS FOR THE NONLINEAR HELMHOLTZ EQUATION IN KERR MEDIUM

. In this paper, we consider a kind of eﬃcient ﬁnite diﬀerence methods (FDMs) for solving the nonlinear Helmholtz equation in the Kerr medium. Firstly, by applying several iteration methods, we linearize the nonlinear Helmholtz equation in several diﬀerent ways. Then, based on the resulted linearized prob- lem at each iterative step, by rearranging the Taylor expansion and using the ADI method, we deduce a kind of new FDMs, which also provide a route to deal with the problem with discontinuous coeﬃcients.Finally, some numerical results are presented to validate the eﬃciency of the proposed schemes, and to show that our schemes perform with much higher accuracy and better convergence compared with the classical ones.


1.
Introduction. The propagation of electromagnetic waves in some materials is usually modeled by the famous Maxwell's equations [4] with various proper medium responses. These significant responses reflect the material's properties, such as the magnetic permeability and electric permittivity with respect to the location and the frequency of the propagating field. When it turns to high intensity radiation situation, not only the medium quantities may depend on the magnitude of the propagating field, but also the response will become nonlinear. In nonlinear optics, one may often focus on the propagation of monochromatic waves, such as continuous high intensity laser beams. In this case, some reasonable assumptions (see [3]) simplify the Maxwell's models to a nonlinear Helmholtz (NLH) equation [5,13] ∆E(X) + ω 2 0 c 2 n 2 (X, |E|)E(X) = f (X), n 2 (X, |E|) = n 2 0 (X) + 2n 0 (X)n 2 (X)|E| 2 , where E(X) denotes the electric field, X = [x 1 , · · · , x D ] is the spatial coordinate (x D is called the longitudinal direction which also denotes as z as that in Fig. 1 and the rest coordinates are regarded as transverse directions), f (X) is a given function, ω 0 is the frequency, c is the speed of light in vacuum, ∆ = ∂ 2 x1 + · · · + ∂ 2 x D is the Ddimensional Laplacian operator, n 0 (X) and n 2 (X) are the linear index of refraction and the Kerr coefficient, respectively. And n 0 , n 2 are always assumed to be real which means that all mediums are transparent. The above NLH equation can be used to govern the propagation of linearly polarized, time-harmonic electromagnetic waves in Kerr-type dielectrics, which can produce some important nonlinear optical effects , such as the optical bistability [6,30] and spatial solitons [20]. Let the Kerr medium be surrounded by the linear homogeneous medium in which n 0 = n ext 0 , n 2 = 0. We introduce the linear wave number k 0 = ω 0 n ext 0 /c and denote the normalized quantities v(X) = n 0 (X) n ext 0 2 , ε(X) = 2n 2 (X)n 0 (X) Then, equation (1) can be rewritten as ∆E(X) + k 2 0 v(X) + ε(X)|E(X)| 2 E(X) = f (X).
Since the Kerr medium coefficient n 2 is discontinuous at the interfaces z = 0 and z = Z max , and n 0 may also be discontinuous at these interfaces, it naturally leads to the discontinuities of the coefficients v and ε in equation (2). Thus, for the homogeneous Kerr medium, the coefficients v and ε are piecewise constants as follows Remark 1. When the Kerr medium is inhomogeneous which means the whole Kerr medium will be cut into pieces by some other mediums, the linear medium for example. In this case, the discontinuities of the coefficients v and ε in equation (2) will be more complex, we refer [2,3] for the detail.
When the electric field E and the material coefficients n 0 and n 2 vary only in one direction z, the model (2) reduces to the 1D nonlinear Helmholtz equation To solve (3) in the interval [0, Z max ], some boundary conditions are needed. It is well-known that, at the interfaces z = 0 and z = Z max , the electronic field E and its first derivative are continuous [3], According to this fact, the so-called two-way boundary conditions which are also used in [3,14,15] are as follows where i = √ −1 is the imaginary unit.
Remark 2. According to [3], the above boundary condition (4) is developed from the inhomogeneous Sommerfeld type relation and its complete form is where E 0 inc and E Zmax inc are the incident waves which shoot into the computational domain from z = −∞ and z = +∞, respectively. In this paper, we assume that the laser beam only shoots in the Kerr medium from z = −∞, so E Zmax inc = 0. And the problem (3) can be rescaled by E = E/E 0 inc ,ε = ε E 0 inc 2 . Without loss of generality , we assume E 0 inc = 1. In 2D case, assuming the linear homogeneous medium in R 2 \Ω 0 is truncated by a finite domain Ω and Ω 0 ⊂ Ω is fully filled with the Kerr medium (see Fig. 2), then (2) is followed by, where n denotes the unit outer normal vector to Γ and g = ∂Einc Many researches are done for the NLH equation. Through transforming it into a phase-amplitude equation, Chen and Mills proposed an approach to obtain the closed form solution of the NLH equation with a single nonlinear layer [7] and multilayered structures [8]. And by using the multidimensional generalization of the nonlinear Schrödinger equation, the exact solution of the NLH equation in special case was also considered in [18]. Moreover, in [11,12], the existence and asymptotic behavior of the real-valued standing wave solution of the NLH equations were analyzed. On the other hand, numerical methods for the NLH equations are also investigated. Fibich and the collaborators studied the NLH equation in [1,2,3,14,15]: In [14], the authors constructed a two-way artificial boundary condition for the NLH equation to ensure that not only the backscattered waves generate no reflection but also the correct value of the incoming wave can be imposed; the NLH equation was also solved by using nonorthogonal expansions in [15]; by coupling with a new technique of the separation of variables, a fourth order finite difference scheme was developed in [1], and the algorithm was also extended to the three-dimensional axially symmetric problem; in [2,3], the authors solved the NLH equations by developing an efficient Newton's iteration method to deal with the strong nonlinearity. In addition, a finite element method which can approximate the discontinuous coefficient problem is constructed in [22,23]. Recently, Wu and Zou proved the existence and uniqueness of the NLH equation, and also analyzed the stability and the error estimate with explicit wave numbers for the finite element approximation in [29]. And the author proposed a robust modified Newton's method in [31].
There are many difficulties when the NLH equation is approximated by using numerical schemes. Firstly, since the NLH equation is a strong nonlinear problem, we need to search a robust iteration method for solving it. Secondly, similar to the Helmholtz equation, the solution of this problem is highly oscillating with a high number (see [19]). Moreover, it usually contains discontinuous coefficients due to the different propagating mediums. There are also other issues to be solved, such as the strong indefinite linear system generated from this equation and so on. But, in this paper, we mainly focus on the case that it admits highly oscillated solutions. It is well-known that, the FDM is one of the most frequently used numerical methods due to its simple structure. Usually, the FDM is constructed through approximating the derivative terms in the original equation with some difference quotients which are obtained by the Taylor's expansion directly [21,32,33], this may lead the FDM to suffer from some disadvantages, low computation accuracy for example. To simulate the Helmholtz equation with high wave numbers in one dimension, a new finite difference scheme was proposed in [25]. Being different from the classical one directly based on the Taylor expansion, the new finite difference method is constructed by a rearranged formula which can contain more regularity information of the solution, and thus more accurate approximations were achieved. Recently, this kind of schemes were applied to the higher dimensional problems in [16,17,26,27,28]. But all of these problems investigated above are linear problems. In this article, we will extend this idea to solve the NLH equation. Through some iteration methods, the NLH equation will be linearized as a linear one at each iterative step firstly. Several iteration methods are considered, including the classical ones and the error correction iteration method [24] in which the original iteration solution was modified by a residual. Then, based on the above resulted linear problem, the new finite difference scheme is constructed, which is naturally suitable for the problem with discontinuous coefficients to match the different propagating mediums (Kerr and linear mediums).
The rest of this paper is organized as follows. In the next section, we apply several iteration methods including the error correction one to linearize the NLH equations. Then, in Section 3, after constructing the new finite difference scheme for the 1D linearized equation, we extend the scheme to the 2D problem by using the ADI technique. To test the efficiency of the numerical schemes, some numerical experiments are performed in Section 4. We finally make conclusions in Section 5.
2. Iteration methods. To solve the nonlinear equation, an iteration method is needed. In this section, we introduce several kinds of iteration methods for solving the NLH equation. For convenience, we rewrite the problem as, where L := ∆ + k 2 0 v denotes the linear operator. Frozen-nonlinearity iteration may be the simplest iteration method in which the nonlinear term is frozen as a known quantity. For example, by replacing |E| 2 with the previous iteration solution, (6) is linearized by where E l donates the lth iteration solution. However, it is well known that the frozen-nonlinearity iteration is only the first order convergent. An alternative iteration method is the Newton's method. Letting where E is the conjugation of E. Assuming E = s + E l (s is small enough) and using the Taylor's expansion, we can approximate (8) as Letting s = E − E l ≈ E l+1 − E l and substituting (8) and (9) into (10), we can get the Newton's iteration method, Furthermore, a modified Newton's method was proposed in [31] by replacing E l+1 with E l in (11) , that is In this paper, we also employ the error correction method in [24] for solving the nonlinear Helmholtz equation. Next, we will show the process. For simplicity, the above three iteration methods can be rewritten as a general formula where in (12), (12).
Assuming that µ l+1 is the error between E l+1 and the exact solution E, i.e., µ l+1 = E − E l+1 , then, subtracting (13) from the original equation (6), we get the error equation as follows Since freezing the terms µ l+1 2 and µ l+1 2 at µ l in the above equation, (14) yields where Furthermore, for the error on the boundary, there holds Obviously, by solving (15), the original solution E l+1 could be modified by a more accurate one E l+1 = E l+1 + µ l+1 which is expected to be useful in decreasing the iteration number. And the iterative procedure with the error correction can be concluded in the following Algorithm 1.
3. Efficient FDMs for linearized equations. After applying the iteration methods introduced in the above section, the NLH equation is linearized to linear problems at each iterative step. To implement these methods, a spatial discretization scheme is needed. From the stability analysis in [29], we know that the solution of the NLH equation satisfying ||E|| H j ≤ c 0 k j−1 0 (j = 0, 1, 2) with c 0 being a positive constant and H j being the classical norm in the Sobolev space. This implies the oscillation of the field E when k 0 is large. In this section, we will develop a family of efficient finite difference schemes for solving (13) and (15) with corresponding boundary conditions, respectively. Since these two equations have the same form, by taking (13) for an example, we will exhibit the construction of the new finite difference scheme.
3.1. 1D problem. Recalling the iteration formula for the 1D NLH equation When the frozen-nonlinearity iteration (7) and the modified Newton's method (12) are used, i.e., M = 0, we can rewrite (17) as where τ = − k 2 0 v + N . (19) has the same form as the linear Helmholtz equation, so we can extend the idea in [25,27] to solve it.
Since the parameters τ and f may be discontinuous at the interfaces of the Kerr medium and the linear medium, we divide the computational domain by a mesh in which these interfaces are specially selected as mesh points. Assume that Taylor's expansion tells that And according to (19), there hold Substituting (22)-(23) into (20) and (21), it yields where Then, eliminating ( E l+1 m ) (1) in (24) and (25), we get the scheme for the interior point as follows For the boundary points z = 0, z = Z max , directly letting m = 1 and m = N in (24) and (25), respectively, we have Then, substituting the boundary condition (18) into the above formulas, we have the numerical schemes for z = 0 and z = Z max as follows Obviously, taking different k in (26)-(28) could obtain different finite difference schemes. For example, letting k = 1, the new finite difference schemes are Remark 3. In fact, a more accurate numerical scheme could be developed. For example, when the frozen-nonlinearity iteration method is used, (19) is the equation with variable coefficient, it yields, . In this case, like (22)-(23), we can get more precise formulas Then, according to the Taylor's series, we get Thus, eliminating the terms E Finally, by approximating τ m with difference quotients in (30), we can get a more accurate scheme. Moreover, the corresponding schemes for boundary points can also be constructed easily like (27)- (28).
When the Newton's iteration method (11) is considered, (17) needs to be separated into the real and imaginary parts due to the existence of E l+1 , which is followed by where Taking the real part equation (31) for an example, following the same process for (19), we have, at any interior point z m , that and where Eliminating R m in (37) and (38), we get where Similarly, setting m = 1 and m = N in (37) and (38), respectively, and according to the boundary condition (32), we get Obviously, by retaining different terms in the right hand side of (39)-(41), we can also get a series of finite difference schemes for the real part equation (31)- (32). For example, taking k = 1 in (39), it yields

By approximating φ
(1) in the above formula, we can get a finite difference scheme for (31) at the interior point z m (m = 2, · · · , N − 1), where Similarly, letting k = 1 in (40)-(41), we have By using the boundary condition (34) and approximating ( f + R,1 ) (1) and respectively, we get a finite difference scheme at the boundary points z = 0 and z = Z max , For the imaginary part (33)-(34), the same procedure can be applied to develop the new finite difference scheme at any interior point z m (m = 2, · · · , N − 1), where And the schemes for boundary points z = 0 and z = Z max are Obviously, (42)-(47) constitute a finite difference scheme for the system of equations (31)-(34). According to the above process, it can be found that, through translating the high order terms E (2n) and E (2n+1) into the lower ones, much more terms in the Taylor's expansion could be included in the new finite difference scheme. Therefore, the new scheme is expected to achieve much better computational accuracy. Moreover, the case of discontinuous coefficients is considered fully and naturally in these schemes.
According to (5) and Fig. 2, the parameters R, I, I , R are also discontinuous at the interface ∂Ω 0 . So, we divide the computational domain Ω into N x × N y parts such that ∂Ω 0 (including its four vertexes) is overlapped with some mesh points. And let (x m , y n )(1 m N x , 1 n N y ) be the mesh points and h − m = x m − x m−1 , h + m = x m+1 − x m , k − n = y n − y n−1 , k + n = y n+1 − y n be the mesh sizes.
It is well-known that the ADI method is used to simulate a high-dimensional problem by solving a series of one-dimensional problems. Based on this, by directly separating the real part equation (48) into two 1D equations in x, y directions at the interior point (x m , y n ), we have ∂ 2 R m,n ∂x 2 = γ x R m,n R m,n + γ x I m,n I m,n + f xR;m,n , ∂ 2 R m,n ∂y 2 = γ y R m,n R m,n + γ y I m,n I m,n + f yR;m,n , where γ x + γ y = 1, f xR + f yR = f R . Similar to (42), the new finite difference schemes for (52) and (53) can be directly got as follows − H − y γ y I + m,n L + y,1 + H + y γ y I − m,n L − y,1 I m,n − H − y γ y I + m,n X + y,1 + H + y γ y I − m,n X − y,1 where Combining (54) and (55), we get the new finite difference scheme for (48) at the interior point (x m , y n )(2 m N x − 1, 2 n N y − 1), where Similar to the interior points, the new finite difference scheme for each boundary point is also constructed by developing two schemes in x and y directions, respectively. For example, when (x m , y n ) ∈ Γ 1 \{Γ 1 ∩ Γ 3 , Γ 1 ∩ Γ 4 }, the scheme (54) with n = 1 can be seen as the new finite difference scheme in x direction, and in y direction, according to (53), we get ∂y .
Remark 4. In fact, for the 2D equation, when the frozen-nonlinearity iteration and the modified Newton's method are used, it is not necessary to separate (5) into real and imaginary parts. In this case, the 2D equation can be divided into two 1D equations directly. Furthermore, since two 1D equations are separated from a 2D equation, we assume that f xR , f yR are not discontinuous, and ε = 0 in f at the interface ∂Ω 0 of the Kerr medium and the linear medium. 4. Numerical experiments. In this section, we will show some numerical tests to verify the efficiency of the scheme proposed in the above section. And we set δ = 10 −8 in Algorithm 1.
To test the accuracy of the proposed scheme, with v = 1, k 0 = 8, Z max = 10 and the initial guess E 0 (z) = e ik0z , we firstly solve the 1D equation (3) with ε = 0.01 and 0.1 by using the new scheme based on the frozen-nonlinearity iteration method with k 0 h = 1. We can see that the numerical solutions can highly match the reference ones (see Figs. 3-4). Then, selecting frozen-nonlinearity iteration, the errors in l ∞ -norm with ε = 0.01 are compared in Table 1 among different numerical schemes: the standard finite difference (SFD) scheme, the compact finite difference (CFD) scheme, the finite volume (FV) method proposed in [2] and two new schemes (26)(k = 1) and (29). Clearly, the newly proposed finite difference schemes can achieve the best accuracy. These all imply that the newly proposed finite difference scheme can approximate the high oscillation solution of the NLH equation effectively. Then, under the same computational environments, but with Z max = 1, the iteration numbers are compared among four iteration methods for varying ε and k 0 in Table 2. It can be found that, among the classical ones, the Newton's iteration method converges fastest, but the error correction has obvious advantage in decreasing the iteration number especially for the cases of large ε and k 0 . The values of k 0 or ε seem have little influence on the iteration number when the error correction method is convergent.
Furthermore, we simulate the optical bistability by using the proposed finite difference scheme. Firstly, letting the transmittance T := E(Z max ), we solve (3)-(4) with varying ε and plot |T | 2 in Fig. 5, the result is very similar to that in [2]. It worth noting that, since the Newton' s method is sensitive to the initial guess, we start the Newton's method with a initial guess which is obtained by solving (3)-(4) with the frozen-nonlinearity iteration method when ε = 0.01, E 0 = 0. The solutions corresponding to the rest ε are obtained by increasing the value of ε step by step through some proper ∆ε and the solution with ε i is selected as the initial guess for ε i+1 = ε i + ∆ε. Then, we select much smaller ∆ε to compute the switchback-type multi-solution near ε = 0.724. In Fig. 6, we choose the solution with ε = 0.722 (point A) and ε = 0.726 (point E) as the initial guess to calculate the solutions corresponding to the lower branch and the upper branch, respectively. It worth to note that, in the lower branch, when we increase the value of ε at point F (ε ≈ 0.72524), the value of |T | 2 will directly jump to the one corresponding to point H and then varies following the route H → E. Similarly, in the upper branch, decreasing the value of ε at point S(ε ≈ 0.72376) leads to the value of |T | 2 follows the route S → G → A. Obviously, |E| 2 has three different values when ε ∈ (0.72376, 0.72524). To simualte the middle branch, we firstly set ε = 0.724226 and solve (3)-(4) to obtain the solution corresponding to point C with the initial guess which is selected as the linear combination of the solution at point B (ε = 0.724226 in the lower branch) and the solution at point D (ε = 0.724226 in the upper branch). Then, the middle branch is captured by selecting the solution corresponding to point C as the initial guess. Successfully approximating the optical bistability also indicates the robustness of the new finite difference method.       Figure 6. Switchback-type non-uniqueness of |T | 2 near ε = 0.724 for the 1D problem.
In Fig. 7, we exhibit the exact solution and the numerical solution obtained by the new finite difference scheme with k 0 = 100, h = 1/400, it can be observed that the solution is highly oscillating, but the numerical solution can match the exact solution well.
The intensities of the incident field and the total field for different cases are shown in Fig. 8 and Fig. 9. When only one incident wave E 1 inc comes into Ω from south (see Fig. 8(A)), it can be found that the incident wave can pass through the Kerr medium almost without any change when ε = k −2 0 (see Fig. 8(B)), while the field is totally different when ε = 0 (see Fig. 8(C)). When there are two solitons with the field E 2 inc come into the medium from south and east respectively, they meet and collide in Ω 0 . In the case of ε = k −2 0 , these two solitons are also nearly unchange when they pass through the Kerr medium (see Fig. 9(B)). It is worth to note that, comparing with Fig. 8(B), the amplitudes near (0, −1/2) and (1/2, 0) have bigger values, which means that the backward scattering becomes stronger when two solitons colliding. Similarly, the total field changes a lot when the same medium is fulfilled in Ω 0 and Ω \ Ω 0 (see Fig. 9(C)). Therefore, from the transmission and collision of the nonparaxial solitons examples, we can conclude that, the Helmholtz equation and the NLH equation are much different despite the nonlinear coefficient is very small (here ε = 10 −4 ). And the scheme studied in this paper can release these clearly.

5.
Conclusions. In this paper, we construct a kind of new finite difference schemes for solving the nonlinear Helmholtz equation based on some iteration methods.
Numerical results indicate that, the proposed scheme not only can approximate the high oscillation solution with better computational accuracy, but also can be used to simulate some physical phenomenons in the Kerr medium, such as the optical bistability and the collision of nonparaxial solitons. Moreover, without any extra consideration, this finite difference scheme also provides a route to deal with the problems with discontinuous coefficients or source terms. Thus, it can be extended to much more complex cases, such as the multi-layered Kerr mediums propagating problem and the nonlinear Maxwell's equations.