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

We consider a system of first order coupled mode equations in $\mathbb{R}^d$ describing the envelopes of wavepackets in 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 bifurcating solution is approximated in terms of a solution of an effective nonlinear Schr\"odinger equation. The proof is based on a decomposition in Fourier variables and a nested Banach fixed point argument. A numerical example in $\mathbb{R}^2$ is provided.


Introduction
In [6] the authors have derived and justified coupled mode equations (CMEs) as modulation equations for the envelopes of a class of wavepackets in the d-dimensional periodic nonlinear Schrödinger equation. These CMEs have the form (1.1) ipB t A j`v pjq g¨∇ A j q`N ÿ r"1 κ jr A r`Nj p Aq " 0, j " 1, . . . , N, where for j, r P t1, . . . , N u N j p Aq :" ÿ pα,β,γqPt1,...,N u 3 γ pα,β,γq j A α A β A γ , γ j P C, κ jr P C, v pjq g P R d , and where the matrix κ " pκ jr q N j,r"1 is Hermitian. The aim of this paper is to prove the existence of localized time harmonic solutions (1.2) Apx, tq " e´i ωt Bpxq, | Bpxq| Ñ 0 as |x| Ñ 8, with ω in a gap of the linear spatial operator of (1.1). Such solutions are often called (standing) gap solitons. We prove their existence in an asymptotic region near a spectral edge point ω 0 . The result can be interpreted as a bifurcation from the zero solution at the spectral edge. The equation for B is Our proof is constructive in that we use an asymptotic approximation of a solution B at ω " ω 0`O pε 2 q, ε Ñ 0. The approximation is a modulation ansatz with a slowly varying envelope modulating the bounded linear solution at the spectral edge ω 0 . The envelope is shown to satisfy a d´dimensional nonlinear Schrödinger equation (NLS) with constant coefficients. We prove that for sufficiently smooth solutions of the NLS there are solutions B of (1.3) at ω " ω 0`O pε 2 q which are close to the asymptotic ansatz. The proof is carried out in Fourier variables in L 1 pR d q. It is based on a decomposition of the solution in Fourier variables according to the eigenvectors of Lpikq P C nˆn and on a nested Banach fixed point argument.
First order CMEs have been studied previously to describe wavepackets in periodic structures [2,11,8,3,7,5]. However, the question of the existence of solitary waves has 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 small contrast media. These are parametrized by velocity v P p´1, 1q (after a rescaling). In [4] a numerical continuation was used to construct gap solitons also in one dimensional CMEs for finite contrast periodic structures. In higher dimensions it was shown in [6] that for (1.1) a spectral gap of Lp∇q exists only for N ě 4. Next, standing gap solitons were computed numerically for d " 2, N " 4. Here we prove the existence of standing GSs of the form (1.2) for ω 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 B and a numerical verification of the convergence of the asymptotic error.

Formal Asymptotics of Gap Solitons
The formal asymptotics of localized solutions of (1.3) were performed already in [6]. We repeat here the calculation for readers' convenience.
The spectrum of Lp∇q can be determined using Fourier variables. We employ the Fourier transform Lpikq η pjq pkq " λ j pkq η pjq pkq for some η pjq pkq P C N . Because Lpikq˚" Lpikq for all k P R d , we have λ j : R d Ñ R. The mapping k Þ Ñ pλ 1 pkq, . . . , λ N pkqq T is the dispersion relation of (1.1).
The central assumptions of our analysis are (A.1) The spectrum σpLp∇qq Ă R has a gap, denoted by pα, βq with α ă β. (A.2) ω 0 P tα, βu and for some j 0 P N, k 0 P B we have ω 0 " λ j pkq if and only if pj, kq " pj 0 , k 0 q. Assumption (A.2) means that the spectral edge ω 0 is defined by one isolated extremum of the eigenvalue λ j 0 and that this is separated at k " k 0 from all other eigenvalues.
In Fourier variables this is Substituting (2.2) and ω " ω 0`ε 2 ω 1 into the Fourier transform of the left hand side of (1.3), we get where G 0 :" 1 2 D 2 λj 0 pk 0 q, Γ :" η pj 0 q pk 0 q˚x N p η pj 0 q pk 0 qq, and where Rpkq :"ε 3˜λ j 0 pkq´ω 0´ˆk´k 0 ε˙T G 0 k´k 0 ε¸ η pj 0 q pk 0 q p Bˆk´k 0 εε 3ˆx N p η pj 0 q pk 0 qq´ η pj 0 q pk 0 q˚x N p η pj 0 q pk 0 qq η pj 0 q pk 0 qq˙p p B˚p B˚p Bqˆk´k 0 εi s small as shown in Sec. 3. A necessary condition for the smallness of the residual corresponding to B app is the vanishing of the square brackets. This is equivalent to 3) is the effective nonlinear Schrödinger equation (NLS) for the envelope C.  3) are satisfied. Let pα, βq Ă R be the spectral gap from (A.1) and let ω 1 P R be such that signpω 1 q " 1 if ω 0 " α and signpω 1 q "´1 if ω 0 " β. If C is a P T -symmetric (i.e. Cp´xq " Cpxq) solution of (2.3) with p C P L 1 4 pR d q, then there are constants c 1 , c 2 , ε 0 ą 0 such that for each ε P p0, ε 0 q there is a solution B of equation

The Bifurcation and Approximation Result
In particular, } B´εCpε¨q η pj 0 q pk 0 qe ik 0¨} C 0 b pR d q ď c 2 ε 9{5 . The constants c 1 and c 2 depend polynomially on } p C} L 1 4 pR d q .
The space L 1 s pR d q for s ě 0 is defined as L 1 s pR d q :" tf P L 1 pR d q : p1`|¨|q s f P L 1 pR d qqu. Clearly, due to B P L 1 pR d q we have the decay Bpxq Ñ 0 as |x| Ñ 8.
The existence of a P T -symmetric solution C is satisfied, e.g., if G : 0 is definite and signpΓq "´signpω 1 q. Due to the extremum of λ j 0 at k " k 0 we have then that D 2 λj 0 pk 0 q is positive/negative definite if Γ is positive/negative 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 signpω 1 q implies We proceed with the proof of Theorem 1. Like in Sec. 2 we work here in Fourier variables.
We decompose for each k P R d the solution p Bpkq P C N into the component proportional to the eigenvector η pj 0 q pkq and the l 2 pC N q´orthogonal complement. We define the projections We aim to construct a solution p B with ψ approximated by the envelope in our ansatz, i.e. by ε 1´d p C`¨´k 0 ε˘. We choose for ψ a decomposition according to the support At the moment r is a free parameter; it will be specified below. We also define For the selected ansatz equation (3.2) becomes Because ω 0 P σpLpik 0 qq, the inverse of the matrix pω 0`ε 2 ω 1 qI´Lpikq is not bounded uniformly in ε. In a neighbourhood of k 0 the norm of the inverse blows up as ε Ñ 0. However, We separate the explicit part of (3.4) by writing where p B Q,1 solves the explicit part, i.e.
The system to solve is thus Our procedure for constructing a solution can be sketched as follows.
(2) For any p D, p R P L 1 pR d q and p B Q,1 from step 1 we solve (3.7) by a fixed point argument for a small p B Q,2 .
(3) For any p D P L 1 pR d q and for p B Q from steps 1 and 2 we solve (3.6) with k P B ε r pk 0 q c for a small p R by a fixed point argument. (4) With the components obtained in the above steps we find a solution p D P L 1 2 pR d q of (3.6) with k P B ε r pk 0 q close to a p C P L 1 2 pR d q (with C a solution of (2.3)) -provided such a C exists. In addition C needs to satisfy a certain symmetry, the P T -symmetry. Also here a fixed point argument is used -roughly speaking for the difference p D´p C.
The estimate for } x N p p B P q} L 1 pR d q follows by Young's inequality for convolutions. l Due to the cubic structure of N we have } v} L 1 pR d q ď αu. The constants c and c 0 depend polynomially on } p D} L 1 and } p R} L 1 . Similary, we obtain the contraction (for ε ą 0 small enough) Hence, for ε ą 0 small enough we have a unique solution For any p D P L 1 pR d q and with the above estimate on p B Q we look for a small p R. The support of p R is B c ε r´1 , whence for k P R d zB ε r pk 0 q we can divide in (3.6) by ω 0`ε 2 ω 1´λj 0 pkq and obtain p Rˆk´k 0 ε˙" ε d´1 νpkq η pj 0 q pkq˚x N p p Bqpkq, k P R d zB ε r pk 0 q with νpkq :" pλ j 0 pkq´ω 0´ε 2 ω 1 q´1.
In order to exploit the localized nature of p B D and the smallness of p B Q , we write for κ :" k´k 0 where (3.10) Finally, we consider the component p B D . For k P suppp p B D q " B ε r pk 0 q we rewrite (3.6) as follows. We add and subtract x N p p B D q like in (3.10), we Taylor expand λ j 0 pkq at k " k 0 , and we use the variable κ " ε´1pk´k 0 q. This leads to δpk 0`ε κq p Dpκq`ε d´1 χ B ε r´1 pκq η pj 0 q pk 0`ε κq˚ˆx N p p Bqpk 0`ε κq´x N p p B D qpk 0`ε κq˙ and δpkq :" λ j 0 pkq´ω 0´1 2 pk´k 0 q T D 2 λ j 0 pk 0 qpk´k 0 q.
Finally, we estimate ρ 3 . Note that x N p p Bqpk 0`ε¨q´x N p p B D qpk 0`ε¨q appears also in (3.10). We have The whole right hand side of (3.12) is thus estimated as q`ε 1´βp1´rq`εαf or any β P r0, 3q. Once again, the constant c depends polynomially on its arguments.
For example, the choice s D " 2, β " 1, r " 1{2 leads to α " 2 and 1´βp1´rq " 1{2, i.e. }ρp p Dq} L 1 pR d q ď cε 1{2 provided p D P L 1 2 pR d q. For s D " 2, β " 1 the largest value of mint1´βp1´rq, αu is 4{5 attained at r " 4{5. Hence, with r " 4{5 we get the estimate (3.14) }ρp p Dq} L 1 pR d q ď cp} p D} L 1 2 qε 4{5 . We return now to equation (3.12), which is a perturbation of the NLS (2.3). Writing We search for a solution D close to a solution C of the NLS, i.e. of f NLS pCq " 0. For that we set and look for a solution with a small p d. In order to obtain a differentiable function (to use the Jacobian of the NLS), we write f NLS in real variables. Writing D " D R`i D I , C " C R`i C I , d " d R`i d I , and ρ " ρ R`i ρ I , we define F NLS pD R , D I q :"ˆR epf NLS pD R`i D I qq Impf NLS pD R`i D I qq˙, the Jacobian (3.16) J :" DF NLS pC R , C I q as well as the Fourier-truncation of the Jacobian With this notation (3.15) reads .
Writing y C R,I pεq ": y C R,I`y γ R,I pεq , it is supppy γ R,I pεq q Ă B c ε r´1 . The difference p p J ε´p J 0 q p d consists of terms that are linear or quadratic in y γ R,I pεq ; for instance terms like Young's inequality for convolutions yields with c depending polynomially on } p C} L 1 , then there is some c ą 0 such that . The loss of 2 in the weight is due to the factor κ T G 0 κ. The entries ω 1`κ T G 0 κ clearly preserve the P T -symmetry. For the convolution terms we have, for instance pεq q ñ C pεq p´xq " C pεq pxq @x. . The term ρ is estimated in (3.14) and dictates the order ε 4{5 . The difference F NLS p p C pεq`p dq´p J ε p d consists of terms quadratic in p d and hence is bounded in L 1 pR d q by c 1 p} p d} 2 is now clear due to the quadratic nature of We conclude that if the solution C of the NLS (2.3) satisfies p C P t p f P L 1 2 pR d q : Imp x C R q " Rep x C I qu, then there is c ą 0 such that for all ε ą 0 small enough the constructed solution D of (3.12) satisfies p D P t p f P L 1 2 pR d q : Imp x f R q "´Rep p f I qu and This allows us to estimate p B D´p B app . We have (3.20), the Lipschitz continuity of η pj 0 q , and the estimate } p for all s C ě 0. This produces at r " 4{5 We can now summarize the error estimate The components p B Q and p B R are estimated in (3.8), (3.9), and (3.11). Having now estimated } p R} L 1 in terms of } p C} L 1 and } p D} L 1 in terms of } p C} L 1 2 , we get for r " 4{5 and s D " 2 , where c 1 and c 2 depend polynomially on } p C} L 1 2 pR d q . Hence, the estimate in Theorem 1 is proved. Remark on the ε-convergence: We expect that the error converges like ε 2 rather than ε 9{5 , see also the numerical test in Sec. 4. In order to have }ρp p Dq} L 1 ď cε, we need to set β " 0 in (3.13), obtaining }ρ 2 p p Dq} L 1 ď cε} p D} L 1 3 . In the fixed point construction of p d we could then work in B 3, sym cε . This, however, requires that p J´1 ε : X sym 0 Ñ X sym 3 , which is not obvious.
Using the numerical Petviashvili iteration [10,1], we also produce a numerical approximation of a solution B at ω " ω 0`ε 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 B app . Note that although B app can be real (if a real solution C of the NLS is chosen), equation (1.3) does not allow real solutions B due to the term i∇ B and due to the realness of α 1 , α 2 , and α 3 and of γ pα,β,γq j . Nevertheless, if B app is real, there must be a solution B with Imp Bq " Opε 9{5 q. Figure 2 shows B app and B for ε " 0.05.
The numerical parameters for the Petviashvili iteration were selected as follows: we compute on the domain x P r´3{ε, 3{εs 2 with the discretization given by 160x160 grid points, i.e. dx 1 " dx 2 " 3{p80εq. Note that because k 0 " 0, the relatively coarse discretization for small values of ε does not matter (there are no oscillations to be resolved). For each ε we evaluate the asymptotic error E :" } B´ B app } C 0 b . Figure 3 shows the convergence of the error in ε. Clearly, Epεq " cε 2 , which suggests that our proof (which produces cε 9{5 ) is suboptimal.