On nonparaxial nonlinear Schr\"odinger-type equations

TL;DR

This paper investigates the mathematical properties and numerical methods for the one-dimensional nonparaxial nonlinear Schrödinger equation, an alternative model to the classical equation when paraxial assumptions are invalid.

Contribution

It provides theoretical analysis of well-posedness, Hamiltonian, and multi-symplectic structures, and develops structure-preserving numerical schemes for the nonparaxial nonlinear Schrödinger equation.

Findings

Established linear well-posedness of the equation

Derived Hamiltonian and multi-symplectic formulations

Developed numerical methods that preserve these structures

Abstract

In this paper the one-dimensional nonparaxial nonlinear Schr\"odinger equation is considered. This was proposed as an alternative to the classical nonlinear Schr\"odinger equation in those situations where the assumption of paraxiality may fail. The paper contributes to the mathematical properties of the equation in a two-fold way. First, some theoretical results on linear well-posedness, Hamiltonian and multi-symplectic formulations are derived. Then we propose to take into account these properties in order to deal with the numerical approximation. In this sense, different numerical procedures that preserve the Hamiltonian and multi-symplectic structures are discussed and illustrated with numerical experiments.

Equations

$$\text{textkappa}{u}_{tt}+i{u}_t+\text{textbeta}{u}_{xx}+f(u)=0,$$

$$\text{textkappa}{u}_{tt}+i{u}_t+\text{textbeta}\Delta u+f(u)=0,$$

$$u(x,t)=A{e}^{i(kx+\text{textomega}t)},$$

$$\text{textkappa}{\text{textomega}}^2+\text{textomega}+\text{textbeta}{k}^2-A^2=0,$$

$$u(x,t)=\rho (x-Vt,\text{texteta},V,\text{textkappa})\exp \left(i\text{texttheta}(x-Vt,\text{texteta},V,\text{textkappa},t)\right),$$

$$u_t=\text{mathds}J\frac{\delta \mathcal{H}}{\delta u}$$

$$\text{mathds}K{\mathbf{z}}_t+\text{mathds}M{\mathbf{z}}_x={\nabla }_{\mathbf{z}}\mathfrak{S}(\mathbf{z}),$$

$$k(U,V)\stackrel{\mathrm{def}}{:=}(\text{mathds}MU{)}^{T}V,\text{textomega}(U,V)\stackrel{\mathrm{def}}{:=}(\text{mathds}KU{)}^{T}VU,V\in {\mathbb{R}}^d.$$

$$\text{mathds}K{Z}_t+\text{mathds}M{Z}_x={\mathfrak{S}}^{\prime \prime }(\mathbf{z})Z,$$

$${\partial }_t\text{textomega}(U,V)+{\partial }_{x}k(U,V)=0,$$

$${\partial }_t\left(\mathrm{d}\mathbf{z}\wedge (\text{mathds}K\cdot \mathrm{d}\mathbf{z})\right)+{\partial }_x\left(\mathrm{d}\mathbf{z}\wedge (\text{mathds}M\cdot \mathrm{d}\mathbf{z})\right)=0,$$

$${\mathfrak{E}}_t+{\mathfrak{F}}_x=0,{\mathfrak{I}}_t+{\mathfrak{M}}_x=0,$$

$$\mathfrak{E}(\mathbf{z})\stackrel{\mathrm{def}}{:=}\mathfrak{S}(\mathbf{z})-\frac{1}{2}k({\mathbf{z}}_x,\mathbf{z}),\mathfrak{F}(\mathbf{z})\stackrel{\mathrm{def}}{:=}\frac{1}{2}k({\mathbf{z}}_t,\mathbf{z}),$$

$$\mathfrak{I}(\mathbf{z})\stackrel{\mathrm{def}}{:=}\frac{1}{2}\text{textomega}({\mathbf{z}}_x,\mathbf{z}),\mathfrak{M}(\mathbf{z})\stackrel{\mathrm{def}}{:=}\mathfrak{S}(\mathbf{z})-\frac{1}{2}\text{textomega}({\mathbf{z}}_t,\mathbf{z}).$$

$$\mathcal{F}f(\xi )=\hat{f}(\xi )\stackrel{\mathrm{def}}{:=}{\int }_{-\infty }^{\infty }{e}^{-ix\xi }f(x)dx,$$

$$f(z)=\frac{\partial V(z)}{\partial \overline{z}},$$

$$\text{textkappa}{u}_{tt}+i{u}_t+\text{textbeta}{u}_{xx}=0,$$

$$u_t=v,\text{textkappa}{v}_t=-iv-\text{textbeta}{u}_{xx}.$$

$$\frac{d}{dt}\left(\begin{array}{c}\hat{u}(\xi ,t)\\ \hat{v}(\xi ,t)\end{array}\right)+A(\xi )\left(\begin{array}{c}\hat{u}(\xi ,t)\\ \hat{v}(\xi ,t)\end{array}\right)=0,A(\xi )=\left(\begin{array}{cc}0& -1\\ \frac{-\text{textbeta}{\xi }^2}{\text{textkappa}}& \frac{i}{\text{textkappa}}\end{array}\right).$$

$$\left(\begin{array}{c}\hat{u}(\xi ,t)\\ \hat{v}(\xi ,t)\end{array}\right)=m(\xi ,t)\left(\begin{array}{c}\widehat{u_0}(\xi )\\ \widehat{v_0}(\xi )\end{array}\right),$$

$$m(\xi ,t)\stackrel{\mathrm{def}}{:=}{e}^{-tA(\xi )}.$$

$${\lambda }_{\pm }(\xi )=\frac{i}{2\text{textkappa}}\pm \frac{1}{2\text{textkappa}}\sqrt{4\text{textbeta}\text{textkappa}{\xi }^2-1}.$$

$$m(\xi ,t)$$

$$M(\xi ,t)$$

$$f⟼m(\xi )\hat{f}(\xi )⟼{\mathcal{F}}^{-1}(m(\xi )\hat{f}(\xi )),$$

$$\underset{\xi \to {\xi }_{\pm }}{\lim }m(\xi ,t)=\left(\begin{array}{cc}(1+t\lambda )e^{-t\lambda }& \lambda t{e}^{-t\lambda }\\ -t\lambda e^{-t\lambda }& (1-t\lambda )e^{-t\lambda }\end{array}\right)=e^{-tA({\xi }_{\pm })}.$$

$$\underset{P(\text{textkappa})}{\underbrace{\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& -1& \text{textkappa}& 0\\ 1& 0& 0& \text{textkappa}\end{array}\right)}}\left(\begin{array}{c}{p}_t\\ q_t\\ {\varphi }_t\\ {\phi }_t\end{array}\right)=\left(\begin{array}{c}\varphi \\ \phi \\ -\text{textbeta}{p}_{xx}-\mathrm{Re}f(p,q)\\ -\text{textbeta}{q}_{xx}-\mathrm{Im}f(p,q)\end{array}\right),$$

$$\left(\begin{array}{c}{p}_t\\ q_t\\ {\varphi }_t\\ {\phi }_t\end{array}\right)$$

$$\text{mathds}J(\text{textkappa})\stackrel{\mathrm{def}}{:=}\left(\begin{array}{cccc}0& 0& -\frac{1}{\text{textkappa}}& 0\\ 0& 0& 0& -\frac{1}{\text{textkappa}}\\ \frac{1}{\text{textkappa}}& 0& 0& -\frac{1}{{\text{textkappa}}^2}\\ 0& \frac{1}{\text{textkappa}}& \frac{1}{{\text{textkappa}}^2}& 0\end{array}\right),$$

$$\mathcal{H}(p,q,\varphi ,\phi )\stackrel{\mathrm{def}}{:=}{\int }_{-\infty }^{\infty }\left(\frac{\text{textbeta}}{2}(p_x^2+q_x^2)-V(p,q)-\frac{\text{textkappa}}{2}({\varphi }^2+{\phi }^2)\right)dx.$$

$$\left(\begin{array}{c}{u}_t\\ v_t\end{array}\right)=\text{mathds}J(\text{textkappa})\delta \mathcal{H}(u,v),$$

Full text

††thanks: Corresponding author

1 Introduction

The present paper is concerned with one-dimensional nonparaxial nonlinear Schrödinger equations (NNLS) of the form

[TABLE]

where $u=u(x,t)$ is a complex-valued function of $x\in\mathfrak{R}$ and $t\geqslant 0$. The parameters \textkappa and \textbeta are positive (with \textkappa, in general, small) and $f$ is a complex-valued function of a complex variable. Equation (1.1) and its two-dimensional version

[TABLE]

(where $\Delta\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$ is the Laplace operator) were proposed in the mathematical modelling of nonlinear optical devises under Kerr-type nonlinear media to generate soliton beams for transmission of information (see e. g. [9] and references therein). In (1.1) $u$ represents the scalar (complex) field envelope of a continuous monochromatic beam in a self-focusing Kerr-type nonlinear medium, governed by $f$, and under a linear diffraction in one transverse direction (with two directions in the case of (1.2)). Some examples of $f$ that appear in the physical applications are:

$f(u)=|u|^{q}u,\;q>0$, [13]. This contains, as particular case, the cubic NNLS equation, [10]. 2. 2.

$f(u)=\text{\textalpha}|u|^{\sigma}u-\gamma|u|^{2\sigma}u,\;\text{\textalpha},\gamma,\sigma>0$. This contains, as particular case, the cubic-quintic NNLS equation, [12]. 3. 3.

$f(u)=\text{\textalpha}_{0}|u|^{2}u+(\text{\textalpha}_{1}+\text{\textalpha}_{2}|u|^{2})H(x)u$, where $H(x)$ is the Heaviside function and $\text{\textalpha}_{j}\geqslant 0,j=1,2,3$, [26]. 4. 4.

$f(u)=\left(\displaystyle\frac{\delta}{4\text{\textkappa}}-\text{\textalpha}|u|^{2}\right)u,\;\text{\textalpha},\delta\geqslant 0$, [26]. 5. 5.

$f(u)=\displaystyle\frac{1}{2}\displaystyle\frac{2+\gamma|u|^{2}}{(1+\gamma|u|^{2})^{2}}|u|^{2}u,\;\gamma>0$, [12].

When $\text{\textkappa}=0$, Equation (1.1) reduces to the classical family of nonlinear Schrödinger equations (NLS), [30]. In the context of nonlinear optics, the NLS equation is used in those experiments under an assumption of paraxiality. This means that the diffraction of the light is allowed to develop structure in only one of the coordinates transverse to the direction of propagation. However, some other observations may go beyond this paraxial approximation and require alternatives of modelling where nonparaxial effects must be taken into account. This is the case of, e. g., the study of ultranarrow or high-intensity beams (required in the miniaturization of Information Technology devices) or in the interaction of individual paraxial soliton beams but propagating in different directions which form a significant angle, [10]. An additional, relevant observation in this sense was made by Feit & Fleck, [15], who pointed out that the unphysical catastrophic collapse of self-focusing beams predicted by the paraxial theory is due to the invalidity of this approximation in the neighborhood of a self-focus, see also [21, 1, 28].

In (1.1) and (1.2) the nonparaxial effects are mathematically represented by the inclusion of the second-order time derivative term and its associated parameter \textkappa. This can be expressed in different ways, depending on the above mentioned nonparaxial situations, [11]: a small and positive value of \textkappa may mean that the optical wavelength is a small but non-negligible magnitude when compared to the width of the beam or that the beam is having some degree of spread with respect to the paraxial propagation.

To our knowledge, mathematical properties of (1.1) are known for particular cases of $f$, like some of those mentioned above. We make now a brief summary of them, see the corresponding references for details. The properties mainly concern the existence of conservation laws and special solutions. In the first case, three quantities, the energy-flow, the momentum and the Hamiltonian, are known to be preserved in time by smooth enough solutions $u$ when $f$ is of the form 1 and 2 in the list above, [10, 12, 13]. (This includes the cubic and the cubic-quintic equations.) On the other hand, for the cubic case, plane wave solutions

[TABLE]

will satisfy the dispersion relation

[TABLE]

defining elliptic curves in the $(k,\text{\textomega})$ plane. Furthermore, equation (1.1) also admits soliton-type solutions. This is known for almost all the cases in the list above. The solutions have the form

[TABLE]

for some real-valued functions $\rho,\text{\texttheta}$ depending on \textkappa, the amplitude (\texteta) and the transverse velocity ($V$) parameters, see [10, 12, 13] for the specific form in the corresponding equation. Contrary to the NLS equation, the NNLS is a two-way model and admits solutions (1.3) propagating backward or forward (thus $V$ may be positive or negative). As mentioned in [11] for the cubic case (see also [12]), recovering the solitons of the NLS requires a multiple limit $\text{\textkappa},\text{\textkappa}\text{\texteta}^{2},\text{\textkappa}V^{2}\rightarrow 0$. Finally, it is worth mentioning that a perturbation theory for analyzing the effect of small terms in the self-focusing, cubic NLS equation in critical dimension, developed by Fibich and Papanicolau in [17] (see also [16]), includes the NNLS equation (1.2), studied here as perturbation of the NLS. The prediction of the modulation theory in that case is the formation of decaying focusing-defocusing oscillations, instead of singular solutions, and is in agreement with the observations of Feit & Fleck and others, [15, 1, 28].

The numerical approximation to (1.1) and (1.2) presented in the literature is focused on the cubic NNLS equation and investigates, by computational means, the dynamics of the nonparaxial model, with special emphasis on the description of the self-focusing of the beam, the elimination of backward wave which accompany the propagation of the beam and the evolution of the nonparaxial solitons. As far as the numerical techniques are concerned, the algorithm used by Feit & Fleck, [15], is based on a split-step approach, in the forward in time direction (see [11] for a modified version). On the other hand, Fibich & Tsynkov, [18], introduce a finite difference, fourth-order method, with nonlocal, two-way absorbing boundary conditions (ABC) in the direction of beam propagation, in order to obtain a direct simulation of self-focusing in the nonparaxial case. An improved version, based on introducing Sommerfield-type local radiation boundary conditions in the discretization, was proposed in [19]. The nonparaxial beam propagation method (NBPM), developed by Chamorro et al., [11], is derived by using finite differences leading to an explicit algorithm for the time evolution. The resulting difference-differential equation is computationally solved in the spectral domain with FFT techniques. An efficient parallel implementation of the NBPM can be seen in [25]. Finally, it is also worth mentioning the split-step methods, based on Padé approximation, proposed in [22] for the forward in time equation from (1.2), with time discretization of Crank-Nicolson type.

The present paper contributes to the mathematical analysis of (1.1) and (1.2) in a two-fold way:

Three new (to our knowledge) theoretical properties are presented. We first prove that the initial-value problem (ivp) of (1.1) is linearly well-posed (in the sense of existence and uniqueness of solution). The second property is the Hamiltonian structure of (1.1), which means that under suitable hypotheses on $f$, the NNLS equation can be written in the form

[TABLE]

on a suitable functional space for $u$, where the symplectic structure is given by some matrix operator $\mathds{J}$ and $\dfrac{\delta\mathscr{H}}{\delta\hskip 0.55603ptu}$ stands for the Fréchet derivative of some Hamiltonian function $\mathscr{H}$. The Hamiltonian structure is a property shared by many partial differential equations (PDEs) which appear in the mathematical modelling, including the NLS equation, [30]. One of the consequences of the Hamiltonian formulation is the time conservation of the Hamiltonian $\mathscr{H}$ by the solutions, and the functional derived in this paper generalizes those obtained for particular cases of $f$, [12]. Additionally, the other two invariants, the energy-flow and the momentum, are also generalized and associated to symmetry groups of (1.1). Finally, the Hamiltonian structure can also be extended to the two-dimensional version (1.2).

A third theoretical property studied in the present paper is the formulation of (1.1) (and its two-dimensional version (1.2)) as multi-symplectic. We recall that a system of PDEs is said to be multi-symplectic (MS) in one dimension if it can be written in the form, [4]

[TABLE]

where $\boldsymbol{z}\,(x,\,t)\,:\ \mathds{R}\times\mathds{R}^{\,+}\ \longmapsto\ \mathds{R}^{\,d}\,,d\geqslant 3$, $\mathds{K}$ and $\mathds{M}$ are real, skew-symmetric $d\times d$ matrices, $\boldsymbol{\nabla}_{\,\boldsymbol{z}}$ is the gradient operator in $\mathds{R}^{\,d}\,$ and the potential $\mathfrak{S}\,(\boldsymbol{z})\,$ is assumed to be a smooth function of $\boldsymbol{z}$. The MS theory generalizes the Hamiltonian formulation in the sense of the presence of a symplectic structure with respect to each of the space and time variables. These structures are respectively defined by the two-forms, [4, 6]

[TABLE]

Then, the multi-symplectic property of (1.4) means that if $U$ and $V$ are solutions of the corresponding variational equation

[TABLE]

(where $\mathfrak{S}^{\prime\prime}(\boldsymbol{z})$ stands for the Hessian matrix of $\mathfrak{S}(\boldsymbol{z})$) then

[TABLE]

The two-forms \textomega and $k$ can be written in terms of the differentials $d\boldsymbol{z},\boldsymbol{z}\in\mathbb{R}^{d}$ in such a way that the MS conservation law (1.5) is alternatively expressed as

[TABLE]

with $\wedge$ being the standard exterior product of differential forms, [29, 23].

The symplecticity must be understood locally, as the forms vary in space and time. This local character also affects the preservation of quantities in MS systems (1.4); specifically, when the function $\mathfrak{S}\,(\boldsymbol{z})$ does not depend explicitly on $x$ or $t\,$, then local energy and momentum conservation laws are satisfied:

[TABLE]

where, [6]

[TABLE]

As in the case of the Hamiltonian structure, the MS formulation also holds in many PDEs used in modelling, including the NLS equation, [24]. We finally note that any MS PDE system is also Lagrangian, [4]. 2. 2.

The second type of contributions of this paper concerns the numerical approximation to (1.1). Compared to the references in the literature on this subject, commented above, here we adopt a different point of view. In accordance with the theoretical structures of the equation, we are interested in the geometric numerical integration, [20]. As is well known, this approach aims at analyzing qualitative properties of the numerical approximation which may improve the accuracy, beyond the quantitative measure given by the classical order of convergence. These properties come typically from emulating, in a discrete sense, geometric structures of the equation under study and which may have influence on the numerical integration, especially for long term simulations.

In this case, the present paper is focused on the Hamiltonian and MS structure to propose geometric numerical methods to approximate the periodic initial-value problem associated to (1.1). By using the method of lines, the discretization in space is first studied. We observe that the use of a symmetric operator to approximate the second partial derivative in space generates a semi-discrete system with a Hamiltonian structure as in the paraxial case, cf. [8]. If this symmetric character is obtained from the approximation to the first partial derivative with a skew-symmetric operator, then the resulting semi-discrete system is also multi-symplectic, in the sense of the preservation of some discrete MS conservation law, cf. [5]. These properties of the spatial discretization enable us to choose a symplectic time integration with the aim of providing the full discretization with a symplectic,[27], an a multi-symplectic, [6], structure.

The paper is structured as follows. In Section 2, the theoretical results, concerning linear well-posedness, Hamiltonian structure and MS formulation are introduced and proved. Additional conserved quantities and the specific forms of the MS conservation law (1.5) and the local conservation laws (1.6), (1.7) are derived. Section 3 is devoted to the description of geometric numerical methods and their properties when approximating (1.1), mainly focused on the preservation of the Hamiltonian and MS structures, as well as the invariants of the problem. The performance of the geometric approximation is illustrated in some numerical experiments by taking a full discretization based on the Fourier pseudospectral collocation method in space along with the symplectic time integration given by the implicit midpoint rule. The experiments involve soliton simulations and the evolution of errors in discrete versions of the conserved quantities. Conclusions are summarized in Section 4.

The following notation will be used throughout the paper. We will alternatively use the complex form (1.1) and its equivalent formulation as a real system for the real and imaginary parts of $u$ and $u_{t}$. On the other hand, $H^{s}=H^{s}(\mathbb{R}),s\geqslant 0$ will stand for the $L^{2}-$based Sobolev space of order $s$, with $H^{0}=L^{2}$. The Fourier transform of an integrable function $f$ is defined as

[TABLE]

with inverse operator denoted by $\mathcal{F}^{-1}$. The transform is extended to $f\in L^{2}$ by using density arguments in the usual way.

Finally, $\langle\cdot,\cdot\rangle_{n}$ will denote the Euclidean inner product in some $\mathbb{R}^{n}$, where $n\geqslant 1$ will take different values throughout Sections 2 and  3.

2 Some theoretical properties of the NNLS equation

In this section we will assume that the function $f$ in (1.1) satisfies

[TABLE]

(where $\overline{z}$ denotes complex conjugate) for some smooth, real-valued potential $V$.

2.1 Linear well-posedness

We first study well-posedness of the linearized equation associated to (1.1), that is

[TABLE]

or, as a first-order system

[TABLE]

By well-posedness we mean existence and uniqueness of solutions and continuous dependence on the initial data in the corresponding spaces. We adopt the strategy considered in, e. g. [3], based on the Fourier transform of (2.2) and the representation of the solutions of the resulting system. Taking the Fourier transform with respect to $x$ in (2.2) we have

[TABLE]

The solution of the ivp for (2.2) with initial data $u(x,0)=u_{0}(x),v(x,0)=v_{0}(x)$ can be written, in the Fourier space, as

[TABLE]

where $\widehat{u_{0}}(\xi),\widehat{v_{0}}(\xi)$ are the Fourier transforms of $u_{0},v_{0}$, respectively, and $m(\xi,t)$ is the Fourier multiplier

[TABLE]

We now study the structure of (2.3). Note that the eigenvalues of $A(\xi)$ are

[TABLE]

It holds that $\lambda_{+}(\xi)\neq\lambda_{-}(\xi)$ except when $4\text{\textbeta}\text{\textkappa}\xi^{2}-1=0$, that is, when $\xi=\xi_{\pm}=\pm 1/2\sqrt{\text{\textbeta}\text{\textkappa}}$ and for which $\lambda\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\lambda_{+}=\lambda_{-}=i/2\text{\textkappa}$. According to the corresponding spectral decomposition of $A(\xi)$, the Fourier multiplier (2.3) can be written, for $\xi\neq\xi_{\pm}$, as

[TABLE]

As observed in [3], if $m$ is unbounded at finite values of $\xi$, then the linear ivp cannot be well-posed in any of the Sobolev spaces $H^{s},s\geqslant 0$ because the operators

[TABLE]

(where $\mathcal{F}^{-1}$ denotes the inverse Fourier transform) are not bounded maps from $H^{s}$ to $L^{2}$. In our case, the only possible poles of $m$ are precisely given by $\xi=\xi_{\pm}$. A tedious but direct computation shows that

[TABLE]

Therefore $m(\xi,t)$ is bounded on bounded intervals and linear well-posedness in $L^{2}$ based Sobolev spaces follows, [3].

2.2 Hamiltonian structure, conserved quantities and symmetry groups

The second theoretical property considered in this section concerns the extension of the conservation laws, derived for some particular equations of (1.1), [12], to the more general case of $f$ satisfying (2.1). Note that, as a first-order real system, (1.1) has the form

[TABLE]

where $u=p+iq,u_{t}=\phi+i\varphi$. For $\text{\textkappa}>0$ we invert the matrix $P(\text{\textkappa})$ to write (2.4) as

[TABLE]

where $\delta\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\left(\frac{\delta}{\delta p},\frac{\delta}{\delta q},\frac{\delta}{\delta\phi},\frac{\delta}{\delta\varphi}\right)$,

[TABLE]

and, using (2.1),

[TABLE]

Then (2.5) gives the Hamiltonian formulation of (1.1) with structure matrix given by (2.6) and Hamiltonian (2.7). In complex form, with $v=u_{t}$, (2.5) is of the form

[TABLE]

with

[TABLE]

The Hamiltonian formulation can be extended to the two-dimensional case with $\mathds{J}$ as in (2.8) and

[TABLE]

Besides the Hamiltonian (2.7), some additional conserved quantites can also be derived from the symmetry groups of (1.1). The co-symplectic matrix (2.6) defines the Poisson bracket of two functionals $F(p,q,\phi,\varphi)$ and $G(p,q,\phi,\varphi)$, [23]

[TABLE]

Note now that the functional

[TABLE]

satisfies $\{\mathcal{I}_{2},\mathscr{H}\}=0$, being $\mathscr{H}$ the Hamiltonian (2.7). This implies, [23], that $\mathcal{I}_{2}$ is invariant by the solutions of (2.4). In complex form, (2.9) reads

[TABLE]

It is not hard to see the connection between (2.9) and the symmetry group of (2.4) consisting of spatial translations

[TABLE]

since the infinitesimal generator of (2.11) is

[TABLE]

The invariants $\mathscr{H}$ and $\mathcal{I}_{2}$ were obtained for particular cases of $f$ in, e. g. [10, 13, 12] and, in this sense, (2.7) and (2.10) extend the existence of the Hamiltonian and momentum to the general case of (1.1) with $f$ satisfying (2.1). When $\text{\textkappa}=0$, (2.7) and (2.9) correspond to the formulas obtained in [14]. For the preservation of an analogous quantity to the mass of the paraxial case (called energy-flow in [12]) additional hypotheses on $f$ are required, [14].

Lemma 2.1

Assume that $f:\mathbb{C}\rightarrow\mathbb{C}$ satisfies

[TABLE]

Then $f(z)=zg(|z|)$ for some real-valued function $g$.

Proof. Note first that (2.12) implies that $f(z)$ is real when $z$ is real. On the other hand, if $z\in\mathbb{C},z\neq 0$, using (2.13) we have

[TABLE]

Therefore

[TABLE]

where $g(|z|)\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\displaystyle\frac{1}{|z|}f\left(|z|\right)$ is real. Finally, (2.13) also implies that $f(-z)=-f(z),z\in\mathbb{C}$. In particular, $f(0)=0$, which completes the proof.$\Box$

Note that, under the hypotheses (2.12), (2.13), the functional

[TABLE]

satisfies $\{\mathcal{I}_{1},\mathscr{H}\}=0$ and, therefore, $\mathcal{I}_{1}$ is another conserved quantity of (1.1). In complex form, (2.9) is

[TABLE]

(For the particular cases of $f$ considered in [10, 12], an equivalent expression for (2.15) is derived.) In this case, the symmetry group of (1.1) associated to (2.14) consists of rotations

[TABLE]

and its infinitesimal generator can be written as

[TABLE]

For the paraxial case $\text{\textkappa}=0$, the corresponding quantity (2.14) can be seen in [14].

Note that, in order for (1.1) to admit (2.16) as symmetry group, only the condition (2.13) is needed; but the connection with (2.14) as conserved quantity additionally requires to assume (2.12).

The relation of the quantities $\mathcal{I}_{1},\mathcal{I}_{2}$ with the symmetry groups of (1.1) motivates the study of solitary-wave solutions for the general case of $f$ satisfying (2.1), (2.12)-(2.13). They can be found as relative equilibria (see [14] and references therein), that is, equilibria $(p_{0},q_{0},\phi_{0},\varphi_{0})$ of the Hamiltonian on fixed level sets of the first two invariants

[TABLE]

for real multipliers $\text{\textmugreek}_{1},\text{\textmugreek}_{2}$ and real $c_{1},c_{2}$ determining the fixed level sets of $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ respectively. After some computations and if $u_{0}=p_{0}+iq_{0}$, the first equation in (2.17) reads

[TABLE]

Once $u_{0}$ is obtained from (2.18), the solitary-wave solution of the ivp of (1.1) with initial conditions $u(x,0)=u_{0}(x),u_{t}(x,0)=i\text{\textmugreek}_{1}u_{0}(x)+\text{\textmugreek}_{2}u_{0}^{\prime}(x)$ is derived from this profile by a coupled rotation and translation determined by the multipliers $\text{\textmugreek}_{1}$ and $\text{\textmugreek}_{2}$ respectively, that is

[TABLE]

The resolution of (2.17) involves to obtain the corresponding relations between the level set values $c_{j}$ and the multipliers $\text{\textmugreek}_{j},j=1,2$. Equation (2.18) can be explicitly solved for some functions $f$ in a similar way, e g., to that of [14] for the NLS equation. This is the case, for example, of the explicit formulas derived in [10, 12, 13]. The existence of solutions $u_{0}$ of (2.18) for a more general term $f$ is, to our knowledge, an open question. In this sense, classical techniques of numerical generation, [31], may serve as a first approach.

2.3 Multi-symplectic structure

In this section we derive the MS structure of (1.1). We define the variables $v,w$ such that $p_{x}=v,q_{x}=w$. Then (1.1) can be written as a system

[TABLE]

Now if $\boldsymbol{z}=(p,q,v,w,\phi,\varphi)^{T}\in\mathds{R}^{6}$, consider the vector field $F:\mathds{R}^{6}\rightarrow\mathds{R}^{6}$ given by the right hand side of (2.19)

[TABLE]

Note that the Jacobian $F^{\prime}(\boldsymbol{z})$ is symmetric for all $\boldsymbol{z}$ and therefore Poincaré’s lemma implies that $F$ is conservative, that is $F=\boldsymbol{\nabla}_{\boldsymbol{z}}\,\mathfrak{S}$ for some potential $\mathfrak{S}\,$. Thus system (2.19) can be written in the form (1.4) in $\mathds{R}^{6}$ with

[TABLE]

leading to the MS formulation of (1.1). A potential is given by

[TABLE]

In this case, the conservation law of multi-symplecticity has the form (1.5) with

[TABLE]

Similarly, the local energy and momentum conservation laws have the form (1.6), (1.7) with

[TABLE]

The MS formulation can be extended to the two-dimensional version (1.2). As before, we write $u=p+iq$ and define the variables $\phi,\varphi,v_{1},w_{1},v_{2},w_{2}$ such that

[TABLE]

Then, in terms of $\boldsymbol{z}=(p,q,v_{1},w_{1},v_{2},w_{2},\phi,\varphi)^{T}\in\mathfrak{R}^{8}$, equation (1.2) admits a MS formulation in 2D (cf. [4, 5])

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

In this case, the MS conservation law has the form

[TABLE]

where

[TABLE]

Corresponding formulas for the local conservation laws can be derived.

3 Numerical discretization of the NNLS equation

In this section we study the numerical approximation to the periodic initial-value problem of (1.1) with $f$ satisfying (2.1), (2.12), (2.13). As mentioned in the introduction, the purpose here is searching for geometric discretizations emulating qualitative properties of the continuous NNLS equations, mainly focused on the Hamiltonian and MS structures.

First it may be worth mentioning the influence of the imposition of periodic boundary conditions on these two formulations. In the case of the MS structure, note that the symplecticity is understood locally, the conservation law (1.5) does not depend on specific boundary conditions. However, by integrating (1.6) on a one-period interval $(a,b)$ (with period $L=b-a$) in the spatial domain, periodic boundary conditions imply the preservation of the global energy and momentum

[TABLE]

where $\mathfrak{E},\mathfrak{I}$ are given by (2.23) and (2.24), respectively.

Note also that the periodic initial-value problem of (1.1) retains a Hamiltonian structure, with co-symplectic matrix (2.6) and a Hamiltonian $\mathscr{H}_{L}\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}-\mathcal{E}$ of (3.1). Additionally, $\mathcal{I}$ in (3.1) defines the corresponding version $\mathcal{I}_{2L}\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}-\mathcal{I}$ of (2.9). Finally, it is not hard to see that when $f$ satisfies (2.12), (2.13), then the functional

[TABLE]

analogous to (2.14), is preserved in time by the solutions. Note finally that by integrating (1.5), with \textomega and $k$ given by (2.22), on $(a,b)$ and using the periodic boundary conditions, then

[TABLE]

is constant in time.

3.1 Semi-discretization in space

By using the method of lines, we firstly discretize in space. For that, we consider an integer $N\geqslant 1$. On a uniform grid of $N$ points $x_{j},j=0,\ldots,N-1$, on an interval $(a,b)$ and stepsize $h=(b-a)/N$, the second partial derivative in space is approximated by some grid operator $A_{h}$ in such a way that the following semi-discrete second-order equation holds

[TABLE]

where $u_{h}(t)$ is a complex, $N$-vector approximating the exact solution $u(\cdot,t)$ at the grid values and $f(u_{h})$ is the vector whose components are obtained by evaluating $f$ at the components of $u_{h}$. With a similar proof to that of Theorem 6.1 in [8], it can be seen that, whenever $A_{h}$ is a symmetric matrix, (3.4) admits a Hamiltonian formulation

[TABLE]

with $\mathds{J}(\text{\textkappa})$ as in (2.6), $u_{h}=p_{h}+iq_{h},u_{h}^{\prime}=\phi_{h}+i\varphi_{h}$ and

[TABLE]

where $1_{N}$ denotes the vector of length $N$ with all its components equal to one and $V(p_{h},q_{h})$ is obtained by evaluating $V$ at each component of $(p_{h},q_{h})$. We also notice that $hH_{h}$ is the natural discretization of the Hamiltonian $\mathscr{H}_{L}=-\mathcal{E}$ in (3.1).

The behaviour of the spatial semi-discretization with respect to the invariants (2.9) and (2.14) is as follows. With a similar proof to that of Theorem 6.2 in [8], it can also be checked that, when $A_{h}$ is symmetric and $f(z)$ as in Lemma 2.1, then

[TABLE]

is an invariant of the semidiscrete system (3.5). Note that $hI_{1,h}$ is the natural discretization of $\mathcal{I}_{1L}$ in (3.2).

On the other hand, the behaviour with respect to (2.9) can also be studied in a similar way to that of the paraxial case, see [8]. If in (3.4) we assume that

[TABLE]

then a natural discrete version of $\mathcal{I}_{2L}$ is $hI_{2,h}$ where

[TABLE]

In a similar way to [8], it can be proved that, when $DB_{h}=-B_{h}D$ for the $N\times N$ matrix $D$ which reverses the order of the components of the vector to which it is applied (i.e. $(Dx)_{j}=x_{N-j-1},j=0,\ldots,N-1$) and the initial conditions are symmetric in the space interval of integration (i.e $Dp_{h}(0)=p_{h}(0)$, $Dq_{h}(0)=q_{h}(0)$, $D\varphi_{h}(0)=\varphi_{h}(0)$, $D\phi_{h}(0)=\phi_{h}(0)$), it happens that $I_{2,h}(t)=0$ for every time $t$. On the other hand, for general initial conditions, if $B_{h}$ is a skew-symmetric matrix, under the same hypotheses of Theorem 5.1 in [8], $hI_{2,h}$, where in (3.8) $B_{h}$ is substituted by the pseudospectral differentiation operator, is a quasiinvariant in the sense described in the same theorem.

The last property of the semi-discretization in space considered here concerns the multi-symplecticity. Let ${\boldsymbol{z}}_{h}(t)=({\boldsymbol{z}}_{h,\,j}(t))_{j=0}^{N-1}\in\mathbb{R}^{6N}$, where ${\boldsymbol{z}}_{h,\,j}(t)$ is an approximation to

[TABLE]

(with ${\boldsymbol{z}}_{h,\,j}(0)=\boldsymbol{z}(x_{j},0)$). If $A_{h}$ satisfies (3.7), then (1.4), (2.20), (2.21) can be discretized in the form

[TABLE]

for $j=0,\ldots,N-1$, where $C_{h}\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}B_{h}\otimes I_{6}$, with $I_{6}$ the $6\times 6$ identity matrix and $\otimes$ standing for the Kronecker product of matrices. In terms of ${\boldsymbol{z}}_{h}(t)$, (3.9) can be written as

[TABLE]

where $I_{N}$ denotes the $N\times N$ identity matrix, $\boldsymbol{\nabla}_{\,{\boldsymbol{z}}}\,\mathfrak{S}\,({\boldsymbol{z}}_{h})$ stands for the vector of length $6N$ matrix with components $\boldsymbol{\nabla}_{\,{\boldsymbol{z}}}\,\mathfrak{S}\,({\boldsymbol{z}}_{h\,j})\in\mathbb{R}^{6},j=0,\ldots,N-1$. By using some properties of the Kronecker product, we have $(I_{N}\otimes\mathds{M})C_{h}=B_{h}\otimes\mathds{M}$ and therefore (3.10) reads

[TABLE]

The approximation is multi-symplectic in the following sense: let $U,V\in\mathbb{R}^{6N}$ be solutions of the variational equation associated to (3.11),

[TABLE]

where $\mathfrak{S}^{\prime\prime}\,({\boldsymbol{z}}_{h})$ is the block diagonal matrix with $6\times 6$ blocks $\mathfrak{S}^{\prime\prime}\,({\boldsymbol{z}}_{h,j}),j=0,\ldots,N-1$ and the $j$-th component of $\mathfrak{S}^{\prime\prime}\,({\boldsymbol{z}}_{h})Z$ is given by $\mathfrak{S}^{\prime\prime}\,({\boldsymbol{z}}_{h,j})Z_{j}$, being $Z=\{Z_{j}\}_{j=0}^{N-1},Z_{j}\in\mathbb{R}^{6}$. Then, we have

[TABLE]

Using the property $(B_{h}\otimes\mathds{M})^{T}=(B_{h}^{T}\otimes\mathds{M}^{T})$ and the skew-symmetry of $\mathds{M}$, then (3.12) reads

[TABLE]

which represents the discrete MS conservation law preserved by (3.11). We finally note that

[TABLE]

Therefore, if $B_{h}$ is skew-symmetric, then (3.13) leads to the conservation of total symplecticity in time,

[TABLE]

cf. [5]. Note that $h\displaystyle\sum_{j=0}^{N-1}\text{\textomega}_{\,j}$ is the natural discretization of $\widetilde{\text{\textomega}}$ in (3.3).

3.2 Full discretization

The Hamiltonian structure of the semi-discrete system (3.5) suggests to use symplectic methods in order to preserve the geometric character of the full discretization. Note first that, although the integrability of (3.5) (and of (1.1)) is, to our knowledge, not known, a good behaviour with respect to the preservation of the discrete Hamiltonian (3.1) is expected, at least when approximating soliton-type solutions, [14]. On the other hand, as far as (3.6) is concerned, observe that this is a quadratic invariant associated to the symmetric matrix

[TABLE]

Then, with a similar proof to that of Theorem 6.3 in [8], we have

[TABLE]

where $(p_{h}^{n},q_{h}^{n},\phi_{h}^{n},\varphi_{h}^{n})$ denotes the approximation to $(p_{h}(t_{n}),q_{h}(t_{n}),\phi_{h}(t_{n}),\varphi_{h}(t_{n}))$ at $t_{n}=n\Delta t$ (with time step $\Delta t$) given by a symplectic Runge-Kutta method. We also notice that, after full discretization, $\mathcal{I}_{1L}$ is also conserved.

In the case of the quadratic quantity (3.8), the associated matrix is

[TABLE]

When $B_{h}$ is skew-symmetric, $S_{2,h}$ is symmetric and the invariant would also be conserved by integrating (3.5) in time with a symplectic Runge-Kutta method. In particular, this holds when $B_{h}$ is the pseudospectral differentiation operator.

Similarly, when (3.9) is discretized in time with a symplectic method, we obtain a fully discrete scheme which preserves by construction a discrete version of (3.13) and, consequently, is multi-symplectic, [6, 24].

3.3 Numerical experiments

In order to illustrate the previous results, the periodic initial-value problem of (1.1) on a long enough interval $(a,b)$ with $\text{\textbeta}=1/2$ and $f(u)=|u|^{2}u$ was numerically integrated to approximate the soliton-type solution of the form (1.3) given by

[TABLE]

The spatial discretization was performed with the Fourier pseudospectral collocation method, so that $A_{h}$ in (3.4) corresponds to the evaluation at the collocation points of the second derivative of the trigonometric interpolant polynomial based on the nodal values of the semidiscrete approximation $u_{h}$. It is well known that $A_{h}$ satisfies (3.7) with

[TABLE]

where $\text{\texttheta}_{N}=e^{-\frac{2\pi i}{N}}$. Therefore, $B_{h}$ is skew-symmetric and $A_{h}$ is symmetric. Using the fact that the Discrete Fourier Transform diagonalizes $A_{h}$, the pseudospectral method is in practice implemented in Fourier space for the discrete Fourier coefficients of $u_{h}$.

For the time integration, we have chosen the implicit midpoint rule, which is a symplectic Runge-Kutta method, and therefore the conservation of quadratic invariants is expected. On the other hand, the formulation of the fully discrete method as multi-symplectic (a property which is a consequence of the study developed in Sections 3.1 and 3.2) can be used to derive a numerical dispersion relation in a more direct way, [2, 6, 7]. Let $\boldsymbol{z}_{i}^{n}$ be the numerical approximation at time $t_{n}$ to the value of $\boldsymbol{z}$ at the $i$-th grid point and define the operators on $\mathbb{R}^{6}$

[TABLE]

and $i=0,\ldots,N-1$. When (3.9) is discretized in time with the implicit midpoint rule we obtain

[TABLE]

where $\mathcal{D}_{t},\mathcal{M}_{t}$ are given by (3.15), $\boldsymbol{z}^{n}=(\boldsymbol{z}_{j}^{n})_{j=0}^{N-1}$ and $z_{j}^{n}$ is an approximation to $\boldsymbol{z}_{h,j}(t_{n})$. Note that since $\mathcal{M}_{t}$ is an operator on $\mathbb{R}^{6}$ and $(I_{N}\otimes\mathcal{M}_{t})(B_{h}\otimes I_{6})=B_{h}\otimes\mathcal{M}_{t}$, then system (3.16) can be solved iteratively for $\boldsymbol{z}^{n+1/2}\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}(I_{N}\otimes\mathcal{M}_{t})\boldsymbol{z}^{n}$ by the fixed point algorithm

[TABLE]

for $\text{\textnu}=0,1,\ldots$, where

[TABLE]

with the linear ($\mathcal{L}$) and the nonlinear ($\mathds{N}$) part of the gradient $\boldsymbol{\nabla}_{\,\boldsymbol{z}}\,\mathfrak{S}\,({\boldsymbol{z}})$ and, for $Z=\{Z_{j}\}_{j=0}^{N-1},Z_{j}\in\mathbb{R}^{6},\left(\mathcal{N}(Z)\right)_{j}=\mathds{N}(Z_{j}),j=0,\ldots,N-1$. For the case of (1.4), (2.20), (2.21) and in terms of the discrete Fourier coefficients of $z^{n+1/2}=(p,q,v,w,\phi,\varphi)$, the fixed-point system (3.16) will have the form

[TABLE]

for $-N/2\leqslant m\leqslant N/2-1$, where $\text{\textmugreek}_{m}\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}im\frac{2\pi}{b-a}$, $\widehat{A}(m),\widehat{B}(m)$ are the $m$-th discrete Fourier coefficients of

[TABLE]

respectively and

[TABLE]

as the approximation at $t_{n+1}$ from the resolution of (3.18) by the corresponding algorithm (3.17). Note that the contribution of the variables $v$ and $w$ is given by

[TABLE]

and this is used to simplify the rest of the equations leading to the final form (3.18). Its complex version is

[TABLE]

for $-N/2\leqslant m\leqslant N/2-1$, where $u=p+iq,\chi=\phi+i\varphi$ and $u^{n+1}=2u-u^{n},\chi^{n+1}=2\chi-\chi^{n}$.

The previous formulation is now used to derive a dispersion relation for the numerical method. Observe that the linear dispersion relation of (1.4), (2.20), (2.21) is

[TABLE]

Note that using the operators defined in (3.15), the scheme (3.16) can be written in Fourier space as

[TABLE]

Multiplying the equations 1 and 2 by $\mathcal{M}_{t}$ and using equations 3 to 6, we obtain

[TABLE]

or, in complex form ($u=p+iq$)

[TABLE]

for $-N/2\leqslant m\leqslant N/2-1$. Then, substituting $u_{j}^{n}=e^{i(j\xi+n\widetilde{\text{\textomega}})}$ into (3.20) leads to the numerical dispersion relation

[TABLE]

where $\mathfrak{D}$ is given by (3.19) and

[TABLE]

Equation (3.21) shows the approximate preservation of the linear dispersion relation (3.19) by (3.20).

For the nonlinear case, assuming that $f$ satisfies (2.12), (2.13), the dispersion relation is

[TABLE]

where $f(z)=zg(z),g$ as in Lemma 2.1. Now, the inclusion of the operator $\mathcal{M}_{t}$ in the nonlinear term of (3.20) leads to the numerical dispersion relation of the form

[TABLE]

The following numerical experiments study the accuracy and geometric properties of the scheme when approximating (3.14). As parameter values, we have taken $\text{\texteta}=1$, $V=1$, $\text{\textkappa}=10^{-4}$ and, as space interval of integration, $[a,b]=[-150,50]$. The final time of integration is $T=100$.

We first check the order of convergence. Figures 1 and 2 show the relative errors in the discrete $L^{2}$-norm of the solution and the derivative respectively when $N=1000$ nodes have been considered in space, and the stepsizes $\Delta t=1/5,1/10,1/20,1/40$ have been taken in time, apart from the tolerance $10^{-13}$ for the fixed-point iteration associated to the midpoint rule. The spectral accuracy of the pseudospectral approximation and the regularity of (3.14) make the error in space negligible and it can be observed that the error is divided by $4$ when the time-stepsize is halved, as it corresponds to the second order of the midpoint rule. Moreover, the growth of error with time is at most linear in both the solution and the derivative, as it also happened when integrating solution solutions of the paraxial NLS equation with the same integrators, [14].

The behaviour of the method with respect to the quantities $I_{1,h}$, $I_{2,h}$ and $H_{h}$ is illustrated in Figures 3, 4 and 5, respectively. They display the evolution of the relative error between the values of the corresponding quantity at the numerical solution and that of the (exact) initial condition. Figures 3 and 4 show that the errors in $I_{1,h}$ and $I_{2,h}$ are of the size of the tolerance and double precision round-off, which corroborates the previous results of preservation when using a symplectic method. As for the Hamiltonian, Figure 5 shows that the relative error keeps small and, up to the final time of integration, does not grow with time. Note also that, when the time-stepsize $\Delta t$ is halved, the error is divided by $16$. This fourth order of the error in the Hamiltonian may be explained, as in the paraxial case (see [14]) by using (2.17). The fact that, at the soliton-type solution (3.14), the variational derivative of the Hamiltonian is a linear combination of the variational derivatives of $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ implies that, in first approximation, the error in the natural discretization of $\mathscr{H}$ is also negligible and just the error corresponding to $O(\Delta t^{4})$ can be observed.

4 Concluding remarks

In the present paper we consider nonlinear Schrödinger equations of nonparaxial type (NNLS), proposed as an alternative to the NLS approach in those models with non-negligible nonparaxial effects. The formulation considered here attempts covering different versions of the NNLS equation presented in the literature, like the cubic, the cubic-quintic and other cases, [10, 12, 13, 26]. The paper contributes to several mathematical properties of the equations. We first establish linear well-posedness results (existence and uniqueness of solution of the linerized initial-value problem in suitable Sobolev spaces). Then the Hamiltonian structure and two additional conservation laws are derived, extending the results obtained in [12], relating the invariants with symmetry groups of the equations and with the existence and generation of solitary-wave solutions as relative equilibria. A third theoretical property proved in this paper is the derivation of a multi-symplectic formulation, meaning the existence of a symplectic structure with respect to both the space and time variables.

The present paper finally studies the preservation of some of these properties in the numerical approximation to the NNLS equations. We establish conditions on the spatial and time discretizations in order for the resulting schemes to preserve discrete versions of the continuous invariants as well as the MS structure. These results are compared with those of the paraxial case and illustrated with some numerical experiments for the cubic NNLS by using a Fourier pseudospectral discretization in space and the implicit midpoint rule as time integrator.

Several features motivate to extend this research to the two-dimensional version of the NNLS equations. The first are suggested in some remarks of the present paper and concern the extension of the theoretical results (mainly the Hamiltonian and the MS structures) to the bi-dimensional case. An additional challenge, specially from the viewpoint of the numerical approximation, is the presence of new phenomena in 2D, like singularity formation in the NLS and its relation with the inclusion of nonparaxiality in the new models.

Acknowledgments

This work has been supported by Ministerio de Ciencia, Innovación y Universidades, FEDER and Junta de Castilla y León through projects MTM 2015-66837-P, VA024P17, VA041P17 and VA105G18.