Bifurcation of Gap Solitons in Coupled Mode Equations in $d$ Dimensions

TL;DR

This paper proves the existence of standing gap solitons in coupled mode equations in multiple dimensions, using bifurcation theory and spectral analysis, supported by a numerical example in two dimensions.

Contribution

It establishes the bifurcation of gap solitons from zero solutions in coupled mode equations under spectral gap conditions, extending previous results to higher dimensions.

Findings

Existence of standing gap solitons in $ ext{d}$-dimensional coupled mode equations.

Reduction to a perturbed nonlinear Schrödinger equation.

Numerical demonstration of gap solitons in 2D.

Abstract

We consider a system of first order coupled mode equations in $\mathbb{R}^d$ describing the envelopes of wavepackets in nonlinear periodic media. Under the assumptions of a spectral gap and a generic assumption on the dispersion relation at the spectral edge, we prove the bifurcation of standing gap solitons of the coupled mode equations from the zero solution. The proof is based on a Lyapunov-Schmidt decomposition in Fourier variables and a nested Banach fixed point argument. The reduced bifurcation equation is a perturbed stationary nonlinear Schr\"odinger equation. The existence of solitary waves follows in a symmetric subspace thanks to a spectral stability result. A numerical example of gap solitons in $\mathbb{R}^2$ is provided.

Equations

$$\mathrm{i}({\partial }_t{A}_j+v_g^{(j)}\cdot \nabla A_j)+\sum _{r=1}^N{\kappa }_{jr}{A}_r+N_j(\vec{A})=0,j=1,\dots ,N,$$

$$N_j(\vec{A})$$

$${\gamma }_j^{(m,n,o)}$$

$$\vec{A}(x,t)=e^{-\mathrm{i}\omega t}\vec{B}(x),|\vec{B}(x)|\to 0\text{ as }|x|\to \infty ,$$

$$\omega \vec{B}-L(\nabla )\vec{B}+\vec{N}(\vec{B})=0,$$

$$L(\nabla )=\left(\begin{array}{ccc}-\mathrm{i}{v}_g^{(1)}\cdot \nabla & & \\ & \ddots & \\ & & -\mathrm{i}{v}_g^{(N)}\cdot \nabla \end{array}\right)-\kappa .$$

$$\hat{f}(k):=(\mathcal{F}f)(k):=(2\pi {)}^{-d/2}{\int }_{{\mathbb{R}}^d}f(x)e^{-\mathrm{i}k\cdot x}\mathrm{d}x$$

$$\sigma (L(\nabla ))={\cup }_{j\in \{1,\dots ,N\}}{\lambda }_j({\mathbb{R}}^d),$$

$$L(\mathrm{i}k){\vec{\eta }}^{(j)}(k)={\lambda }_j(k){\vec{\eta }}^{(j)}(k)$$

$${\omega }_0={\lambda }_j(k)\text{if and only if }(j,k)=(j_0,k_0).$$

$${\vec{B}}_{\text{app}}(x):=\epsilon C(\epsilon x)e^{\mathrm{i}{k}_0\cdot x}{\vec{\eta }}^{(j_0)}(k_0).$$

$${\hat{\vec{B}}}_{\text{app}}(k)={\epsilon }^{1-d}\hat{C}\left(\frac{k-k_0}{\epsilon }\right){\vec{\eta }}^{(j_0)}(k_0).$$

$$({\omega }_0+{\epsilon }^2{\omega }_1){\hat{\vec{B}}}_{\text{app}}-L(\mathrm{i}k){\hat{\vec{B}}}_{\text{app}}+\hat{\vec{N}}({\vec{B}}_{\text{app}})$$

$$=$$

$$G_0:=\frac{1}{2}{D}^2{\lambda }_{j_0}(k_0),\Gamma :={\vec{\eta }}^{(j_0)}(k_0{)}^{\ast }\hat{\vec{N}}({\vec{\eta }}^{(j_0)}(k_0)),$$

$$\vec{R}(k):=$$

$$+$$

$${\omega }_{1}C+{\nabla }^T(G_0\nabla C)+\Gamma |C{|}^{2}C=0$$

$$L_s^1({\mathbb{R}}^d):=\{f\in L^1({\mathbb{R}}^d):(1+|\cdot |{)}^{s}f\in L^1({\mathbb{R}}^d))\}.$$

$${\Vert \hat{\vec{B}}-{\epsilon }^{1-d}\hat{C}\left(\frac{\cdot -k_0}{\epsilon }\right){\vec{\eta }}^{(j_0)}(k_0)\Vert }_{L^1({\mathbb{R}}^d)}\le c_1{\epsilon }^2.$$

$$\Vert \vec{B}-\epsilon C(\epsilon \cdot ){\vec{\eta }}^{(j_0)}(k_0)e^{\mathrm{i}{k}_0\cdot }{\Vert }_{C_0({\mathbb{R}}^d)}\le c_2{\epsilon }^2.$$

$$\omega ={\omega }_0+{\epsilon }^2{\omega }_1\notin \sigma (L(\nabla )).$$

$$P_k:{\mathbb{C}}^N\to \text{span}{\vec{\eta }}^{(j_0)}(k),\vec{v}\mapsto ({\vec{\eta }}^{(j_0)}(k{)}^{\ast }\vec{v}){\vec{\eta }}^{(j_0)}(k)$$

$$Q_k:=I-P_k.$$

$$\hat{\vec{B}}(k)={\hat{\vec{B}}}_P(k)+{\hat{\vec{B}}}_Q(k),$$

$${\hat{\vec{B}}}_P(k)=P_k\hat{\vec{B}}(k)=:\psi (k){\vec{\eta }}^{(j_0)}(k),{\hat{\vec{B}}}_Q(k)=Q_k\hat{\vec{B}}(k).$$

$$\psi (k)={\epsilon }^{1-d}\hat{D}\left(\frac{k-k_0}{\epsilon }\right)+{\epsilon }^{1-d}\hat{R}\left(\frac{k-k_0}{\epsilon }\right),$$

$$\text{supp}\hat{D}\subset B_{{\epsilon }^{r-1}}:=\{k\in {\mathbb{R}}^d:|k|<{\epsilon }^{r-1}\},\text{supp}\hat{R}\subset B_{{\epsilon }^{r-1}}^c:={\mathbb{R}}^d\setminus B_{{\epsilon }^{r-1}}$$

$${\hat{\vec{B}}}_D(k):={\epsilon }^{1-d}\hat{D}\left(\frac{k-k_0}{\epsilon }\right){\vec{\eta }}^{(j_0)}(k),{\hat{\vec{B}}}_R(k):={\epsilon }^{1-d}\hat{R}\left(\frac{k-k_0}{\epsilon }\right){\vec{\eta }}^{(j_0)}(k)$$

$${\hat{\vec{B}}}_P={\hat{\vec{B}}}_D+{\hat{\vec{B}}}_R.$$

$$\omega \hat{\vec{B}}(k)-L(\mathrm{i}k)\hat{\vec{B}}(k)+\hat{\vec{N}}(\hat{\vec{B}})(k)=0,k\in {\mathbb{R}}^d.$$

Full text

Key-words: coupled mode equations, gap soliton, bifurcation, fixed point

MSC: 35Q51, 35Q55, 35Q60, 35L60, 47H10

1. Introduction

First order coupled mode equations (CMEs) are used to describe a class of wavepackets in periodic structures [2, 15, 11, 3, 10, 5, 9]. They are modulation equations for the envelopes of asymptotically broad and small wavepackets. They were first studied in nonlinear optical fiber gratings, see e.g. [2, 3]. A rigorous justification of such an approximation was performed in [15] and [11] for the one dimensional cubic nonlinear wave equation. In [9] the authors derived and justified CMEs as modulation equations for the $d$-dimensional periodic Gross-Pitaevskii equation. These CMEs have the form

[TABLE]

where for $j,r\in\{1,\dots,N\}$

[TABLE]

and where the matrix $\kappa=(\kappa_{jr})_{j,r=1}^{N}$ is Hermitian. This system (although only for the setting with $\kappa=0$) was first derived in [10]. Like with all modulation equations, the application of CMEs is not limited to the Gross-Pitaevskii equation. It applies to wavepackets centered around $N$ Bloch waves ($N$-wave mixing) with nonzero group velocities in models with nonlinearities that are cubic at lowest order.

The aim of this paper is to prove the existence of localized time harmonic solutions

[TABLE]

with $\omega$ in a gap of the linear spatial operator of (1.1). Such solutions are often called (standing) gap solitons. A necessary condition for the existence of a spectral gap is $\kappa\neq 0$. Hence, gap solitons cannot be obtained in the setting of [10]. We prove the existence of gap solitons in an asymptotic region near a spectral edge point $\omega_{0}$. The result can be interpreted as a bifurcation from the zero solution at the spectral edge.

The equation for $\vec{B}$ is

[TABLE]

where

[TABLE]

Our proof is constructive in that we use an asymptotic approximation of a solution $\vec{B}$ at $\omega=\omega_{0}+O(\varepsilon^{2}),\varepsilon\to 0$. The approximation is a modulation ansatz with a slowly varying envelope modulating the bounded linear solution at the spectral edge $\omega_{0}$. The envelope is shown to satisfy a $d-$dimensional nonlinear Schrödinger equation (NLS) with constant coefficients. We prove that for sufficiently smooth PT symmetric (parity time symmetric) solutions of the NLS there are solutions $\vec{B}$ of (1.3) at $\omega=\omega_{0}+O(\varepsilon^{2})$ which are close to the asymptotic ansatz. The proof is carried out in Fourier variables in $L^{1}(\mathbb{R}^{d})$. It is based on a decomposition of the solution in Fourier variables according to the eigenvectors of $L({\rm i}k)\in\mathbb{C}^{n\times n}$ and on a nested Banach fixed point argument. The reduced bifurcation equation is a perturbed stationary nonlinear Schrödinger equation (NLS). Solitary waves are then found via a persistence argument starting from solitary waves of the unperturbed NLS. The persistence holds in a symmetric subspace thanks to a spectral stability result of Kato.

The chosen approach is similar to that used in [7, 8, 6]. Unlike in these papers, where $L^{2}$-based spaces were used, we work here in $L^{1}$ in Fourier variables. This avoids the unfavorable scaling property of the $L^{2}$ norm of functions with an asymptotically slow dependence on $x$, namely $\|f(\varepsilon\cdot)\|_{L^{2}(\mathbb{R}^{d})}=\varepsilon^{-d/2}\|f\|_{L^{2}(\mathbb{R}^{d})}$. The $L^{1}-$approach was first used in this context for the bifurcation of time harmonic gap solitons in the one dimensional wave equation in [13].

The question of the existence of solitary waves of CMEs has previously been addressed only in one dimension in [2], where an explicit family of gap solitons was found for CMEs describing the asymptotics of wavepackets in media with infinitesimally small contrast. These gap solitons are parametrized by the velocity $v\in(-1,1)$ (after a rescaling). In [4] a numerical continuation was used to construct gap solitons also in one dimensional CMEs for finite contrast periodic structures. It was shown in [9] that for (1.1) in dimensions $d>1$ a spectral gap of $L(\nabla)$ does not exist for $N$ ungerade and for $N=2$. Next, a gap was found in a special case of (1.1) with $d=2,N=4$ and standing gap solitons were computed numerically for this case. Here we assume the presence of a spectral gap and prove the existence of standing gap solitons of the form (1.2) for $\omega$ asymptotically close to the spectrum under the condition that the spectral edge is given by an isolated extremum of the dispersion relation.

The rest of the paper consists firstly of a formal derivation of the effective NLS equation for the modulation ansatz in Section 2. Next, in Section 3 we state and prove the main approximation result. Finally, Section 4 presents a numerical example of a solution $\vec{B}$ and a numerical verification of the convergence of the asymptotic error.

2. Formal Asymptotics of Gap Solitons

The formal asymptotics of localized solutions of (1.3) were performed already in [9]. We repeat here the calculation for readers’ convenience.

The spectrum of $L(\nabla)$ can be determined using Fourier variables. We employ the Fourier transform

[TABLE]

with the inverse formula $f(x)=(\mathcal{F}^{-1}\widehat{f})(x)=(2\pi)^{-d/2}\int_{\mathbb{R}^{d}}\widehat{f}(k)e^{{\rm i}k\cdot x}\,{\rm d}k$. The spectrum of $L(\nabla)$ is

[TABLE]

where $\lambda_{j}(k)$ is the eigenvalue of $L({\rm i}k)\in\mathbb{C}^{N\times N}$ for each $k\in\mathbb{R}^{d}$, i.e.

[TABLE]

for some $\vec{\eta}^{(j)}(k)\in\mathbb{C}^{N}\setminus\{0\}$. Because $L({\rm i}k)^{*}=L({\rm i}k)$ for all $k\in\mathbb{R}^{d}$, we have $\lambda_{j}:\mathbb{R}^{d}\to\mathbb{R}$. The mapping $k\mapsto(\lambda_{1}(k),\dots,\lambda_{N}(k))^{T}$ is the dispersion relation of (1.1).

The central assumptions of our analysis are

• (A.1)

The spectrum $\sigma(L(\nabla))\subset\mathbb{R}$ has a gap, denoted by $(\alpha,\beta)$ with $\alpha<\beta$.

• (A.2)

$\omega_{0}\in\{\alpha,\beta\}$ and for some $j_{0}\in\mathbb{N},k_{0}\in\mathbb{B}$ we have

[TABLE]

As mentioned in the introduction, assumption (A.1) implies that if $d\geq 2$, then $N\geq 4$ and $N$ even. Assumption (A.2) means that the spectral edge $\omega_{0}$ is defined by one isolated extremum of the eigenvalue $\lambda_{j_{0}}$ and that this is separated at $k=k_{0}$ from all other eigenvalues.

We make the following asymptotic ansatz for a gap soliton at $\omega=\omega_{0}+\varepsilon^{2}\omega_{1}\notin\sigma(L(\nabla))$, where $\varepsilon>0$ is a small parameter and $\omega_{1}=O(1)$ (as $\varepsilon\to 0$),

[TABLE]

In Fourier variables this is

[TABLE]

Substituting (2.2) and $\omega=\omega_{0}+\varepsilon^{2}\omega_{1}$ into the Fourier transform of the left hand side of (1.3), we get

[TABLE]

where

[TABLE]

(with $v^{*}$ being the Hermitian transpose of a vector $v\in\mathbb{C}^{N}$) and where

[TABLE]

is small as shown in Sec. 3.

A necessary condition for the smallness of the residual corresponding to $\vec{B}_{\text{app}}$ is the vanishing of the square brackets. This is equivalent to

[TABLE]

for $C:\mathbb{R}^{d}\to\mathbb{C}$. Equation (2.3) is the effective nonlinear Schrödinger equation (NLS) for the envelope $C$.

3. The Bifurcation and Approximation Result

Under assumptions (A.1-A.2) and the following assumption (A.3) we prove the bifurcation result below.

• (A.3)

The kernel of the Jacobian $J$ corresponding to the NLS equation, as defined in (3.22), is $(n+1)$-dimensional (i.e. generated only by the continuous invariances of the NLS).

We define next the space $L^{1}_{s}(\mathbb{R}^{d})$ for $s\geq 0$ as

[TABLE]

For vector valued functions $f:\mathbb{R}^{d}\to\mathbb{C}^{N}$ we write $f\in L_{s}^{1}(\mathbb{R}^{d})$ if $f_{j}\in L_{s}^{1}(\mathbb{R}^{d})$ for each $j=1,\dots,N$.

The space of continuous functions $f:\mathbb{R}^{d}\to\mathbb{C}$ satisfying the asymptotics $f(x)\to 0$ as $|x|\to\infty$ is denoted by $C_{0}(\mathbb{R}^{d})$. We equip the space with the supremum norm.

Theorem 1**.**

Choose $\omega_{0}$ such that (A.1) and (A.2) are satisfied. Let $(\alpha,\beta)\subset\mathbb{R}$ be the spectral gap from (A.1) and let $\omega_{1}\in\mathbb{R}$ be such that $\text{sign}(\omega_{1})=1$ if $\omega_{0}=\alpha$ and $\text{sign}(\omega_{1})=-1$ if $\omega_{0}=\beta$. If $C$ is a $PT$-symmetric (i.e. $C(-x)=\overline{C(x)}$) solution of (2.3) with $\widehat{C}\in L^{1}_{4}(\mathbb{R}^{d})$ and such that (A.3) holds, then there are constants $c_{1},c_{2},\varepsilon_{0}>0$ such that for each $\varepsilon\in(0,\varepsilon_{0})$ there is a solution $\vec{B}$ of equation (1.3) with $\omega=\omega_{0}+\varepsilon^{2}\omega_{1}$ which satisfies $\widehat{\vec{B}}\in L^{1}_{2}(\mathbb{R}^{d})$ and

[TABLE]

In particular,

[TABLE]

The constants $c_{1}$ and $c_{2}$ depend polynomially on $\|\widehat{C}\|_{L^{1}_{4}(\mathbb{R}^{d})}$.

Clearly, due to $\widehat{\vec{B}}\in L^{1}(\mathbb{R}^{d})$ the lemma of Riemann-Lebesgue implies the decay $\vec{B}(x)\to 0$ as $|x|\to\infty$.

The existence of a $PT$-symmetric solution $C$ is satisfied, e.g., if $G_{0}$ is definite and $\text{sign}(\Gamma)=-\operatorname{sign}(\omega_{1})$. Due to the extremum of $\lambda_{j_{0}}$ at $k=k_{0}$ we have then that $\Gamma$ is positive/negative if $D^{2}\lambda{j_{0}}(k_{0})$ is positive/negative definite respectively. Hence, the NLS is of focusing type and after a rescaling of the $x$ variables it supports a real, positive, radially symmetric solution with exponential decay at infinity (Townes soliton).

Note that the condition on $\text{sign}(\omega_{1})$ implies

[TABLE]

We proceed with the proof of Theorem 1. Like in Sec. 2 we work here in Fourier variables. We employ a Lyapunov-Schmidt-like decomposition. For each $k\in\mathbb{R}^{d}$ we split the solution $\widehat{\vec{B}}(k)\in\mathbb{C}^{N}$ into the component proportional to the eigenvector $\vec{\eta}^{(j_{0})}(k)$ and the $l^{2}(\mathbb{C}^{N})-$orthogonal complement. We define the projections

[TABLE]

and

[TABLE]

Then

[TABLE]

where

[TABLE]

We aim to construct a solution $\widehat{\vec{B}}$ with $\psi$ approximated by the envelope in our ansatz, i.e. by $\varepsilon^{1-d}\widehat{C}\left(\tfrac{\cdot-k_{0}}{\varepsilon}\right)$. We choose for $\psi$ a decomposition according to the support

[TABLE]

where

[TABLE]

with $r\in(0,1)$. At the moment $r$ is a free parameter; it will be specified below. We also define

[TABLE]

such that

[TABLE]

In the ansatz in (3.1) we wish to find the component $\widehat{D}$ close to $\chi_{B_{\varepsilon^{r-1}}}\widehat{C}$ and the component $\widehat{R}$ small. The component $\widehat{\vec{B}}_{D}$ is then approximated by $\widehat{\vec{B}}_{\text{app}}$ in (2.2). If also $\widehat{\vec{B}}_{Q}$ is small, then the whole constructed solution $\widehat{\vec{B}}$ is close to $\widehat{\vec{B}}_{\text{app}}$.

The Fourier transform of (1.1) is

[TABLE]

For the selected ansatz equation (3.2) becomes

[TABLE]

Because $\omega_{0}\in\sigma(L({\rm i}k_{0}))$, the inverse of the matrix $(\omega_{0}+\varepsilon^{2}\omega_{1})I-L({\rm i}k)$ is not bounded uniformly in $\varepsilon$. In a neighbourhood of $k_{0}$ the norm of the inverse blows up as $\varepsilon\to 0$. However,

[TABLE]

is invertible uniformly in $\varepsilon$ due to assumption (A.2).

We separate the explicit part of (3.4) by writing

[TABLE]

and

[TABLE]

where $\widehat{\vec{B}}_{Q,1}$ solves the explicit part, i.e.

[TABLE]

The system to solve is thus

[TABLE]

Our procedure for constructing a solution can be sketched as follows.

1. (1)

For any $\widehat{D},\widehat{R}\in L^{1}(\mathbb{R}^{d})$ equation (3.5) produces a small $\widehat{\vec{B}}_{Q,1}$ (because $\|\widehat{\vec{B}}_{P}\|_{L^{1}}=O(\varepsilon)$). 2. (2)

For any $\widehat{D},\widehat{R}\in L^{1}(\mathbb{R}^{d})$ and $\widehat{\vec{B}}_{Q,1}$ from step 1 we solve (3.7) by a fixed point argument for a small $\widehat{\vec{B}}_{Q,2}$. 3. (3)

For any $\widehat{D}\in L^{1}(\mathbb{R}^{d})$ and for $\widehat{\vec{B}}_{Q}$ from steps 1 and 2 we solve (3.6) with $k\in B_{\varepsilon^{r}}(k_{0})^{c}$ for a small $\widehat{R}$ by a fixed point argument. 4. (4)

With the components obtained in the above steps we find a solution $\widehat{D}\in L_{2}^{1}(\mathbb{R}^{d})$ of (3.6) with $k\in B_{\varepsilon^{r}}(k_{0})$ close to a $\widehat{C}\in L^{1}_{2}(\mathbb{R}^{d})$ (with $C$ a solution of (2.3)) - provided such a $C$ exists. In addition $C$ needs to satisfy a certain symmetry, the $PT$-symmetry. Also here a fixed point argument is used - roughly speaking for the difference $\widehat{D}-\widehat{C}$. 5. (5)

The error $\|\widehat{\vec{B}}-\widehat{\vec{B}}_{\text{app}}\|_{L^{1}(\mathbb{R}^{d})}$ is $O(\varepsilon^{2})$ if $\widehat{C}$ decays fast enough, namely if $\widehat{C}\in L^{1}_{4}(\mathbb{R}^{d})$.

Lemma 1**.**

If $\widehat{D}\in L^{1}(\mathbb{R}^{d})$, $\widehat{R}\in L^{1}_{s_{R}}(\mathbb{R}^{d})$ for some $s_{R}\geq 0$, and $\text{supp}\widehat{D}\subset B_{\varepsilon^{r-1}},$ $\text{supp}\widehat{R}\subset B_{\varepsilon^{r-1}}^{c}$, then there are constant $c_{1},c_{2}>0$ such that for all $\varepsilon>0$ small enough

[TABLE]

and

[TABLE]

Proof.

Because $\|\widehat{f}(\varepsilon^{-1}(\cdot-k_{0}))\|_{L^{1}(\mathbb{R}^{d})}=\varepsilon^{d}\|\widehat{f}\|_{L^{1}(\mathbb{R}^{d})}$, we have

[TABLE]

The estimate for $\|\widehat{\vec{N}}(\widehat{\vec{B}}_{P})\|_{L^{1}(\mathbb{R}^{d})}$ follows by Young’s inequality for convolutions. $\Box$

• (1)

Component $\widehat{\vec{B}}_{Q,1}$

Because $\widehat{\vec{B}}_{Q,1}=-M_{k}^{-1}Q_{k}\widehat{\vec{N}}(\widehat{\vec{B}}_{P})(k)$, we get from Lemma 1 the estimate

[TABLE]

• (2)

Component $\widehat{\vec{B}}_{Q,2}$

With $\widehat{\vec{B}}_{Q,1}$ from above (and $\widehat{D},\widehat{R}$ given) component $\widehat{\vec{B}}_{Q,2}$ satisfies

[TABLE]

Due to the cubic structure of $\vec{N}$ we have

[TABLE]

For $\widehat{\vec{B}}_{Q,2}\in B_{\varepsilon^{\eta}},\eta>0$ we have

[TABLE]

Hence, $\vec{G}:B^{(L^{1})}_{c_{0}\varepsilon^{5}}\to B^{(L^{1})}_{c_{0}\varepsilon^{5}}$ for some $c_{0}(\|\widehat{D}\|_{L^{1}},\|\widehat{R}\|_{L^{1}})>0$, where $B^{(L^{1})}_{\alpha}:=\{\vec{v}\in L^{1}(\mathbb{R}^{d}):\|\vec{v}\|_{L^{1}(\mathbb{R}^{d})}\leq\alpha\}.$ The constants $c$ and $c_{0}$ depend polynomially on $\|\widehat{D}\|_{L^{1}}$ and $\|\widehat{R}\|_{L^{1}}$.

Similary, we obtain the contraction (for $\varepsilon>0$ small enough)

[TABLE]

if $\widehat{\vec{B}}^{(1)}_{Q,R},\widehat{\vec{B}}^{(2)}_{Q,R}\in B^{(L^{1})}_{c_{0}\varepsilon^{5}}$. Hence, for $\varepsilon>0$ small enough we have a unique solution $\widehat{\vec{B}}_{Q,2}\in B^{(L^{1})}_{c_{0}\varepsilon^{5}}$ of $\widehat{\vec{B}}_{Q,2}=\vec{G}(\widehat{\vec{B}}_{Q,2})$, i.e.

[TABLE]

• (3)

Component $\widehat{\vec{B}}_{R}$

For any $\widehat{D}\in L^{1}(\mathbb{R}^{d})$ and with the above estimate on $\widehat{\vec{B}}_{Q}$ we look for a small $\widehat{R}$. The support of $\widehat{R}$ is $B_{\varepsilon^{r-1}}^{c}$, whence for $k\in\mathbb{R}^{d}\setminus B_{\varepsilon^{r}}(k_{0})$ we can divide in (3.6) by $\omega_{0}+\varepsilon^{2}\omega_{1}-\lambda_{j_{0}}(k)$ and obtain

[TABLE]

with

[TABLE]

Since $\nabla\lambda_{j_{0}}(k_{0})=0$, we have $|\lambda_{j_{0}}(k)-\omega_{0}|>c\varepsilon^{2r}$ for all $k\in\mathbb{R}^{d}\setminus B_{\varepsilon^{r}}(k_{0})$ and hence $|\nu(k)|\leq c\varepsilon^{-2r}$. In order to exploit the localized nature of $\widehat{\vec{B}}_{D}$ and the smallness of $\widehat{\vec{B}}_{Q}$, we write for $\kappa:=\tfrac{k-k_{0}}{\varepsilon}$

[TABLE]

where

[TABLE]

with $h(k):=\left(1+\tfrac{|k-k_{0}|}{\varepsilon}\right)^{s_{D}}$. Using $\sup_{k\in\text{supp}\widehat{\vec{B}}_{R}}|h(k)^{-1}|\leq c\varepsilon^{(1-r)s_{D}}$, the cubic form of $\widehat{\vec{N}}$ and the fact that $\left(\widehat{D}\left(\tfrac{\cdot-k_{0}}{\varepsilon}\right)*\widehat{\overline{D}}\left(\tfrac{\cdot+k_{0}}{\varepsilon}\right)*\widehat{D}\left(\tfrac{\cdot-k_{0}}{\varepsilon}\right)\right)(k_{0}+\varepsilon\kappa)=\varepsilon^{2d}(\widehat{D}*\widehat{\overline{D}}*\widehat{D})(\kappa)$, we get

[TABLE]

Since $\|\widehat{\vec{B}}_{R}\|_{L^{1}(\mathbb{R}^{d})}\leq c\varepsilon\|\widehat{R}\|_{L^{1}(\mathbb{R}^{d})}$, $\|\widehat{\vec{B}}_{D}\|_{L^{1}(\mathbb{R}^{d})}\leq c\varepsilon\|\widehat{D}\|_{L^{1}(\mathbb{R}^{d})}$, and $\|\widehat{\vec{B}}_{Q}\|_{L^{1}(\mathbb{R}^{d})}\leq c(\|\widehat{D}\|_{L^{1}},\|\widehat{R}\|_{L^{1}})\varepsilon^{3}$ (with a polynomial $c$), we have

[TABLE]

with $c$ depending polynomially on $\|\widehat{D}\|_{L^{1}}$ and $\|\widehat{R}\|_{L^{1}}$.

If $\widehat{D}\in L^{1}_{s_{D}}(\mathbb{R}^{d})$, $s_{D}\geq 0$, then there is a constant $c$ depending polynomially on $\|\widehat{D}\|_{L^{1}_{s_{D}}(\mathbb{R}^{d})}$ such that

[TABLE]

for all $R$ with $\|\widehat{R}\|_{L^{1}(\mathbb{R}^{d})}\leq c\varepsilon^{\alpha}$ and all $\varepsilon>0$ small enough. Hence $H:B^{(L^{1})}_{c\varepsilon^{\alpha}}\to B^{(L^{1})}_{c\varepsilon^{\alpha}}$ for some $c>0$ if $\varepsilon>0$ is small enough and if $\widehat{D}\in L^{1}_{s_{D}}(\mathbb{R}^{d})$.

Similarly, we get the contraction property of $H$ on $B^{(L^{1})}_{c\varepsilon^{\alpha}}$ for $\varepsilon>0$ small enough. The constructed fixed point $\widehat{R}\in B^{(L^{1})}_{c\varepsilon^{\alpha}}$ yields for any $s_{D}\geq 0$

[TABLE]

• (4)

Component $\widehat{\vec{B}}_{D}$

Finally, we consider the component $\widehat{\vec{B}}_{D}$. For $k\in\text{supp}(\widehat{\vec{B}}_{D})=B_{\varepsilon^{r}}(k_{0})$ we rewrite (3.6) as follows. We add and subtract $\widehat{\vec{N}}(\widehat{\vec{B}}_{D})$ like in (3.10), we Taylor expand $\lambda_{j_{0}}(k)$ at $k=k_{0}$, and we use the variable $\kappa=\varepsilon^{-1}(k-k_{0})$. This leads to

[TABLE]

where

[TABLE]

and $\delta(k):=\lambda_{j_{0}}(k)-\omega_{0}-\frac{1}{2}(k-k_{0})^{T}D^{2}\lambda_{j_{0}}(k_{0})(k-k_{0})$.

Next, we write $\rho=\rho_{1}+\rho_{2}+\rho_{3}$, where

[TABLE]

For $\rho_{1}$ we get

[TABLE]

Because $\Gamma=\sum_{m,n,o,j\in\{1,\dots,N\}}\gamma_{j}^{(m,n,o)}\overline{\eta}_{j}^{(j_{0})}(k_{0})\eta_{m}^{(j_{0})}(k_{0})\overline{\eta}_{n}^{(j_{0})}(k_{0})\eta_{o}^{(j_{0})}(k_{0})$, we get

[TABLE]

due to the Lipschitz continuity of $k\mapsto\eta^{(j_{0})}(k)$. Hence, by Young’s inequality for convolutions,

[TABLE]

For $\rho_{2}$ we note that $|\delta(k_{0}+\varepsilon\kappa)|\leq c\varepsilon^{3}|\kappa|^{3}$ for $\kappa\in B_{\varepsilon^{r-1}}$ and $\varepsilon>0$ small enough. We estimate

[TABLE]

for any $\beta\in[0,3)$.

Finally, we estimate $\rho_{3}$. Note that $\widehat{\vec{N}}(\widehat{\vec{B}})(k_{0}+\varepsilon\cdot)-\widehat{\vec{N}}(\widehat{\vec{B}}_{D})(k_{0}+\varepsilon\cdot)$ appears also in (3.10). We have

[TABLE]

for $\varepsilon>0$ small enough, where we have made use of (3.8), (3.9), (3.12), the fact that $\|\widehat{\vec{B}}_{D}\|_{L^{1}}\leq\varepsilon\|\widehat{D}\|_{L^{1}}$ and the estimate $\|\widehat{R}\|_{L^{1}}\leq c(\|\widehat{D}\|_{L^{1}_{s_{D}}})\varepsilon^{\alpha}$. The dependence of $c_{1}$ and $c_{2}$ on $\|\widehat{D}\|_{L^{1}_{s_{D}}}$ is polynomial. As a result

[TABLE]

The whole right hand side of (3.13) is thus estimated as

[TABLE]

for any $\beta\in[0,3)$. Once again, the constant $c$ depends polynomially on its arguments.

In order to solve (3.13) for $\widehat{D}$ below (using a fixed point argument), we need to consider $\rho$ as an inhomogeneity. The linearized operator to be inverted in the iteration is of second order such that in Fourier space it acts from $L^{1}_{2}(\mathbb{R}^{d})$ to $L^{1}(\mathbb{R}^{d})$. For that reason we need to choose $s_{D}=2$ and $\beta\geq 1$ above. The choice $s_{D}=2,\beta=1,r=1/2$ leads to $\alpha=2$ and $1-\beta(1-r)=1/2$, i.e. $\|\rho(\widehat{D})\|_{L^{1}(\mathbb{R}^{d})}\leq c\varepsilon^{1/2}$ provided $\widehat{D}\in L^{1}_{2}(\mathbb{R}^{d})$. For $s_{D}=2,\beta=1$ the largest value of $\min\{1-\beta(1-r),\alpha\}$ is $4/5$ attained at $r=4/5$. Hence, with $r=4/5$ we get the best possible estimate

[TABLE]

This order determines the accuracy of the approximation and turns out to be insufficient. It leads to $\|\widehat{\vec{B}}-\widehat{\vec{B}}_{\text{app}}\|_{L^{1}(\mathbb{R}^{d})}\leq c\varepsilon^{9/5}$ instead of $c\varepsilon^{2}$.

Clearly, the leading order term in the residual $\rho$ is caused by the error from the Taylor expansion of $\lambda_{j_{0}}$. We introduce a refined ansatz for $\widehat{D}$ in order to make this error of higher order. Note that equation (3.13) is a perturbation of the NLS (2.3). Writing

[TABLE]

equation (3.13) is

[TABLE]

We search for $D$ close to a solution $C$ of the NLS, i.e. of $f_{\text{NLS}}(C)=0$. For that we need to look for $\widehat{D}$ in a vicinity of

[TABLE]

We choose the following ansatz

[TABLE]

where $D^{3}\lambda_{j_{0}}(k_{0})(\kappa,\kappa,\kappa)$ is the third order term in the Taylor expansion of $\varepsilon^{-3}\lambda_{j_{0}}(k_{0}+\varepsilon\kappa)$ and where $\text{supp}(\widehat{d})\subset B_{\varepsilon^{r-1}}$. We look for a solution with a small $\widehat{d}$.

We also define

[TABLE]

such that $\widehat{D}=\widehat{D}_{0}+\varepsilon\nu\widehat{C}^{(\varepsilon)}$. With this notation equation (3.17) reads

[TABLE]

Compared to $\rho$ the right hand side $\tilde{\rho}$ is smaller as we show next. It is

[TABLE]

where

[TABLE]

To estimate $\tilde{\rho}_{2}$ note that the first Taylor expansion error term (i.e. the first line) can be estimated in $L^{1}(\mathbb{R}^{d})$ by $c\varepsilon^{2}\|\widehat{C}^{(\varepsilon)}\|_{L^{1}_{4}(\mathbb{R}^{d})}\leq c\varepsilon^{2}\|\widehat{C}\|_{L^{1}_{4}(\mathbb{R}^{d})}$. For the second Taylor expansion error note first that $|\nu(\kappa)|\leq c(1+|\kappa|)$ for all $\kappa\in\mathbb{R}^{d}$ such that

[TABLE]

Just like in (3.15) we get

[TABLE]

for any $\beta\in[0,3)$. Once again, we need to choose $\beta\geq 1$ such that $\tilde{\rho}_{2}:\widehat{d}\mapsto\tilde{\rho}_{2}(\widehat{d})$ is a mapping from $L^{1}_{2}(\mathbb{R}^{d})$ to $L^{1}(\mathbb{R}^{d})$. With $\beta=1$ and $r=1/2$ we get $1-\beta(1-r)=1/2$.

The last term in $\tilde{\rho}_{2}$ is

[TABLE]

In the $L^{1}$-norm this can be estimated using Young’s inequality by

[TABLE]

In summary, with the above choice of $\beta$ and $r$ we get

[TABLE]

The terms $\rho_{1}$ and $\rho_{3}$ are estimated in (3.14) and (3.16). As discussed above, we seek $\widehat{D}$ in $L^{1}_{2}(\mathbb{R}^{d})$ and hence we set $s_{D}=2$ in (3.11). This yields $\alpha=2$ and

[TABLE]

In summary, for any $\widehat{C}\in L^{1}_{4}(\mathbb{R}^{d})$ fixed (with $C$ being a solution of the NLS), the right hand side of (3.19) is estimated as

[TABLE]

We proceed with a fixed point argument for the correction $d$. In order to obtain a differentiable function (to use the Jacobian of the NLS), we write $f_{\text{NLS}}$ in real variables. Writing $D=D_{R}+{\rm i}D_{I},C=C_{R}+{\rm i}C_{I}$, $d=d_{R}+{\rm i}d_{I}$, and $\tilde{\rho}=\tilde{\rho}_{R}+{\rm i}\tilde{\rho}_{I}$, we define

[TABLE]

the Jacobian

[TABLE]

as well as the Fourier-truncation of the Jacobian

[TABLE]

where $\widehat{C_{R}}^{(\varepsilon)}:=\chi_{B_{\varepsilon^{r-1}}}\widehat{C_{R}},C_{R}^{(\varepsilon)}:={\mathcal{F}}^{-1}(\widehat{C_{R}}^{(\varepsilon)})$ and $\widehat{C_{I}}^{(\varepsilon)}:=\chi_{B_{\varepsilon^{r-1}}}\widehat{C_{I}}$, $C_{I}^{(\varepsilon)}:={\mathcal{F}}^{-1}(\widehat{C_{I}}^{(\varepsilon)})$. $\widehat{J}_{\varepsilon}$ has the form

[TABLE]

With this notation (3.19) reads

[TABLE]

where

[TABLE]

with $\widehat{\vec{C}}^{(\varepsilon)}:=(\widehat{C_{R}}^{(\varepsilon)},\widehat{C_{I}}^{(\varepsilon)})^{T}$ and $\vec{\tilde{\rho}}:=(\tilde{\rho}_{R},\tilde{\rho}_{I})^{T}$.

The aim is to construct a small fixed point $\widehat{\vec{d}}\in L^{1}_{s_{D}}(\mathbb{R}^{d})$ of $\widehat{J}_{\varepsilon}^{-1}W$. The difficulty is that the inverse of $\widehat{J}_{\varepsilon}$ is not bounded uniformly in $\varepsilon$. This is due to the presence of the $d+1$ zero eigenvalues of $J$ caused by the $d$ spatial shift invariances and the phase invariance of the NLS, see assumption (A.3). The distance of the essential spectrum of $\widehat{J}_{\varepsilon}$ from zero is $|\omega_{1}|$ due to the choice of $\text{sign}(\omega_{1})$. To eliminate the zero eigenvalues, we work in a symmetric subspace of $L^{1}_{s_{D}}(\mathbb{R}^{d})$ in which the invariances do not hold. A natural symmetry is the $PT$-symmetry. Hence, we consider the fixed point problem

[TABLE]

in the space

[TABLE]

Note that $\widehat{J}_{0},\widehat{J}_{\varepsilon}:L^{1}_{q}(\mathbb{R}^{d})\to L^{1}_{q-2}(\mathbb{R}^{d})$ for any $q\geq 2$. Because of assumption (A.3) $\widehat{J}_{0}^{-1}$ is bounded in $X_{s_{D}}^{\text{sym}}$ for any $s_{D}\geq 2$. Since $\widehat{J}_{\varepsilon}$ is a perturbation of $\widehat{J}_{0}$, we still need to ensure that [math] is not an eigenvalue of $\widehat{J}_{\varepsilon}$. For that we use a spectral stability result of Kato, see [12, Theorem IV.3.17]. Applied to our problem in the Banach space $X_{s_{D}-2}^{\text{sym}},s_{D}\geq 2$ with the domain of $\widehat{J}_{0}$ being $D(\widehat{J}_{0})=X_{s_{D}}^{\text{sym}}$, it reads:

Assume $\widehat{J}_{\varepsilon}-\widehat{J}_{0}$ is $\widehat{J}_{0}$-bounded, i.e. $\mbox{domain}(\widehat{J}_{0})\subset\mbox{domain}(\widehat{J}_{\varepsilon}-\widehat{J}_{0})$ and for some $a,b\geq 0$ is

[TABLE]

If for some $\zeta\in\rho(\widehat{J}_{0})$

[TABLE]

then $\zeta\in\rho(\widehat{J}_{\varepsilon})$ and

[TABLE]

We check now (3.23) and (3.24) for $\zeta=0$. For $\widehat{\vec{d}}\in X_{s_{D}}^{\text{sym}}$ one has

[TABLE]

Writing $\widehat{C_{R,I}}^{(\varepsilon)}=:\widehat{C_{R,I}}+\widehat{\gamma_{R,I}}^{(\varepsilon)}$, it is $\text{supp}(\widehat{\gamma_{R,I}}^{(\varepsilon)})\subset B_{\varepsilon^{r-1}}^{c}$. The difference $(\widehat{J}_{\varepsilon}-\widehat{J}_{0})\widehat{\vec{d}}$ consists of terms that are linear or quadratic in $\widehat{\gamma_{R,I}}^{(\varepsilon)}$; for instance terms like $\widehat{\gamma_{R}}^{(\varepsilon)}*\widehat{C_{R}}^{(\varepsilon)}*\widehat{d_{R}}$ or $\widehat{\gamma_{I}}^{(\varepsilon)}*\widehat{C_{R}}^{(\varepsilon)}*\widehat{d_{I}\hskip 1.42271pt}$. Because

[TABLE]

Young’s inequality for convolutions yields

[TABLE]

with $c$ depending polynomially on $\|\widehat{C}\|_{L^{1}_{s_{D}}}$. Conditions (3.23) and (3.24) for $\zeta=0$ are thus satisfied with $a=c\varepsilon^{2(1-r)}$ and $b=0$ if $\varepsilon>0$ is small enough and if $\widehat{C}\in L^{1}_{s_{D}}(\mathbb{R}^{d})$.

For the fixed point problem we use $s_{D}=2$ and $r=1/2$ and show firstly that if $\widehat{C}\in X_{4}^{\text{sym}}$, then there is some $c>0$ such that for all $\varepsilon>0$ small enough

[TABLE]

where $B_{c\varepsilon}^{2,\text{ sym}}$ is the $c\varepsilon$-ball in $X_{2}^{\text{sym}}$, i.e.

[TABLE]

Note that the requirement $\widehat{C}\in L^{1}_{4}$ is dictated by (3.20).

Secondly, we prove that $\widehat{J}_{\varepsilon}^{-1}W$ is contractive provided $\widehat{C}\in\{\widehat{f}\in L^{1}_{4}(\mathbb{R}^{d}):\operatorname{Im}(\widehat{f_{R}})=-\operatorname{Re}(\widehat{f_{I}})\}$. We start by showing $\widehat{J}_{\varepsilon}:X_{2}^{\text{sym}}\to X_{0}^{\text{sym}}$. The loss of $2$ in the weight is due to the second order nature of the operator $J$, i.e. due to the factor $\kappa^{T}G_{0}\kappa$ in Fourier variables. The entries $\omega_{1}+\kappa^{T}G_{0}\kappa$ clearly preserve the $PT$-symmetry. For the convolution terms we have, for instance

[TABLE]

Hence $C_{R}^{(\varepsilon)^{2}}d_{R},C_{I}^{(\varepsilon)^{2}}d_{R}$, and $C_{R}^{(\varepsilon)}C_{I}^{(\varepsilon)}d_{I}$ are even and $C_{I}^{(\varepsilon)^{2}}d_{I},C_{R}^{(\varepsilon)^{2}}d_{I}$, and $C_{R}^{(\varepsilon)}C_{I}^{(\varepsilon)}d_{R}$ are odd such that the $PT-$symmetry is preserved also by the convolution terms. In end effect, $\operatorname{Im}((\widehat{J}_{\varepsilon}\widehat{\vec{d}})_{1})=-\operatorname{Re}((\widehat{J}_{\varepsilon}\widehat{\vec{d}})_{2})$. Hence, $\widehat{J}_{\varepsilon}:X_{2}^{\text{sym}}\to X_{0}^{\text{sym}}$ and for $\widehat{J}_{\varepsilon}^{-1}$ we get $\widehat{J}_{\varepsilon}^{-1}:X_{0}^{\text{sym}}\to X_{2}^{\text{sym}}$.

Next, we show that $W:B_{c\varepsilon}^{2,\text{ sym}}\to B_{c\varepsilon}^{0,\text{sym}}$ if $\widehat{C}\in X_{4}^{\text{sym}}$. The term $\rho$ is estimated in (3.21) and dictates the order $\varepsilon^{1}$.

The difference $F_{\text{NLS}}(\widehat{\vec{C}}^{(\varepsilon)}+\widehat{\vec{d}})-\widehat{J}_{\varepsilon}\widehat{\vec{d}}$ consists of terms quadratic in $\widehat{\vec{d}}$ and hence is bounded in $L^{1}(\mathbb{R}^{d})$ by $c_{1}(\|\widehat{\vec{d}}\|_{L^{1}(\mathbb{R}^{d})}^{2}+\|\widehat{\vec{d}}\|_{L^{1}(\mathbb{R}^{d})}^{3})$. In summary,

[TABLE]

if $\|\widehat{\vec{d}}\|_{L^{1}_{2}(\mathbb{R}^{d})}\leq c\varepsilon$ and $\varepsilon>0$ is small enough. Due to the boundedness of $\widehat{J}_{\varepsilon}^{-1}:X_{0}^{\text{sym}}\to X_{2}^{\text{sym}}$ we thus have (3.25).

The contractive property of $\widehat{J}_{\varepsilon}^{-1}W$ in $B_{c\varepsilon}^{2,\text{ sym}}$ is now clear due to the quadratic nature of $F_{\text{NLS}}(\widehat{\vec{C}}^{(\varepsilon)}+\widehat{\vec{d}})-\widehat{J}_{\varepsilon}\widehat{\vec{d}}$.

We conclude that if the solution $C$ of the NLS (2.3) satisfies $\widehat{C}\in\{\widehat{f}\in L^{1}_{4}(\mathbb{R}^{d}):\operatorname{Im}(\widehat{C_{R}})=-\operatorname{Re}(\widehat{C_{I}})\}$, then there is $c>0$ such that for all $\varepsilon>0$ small enough the constructed solution $D$ of (3.13) satisfies

[TABLE]

and due to (3.18)

[TABLE]

Here we have also used $\|\varepsilon\nu\chi_{B_{\varepsilon^{-1/2}}}\widehat{C}\|_{L^{1}_{2}}\leq c\varepsilon\|\widehat{C}\|_{L^{1}_{3}}$.

This allows us to estimate $\widehat{\vec{B}}_{D}-\widehat{\vec{B}}_{\text{app}}$. We have

[TABLE]

Next we use (3.26), the Lipschitz continuity of $\vec{\eta}^{(j_{0})}$, and the estimate $\|\widehat{C}\|_{L^{1}(B^{c}_{\varepsilon^{r-1}})}\leq\varepsilon^{(1-r)s_{C}}\|\widehat{C}\|_{L^{1}_{s_{C}}(\mathbb{R}^{d})}$ for all $s_{C}\geq 0$. This produces at $r=1/2$

[TABLE]

if $s_{C}=4$.

We can now summarize the error estimate

[TABLE]

The components $\widehat{\vec{B}}_{Q}$ and $\widehat{\vec{B}}_{R}$ are estimated in (3.8), (3.9), and (3.12). Having now estimated $\|\widehat{R}\|_{L^{1}}$ in terms of $\|\widehat{C}\|_{L^{1}}$ and $\|\widehat{D}\|_{L^{1}}$ in terms of $\|\widehat{C}\|_{L^{1}_{2}}$, we get for $r=1/2$ and $s_{D}=2$

[TABLE]

where $c_{1}$ and $c_{2}$ depend polynomially on $\|\widehat{C}\|_{L^{1}_{2}(\mathbb{R}^{d})}$. Hence, the estimate in Theorem 1 is proved.

4. Numerical Example of Bifurcating Gap Solitons for $d=2$

In [9] it is shown that assumption (A.1), i.e. the existence of a spectral gap is satisfied for $N=4$ in the symmetric case

[TABLE]

provided $|\alpha_{1}|^{2}>2(|\alpha_{2}|^{2}+|\alpha_{3}|^{2})$. In the following example we choose $v=(0,1)^{T},w=(1,0)^{T}$, $\alpha_{1}=2,$ and $\alpha_{2}=\alpha_{3}=1$. The dispersion relation $\omega_{j}:\mathbb{R}^{2}\to\mathbb{R},j=1,\dots,4$ of (1.1) is plotted in Fig. 1. The gap appears even though the sufficient condition $|\alpha_{1}|^{2}>2(|\alpha_{2}|^{2}+|\alpha_{3}|^{2})$ is not satisfied.

We see that the second eigenvalue $\lambda_{2}$ has an isolated maximum at $k=k_{0}:=0$. The corresponding frequency is $\omega_{0}:=\lambda_{2}(0)=0$. The eigenvector corresponding to $\lambda_{2}(0)$ is $\vec{\eta}_{j_{0}}(0)=\tfrac{1}{\sqrt{2}}(1,1,-1,-1)^{T}$.

We use the following special case of the coefficients $\gamma_{j}^{(m,n,o)}$ in (1.3)

[TABLE]

Clearly, coefficients (4.1) and (4.2) allow symmetric solutions with $B_{2}=\overline{B_{1}}$ and $B_{4}=\overline{B_{3}}$. We do not make a direct use of this symmetry in our computations. We construct the approximation $\vec{B}_{\text{app}}$ of a solution of (1.3) at $\omega=\omega_{0}+\varepsilon^{2}\omega_{1}$ for six values of $\varepsilon$: $0.2,0.1,0.05,0.025,0.0125$, and $0.00625$. The coefficients of the effective NLS (2.3) are $G_{0}=-0.25~{}I_{2x2}$ and $\Gamma=2.25$ and we choose $\omega_{1}=1$. A real $C$ radially symmetric was chosen in this example. It was computed using the shooting method for the NLS in polar variables.

Using the numerical Petviashvili iteration [14, 1], we also produce a numerical approximation of a solution $\vec{B}$ at $\omega=\omega_{0}+\varepsilon^{2}$. The Petviashvili iteration is a fixed point iteration in Fourier variables with a stabilizing normalization factor. The initial guess of the iteration was chosen as $\vec{B}_{\text{app}}$. Note that although $\vec{B}_{\text{app}}$ can be real (if a real solution $C$ of the NLS is chosen), equation (1.3) does not allow real solutions $\vec{B}$ due to the term ${\rm i}\nabla\vec{B}$ and due to the realness of $\alpha_{1},\alpha_{2},$ and $\alpha_{3}$ and of $\gamma_{j}^{(m,n,o)}$. Nevertheless, if $\vec{B}_{\text{app}}$ is real, there must be a solution $\vec{B}$ with $\text{Im}(\vec{B})=O(\varepsilon^{2})$. Figure 2 shows $\vec{B}_{\text{app}}$ and $\vec{B}$ for $\varepsilon=0.05$.

The numerical parameters for the Petviashvili iteration were selected as follows: we compute on the domain $x\in[-3/\varepsilon,3/\varepsilon]^{2}$ with the discretization given by 160x160 grid points, i.e. $dx_{1}=dx_{2}=3/(80\varepsilon)$. Note that because $k_{0}=0$, the relatively coarse discretization for small values of $\varepsilon$ does not matter (there are no oscillations to be resolved).

For each $\varepsilon$ we evaluate the asymptotic error $E:=\|\vec{B}-\vec{B}_{\text{app}}\|_{C_{0}}$. Figure 3 shows the convergence of the error in $\varepsilon$. Clearly, $E(\varepsilon)\sim c\varepsilon^{2}$, which confirms the convergence rate proved in our theorem.

Acknowledgements

This research is supported by the German Research Foundation, DFG grant No. DO1467/3-1.