Numerical solution of the nonlinear Helmholtz equation using nonorthogonal expansions (Q2566671)

The paper aims to develop an advanced simulation algorithm for nonlinear wave transmission within the framework of a model based on the nonlinear Helmholtz (NLH) equation for a complex amplitude \(E(z,x,y)\), \[ (\partial^2_{zz}+\Delta_{\perp})E + k_0^2(1 + \varepsilon | E| ^{2\sigma})E=0, \] where \(\Delta_{\perp} = \partial_{xx}^2 + \partial_{yy}^2\) is the transverse diffraction operator, while the longitudinal coordinate \(z\) plays the role of the evolutional variable (the wave propagation is considered in the spatial domain, without any time dependence). The difference from the usually employed simplified propagation equation (in the spatial domain) of the nonlinear Schrödinger type is that the NLH equation makes it possible to explicitly consider the effect of nonlinear backscattering on nonuniformities of the field density, \(| E| ^2\), and does not rely upon the paraxial approximation (weak diffraction for very broad beams). The equation is to be solved within a finite interval of \(z\), which corresponds to the physically meaningful problem of the transmission in a finite-thickness slab. The physically relevant value of the nonlinear index is \(\sigma =1\), which corresponds to the Kerr effect in optics. The paper treats this case, as well as another one, with \(\sigma=2\) (quintic nonlinearity), with one transverse coordinate \(x\), rather than \(x\) and \(y\) (i.e., \(\Delta_{\perp}=\partial^2_x\)). The latter case, although not directly relevant to nonlinear optics, is interesting, as it corresponds to the critical situation (in the sense of the possibility of weak collapse) in the one-dimensional nonlinear-propagation problem. The former case relates to the critical situation in the propagation setting with two transverse dimensions. In simulations of the NLH equation, a crucially important role is played by lateral boundary conditions set at fixed values of \(x\), a significant issue being correct handling of outcoming and incoming waves (which must avoid unphysical reflections). In a previous paper of the same authors [J. Comput. Phys. 171, 632--677 (2001; Zbl 0991.78016)], a method of artificial boundary conditions was developed for this purpose. In that paper, zero Dirichlet conditions were fixed at the artificial boundaries. In the present paper, a major improvement is replacement of the latter conditions by the Sommerfeld radiation conditions. The implementation of the latter type of the boundary conditions amounts to computation of eigenvalues and eigenvectors of a non-Hermitian matrix. Accordingly, the separation of variables, which is a necessary ingredient of the numerical integration code, must be performed in terms of an expansion over a basis of non-orthogonal eigenvectors. The algorithm is applied to the simulations of both critical field configurations (ones at the collapse threshold) and propagation of subcritical pulses (solitons).