Efficient finite difference methods for the nonlinear Helmholtz equation in Kerr medium

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]

where $E(\mathbf{X})$ denotes the electric field, $\mathbf{X} = [x_1,\cdots,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(\mathrm{\mathbf{X}})$ is a given function, $\omega_{0}$ is the frequency, $c$ is the speed of light in vacuum, $\Delta = \partial_{x_{1}}^{2}+\cdots+\partial_{x_{D}}^{2}$ is the $D$-dimensional Laplacian operator, $n_0(\mathbf{X})$ and $n_2(\mathbf{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_0^{\text{ext}}, n_2 = 0$. We introduce the linear wave number $k_{0} = \omega_{0} n_{0}^{\text{ext}} /c$ and denote the normalized quantities

Then, equation (1) can be rewritten as

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 $\varepsilon$ in equation (2). Thus, for the homogeneous Kerr medium, the coefficients $v$ and $\varepsilon$ 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 $\varepsilon$ 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 = \sqrt{-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_{\mathrm{inc}}^{0}$ and $E_{\mathrm{inc}}^{Z_{\max }}$ are the incident waves which shoot into the computational domain from $z = -\infty$ and $z = +\infty$, respectively. In this paper, we assume that the laser beam only shoots in the Kerr medium from $z = -\infty$, so $E_{\mathrm{inc}}^{Z_{\max }} = 0$. And the problem (3) can be rescaled by $\widetilde{E} = E / E_{\mathrm{inc}}^{0}, \tilde{\varepsilon} = \varepsilon\left|E_{\mathrm{inc}}^{0}\right|^{2}$. Without loss of generality, we assume $E_{\mathrm{inc}}^{0} = 1$.

In 2D case, assuming the linear homogeneous medium in $\mathbb{R}^2\backslash \Omega_0$ is truncated by a finite domain $\Omega$ and $\Omega_0\subset\Omega$ is fully filled with the Kerr medium (see Fig. 2), then (2) is followed by,

where $n$ denotes the unit outer normal vector to $\Gamma$ and $g = \frac{\partial E_{\mathrm{i n c}}}{\partial n}+i k_{0} E_{\mathrm{i n c}}$, $\varepsilon\equiv 0$ in $\Omega\setminus\Omega_{0}$, $\varepsilon\neq 0$ in $\Omega_{0}$.

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 $\mathcal{L}: = \Delta+k_{0}^2v$ 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 $l$th 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

then there hold,

where $\overline{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}\thickapprox 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 $\overline{E^{l+1}}$ with $\overline{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

Assuming that $\mu^{l+1}$ is the error between $\widetilde{E}^{l+1}$ and the exact solution $E$, i.e., $\mu^{l+1} = E-\widetilde{E}^{l+1}$, then, subtracting (13) from the original equation (6), we get the error equation as follows

Since

freezing the terms $\left(\mu^{l+1}\right)^{2}$ and $\left|\mu^{l+1}\right|^{2}$ at $\mu^{l}$ in the above equation, (14) yields

where

Furthermore, for the error on the boundary, there holds

Obviously, by solving (15), the original solution $\widetilde{E}^{l+1}$ could be modified by a more accurate one $E^{l+1} = \widetilde{E}^{l+1}+\mu^{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}\leq c_0k_0^{j-1}\; (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., $\mathcal{M} = 0$, we can rewrite (17) as

where $\tau = -\left(k_{0}^{2} v+\mathcal{N}\right)$. (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 $\tau$ and $\widetilde{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 $\{z_m\}(1\leqslant m\leqslant N)$ with $z_1 = 0,z_N = Z_{\text{max}}$ are the mesh points. And for any interior point $z_m$, setting $h_{m}^{-} = z_{m}-z_{m-1}, h_{m}^{+} = z_{m+1}-z_{m}$ and

Taylor's expansion tells that

And according to (19), there hold

Substituting (22)-(23) into (20) and (21), it yields

where

with

Then, eliminating $(\widetilde{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_{\text{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_{\text{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,

where $\tau(z) = -k_0^2\left[v+\varepsilon E^{l}(z)\overline{E^{l}}(z)\right]$.

In this case, like (22)-(23), we can get more precise formulas

where $C_{2 n-1}^{2 n-2} = (2n-1)!/(2n-2)!$.

Then, according to the Taylor's series, we get

where

Thus, eliminating the terms $\widetilde{E}_{m}^{(1)}$, we can obtain

Finally, by approximating $\tau_{m}^{(1)},\widetilde{f}_{m}^{(1)}$ 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 $\overline{\widetilde{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^{(1)}$ 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 $\phi_{m}^{(1)} (\phi = I,\widetilde{f}_R)$ with $\frac{\phi_{m+1}-\phi_{m-1}}{h_{m}^{+}+h_{m}^{-}}$ in the above formula, we can get a finite difference scheme for (31) at the interior point $z_m (m = 2,\cdots,N-1)$,

where

with

Similarly, letting $k = 1$ in (40)-(41), we have

By using the boundary condition (34) and approximating $(\widetilde{f}_{R, 1}^{+})^{(1)}$ and $(\widetilde{f}_{R, N}^{-})^{(1)}$ with $\frac{\widetilde{f}_{R, 2}^{+}-\widetilde{f}_{R, 1}^{+}}{h_1^+}$ and $\frac{\widetilde{f}_{R, N}^{-}-\widetilde{f}_{R, N-1}^{-}}{h_N^-}$, respectively, we get a finite difference scheme at the boundary points $z = 0$ and $z = Z_{\text{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,\cdots,N-1)$,

where

with

And the schemes for boundary points $z = 0$ and $z = Z_{\text{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.

3.2.   2D problem

In this section, we extend the new finite difference scheme to the 2D problem by applying the ADI method [9,10]. Similar to (31)-(34), the 2D equation (5) need to be divided into real and imaginary parts like

where $g_{l R} = \operatorname{real}\left(g_{l}\right), g_{l I} = \operatorname{imag}\left(g_{l}\right)(l = 1,2,3,4)$.

According to (5) and Fig. 2, the parameters $\widehat{\mathcal{R}},\widehat{\mathcal{I}},\overset{\frown}{\mathcal{I}},\overset{\frown}{\mathcal{R}}$ are also discontinuous at the interface $\partial\Omega_0$. So, we divide the computational domain $\Omega$ into $N_x\times N_y$ parts such that $\partial\Omega_0$ (including its four vertexes) is overlapped with some mesh points. And let $(x_m,y_n)(1 \leqslant m \leqslant N_{x}, 1 \leqslant n \leqslant 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

where $\gamma_{x}+\gamma_{y} = 1, \widetilde{f}_{xR}+\widetilde{f}_{yR} = \widetilde{f}_R$.

Similar to (42), the new finite difference schemes for (52) and (53) can be directly got as follows

where

Combining (54) and (55), we get the new finite difference scheme for (48) at the interior point $(x_m,y_n) (2\leqslant m\leqslant N_x-1, 2\leqslant n\leqslant N_y-1)$,

where

with

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)\in\Gamma_1\backslash\{\Gamma_1\cap\Gamma_3,\Gamma_1\cap\Gamma_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

Then, substituting the corresponding boundary condition into the above formula, we have

So, combining (54)($n = 1$) and (57), the new finite difference scheme for the boundary points on $\Gamma_1\backslash\{\Gamma_1\cap\Gamma_3,\Gamma_1\cap\Gamma_4\}$ is

where

Similar to (57), on other three boundaries (excluding the vertexes), we also have: for the boundary points on $\Gamma_2\backslash\{\Gamma_2\cap\Gamma_3,\Gamma_2\cap\Gamma_4\}$ ($y$ direction):

for the boundary points on $\Gamma_3\backslash\{\Gamma_3\cap\Gamma_1,\Gamma_3\cap\Gamma_2\}$($x$ direction):

for the boundary points on $\Gamma_4\backslash\{\Gamma_4\cap\Gamma_1,\Gamma_4\cap\Gamma_2\}$($x$ direction):

Thus, applying the same process, the new finite difference schemes for these boundary points can also be written as (58) with different $\mathbf{A_b,R_b}$ and $F_b$. The details are: for the points on $\Gamma_2\backslash\{\Gamma_2\cap\Gamma_3,\Gamma_2\cap\Gamma_4\}$, there hold

and for the points on $\Gamma_3\backslash\{\Gamma_3\cap\Gamma_1,\Gamma_3\cap\Gamma_2\}$, there hold

and for the points on $\Gamma_4\backslash\{\Gamma_4\cap\Gamma_1,\Gamma_4\cap\Gamma_2\}$, there hold

According to (57), (59)-(61), the new finite schemes for four vertexes can also be obtained. For example setting $m = 1$ and $n = 1$ in (57) and (60), respectively, and combining them together, we get the finite difference scheme at the vertex $\Gamma_1\cap\Gamma_3$,

where

Similarly, for the rest three vertexes, their new finite difference schemes can be also concluded in (62) with different $\mathbf{A_v,R_v}$ and $F_v$: for the vertex $\Gamma_1\cap\Gamma_4$, there hold

and for the vertex $\Gamma_2\cap\Gamma_3$, there hold

and for the vertex $\Gamma_2\cap\Gamma_4$, there hold

For the imaginary part (50)-(51), the new finite difference scheme can be produced in the same way. And the details are omitted here.

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 $\widetilde{f}_{xR},\widetilde{f}_{yR}$ are not discontinuous, and $\varepsilon\neq0$ in $\widetilde{f}$ at the interface $\partial\Omega_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 $\delta = 10^{-8}$ in Algorithm 1.

4.1.   1D problem

Firstly, let $k_0,\; v$ and $\varepsilon$ be given, and set the exact solution of the 1D equation to be (see [7])

where $\lambda(z)$ and $\varphi(z)$ are both real and satisfy

with

and $\text{cn}[\cdot|\cdot]$ being the Jacobi elliptic function and $W\in[0,k_0]$ being a parameter need to be determined.

Moreover, at $z = 0$, there hold

with

Thus, putting (66) with $z = 0$ into (70) and solving the nonlinear equation, we can determine $W$. Then, we obtain $\beta(z_m)(1\leqslant m\leqslant N)$ by (66), $A$ by $A = k_{0}\left(v+1\right) W+\frac{ \varepsilon}{2} W^{2}$, $\lambda'(0)$ by (72), $\varphi(0)$ by (71), and $\varphi(z_m)$ by (64) one by one. Finally, putting $\lambda(z_m) = \sqrt{\beta(z_m)}$ and $\varphi(z_m)$ into (63), the exact solution is determined.

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^{i k_{0} z}$, we firstly solve the 1D equation (3) with $\varepsilon = 0.01$ and 0.1 by using the new scheme based on the frozen-nonlinearity iteration method with $k_0h = 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^{\infty}$-norm with $\varepsilon = 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 $\varepsilon$ 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 $\varepsilon$ and $k_0$. The values of $k_0$ or $\varepsilon$ 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_{\text{max}})$, we solve (3)-(4) with varying $\varepsilon$ 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 $\varepsilon = 0.01,E^0 = 0$. The solutions corresponding to the rest $\varepsilon$ are obtained by increasing the value of $\varepsilon$ step by step through some proper $\Delta\varepsilon$ and the solution with $\varepsilon_i$ is selected as the initial guess for $\varepsilon_{i+1} = \varepsilon_i+\Delta\varepsilon$. Then, we select much smaller $\Delta\varepsilon$ to compute the switchback-type multi-solution near $\varepsilon = 0.724$. In Fig. 6, we choose the solution with $\varepsilon = 0.722$ (point A) and $\varepsilon = 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 $\varepsilon$ at point $F(\varepsilon\thickapprox0.72524)$, the value of $|T|^2$ will directly jump to the one corresponding to point $H$ and then varies following the route $H\rightarrow E$. Similarly, in the upper branch, decreasing the value of $\varepsilon$ at point $S(\varepsilon\thickapprox0.72376)$ leads to the value of $|T|^2$ follows the route $S\rightarrow G\rightarrow A$. Obviously, $|E|^2$ has three different values when $\varepsilon\in(0.72376,0.72524)$. To simualte the middle branch, we firstly set $\varepsilon = 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$ ($\varepsilon = 0.724226$ in the lower branch) and the solution at point $D$ ($\varepsilon = 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.

4.2.   2D problem

Now, we turn to a 2D problem. Setting $\Omega = [-1 / 2,1 / 2]\times[-1 / 2,1 / 2], \Omega_{0} = [-1 / 4,1 / 4]\times[-1 / 4,1 / 4]$ and $v = 1$ in (6). In this case, we set the parameter $\gamma_x = \gamma_y = 1/2$ in (51)-(52) in our new finite difference scheme. Firstly, we examine the accuracy of the proposed finite difference scheme by taking $\varepsilon = k_0^{-2}$ and the exact solution [29]

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.

To simulate the transmission and collision of the nonparaxial solitons which are also considered in [3,29], we solve the NLH equation (5) with two different incident waves

And the source term is set as

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_{\text {inc }}^1$ comes into $\Omega$ 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 $\varepsilon = k_0^{-2}$ (see Fig. 8(B)), while the field is totally different when $\varepsilon = 0$ (see Fig. 8(C)). When there are two solitons with the field $E_{\text {inc }}^2$ come into the medium from south and east respectively, they meet and collide in $\Omega_0$. In the case of $\varepsilon = k_0^{-2}$, 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 $\Omega_0$ and $\Omega\setminus\Omega_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 $\varepsilon = 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.