# Bifurcation of Nonlinear Bloch Waves from the Spectrum in the Gross-Pitaevskii Equation

## Abstract

We rigorously analyze the bifurcation of stationary so called nonlinear Bloch waves (NLBs) from the spectrum in the Gross-Pitaevskii (GP) equation with a periodic potential, in arbitrary space dimensions. These are solutions which can be expressed as finite sums of quasi-periodic functions, and which in a formal asymptotic expansion are obtained from solutions of the so called algebraic coupled mode equations. Here we justify this expansion by proving the existence of NLBs and estimating the error of the formal asymptotics. The analysis is illustrated by numerical bifurcation diagrams, mostly in 2D. In addition, we illustrate some relations of NLBs to other classes of solutions of the GP equation, in particular to so called out–of–gap solitons and truncated NLBs, and present some numerical experiments concerning the stability of these solutions.

tomas.dohnal@math.tu-dortmund.de

hannes.uecker@uni-oldenburg.de

## 1Introduction

The nonlinear Schrödinger/Gross–Pitaevskii (GP) equation in dimensions,

with a real potential is a canonical model in mathematics and physics. It appears in various contexts, e.g., nonlinear optics [33], and Bose–Einstein condensation [26]. See also, e.g., [34] for mathematical and modeling background. Plugging into (Equation 1), where is the frequency of time–harmonic waves in nonlinear optics, and where is called the chemical potential in Bose–Einstein condensation, we obtain the stationary problem

Here we consider the case that the potential is real and periodic. For simplicity, we let be periodic in each coordinate direction, i.e.,

where denotes the -th Euclidean unit vector in . In other words, we consider the periodic lattice . We make the basic assumption that for some , where . This smoothness assumption on ensures -smoothness of linear Bloch waves, i.e., solutions of with . See §Section 1.1 for a review of spectral properties of

and linear Bloch waves. For suitable the spectrum of shows so called spectral gaps and in recent years a focus has been on the bifurcation of so called gap solitons from the the zero solution at band edges into the gaps. These are localized solutions, which, in the near edge asymptotics have small amplitude and long wave modulated shape. In detail, the asymptotic expansion at with is

where , are Bloch waves at the edge , and the are localized solutions of a system of (spatially homogeneous) nonlinear Schrödinger equations. See, for instance, [32], and the references therein.

Here we seek solutions of which can be expressed as finite a sum of quasi-periodic functions and call such solutions *nonlinear Bloch waves* (NLBs), with quasi–periodicities determined from a selected finite subset of the Bloch waves at . NLBs have been studied in, for instance, [12], where in [38] the approaches are numerical and formal. They have been observed even experimentally, see e.g. [10] for experiments in Bose-Einstein condensates. In [12] the special case of a bifurcation of NLBs into an asymptotically small spectral gap for a separable periodic potential in two dimensions is studied rigorously. In [9] the bifurcation of single component () NLBs in one dimension is proved, including results on secondary bifurcations and exchange of stability. Similarly to Bloch waves in linear lattices NLBs can be understood as the fundamental bounded oscillatory states of the nonlinear system. From the applied point of view one motivation for studying NLBs is the continuation of gap–solitons to “out-of-gap” solitons, i.e., the continuation of localized solutions from one band edge across the gap and into the spectrum on the other side of the gap, where their tails start interacting with the NLBs. For this reason, the study of bifurcation of NLBs from the zero–solution has been mostly restricted to band edges. Here we show that nonlinear Bloch waves bifurcate in from *generic points in the spectrum* of , and give their asymptotic expansions in terms of solutions of the so called algebraic coupled mode equations (ACME), together with error estimates.

In addition to the rigorous analysis we illustrate our results numerically. For this we focus on 2D, as this is much richer than 1D, and use the same potential as in [15], i.e.

with

This represents a square geometry with smoothed-out edges. The function in is extended periodically to to obtain . The numerical band structure of over the Brillouin zone , and also along the boundary of the irreducible Brillouin zone, is plotted in Fig. ?(a),(b), respectively. We denote the so called high symmetry points in for by

(a) | (b) |

As already said, one motivation for studying NLBs are the intriguing properties of their interaction with localized solutions, which we illustrate numerically in §Section 8. For instance, when a gap soliton is continued from the gap into the spectrum, we get a so called “out–of–gap” soliton (OGS) with oscillating (delocalized) tails, see also [41]. In 1D, numerically these OGS can be seen to be homoclinic orbits approaching NLBs, and essentially the same happens in 2D. Moreover, the NLB can form building blocks of so called truncated NLBs (tNLBs), see also [4]. These are localized solutions for in a gap which are close to a NLB on some finite interval but approach as . Then, continuing a branch of tNLBs from the gap into the spectrum, the same interaction scenario as for GS happens, i.e., the tails of the tNLBs pick up NLBs bifurcating from the gap edge into the spectrum, and the tNLBs become delocalized, for which we use the acronym dtNLB. Note that both gap solitons and tNLBs have been observed experimentally, see e.g. [5] for experiments in optical lattices. However, even in 1D at present it is unclear how to analyze OGS, tNLBs and dtNLBs rigorously, i.e., so far there only exist heuristic asymptotics, see §Section 8 for further comments.

Stability of most of these solutions is an open problem. Thus, at the end of the paper we also give a numerical outlook on this, and obtain stability of some NLBs in 1D, and, consequently, stability of some tNLBs and some OGS and dtNLBs. In 2D, we did not find stable NLBs for the potential from (Equation 4), and we close with summarizing the open questions. The broad spectrum of applications of NLBs clearly motivates our rigorous bifurcation analysis.

In the remainder of this introduction we explain the linear band structure, a simple analytical bifurcation result, formulate the main theorem, and describe the structure of the paper in more detail.

### 1.1Linear Bloch waves

For in the Brillouin zone, , consider the Bloch eigenvalue problem

where is the th Euclidean unit vector in . The spectrum of is continuous and is given by the union of the bands defined by the band structure , i.e.

see Theorem 6.5.1 in [16]. The functions are called band functions. The Bloch waves have the form with for all and all . Clearly, both and are periodic in each component of . We assume the normalization

For a given point in the band structure, i.e. with for some , also the point lies in the band structure, which follows from the symmetry

This symmetry is due to the equivalence of complex conjugation and replacing in the eigenvalue problem . Hence, we also have the conjugation symmetry of the Bloch waves, namely

For we have and the point must be understood as the -periodic image within . When is one of the so called *high symmetry points*, i.e. for all , then and are identified via this periodicity. Equation then implies that is real. This can be seen directly from the eigenvalue problem , where for all implies that the boundary condition is real such that a real eigenfunction must exist. Note that the above symmetries rely only on the realness of .

### 1.2The Bifurcation Problem

Our aim is to prove a general bifurcation theorem for NLBs, and, moreover, to derive and justify an effective asymptotic model related to the Lyapunov-Schmidt reduction of the bifurcation problem including an estimate on the asymptotic error. In our approach we select a frequency in the spectrum and choose points in the level set of the band structure at , such that for each we have for some . Our method requires that each of the points is either one of the high symmetry points or belongs to a pair with . See (H1)–(H6) on page for a summary of our assumptions. We seek NLBs bifurcating from and having the asymptotic form

at . The coefficients , i.e. the (complex) amplitudes of the waves, are given by solving the ACME as an effective algebraic system of equations. Generally a sum of quasiperiodic functions with the quasiperiodicity of each given by one of the vectors cannot be an exact solution of as the nonlinearity generates functions with other quasiperiodicities. Our ansatz for the exact solution is thus

with and for . Importantly, this set (defined in ) can be a proper subset of the level set. The subset may be finite even if the level set is, for instance, uncountable. In fact, our assumption (H4) ensures the finiteness. Besides, the subset can be much smaller than and hence the effective ACME-system can be rather small. The subset has to satisfy only (H2-H6).

The major assumptions of our analysis are rationality (assumption (H4)) and certain non-resonance conditions (H5) on the -vectors . In addition, the solutions of the coupled mode equations need to satisfy certain symmetry (“reversibility”) and non-degeneracy conditions, see Definitions ? and ?, in order for us to guarantee that approximates a solution of .The main result is the following

There are three relatively straightforward generalizations of the result. Firstly, the periodic lattice can be replaced by any lattice with linearly independent vectors . Of course, the periodicity cell and the Brillouin zone have to be redefined accordingly. Except for the examples in §Section 6 the results, in particular Theorem ?, hold for a general lattice. Secondly, the nonlinearity can be replaced by other locally Lipschitz nonlinearities which are phase invariant and satisfy for . This may, however, change the powers of in the expansion and the error estimate. Also, the linear operator can be generalized to self adjoint second order differential operators with periodic coefficients such that the asymptotic distribution of eigenvalues remains that in .

### 1.3The Structure of the Paper

In §Section 2 we present a formal asymptotic approximation of nonlinear Bloch waves and a derivation of the ACMEs as effective amplitude equations. In §Section 3 we pose conditions on the solution ansatz and the band structure which are necessary for our analysis, and apply the Lyapunov-Schmidt decomposition to the bifurcation problem. The invertible part of the decomposition is estimated in §Section 4. The singular part and its relation to the ACMEs is described in §Section 5, where also the proof of the main theorem is completed. In §Section 6 we present the ACMEs and their solutions in the scalar case () and in the case of two equations (). Section 7 presents numerical computations of nonlinear Bloch waves in two dimensions for and . The convergence rate of the approximation error is confirmed by numerical tests. Finally, in §Section 8 we give a numerical outlook on the interaction of localized solutions with NLBs, first for some 1D and 2D GS, and second for tNLBs, and we report numerical experiments on stability of NLBs and other solutions.

## 2Formal Asymptotics

Let and choose vectors in the level set of the band structure at . For the asymptotics of nonlinear Bloch waves near we make an analogous ansatz to that used in [12] for nonlinear Bloch waves near band edges in (Equation 2) with a *separable* periodic potential. Formally we write

where the amplitudes are to be determined and where satisfies the quasiperiodicity given by the vector .

Substituting in we get at for each

where

The condition in the sum ensures that the nonlinear terms have the same quasi-periodicity as . Nonlinear terms generated by the ansatz and having other quasi-periodicity than one of those defined by have been ignored in this formal calculation.

Imposing the solvability condition, i.e. making the right hand side -orthogonal to , we get the *algebraic coupled mode equations* (ACMEs)

To make the approximation rigorous, we must account for the nonlinear terms left out above and provide an estimate on the correction .

## 3Solution Ansatz, Assumptions, Lyapunov-Schmidt Decomposition

As mentioned above, one of the difficulties of the analysis is that for a sum of functions with distinct quasi-periodic conditions the nonlinearity can generate functions with a new quasi-periodicity. If the -points defining these new quasi-periodic boundary conditions lie in the -level set of the band structure, then a resonance with the kernel of the linear operator occurs. Also, if the points generated by a repeated iteration of the nonlinearity merely converge to the level set, our techniques fail because a lower bound on the inverse of the linear operator cannot be obtained. These obstacles are avoided if for a selected assumptions (H4) and (H5) below hold.

We select points in the -level set of the band structure. Suppose we seek solutions of with given by the sum of quasiperiodic functions. The ansatz with quasiperiodic such that for all can be a solution of only if each term generated by the nonlinearity applied to this sum has quasiperiodicity defined by one of the vectors in , i.e. if the *consistency condition*

where

is satisfied. In other words the consistency condition says that all combinations for must lie in , with from (Equation 10).

An example of a consistent ansatz for is with , like e.g. for in [15]. On the other hand, for , where with , see Sec. 3.2.2.5 in [15], the ansatz is inconsistent. It is also inconsistent for typical in the interior of with generic in the level set. Therefore, we drop the consistency condition and pursue the more general case where the nonlinearity generates quasiperiodic functions with quasi-periodicity vectors not necessarily contained in .

Hence, we define the set of -points generated by iterations of the nonlinear operator

and write, with ,

At this point is possible but as explained below, our assumption (H4) ensures , i.e. only finitely many new vectors are generated. Thus we can search for a solution in the form of the sum of finitely many quasiperiodic functions

with . The choice of the function space for is made clear below.

We make the following assumptions:

for some , where ;

and are points in the -level set of the band structure, i.e., there are such that

each point is repeated according to the multiplicity of at . In detail, if band functions touch at , then points in equal and ;

the points have rational coordinates, i.e.

the intersection of the set with the level set of the band structure at is exactly the set , i.e.

where

for each the reflection w.r.t. the origin lies in the set too, i.e.

where and denotes congruence with respect to the -periodicity in each component.

With (H3) the bifurcation from multiple Bloch eigenvalues is allowed. In one dimension () multiplicity is at most two, which occurs for so called finite band potentials, see e.g. [27], and only at or . In higher dimensions () crossing or touching of band functions is abundant in generic geometries (potentials ). Our results thus apply also to Dirac points in two dimensions studied, e.g., in [18].

Due to the rationality condition (H4) the sought solution is, in fact, periodic. (H4) also ensures that the set is finite (). Indeed, iterating the operator on a set of points with rational coordinates on a -dimensional torus generates a periodic orbit, i.e. only finitely many distinct points are generated, and the number depends solely on . Condition (H4) is satisfied, e.g. if is a subset of the high symmetry points of , i.e. for all . This is frequently the case for the locations of extrema defining a spectral edge. In general, (H4) is, however, a serious limitation, and removing this assumption would be a major improvement.

The non-resonance condition (H5) is satisfied, for instance, if , i.e. for at one of the band edges, and are all the extremal points of the band structure at which the edge is attained.

The symmetry condition (H6) is needed in the persistence step of the proof, see §Section 5. Note that if , then also by and the periodicity in . For , clearly, . For is (e.g. for with ) and then is the -periodic image of within (for the example ). Moreover, also (H6) is automatically satisfied if is a subset of the high symmetry points of .

Note again that as well as may be proper subsets of . This is, for instance, the case in 2D if , where we could choose (if is simple), and or , which yields two decoupled scalar bifurcation problems; see §Section 6.2 for further discussion.

The remaining two assumptions in Theorem ? are non-degeneracy and reversibility of , defined as follows.

Note that due to the phase invariance of the Jacobian is singular.

Reversibility is a symmetry of the solution. The motivation for restricting to reversible non-degenerate solutions is to ensure the invertibility of in the fixed point iteration for the singular part of the Lyapunov-Schmidt decomposition in §Section 5. Within the phase invariance is, indeed, no longer present. The choice of in Definition ? is natural and based on the intrinsic symmetry of the Bloch eigenfunctions which ensures the complex conjugation symmetry among the coefficients in and, hence, the possibility of reversible solutions. Note that follows directly from .

Next, we assume (H1-H6) and use the Lyapunov-Schmidt decomposition in Bloch variables together with the Banach fixed point theorem to prove the main result, i.e., Theorem ?, which justifies the formal asymptotics for solutions at .

### 3.1Lyapunov-Schmidt Decomposition

Due to the completeness of the Bloch waves in we can expand

As the following lemma shows, working with in the space is equivalent to working with , where

The proof is analogous to that of Lemma 3.3 in [8], see also [15]. The main ingredients are firstly the fact that for large enough (such that for all and , e.g. ) the squared norm is equivalent to

Secondly, one uses the asymptotic distribution of bands in dimensions, see [21]: there are constants such that

For the subsequent analysis we define for each the set of indices producing through the nonlinearity analogously to the definition of , i.e.

For the ansatz , equation is equivalent to the algebraic system

where

Due to the kernel of the linear multiplication operator at in we use a *Lyapunov-Schmidt decomposition* in order to characterize the bifurcation from (i.e. from ). For we let

and write

with , and . In other words we set

where is the -th Euclidean unit vector in . Analogously to we also define

This decomposition splits problem into

The following program is analogous to that in [12]. Namely, for given, we first show the existence of a small solution of the regular part and then prove a persistence result relating certain (reversible and non-degenerate) solutions of to solutions of the singular part including an estimate on their difference, and finally provide an estimate of .

## 4Regular Part of the Lyapunov-Schmidt Decomposition

We define the following spaces and norms

Note that the condition in the definition of can be replaced by because each can be written as , where and . Also note that is a sequence of sequences. Similarly we denote

Clearly, the ansatz satisfies if and only if for all . Therefore, for the problem at hand, where the solution consists of components , the spaces and could be defined with finite sums over . However, since the use of infinite sums in the definitions does not increase the complexity and since it may prove useful in future work on the case of irrational coordinates of , we keep these general definitions.

We will need the following two lemmas, the first following directly from Lemma ?.

*Proof.* We define the sets and of points, which determine the quasiperiodicity of the functions and , i.e.

We have

where the inequality follows by the triangle inequality and by the algebra property of the norm

see Theorem 5.23 in [2].

Our result on the regular part of the Lyapunov-Schmidt decomposition is the following

*Proof.* Writing in the fixed point formulation

we seek a fixed point with . Lemma ? allows us to work interchangeably in in the physical variables. We show the contraction property of within

for some .

The nonlinearity is