q-breathers in Discrete Nonlinear Schroedinger lattices

$q$-breathers are exact time-periodic solutions of extended nonlinear systems continued from the normal modes of the corresponding linearized system. They are localized in the space of normal modes. The existence of these solutions in a weakly anharmonic atomic chain explained essential features of the Fermi-Pasta-Ulam (FPU) paradox. We study $q$-breathers in one- two- and three-dimensional discrete nonlinear Sch\"{o}dinger (DNLS) lattices -- theoretical playgrounds for light propagation in nonlinear optical waveguide networks, and the dynamics of cold atoms in optical lattices. We prove the existence of these solutions for weak nonlinearity. We find that the localization of $q$-breathers is controlled by a single parameter which depends on the norm density, nonlinearity strength and seed wave vector. At a critical value of that parameter $q$-breathers delocalize via resonances, signaling a breakdown of the normal mode picture and a transition into strong mode-mode interaction regime. In particular this breakdown takes place at one of the edges of the normal mode spectrum, and in a singular way also in the center of that spectrum. A stability analysis of $q$-breathers supplements these findings. For three-dimensional lattices, we find $q$-breather vortices, which violate time reversal symmetry and generate a vortex ring flow of energy in normal mode space.


A. Background
In 1955, Fermi, Pasta and Ulam (FPU) published their seminal report on the absence of thermalization in arrays of particles connected by weakly nonlinear springs [1]. In particular they observed that energy, initially seeded in a low-frequency normal mode of the linear problem with a frequency ω q and corresponding normal mode number q, stayed almost completely locked in a few neighbor modes in frequency space, instead of being distributed quickly among all modes of the system. The latter expectation was due to the fact, that nonlinearity does induce a long-range network of interactions among the normal modes. ¿From the present perspective, the FPU problem appears to consist of the following parts: (i) for certain parameter ranges (energy, system size, nonlinearity strength, seed mode number), excitations appear to stay exponentially localized in q-space of the normal modes for long but finite times [2]; (ii) this intermediate localized state is reached on a fast time scale τ 1 , and equipartition is reached on a second time scale τ 2 , which can be many orders of magnitude larger than τ 1 [3]; (iii) tuning the control parameters may lead to a drastical shortening of τ 2 , until both time scales merge, and the metastable regime of localization is replaced by a fast relaxation towards thermal equilibrium [3], which is related to the nonlinear resonance overlap estimate by Izrailev and Chirikov [4]. FPU happened to compute cases with values of τ 2 which were inaccessible within their computation time. Among many other intriguing details of the evolution of the intermediate localized state, we mention the possibility of having an almost regular dynamics of the few strongly excited modes, leading to the familiar phenomenon of beating. That beating manifests as a recurrence of almost all energy into the originally excited mode (similar to the beating among two weakly coupled harmonic oscillators), and was explained by Zabusky and Kruskal with the help of real space solitons of the Kordeweg-de-Vries (KdV) equation, for a particular case of long wavelength seed modes, and periodic boundary conditions [5]. The regular dynamics may be replaced by a weak chaotic one upon crossing yet another threshold in control parameters [6]. Re-markably, that weak chaos is still confined to the core of the localized excitation, leaving the exponential localization almosts unchanged, and perhaps influencing most strongly the parameter dependence of τ 2 . The interested reader may also consult Ref. [7].

B. Motivation
Single-site excitations in translationally invariant lattices of interacting anharmonic oscillators are known to show a similar behaviour of trapping the excitation on a few lattice sites around the originally excited one [8]. As conjectured almost 40 years ago [9], exact time-periodic and spatially localized orbitsdiscrete breathers (also intrinsic localized modes, discrete solitons) persist in such lattices, and the dynamics on such an orbit, as well as the dynamics of the nearby phase space flow, account for many observations, and have found their applications in many different areas of physics [10]. Recently, it was shown, that the regime of localization in normal mode space for the FPU problem can be equally explained by obtaining q-breathers (QB) -time-periodic and modal-space-localized orbits -which persist in the FPU model for nonzero nonlinearity [11]. The originally computed FPU trajectory stays close to such q-breather solutions for long times, and therefore many features of its short-and medium-time dynamics have been shown to be captured by the q-breather solution and small phase-space fluctuations around it. Delocalization thresholds for q-breathers are related to resonances, and corresponding overlap criteria [11]. The method of constructing q-breather solutions was generalized to two-and three-dimensional FPU lattices [12]. Most importantly, a scaling theory was developed, which allowed to construct q-breathers for arbitrarily large system sizes, and to obtain analytical estimates on the degree of localization of these solutions for macroscopic system size [13]. Expectations were formulated, that the existence of q-breathers should be generic to many nonlinear spatially extended systems.

C. Aim
The spectrum of normal mode frequencies of the linear part of the FPU model is acoustic, contains zero, reflecting the fact, that the model conserves total mechanical momentum. In addition it is bounded by a lattice-induced cutoff, the analogue of the Debye cutoff in solid state physics. We may think of two pathways of extending the above discussed results to other system classes. First, we could consider a spatially continuous system, instead of a lattice. However, the initial FPU studies showed confinement to a few long wavelength modes, and the results of Zabusky and Kruskal [5], confirm, that indeed in the spatially continuous system (KdV), trajectories similar to the FPU trajectory persist.
As discussed above, the persistence, localization, and delocalization of q-breathers is due to a proper avoiding of, or harvesting on, resonances between normal modes. It therefore appears to be of substantial interest to extend the previous results to systems with a qualitatively different normal mode spectrum. Such a spectrum is one which does not contain zero, and is called optical. It may correspond to the excitation of certain degrees of freedom in a complex system, or to an FPU type model which is deposited on some substrate, or otherwise exposed to some external fields. Together with the acoustic band, these two band structures, when combined, describe almost any type of normal mode spectrum in a linear system with spatially periodic modulated characteristics.
Adding weak nonlinearities to such a system, taking local action angle representations, and performing various types of multiple scale analysis, such models are mapped on discrete nonlinear Schrödinger models (DNLS) [14]. These models enjoy a gauge invariance, and a conservation of the sum of the local actions (norms), i.e. some global norm. A particular version of these equations is known as the discrete Gross-Pitaevsky equation, and is derived on mean field grounds for Bose-Einstein condensates of ultracold bosonic atoms in optical lattices [15,16]. The norm is simply the conserved number of atoms in this case, and the nonlinearity derives from the atom-atom interaction. The propagation of light in (spatially modulated) optical media is another research area where DNLS models are used [17]. In that case the norm conservation derives from the conservation of electromagnetic wave energy along the propagation distance (in the assumed absence of dissipative mechanisms). But DNLS models served equally well in many other areas of physics, where weak interactions start to play a role, e.g. in the theory of polaron formation due to electron-phonon interaction in solids.
Releasing the norm conservation (e.g. by allowing the number of atoms in a condensate to fluctuate) will reduce the symmetries of the corresponding model, but most importantly, it will lead to possible new resonances of higher harmonics. We will discuss these limitations, and possible effects of releasing these limitations on q-breathers, in the discussion section.
A particular consequence of norm conservation is a corresponding symmetry of the normal mode spectrum, which relates both (upper and lower) band edges, and makes the band center a symmetry point as well.
The paper is structured as follows. Section II is introducing the model and the main equations of motion. Section III gives an existence proof for q-breathers. In section IV we derive analytical expressions for the QB profiles using perturbation theory. These are compared with the numerical results for d = 1 in section V, which also contains a stability analysis of the obtained periodic orbits. Section VI gives results on periodic boundary conditions (as opposed to fixed boundaries). We generalize to d = 2, 3 in section VII, and discuss all results in section VIII.

II. MODEL
We consider a discrete nonlinear Schrödinger (DNLS) equation on a d-dimensional hypercubic lattice of linear size N: Here ψ is a complex scalar which may describe e.g. the probability amplitude of an atomic cloud on an optical lattice site [16], or relates to the amplitudes of a propagating electromagnetic wave in an optical waveguide [17]. The lattice vectors n, m have integer components, and D(n) is the set of nearest neighbors for the lattice site n. If not noted otherwise, we consider fixed boundary conditions: ψ n = 0 if n l = 0 or n l = N + 1 for any of the components of n. Equation (1) is derived from the Hamiltonian using the equations of motion iψ n = ∂H/∂ψ * n . In addition to energy, equation (1) conserves the norm B = n |ψ n | 2 . Note that the change of the nonlinearity parameter µ in (1) is strictly equivalent to changing the norm B. Here we will keep the norm B (alternatively the norm density) fixed, and vary µ.
We perform a canonical transformation to the reciprocal space of normal modes (q-space Together with (1) we obtain the following equations of motion in normal mode space: where cos πq i N +1 are the normal mode frequencies for the linearized system (1) with µ = 0. Nonlinearity introduces a network of interactions among the normal mode oscillators with the following coupling coefficients:

III. PROOF OF EXISTENCE OF t-REVERSIBLE q-BREATHER SOLUTIONS FOR WEAK NONLINEARITY
We look for exact time-periodic solutions, which are stationary solutions of the DNLS equation (1), and which are localized in normal mode space: ψ n (t) = φ n exp(iΩt) with frequency Ω and time-independent amplitudes φ n . In the space of normal modes these stationary solutions have the form Q q (t) = A q exp(iΩt), where the amplitudes of modes A q are time-independent and related to the real-space amplitudes by the transformation (3). At a given norm B they satisfy a system of algebraic equations: We are focusing here (and throughout almost all of the paper) on t-reversible periodic orbits. Therefore we may consider all A q to be real numbers. In this case system (6) contains N +1 equations for N + 1 variables. This system can be condensed into an equation for a vector function: with X = {. . . , A q , . . . , Ω}. The components of F are the left hand sides of (6), while µ, B are parameters. For µ = 0 the normal modes in (4) are decoupled and each oscillator conserves its norm in time: B q (t) = A 2 q . Let us consider the excitation of only one of the oscillators with the seed mode number q 0 : B q = B q 0 δ q,q 0 . The excited normal mode is a time-periodic solution of (4), and is localized in q-space. According to the Implicit Function Theorem [18], the corresponding solution of (6) can be continued into the nonlinear case (µ = 0), if the Jacoby matrix of the linear solution (∂F /∂X) µ=0,B is invertible, i.e. ||∂F /∂X|| µ=0,B = 0. The Jacobian Invertibility of the Jacoby matrix requires a non-degenerated spectrum of normal mode frequencies. Therefore, the continuation of a linear mode with the seed mode number q 0 will be possible if ω q = ω q 0 , for all q = q 0 . That condition holds for the case d = 1, and thus q-breathers exist at least for suitably small values of µ there. For higher dimensions, degeneracies of normal mode frequencies appear. These degeneracies are not an obstacle for numerical continuation of q-breathers, but a formal persistence proof has to deal with them accordingly. Analogous results for two-and three-dimensional β-FPU lattices have been obtained in Ref. [12].

IV. PERTURBATION THEORY FOR q-BREATHER PROFILES
To analyze the localization properties of q-breathers in q-space we use a perturbation theory approach similar to [11], [12]. We consider the general case of a d-dimensional DNLS lattice (1). Taking the solution for a linear normal mode with number q 0 = (q 1,0 , . . . , q d,0 ) as a zero-order approximation, an asymptotic expansion of the solution to the first N equations of (6) in powers of the small parameter σ = µ2 (d−2) /(N + 1) d is implemented. Note, that in the same way as it was done in [11], [12], we fix the amplitude of the seed mode A q 0 . Later on we will apply the norm conservation law (the last equation of (6)) to express the norm B q 0 = |A q 0 | 2 via the total norm B for the case of d = 1. Analytical estimations presented below describe amplitudes of modes located along the directions of the lattice axes starting from the mode q 0 . Mode amplitudes of q-breathers have the slowest decay along these directions. Studying q-breather localization along a chosen dimension i, for the sake of compactness we will use a scalar mode number to denote the i-th component of q, assuming all other components be the same as in q 0 : q j =i = q j,0 .

A. Close to the band edges
Let us start with q-breathers localized in the low-frequency mode domain (q i,0 << N, i = 1, d). According to the selection rules, if the seed mode number q i,0 is even (odd), then only even (odd) modes are excited along i-th dimension. The n-th order of the asymptotic expansion is the leading one for the mode q i,n = (2n + 1)q i,0 . When q i reaches the upper band edge in the i-th direction, a reflection at the edge in q-space takes place if N + 1 is not divisible by q i,0 (cf. Fig. 1d). If N + 1 is divisible by q i,0 then only modes q i,n = (2n + 1)q i,0 < N are excited. In the analytical estimates we assume a large enough lattice size. Then the effect of band edge reflections is appearing in higher orders of the perturbation, which will not be considered. In this case the approximate solution is The corresponding exponential decay of mode norms is < π is the i-th component of the seed wave vector k 0 , and b k 0 = B k 0 /(N + 1) d is the norm density of the seed mode. When using intensive quantities -wave number k i,0 and norm density b k 0 -only, λ does not depend on the system size.
Equations (4) which changes the sign of nonlinearity, and maps modes from one band edge to the other. The replacement q i,0 → N + 1 − q i,0 for all i = 1, d in (6) is equivalent to substitutions µ → −µ and Ω → −Ω. Using this symmetry, we can easily apply the above results to the case of q-breathers localized near the upper band edge, by counting mode indices from the upper edge: q i = N + 1 − q i . We neglect reflections from the lower band edge, such that only modes with numbers q i,n = (2n + 1) q i,0 (n is integer, q i,0 << N) are assumed to be excited. It follows where n is an integer and The analytical expression for the frequency of the q-breather solution turns out to be the same as for the case of small seed mode numbers.

B. Close to the band center
We implement the perturbation theory approach for seed modes close to the band center. Let us first consider the case of odd N. We introduce the new index p i , which is the number of a mode counted from the middle of the spectrum along the i-th dimension: p i = q i − (N + 1)/2. We choose the seed mode with |p i,0 | << (N + 1)/2. In the n-th order of perturbation theory the newly excited mode along the i-th dimension has the number p i,n = (−1) n (2n + 1)p i,0 , and When using intensive quantities, λ again does not depend on the system size.
For the case of even N we use p i = q i −N/2 and assume that the seed mode index p i,0 > 0, i = 1, d. The set of consecutively perturbed modes becomes more complicated: For an even number of perturbation steps the new excited mode along the i-th dimension has index p i,2n = (4n + 1)p i,0 − 2n, and for an odd number of steps it has index p i,2n+1 = −(4n + 3)p i,0 + 2n + 2, n = 0, 1, 2, . . .. The amplitudes satisfy (13), but with For large enough N equation (16) approaches the expression (14), therefore we will use (14) for large N only.

A. Close to the band edge
The key property of q-breathers is that they are localized in the space of linear normal modes. Note that some q-breather solutions may be compact in q-space and contain only one seed mode q 0 [19], or a few modes additionally, due to symmetries of the interaction network spanned by the nonlinear terms [20].  Let us consider the one-dimensional case. We compute q-breathers as the stationary solutions of the nonlinear equations (6) using the single-mode solution for µ = 0 as an initial approximation. Fig. 1a,b shows the distribution of mode norms for q-breathers in the space of wave numbers k = πq/(N + 1) for different chain sizes and different seed wave numbers, located near the lower (k 0 << π, Fig. 1a) and upper (π − k 0 << π, Fig. 1b) edges of the linear mode spectrum.
We find exponential localization of q-breathers, with a localization length which depends strongly on the chosen parameters. The obtained analytical estimations for q-breathers localized at the lower band edge (9),(10) respectively upper band edge (11), (12), are in quantitative agreement with numerical results (see dashed lines in Fig. 1a,b).
The exponential decay in (10) is depending on the seed mode norm density. Many applications (e.g. cold atoms in a condensate) rather fix the total norm, or total norm density. For the obtained q-breather solutions with q 0 << N, the relation between these quantities can be estimated by the sum of infinite geometric series: where λ ≡ λ 1 . From (17) it follows that b k 0 = (1 − λ)b, where b = B/(N + 1) is the total norm density. Substituting the expression for b k 0 into (10) and solving the equation for √ λ we obtain: The same dependence of √ λ on k 0 was obtained for low frequency q-breathers in the β-FPU model [13]. The exponential decay of mode norms in the space of wave numbers, can be now written for q-breathers localized in low-frequency modes: To characterize the degree of localization in k-space, we use the slope of the profile of mode norms in log-normal plots (19) -S [13], where the absolute value of S is equivalent to the inverse localization length ξ: Substituting the expression for √ λ (18) into (20) we obtain: Therefore the slope (inverse localization length) is a function of the rescaled wavenumber z. It therefore parametrically depends on just one effective nonlinearity parameter ν, which is given by the product of the total norm density and the absolute value of nonlinearity strength. S vanishes for z → 0, and it has an extremum min(S) ≈ −0.7432/ν at z min ≈ 2.577. The wave number k 0 = k min ≡ z min ν corresponds to the strongest localization of a q-breather with fixed effective nonlinear parameter ν. With increasing ν the localization length increases. Most importantly, the localization length diverges for small k 0 ≪ k min since |S| ≈ z/(2ν) in that case. This delocalization is due to resonances with nearby normal modes close to the band edge. Note, that our analytical estimations do not depend on the sign of the nonlinearity parameter µ. In Fig. 1a,b, due to the small sizes of the chain, all values of k 0 are greater than k min , so a monotonous dependence of the slope on k 0 is observed. In Fig. 2 we plot theoretical and numerically obtained dependencies S(k 0 ) for different values of the nonlinear parameter µ and different system sizes N (we use large enough N to resolve the theoretically predicted extremum of S). In all cases the numerical results show, that the slopes indeed are characterized by intensive quantities only, and the above derived scaling laws hold.
For positive values of µ and seed wave numbers close to the lower band edge, the extremum in S is reproduced in the numerical data, though the numerical curves deviate from theoretical curves for small k 0 . We also find that k min increases and |S(k min )| decreases with increasing µ as it is predicted by our analytical results (increase of norm density b gives the same effect).
For negative values of µ and seed wave numbers close to the lower band edge, we do not observe an extremum for S. In the region of small k 0 the numerically obtained curves for positive and negative values of µ differ, while the theoretically predicted slopes do not depend on the sign of µ. Fig. 3 shows the mode norm profiles of q-breather solutions with small k 0 obtained for positive and negative values of µ. The reason for the discrepancy between theoretical prediction and numerical results for negative values of µ must be strong contributions from higher order terms in the perturbation expansion. The standard argument is, that the perturbation theory is valid if the localization length is small, i.e. |λ| ≪ 1 as well. When |λ| becomes of the order of one, higher order terms in perturbation theory have to be taken into account, and it would be tempting to conclude that delocalization will take place. That should be true especially when all higher order terms in the series carry the same sign. That is the case for positive nonlinearity here, but for negative µ we obtain alternating signs of higher order terms. These alternating signs therefore effectively  cancel most of the terms in the series, and are responsible for an increasing localization of q-breathers with negative µ in the limit of small wavenumbers. The obtained results for q-breathers with seed wave numbers close to the lower band edge (k 0 << π) are valid for the case of k 0 close to the upper band edge (π − k 0 << π) if we change µ → −µ. Thus, there is a strong asymmetry in the localization properties of q-breather solutions with k 0 < π/2 and k 0 > π/2 for a fixed sign of nonlinearity.
The numerical results in Fig. 2 show that slope values calculated for different system sizes lie on the same curves with corresponding µ even in the region of small k 0 , where higher order corrections to our analytical estimates have to be taken into account. This result is in agreement with the exact scaling of q-breather solutions described in [13]. We plot in Fig. 4 the master slope function S m (z) = νS, which depends on a single variable z [13]. It implies that knowing this single master slope function is sufficient to predict the localization property of a q-breather at any seed wave number k 0 ≪ π, at any energy etc. Numerically obtained slopes presented in Fig. 2 are rescaled and plotted in Fig. 4. We see, that all results corresponding to the same sign of µ condense on a single curve even (and especially) for small k 0 , though these numerical results differ from the analytical estimation.
B. Close to the center of the band Fig. 1c shows the distribution of mode norms for q-breathers in the space of wave numbers k = πq/(N + 1) for different chain sizes and different seed wave numbers, located close to the center of the spectrum (|k 0 − π/2| << π). The analytical estimation for the amplitudes of these q-breathers are in good quantitative agreement with the numerical results, for small enough parameters µ and b, cf. Fig. 1c.
We express again the norm density of the seed mode via the total norm density: 1 ). Substituting this expression into (14) we find: The slope of the mode norm profile in k-space is given by where z = |k 0 − π/2|/(2ν). For z << 1: S = (−1 + z 2 /6 + O(z 4 ))/(2ν), for z >> 1: In the limit z → 0, the slope S → −1/(2ν). Thus, the strongest localization of a q-breather with fixed effective nonlinear parameter ν should be obtained for wave numbers k 0 ≈ π/2. The increase of the effective nonlinearity parameter ν ∼ |µ|b leads to a weaker localization of q-breathers in k-space, since the absolute value of the slope S decreases. Note, that there is a special point for odd N: k 0 = π/2, q 0 = (N + 1)/2. For this seed mode, the q-breather is compact in k-space [19]. In Fig. 5 the theoretical and numerically obtained dependencies S(k 0 ) for different values of µ and different system sizes N are plotted. For small values of µ we observe good agreement between analytical and numerical results. But the increase of the nonlinearity leads to a deviation between the theoretical and numerical curves in the region of k 0 close to π/2. These corrections, as it was for the case of q-breathers localized near the band edges, depend on the location of k 0 (k 0 > π/2 or k 0 < π/2) and on the sign of nonlinearity. Therefore the curves of S(k 0 ) for strong nonlinearity (µ = 30) in Fig. 5 are non-symmetric around the point k 0 = π/2: in contrast to k 0 < π/2, for k 0 > π/2 a local minimum of S is observed. Still, the predicted scaling properties of q-breathers remain correct even for strong nonlinearity: the values of S, computed for different system sizes N, lie on the same curves for fixed µ.

C. Stability of q-breathers
We analyze the linear stability of q-breathers as stationary solutions of DNLS model considering the evolution of small perturbations ε n in the rotating frame of the periodic solution [22]: ψ n (t) = (φ 0 n + ε n (t)) exp(iΩt), where φ 0 n are the non-perturbed time-independent amplitudes, ε n = α n + iβ n . A periodic orbit is stable when all perturbations do not grow in time. Solving the linearized equations for the perturbation, that condition translates into the request, that all eigenvalues s m (m = 1, 2N d ) of the linearized equations must be purely imaginary. Otherwise the orbit is unstable. In Fig. 6  stability analysis. We smoothly continue q-breather solutions for each seed mode number q 0 by increasing the nonlinearity parameter µ and check the stability. If the maximum absolute value of the real parts of all eigenvalues is smaller than 10 −6 , the q-breather is considered as stable (solid line), otherwise it is unstable (dashed line), crosses mark the change of stability. For small values of nonlinearity all q-breathers are stable. Qualitatively different threshold values and dependencies on q 0 for q-breathers with seed modes from different parts of the spectrum q 0 < N/2 and q 0 > N/2 are observed. This is in agreement with the results presented in [21], where stability properties of nonlinear standing waves are studied. Note, that due to the symmetry of the equations, the stability properties of a q-breather with seed mode q 0 for some negative value of the nonlinearity parameter µ is the same as the stability property of the q-breather with seed mode q 0 = N + 1 − q for the nonlinearity parameter µ = −µ. Let us discuss the possible link between linear stability and localization. If a q-breather becomes delocalized, that happens because of resonances between different mode frequencies. Therefore we can expect, that the same resonances will also drive the state unstable. Indeed, these correlations can be clearly observed from the numerical data. However, if a q-breather is well localized, it does not follow that it will be stable as well, since instability can arise due to resonant interaction of modes in the breather core alone.

VI. PERIODIC BOUNDARY CONDITIONS
In the case of periodic boundary conditions, we have used the following transformations between real space and the reciprocal space of normal modes (for even N): ψ n (t) = 1 (24) The DNLS model (1) with periodic boundary conditions has exact solutions for nonlinear traveling waves, which can be written, for instance in case d = 1, as: ψ n (t) = φ 0 exp(iΩt − ik 0 n), where Ω = −2 cos k 0 −µφ 2 0 , k 0 = 2πq 0 /N, q 0 ∈ [−N/2, N/2]. These types of solutions can be considered as compact q-breathers which contain only one mode q 0 . Traveling modes are also not invariant under time reversal. The continuation of a linear standing wave, consisting of two traveling waves with the same norms and wave numbers: k 0 and −k 0 , into the nonlinear regime leads to a time-reversible q-breather solution (see Fig. 7), which is not compact, and its localization properties are similar to the properties of q-breathers in the case of fixed boundary conditions. Here we present the result for decay of mode norms λ (i) d for the case |k i,0 | << π, which differs from (10) by a prefactor: where b k 0 = B k 0 /N d , k i,0 = 2πq i,0 /N. Fig. 7 illustrates good quantitative agreement between analytical and numerical results obtained for small enough values of norm and nonlinearity.

VII. q-BREATHERS IN TWO-AND THREE-DIMENSIONAL LATTICES
For two-and three-dimensional symmetric DNLS lattices (1) with fixed boundary conditions, only linear modes with mode numbers q on the main diagonal have non-degenerate frequencies. Using the Implicit Function Theorem [18], it follows that these modes are continued into the nonlinear regime. However, we also successfully performed numerical continuations of q-breathers with seed mode numbers off the main diagonal as it was done in [12] for the FPU model. For the two-dimensional DNLS model the q-breather, continued from the single linear mode q 0 = (2, 3), is presented in Fig. 8 in real space (a) and in q-space (b). The q-breather with the seed mode q 0 = (3, 2) exists as well and has the same frequency. The slowest decay of the mode norms happens to be along the direction of the main axes. The decay is exponential and in good agreement with the analytical estimation (10) for d = 2. In addition to such asymmetric single mode q-breathers it is possible to construct multi-mode q-breathers continued from a pair of degenerate linear normal modes (q 0 , p 0 ) and (p 0 , q 0 ) (q 0 = p 0 ) with the same norms in both modes and in-phase (Fig. 8c) or antiphase (Fig. 8d) oscillations. It is important to note that the problem of degenerate frequencies is avoided for these solutions. Indeed, system (4) has two invariant manifolds Q q 1 ,q 2 = ±Q q 2 ,q 1 . Looking for a solution on a manifold, the number of independent variables of state is reduced from N 2 to the dimensionality of the manifold, which equals (N 2 + N)/2 for the symmetric manifold and (N 2 −N)/2 for the antisymmetric one. The reduced system of equations contains only modes with non-degenerate frequencies.
For d = 3 we have also verified, that the analytical estimations of mode norm decay (10), (14) agree well with the results of numerical calculations for single-mode q-breathers.
In addition to various time-reversible q-breather solutions, which are constructed in the same way as for d = 2, the three-dimensional DNLS model allows also for non-timereversible ("vortex") multi-mode solutions. Let us consider q-breather solutions on an invariant manifold of the system (4): Q q 1 ,q 2 ,q 3 = exp(2πi/3)Q q 3 ,q 1 ,q 2 = exp(4πi/3)Q q 2 ,q 3 ,q 1 (note, that this manifold has a counterpart with the opposite sign of the phase shifts). On the manifold the number of variables of state is reduced to (N 3 −N)/3. We have constructed numerically a vortex q-breather solution continued from a degenerate triplet of seed modes which have the same norm and a relative phase shift 2π/3: exp(4πi/3)Q q 0 ,p 0 ,q 0 , q 0 = p 0 . The frequency of such a triplet becomes non-degenerate in the reduced system on the manifold. The energy flows in a vortex-like manner in q-space for such excitations (Fig.9).

VIII. DISCUSSION
Comparing the aims we formulated in the introduction with the above results, we can conclude, that indeed the q-breather concept turns out to be generic, and time-periodic orbits, which are localized in normal mode space, are probably as generic as discrete breathers (although perhaps in a different way). Along with the study of the details of q-breathers properties in DNLS models, we found some particularities, which were not (yet) obtained in acoustic (FPU) systems. Let us discuss some of these and other findings from above. A. Scales, delocalization thresholds, and healing length of a BEC Let us consider a finite DNLS chain with N sites. Let us fix the total norm density b and the nonlinearity µ. Therefore we fix the effective nonlinearity parameter ν 2 = |µ|b/16 (18). If the system size was large enough, we will resolve the extremum in S, and therefore the QB with seed wave number k min ≈ 2.577ν is the strongest localized one. It evidently sets an inverse length scale ν. This length scale is known in the GP equation for a BEC in a trap with (similar) fixed boundary conditions. One has a condensate, whose amplitude vanishes at some boundary, yet getting back to some mean value away from the boundary -exactly at the healing length ξ h = (4πan) −1/2 , where n is the condensate density, and a is the scattering length which is proportional to the atom-atom interaction, and therefore to the nonlinearity strength µ within the mean field GP equation [16]. Therefore, the inverse healing length corresponds to the wave number scale, on which the most strongly localized QB is observed.
We now increase the system size further, and compute the localization length (or respectively its negative inverse -the slope S) of the longest wavelength mode. With increasing N the grid of allowed k-values becomes denser, and at some critical N c (ν) the localization length will reach the finite size of the normal mode space, and the q-breather delocalizes. With a little algebra it follows For even larger (and finally infinitely large) lattices the critical value k (c) 0 is marking a border in k-space: at the given norm density ν, all QBs with seed wave numbers k ≫ k (c) 0 are localized, while we obtain delocalization for k ≤ k (c) 0 . So there is a layer of delocalized QBs at the band edge, whose width grows according to (26) with growing ν. Modes launched inside this layer (with the given norm density) will quickly spread their energy among many other modes -they will quickly relax. Modes launched outside this layer will stay localized in normal mode space, at least for sufficiently long times. k (c) 0 is therefore separating a layer of strongly interacting modes from weakly interacting ones. For large enough effective nonlinearity parameter ν ∼ 1 the whole wave number space is filled with strongly interacting normal modes, and the normal mode picture breaks down completely. As long as ν is smaller, the value of k (c) 0 sets a new length scale 2π/k (c) 0 , which is proportional to the squared healing length ξ 2 h . On that new length scale, the normal mode picture breaks down.

B. Delocalization thresholds and modulational instability
It is instructive to remember, that delocalization thresholds of QBs are related to resonances. Indeed, in the present case, the density of states at the band edge diverges in the limit of large system size, and many modes have almost the same frequencies. It is these small differences, which tend to zero as 1/N 2 , and which are responsible for the resonant mode-mode interaction. Notably the analysis of stability of band edge modes [23] (note: for periodic boundary conditions) shows some interesting correlations. The analyzed band edge modes are compact in normal mode space, yet they undergo a tangent bifurcation at amplitudes, which are the smaller, the larger the system size. These instabilities appear however only for a certain sign of the nonlinearity, which exactly corresponds to the actual observed delocalization of QBs. Therefore we expect, that the (so far not studied) case of a FPU chain with negative quartic nonlinearity, which is known not to yield an instability for the band edge mode, will not show delocalization of QBs close to the (upper) band edge. It is furthermore instructive, that if a tangent bifurcation of a band edge mode takes place, the simulation of a lattice shows the onset of modulational instability, which leads to a collection of energy in smaller system volumes, and finally to the formation of discrete breathers -i.e. to localization in real space. Thus we may expect, that the resonant layer of strongly interacting modes may lead to the formation of localized states in real space, while the rest of the normal mode space is evolving in the regime of localization in normal mode space. Thus, we may expect to observe in an actual simulation localization both in normal mode space and in real space.

C. Comparing analytical and numerical results
While the delocalization at one band edge, as predicted by perturbation theory, is reproduced by numerical data, that does not happen at the second band edge (both band edges change their places, when the nonlinearity inverts sign). The breakdown of perturbation theory follows from comparing the leading and next-to-leading order terms in the expansion. When terms are of the same order, we conclude, that the series will diverge, and therefore the QB solutions will delocalize. This is true for the case when all terms in the series have the same sign. However, when the terms have alternating sign, the above conclusion must not be correct. And indeed the numerical data show, that these sign alternations lead to an effective cancelation, and a final strong localization of a QB solution. A similar breakdown of perturbation theory happens for QB solutions localized at the band center for large nonlinearity (or norm). The perturbation theory tells, that QB states are well localized from both sides of the center. Numerical data however show, that this is true only on one side from the center, while the other side shows a tendency towards delocalization.

D. Open questions
Below we list some potentially interesting and important open questions. First, the influence of different boundary conditions has not been systematically studied. QBs are extended states in real space, therefore their spectrum, and the way they interact, may to some extend be sensitive to the choice of boundary conditions. Second, it remains completely open, what kind of excitation we catch by considering a vortex state in normal mode space.
Finally, it would be interesting to see the changes in the QB properties, when the norm conservation is lifted. In general we expect, that this will lead to the generation of higher harmonics, and new types of resonances. In fact it is exactly resonances due to higher harmonics, which give the leading order contribution to the FPU problem. However, for sufficiently narrow optical bands, higher harmonics will be located outside the optical band. These resonances will be therefore important, when the optical band is wide, or when one considers a complex and broader band structure.