Bifurcations of relative periodic orbits in NLS/GP with a triple-well potential
BBifurcations of relative periodic orbits in NLS/GP with a triple-wellpotential
Roy H. GoodmanDepartment of Mathematical SciencesNew Jersey Institute of TechnologyUniversity HeightsNewark, NJ 07102 ∗ November 8, 2018
Abstract
The nonlinear Schr¨odinger/Gross-Pitaevskii (NLS/GP) equation is considered in the pres-ence of three equally-spaced potentials. The problem is reduced to a finite-dimensional Hamil-tonian system by a Galerkin truncation. Families of oscillatory orbits are sought in the neigh-borhoods of the system’s nine branches of standing wave solutions. Normal forms are computedin the neighborhood of these branches’ various Hamiltonian Hopf and saddle-node bifurcations,showing how the oscillatory orbits change as a parameter is increased. Numerical experimentsshow agreement between normal form theory and numerical solutions to the reduced systemand NLS/GP near the Hamiltonian Hopf bifurcations and some subtle disagreements near thesaddle-node bifurcations due to exponentially small terms in the asymptotics.
This paper considers the nonlinear Schr¨odinger/Gross-Pitaevskii equation (NLS/GP) with a triple-well potential as a model problem for investigating the occurrence and bifurcations of approximaterelative periodic orbits (RPOs) in the setting of a Hamiltonian partial differential equation, i.e.solutions that are periodic modulo a symmetry of the system. The largest part of the paper isdedicated to understanding how these orbits bifurcate when the system undergoes a HamiltonianHopf (HH) bifurcation, but will discuss a number of other families of RPOs in this system thatare unrelated to this bifurcation. This paper will analyze, primarily, a finite-dimensional system ofHamiltonian ODE derived as a Galerkin truncation of NLS/GP with the given potential.The HH bifurcation is a mechanism that gives rise to oscillatory instabilities in conservativesystems including nonlinear dispersive wave equations. It occurs when varying some parameter α in the system across a threshold α HH causes pairs of eigenvalues ± iω and ± iω on the imaginaryaxis to collide and split to form so-called Krein quartets of eigenvalues λ = ± µ ± iω . The motionwith α on one side of α HH is oscillatory due to the purely imaginary eigenvalues, while on the otherside, the real parts lead to exponential growth or decay, accompanied by oscillation due to theimaginary parts. ∗ [email protected] a r X i v : . [ n li n . PS ] D ec ur object of study is the cubic nonlinear Schr¨odinger/Gross-Pitaevskii equation, i∂ t u = − ∂ x u + V ( x ) u − ∣ u ∣ u (1.1)with a particular potential. This equation is ubiquitous in mathematical physics, arising due tothe balance between nonlinearity and dispersion in nearly monochromatic wavepackets. In theoptics context, the equation models the propagation of a nearly monochromatic electrical fieldalong a waveguide where t represents the longitudinal distance and V ( x ) describes the geometryof the waveguide in the transverse direction. The sign of the cubic Kerr nonlinearity impliesthat in regions of higher intensity, the refractive index is increased [1]. This system also arisesin describing the evolution of a Bose-Einstein condensate (BEC), a state of matter that occursat extremely low temperatures at which a gas of identical particles obeying Bose statistics losetheir individual identities and share a common wavefunction. When the BEC is confined by strongmagnetic or optical fields to a quasi-one-dimensional cigar-shaped region, its wavefunction satisfiesNLS/GP with the potential V ( x ) describing the weaker potential along the axis of the cigar [2–4]. Here we focus on a particular potential V ( x ) given as the sum of three identical, symmetric, stronglylocalized potentials V ( x ) , hereafter the “triple-well” potential V ( x ) = V ( x + L ) + V ( x ) + V ( x − L ) . (1.2)We assume that the single potential V ( x ) supports a single discrete eigenfunction for the linearsystem (− ∂ x + V ( x )) U ( x ) = Ω U ( x ) , and that V ( x ) satisfies even symmetry V ( x ) = V (− x ) . An example to keep in mind is V ( x ) = − x, (1.3)for which Ω = − U ( x ) = sech ( x ) . The potential V ( x ) and its three eigenfunctions aredisplayed in Figure 1. One can show that in the large- L limit the eigenvalues take the form ( Ω , Ω , Ω ) = ( Ω − ∆ + (cid:15), Ω , Ω + ∆ + (cid:15) ) , (1.4)where, exponentially as L → ∞ , Ω → Ω , ∆ → , (cid:15) → , and (cid:15) ≪ ∆and ∆ itself approaches zero exponentially as a function of L . The frequency ∆ is by definitionpositive, while (cid:15) can be found by asymptotics to be positive. This arrangement of eigenvaluespredisposes a nonlinear mode associated with the frequency Ω to HH bifurcations, as we shall see. Remark 1.1.
In all the numerical calculations, we take V as in equation (1.3) and L = , whichgives approximate parameters ∆ = . × − and (cid:15) = . × − , so that ∆ ≈ (cid:15) , ensuring alarge separation of scales.
10 -5 0 5 10-2-1.5-1-0.50 V ( x ) -10 0 10 x -0.500.5 Ω = -1.0392 -10 0 10 x -0.500.5 Ω = -1.0027 -10 0 10 x -0.500.5 Ω = -0.9632 Figure 1: The potential V ( x ) of equations (1.2) and (1.3) with L = U which is positive everywhere, the antisymmetric first excited state or dipolemode U , and the second excited state U .System (1.1) is Hamiltonian and can be written i∂ t u = δ H δu ∗ , where H is the conserved energy functional H = ∫ (∣ u x ∣ + V ( x ) ∣ u ∣ − ∣ u ∣ ) dx. The system conserves the squared L norm, also known as the optical power or, in the BEC context,the particle number: N = ∫ ∣ u ∣ dx. (1.5)We are interested in the behavior near nonlinear bound states , or standing waves , solutions toequation (1.1) of the form u ( x, t ) = U ( x ) e − i Ω t , consisting of a real-valued function U ( x ; N ) ∈ H ( R ) and a frequency Ω (N ) satisfyingΩ U = −U ′′ + V ( x )U − U . (1.6)These generally exist as one-parameter families where we think of U and Ω as being functions ofthe optical power N .In absence of the cubic term, equation (1.1), the familiar Schr¨odinger equation of quantummechanics, possesses a set of discrete eigenfunctions and frequencies U j ( x ) and ω j . By the implicit3unction theorem, as applied in [5], for each eigenpair, there exists for small N a unique nonlinearnormal mode (NNM), a nonlinear bound state (U j ( x ; N ) , Ω j (N )) satisfying U j ( x, N ) ∼ √N ( U j ( x ) + O (N )) ; ∫ U j ( x ) dx =
1; Ω j (N ) → ω j as N → . (1.7)Potentials of the form (1.2) have three eigenfunctions, which in the limit of large L approachthe form U ( x ) = ( U ( x + L ) + √ U ( x ) + U ( x − L )) U ( x ) = √ ( U ( x + L ) − U ( x − L )) U ( x ) = (− U ( x + L ) + √ U ( x ) − U ( x − L )) . Numerically computed eigenfunctions are shown in Figure 1. The modes U and U are known asthe ground state and the dipole mode , respectively. The modes U and U are symmetric, while U is antisymmetric. In addition U and U are known, respectively, as the first and second excitedstates. Remark 1.2.
We adopt a naming convention from [6] for the nonlinear standing waves. A modeis labeled by a three-character string composed of , + , and − if, in the limit of large N , the solutionhas, respectively, 0, U (⋅) , or − U (⋅) in the indicated location. Thus the continuations of U and U are labeled, respectively, as U − and U − + − . Less obviously, the continuation of U is labelled U because for large amplitudes, the concentration in the first and third wells decreases to zero, asdiscussed in Section 4.2. Interchanging the symbols “-” and “+” leaves the mode unchanged, whilereversing the order of the labels gives a mode’s mirror image. The study of NLS/GP with multi-well potentials generally falls into two categories. The firstis analysis of existence, bifurcation, and stability of stationary solutions to equation (1.6). Aphenomenon of particular interest here is localization or self-trapping , the concentration of energyin one (or at least a small subset) of the wells as the total amplitude is increased. The second isinvestigation of the time-dependent dynamics. Here, we speak of tunneling , as the solution movesbetween the various wells. This can be periodic or chaotic.The first question is addressed, for small N , for the double -well potential V ( x ) by Kirr et al. [5].This system has two NNMs, a ground state and a dipole mode. As the amplitude parameter N isincreased past a critical value, the ground state becomes unstable and a pair of stable asymmetricstanding waves emerges in a symmetry-breaking (supercritical pitchfork) bifurcation, with U ( x ) concentrated in either the left or right well. The behavior in the limit of large N is shown in [7, 8].Here there exist localized solutions that are centered not at the minima of the potential but atpoints in between, including at the maxima.Stationary solutions of the triple -well problem are enumerated by Kapitula et al [6]. In additionto the three NNM’s, they found six additional solutions, which arise in three saddle-node bifurca-tions as N is increased. They show analytically that the symmetry-breaking seen in the double-wellpotential is not possible with a triple well. They numerically calculated stability of each solution4s N is increased, showing that the NNM loses stability and then regains stability in a pair of HHbifurcations. Localization in either of the outer wells occurs due to a saddle-node bifurcation andlocalization in the central well occurs continuously, without any bifurcations, as N is raised. Ourrecent study [9] estimates the critical values of N for the HH bifurcations. Sacchetti [10] uses semi-classical methods to prove the correctness of a finite-dimensional truncation for the general n -wellsystem, and applies this result to the quadruple well, finding that in this case, too, the ground stateundergoes a pitchfork bifurcation, which may be sub- or supercritical, depending on the exponentof the nonlinearity.The second question, regarding time-dependent dynamics and tunneling, has been addressedin several studies of the double-well problem [11–13]. These show the dynamics of supercriticalHamiltonian pitchfork bifurcations: below the critical amplitude, the dynamics resemble a particleoscillating in a single-well potential, while above this amplitude, the dynamics resemble those of adouble-well potential. For other nonlinearities, the pitchfork may be supercritical or subcritical [13,14]. The finite-dimensional models derived in these papers have a two-dimensional phase space.This allows the dynamics to be understood completely by plotting the level sets of the conservedenergy. These studies confirm that solutions to the reduced model that are sufficiently boundedaway from the separatrix are shadowed over long times by PDE solutions but do not addressthe question of solutions on or near the separatrix, nor the interesting question of how radiativedissipation in the PDE could cause the solution to cross the separatrix.The triple well potential V ( x ) leads to significantly more complicated dynamics. Our previousstudy [9] estimates the critical values of N for the HH bifurcation and begins a study of thefour-dimensional phase space of a finite-dimensional model. In particular, it explores some of thedynamics that arise along with the HH bifurcation, and numerically demonstrates chaotic dynamics.It studies an averaged form of the equations which gives some insight to the geometry underlyingthe observed behavior.These phenomena have also been reported in laboratory studies. Notably, Albiez et al. have re-ported both self-trapping and tunneling in a cigar-shaped BEC in a double-well potential [15].Kapitula et al. perform experiments in an optically-induced optical waveguide [6]. They findsymmetry-breaking in an analogue of the double well problem, and verify certain mathematicaldetails of the triple-well problem.The current study aims to extend previous analyses and demonstrate other structures presentin the dynamics. While it was possible and logical in [12] both to describe the relative periodic orbitsin a reduced ODE model and to prove their validity in the full double-well potential in the samepaper, the dynamics the finite-dimensional reduction of the triple well are sufficiently complicatedto require their own paper, and to leave the proof of their validity for later.
A similar analysis to ours, with overlapping but not identical results, appears in the recentwork of Yang [16], who uses the method of multiple scales to derive a normal form for the HHbifurcation directly from the NLS/GP equation, without assuming a triple-well potential. Mostinterestingly, he finds that in the case of a resonance between the continuous spectrum and thediscrete, some coefficients in the normal form equations become complex-valued, which leads to ablow-up behavior similar to one seen by Kevrekidis et al. [17].5 .3 The discrete self-trapping equations
It should be further noted that mathematically and physically, the system under study has muchin common with the discrete self-trapping (DST) equations i ddt Ψ = − D Ψ − F ( Ψ ) ; Ψ ∈ R n , (1.8)where F ( Ψ ) = (∣ Ψ ∣ Ψ , ∣ Ψ ∣ Ψ , . . . , ∣ Ψ n ∣ Ψ n ) T and D is the discrete second derivative or Lapla-cian. This system, also referred to as discrete NLS, was introduced by Eilbeck et al. to model,among other things, the oscillations of small molecules [18]. Eilbeck and collaborators were thefirst to enumerate the stationary solutions for small values of n and to calculate their stabilitynumerically [18, 19]. Susanto has recently surveyed some important results for this system [20].Exact periodic orbits in the dimer case Ψ ∈ R have been known since the 1980’s [21, 22].The trimer has been studied extensively as a model of BEC in a triple-well potential [23, 24].While many groups have noticed that the system can be reduced from three to two degrees offreedom, most of the papers have focused on stationary solutions and numerical computationsusing tools like Poincar´e sections. Johansson computes some RPOs in this system and identifiestwo HH bifurcations of one stationary solution [25], noting a qualitative difference between thetwo, equivalent to the two different normal forms discussed here. Panayotaros has also studied thetrimer, using topological methods to investigate the stability of a different stationary solution [26].Basarab has done some work to apply the methods discussed in the present paper to the DSTtrimer [27]. In this paper, we study the dynamics of a finite-dimensional ODE system that model the dynamicsof solutions to NLS/GP (1.1) whose initial conditions are a small-amplitude linear combination ofthe U j . By exploiting a phase invariance, this ODE system can be reduced to a two-degree-of-freedom Hamiltonian system which, at small amplitudes, is close to the semisimple -1:1 resonance(to be defined below). The phase space is four-dimensional, so it is not possible to simply draw aphase plane, but by reduction techniques, we are able to determine some important features. Inparticular, whereas [6] enumerated nine families of standing wave solutions, we enumerate additionalfamilies of relative periodic orbits in which energy moves periodically between the wells and describehow these families of orbits change as the optical power N is increased.A primary tool used here is canonical normal forms , originally developed in the study of celestialmechanics. One of our goals is simply to understand what occurs in nonlinear waves that undergoHH bifurcations. To this end, we apply two separate types of normal form calculations to thefinite-dimensional model system. The first of these normal forms, which applies small amplitudes,gives more global information than we expected and describes three families of relative periodicorbits. The second, which is more standard, describes how the topology of relative periodic orbitschanges as two HH bifurcation values are crossed. Three saddle-node (SN) bifurcations are alsostudied, using both normal forms and numerical continuations.The remainder of the paper is organized as follows. Section 2 begins with an overview ofthe notation used and continues with a brief review of the necessary concepts from Hamiltonianmechanics. Section 3 describes the finite-dimensional system of ODEs used to model the behaviorof solutions to NLS and a further reduced system which we will study in greater detail. It contains6symptotic expressions for the nonlinear normal mode solutions in the weakly nonlinear limit andalso describes the reduction from three degrees of freedom to two. In Section 4, we briefly enumerateall the branches of standing wave solutions and their bifurcations for reference in later sections.This necessarily repeats some of the calculations of [6]. Section 5 describes the linearization aboutthe dipole mode U − , which undergoes two separate HH bifurcations, and derives asymptoticexpressions for the value of the optical power N at which the bifurcation occurs. Section 6 describesthe normal form of Chow and Kim [28], which applies in the neighborhood to the semisimple -1:1resonance. We apply this normal form in a neighborhood of U − for optical power which is on theorder of the small parameter (cid:15) defined in equation (1.4). The normal form describes the topologicalchanges that occur in branches of relative periodic orbits in which energy oscillates periodicallyamong the three eigenmodes. In a neighborhood of the HH bifurcation, a more standard normalform is applicable. This comes in two flavors that could be called supercritical and subcritical. InSection 7 we examine the dynamics near the mode U − at its two HH bifurcations and find thatthey are of opposite type. Section 8 constructs a normal form that applies in the neighborhood oftwo saddle-node bifurcations. This bifurcation is shown to have two subtypes, similar to what isseen for the HH bifurcation, as well as three different families of small-amplitude periodic orbits.Section 9 contains numerical studies that show the normal form equations accurately describethe bifurcations of the branches of relative periodic orbits that occur in the finite-dimensionalHamiltonian reduction. In the neighborhood of the first HH bifurcation, we find long-lived nearlyperiodic orbits of the PDE that undergo this same bifurcation. It shows similar agreement near thesecond HH bifurcation. Numerical investigation of the branches in a neighborhood of two saddle-node bifurcations is found to be a bit more subtle than the normal form calculation would indicate,owing to terms beyond all orders in the asymptotics. These numerical results are discussed inrelation to the existing literature. Finally, Section 10 concludes with some discussion summarizingthe findings, putting them in context of other research, and describing possible future directions. • An overbar, ¯ z represents the complex conjugate of z . • Boldface quantities, e.g. p , q , are vectors. • Angle brackets ⟨⋅ , ⋅⟩ represent various types of inner products. • P j ( R n ) is the space of homogeneous polynomials of degree j , i.e. polynomials p ( y ) satisfying p ( c y ) = c j p ( y ) . Monomials in P j are written in multi-index notation y α = y α . . . y α n n where α ∈ Z n + and ∣ α ∣ = α + ⋯ + α n = j . The natural inner product on this space is ⟨ F, G ⟩ P j = F ( ∂ y ) G ( y )∣ y = . (2.1)In particular, when applied to monomials, ⟨ y α , y β ⟩ P j = α ! δ α , β , where the multi-index factorialis α ! = ∏ Nj = α j !. Analogous definitions hold for C n . • I n , or just I , is the n × n identity matrix. J n , or just J , is the 2 n × n matrix ( I n − I n ) .7 The canonical Poisson bracket with respect to ( p , q ) coordinates is { F, G } = n ∑ j = ( ∂F∂q j ∂G∂p j − ∂F∂p j ∂G∂q j ) . • The adjoint operator of a quadratic Hamiltonian H ( p , q ) ∈ P ( R n ) is defined asad H F = { F, H } . • Noting that ad H ∶ P j → P j , define ad ( j ) H to be the restriction of ad H to P j . • The kernel and range of a linear operator L are denoted by ker L and ran L , respectively. Familiarity with the basics of Hamiltonian mechanics will be assumed, but we highlight here themain concepts applied in the paper. Most of the ideas used in this paper can be learned from [29]. An n -degree-of-freedom Hamiltonian system consists of a set of 2 n ordinary differential equations˙ q j = ∂ p j H ( q , p , t ) ; ˙ p j = − ∂ q j H ( q , p , t ) , (2.2)where p , q ∈ R n or y = ( q , p ) lies in a 2 n -dimensional manifold, and H ( q , p , t ) is a C function H ∶ R n × R → R . The variables q and p known respectively as the position and momentum vectors.In y coordinates, ˙ y = J ∇ H ( y ; t ) .If H ( y ) ∈ P , then y solves a linear equation˙ y = J S y (2.3)where S is a symmetric matrix. If λ is an eigenvalue of J S , then − λ , ¯ λ , and − ¯ λ are also eigenvalues,and of the same multiplicity. Thus there are four types of eigenvalues: real pairs ± λ , imaginarypairs ± iω , complex Krein quartets λ = ± µ ± iω , and zeros of even multiplicity. These restrictionsimply that the stability of a fixed point can change in only a few ways as parameters in an equationvary, as shown in Figure 2: in a pitchfork bifurcation, imaginary eigenvalues collide at the originand emerge as a pair of real eigenvalues, transitioning from arrangement (a) to arrangement (b)or from (b) to (d), and in a Hamiltonian Hopf bifurcation two pairs collide on the imaginary axisand become a quartet, transitioning from arrangement (a) to (c). Zero eigenvalues that persistas parameters are varied generally indicate the presence of symmetries. The ODE system studiedin this paper has such a symmetry, which we use to reduce the dimension of the problem. Thereduced system has eigenvalues equal to zero at only a few bifurcation points in parameter space. A fixed point in a symmetry-reduced system corresponds to a type of periodic orbit in the non-reduced system that is known as a relative fixed point (RFP): it is a fixed point modulo the actionof a symmetry group, in our case the group S corresponding to gauge invariance. The standing8 e λ I m λ (a) Re λ I m λ (b) Re λ I m λ (c) Re λ I m λ (d) Figure 2: The four possible non-zero arrangements of eigenvalues in a two-degree-of-freedom Hamil-tonian system.wave (1.6) is an RFP of the NLS system (1.1). A periodic orbit in the symmetry reduced systemcorresponds generically to a quasiperiodic orbit in the original system known also as a relativeperiodic orbit (RPO).While our main concern will be RFPs and RFO’s, in the reduced system we can consider simplerfixed points and periodic orbits. Periodic orbits in a neighborhood of a fixed point y ∗ in a finite-dimensional Hamiltonian system are described by the Lyapunov Center Theorem [29, Ch. 9.2],which states that if the linearization about y ∗ has spectrum {± iω, λ , λ , . . . , λ n } , where iω ≠ λ j / iω is never an integer for j = , . . . n , then there exists a one-parameterfamily of periodic orbits emanating from y ∗ . Further, the period of these orbits approaches 2 π / ω as the periodic orbits approach y ∗ . We refer to these as Lyapunov families of periodic orbits.The Lyapunov Center theorem applies in two situations we will encounter when consideringour two degree-of-freedom reduced system. If the spectrum about y ∗ is of the form {± iω, ± λ } withtwo real and two imaginary eigenvalues, then there is a single Lyapunov family of periodic orbitsconnected to y ∗ . If the spectrum consists of two pairs of imaginary eigenvalues {± iω , ± iω } , with0 < ω < ω , then there is always a Lyapunov family with limiting frequency ω . If ω / ω ∉ Z , thenthere is a second Lyapunov family with limiting frequency ω . Suppose that the Hamiltonian is quadratic and can be written H = ∑ nj = ω j ( p j + q j ) / . Then if k ∈ Z n is a nonzero vector of integers such that ⟨ k , ω ⟩ =
0, the Hamiltonian is resonant . The orderof the resonance ∣ k ∣ = ∑ ∣ k i ∣ . If a system has multiple eigenvalues, the leading-order Hamiltonianmay not be expressible in this form, but the resonance is defined similarly.We consider the semisimple -1:1 resonance in Section 6. This is defined by a Hamiltonianwhose leading-order linear part is of the above form with ω = − ω and which has resonancevector k = ( , ) of order 2. The associated linear system has eigenvalues ± i Ω, each of multiplicitytwo. Since the matrix for this system is diagonalizable, it is referred to as the semisimple form ofthe resonance. A linear Hamiltonian system with the same eigenvalues but whose matrix is notdiagonalizable corresponds to a non-semisimple -1:1 resonance which is the generic case and whicharises in Section 7 below. A zero eigenvalue is always resonant. The Hamiltonian 0 iω resonance,with a pair of imaginary eigenvalues and a double zero, is considered in Section 8 The normal form of a given differential equation is another differential equation, obtained by anear-identity change of variables, that is in the “simplest” form possible in a given region of phase9pace. The normal form transformation may lose information. Consider a Hamiltonian of the form H ( y , (cid:15) ) = H ( y ) + (cid:15)H ( y ) + H ( y ) where H ∈ P , H ∈ P , H ∈ P . Assuming that y = O ( (cid:15) / ) makes the term H = O ( (cid:15) ) and theother terms both O ( (cid:15) ) .Simplifying the system has two parts. First, choosing coordinates in which H is in a standardform, and second, constructing a near-identity change of variables that removes as many terms aspossible from (cid:15)H + H , leaving only the resonant part (cid:15)H N2 + H N4 . This transformation introducesterms of higher order in y and (cid:15) that can themselves be simplified by the same procedure, and soon. The sequence of transformations generally does not converge, but useful information can begleaned from truncating the sequence after a finite number of steps. In particular, if the truncatednormal-form system possesses periodic orbits, then, by the implicit function, so does the originalsystem. In our case, we only need to calculate the results of the first change of variables, and wecan do this without explicitly determining the change of variables.An important step involves determining which terms can be removed and which cannot. A termmay be removed from H j if it is in the range of ad ( j ) H . Therefore the normal form of H j is just theprojection of H j onto some complement of the range. Given the inner product (2.1), the Fredholmalternative decomposes P j as P j = ran ( ad ( j ) H ) ⊕ ker ( ad ( j ) T H ) . In fact there exists an easily constructed quadratic Hamiltonian function H T0 such that ( ad ( j ) H ) T = ad ( j ) H T0 . Thus H N j = proj ad ( j ) H T0 H j = ∑ v ∈ B ⟨ v , H j ⟩ P j ⟨ v , v ⟩ P j v , (2.4)where B is an orthogonal basis for ker ( ad ( j ) H T0 ) .It should be further pointed out that the transpose of ad H is defined in terms of an innerproduct on spaces of polynomials. While there are many ways to define such an inner product, thechoice (2.1) is especially elegant and leads to simple formulas. We consider solutions to equation (1.1) with suitable initial data that we assume can be wellapproximated by a time-dependent linear combination of the eigenfunctions u ( x, t ) = n ∑ j = c j ( t ) U j ( x ) , (3.1)where the U j are the eigenfunctions of the linearization of equation (1.1) and the coefficients c j satisfy a system of ordinary differential equations given below.The coefficients c j ( t ) satisfy a Hamiltonian system of equations i ˙ c j = ∂∂ ¯ c j ˜ H ( c, ¯ c ) H ( c, ¯ c ) = Ω ∣ c ∣ + Ω ∣ c ∣ + Ω ∣ c ∣ − ˜ a ∣ c ∣ − ˜ a ∣ c ∣ ( c ¯ c + ¯ c c )− ˜ a ( c ¯ c + ∣ c ∣ ∣ c ∣ + ¯ c c ) − ˜ a ( c ¯ c + ∣ c ∣ ∣ c ∣ + ¯ c c )− ˜ a ( ∣ c ∣ ( c ¯ c + ¯ c c ) + c ¯ c c + ¯ c c ¯ c ) − ˜ a ∣ c ∣ ( c ¯ c + ¯ c c )− ˜ a ∣ c ∣ − ˜ a ( c ¯ c + ∣ c ∣ ∣ c ∣ + ¯ c c ) − ˜ a ∣ c ∣ . (3.2)where the coefficients are defined by the integrals˜ a jklm = ∫ ∞−∞ U j ( x ) U k ( x ) U l ( x ) U m ( x ) dx, which we note are identically zero if ( j + k + l + m ) is odd. A fuller derivation of this system is givenin [9], along with a careful accounting of the terms ignored in the approximation. Because the theHamiltonian satisfies ˜ H ( e iφ c ) = ˜ H ( c ) , Noether’s theorem implies the system conserves a discrete version of the power N in equation (1.5), N = ∣ c ∣ + ∣ c ∣ + ∣ c ∣ . For potentials of the form (1.2), the coefficients ˜ a jklm approach limiting values a jklm as L → ∞ ,up to exponentially small errors, ( a , a , a , a , a , a , a , a , a ) = ( , , , , − , , , , ) ⋅ A, where A = ∫ ∞−∞ U ( x ) dx. When the potential is given by (1.3), A = . Remark 3.1.
Using the approximate values a jklm instead of ˜ a jklm simplifies the analysis (and evenmore the typesetting and reading) of this paper tremendously. The effect of this simplification isnot completely trivial. The approximation introduces additional symmetry which is reflected in thestructure of its solutions, for example in Figure 6. Under these assumptions on the coefficients a jklm , the Hamiltonian (3.2) becomes H ( c , ¯ c ) = Ω ∣ c ∣ + Ω ∣ c ∣ + Ω ∣ c ∣ − A [ (∣ c ∣ + ∣ c ∣ ) + ( c ¯ c + ¯ c c ) + (∣ c ∣ + ∣ c ∣ ) ( c ¯ c + ¯ c c ) + ∣ c ∣ + ∣ c ∣ ∣ c − c ∣ + ( c − c ) ¯ c + ( ¯ c − ¯ c ) c ] (3.3)The symmetry of the potential V ( x ) and of its eigenfunctions is reflected in the Hamilto-nian (3.3) and the equations of motion. The odd and even subspaces, B odd = span { ˆ e } and B even = span { ˆ e , ˆ e } , where ˆ e j is the unit vector in coordinate j , are both invariant under the flow of H ( c , ¯ c ) .11 .1 Nonlinear Normal Modes The normal modes of this system in the linear limit are simply c ( t ) = e − i Ω j t ˆ e j . These are continued into the nonlinear regime via the linear combinations U − = ⎛⎜⎝ √ N ⎞⎟⎠ e i ( Ω − AN ) t , (3.4a) U = ⎛⎜⎝√ N cos θ + √ N sin θ + ⎞⎟⎠ e i ( Ω + ω + ( N )) t , U − + − = ⎛⎜⎝√ N sin θ − √ N cos θ − ⎞⎟⎠ e i ( Ω + ω − ( N )) t , (3.4b)where θ ± ( N ) = ± AN + O ( N ) and ω ± = − AN + O ( N ) . In addition to the NNM solutions thereare ten other solutions (six when symmetries are taken into account) that arise in saddle-nodebifurcations as discussed in Section 4; see [6]. c ≠ We make the canonical change of variables from ( c , c , c , i ¯ c , i ¯ c , i ¯ c ) to ( z , φ, z , i ¯ z , N, i ¯ z ) givenby c = z e iφ ; c = √ N − ∣ z ∣ − ∣ z ∣ e iφ ; c = z e iφ , (3.5)which is valid when c ≠
0. This can be accomplished in steps by, first, converting to canonicalpolar coordinates c j → √ ρ j e iφ j , then noticing that the Hamiltonian depends on the angles onlyin the combinations φ − φ and φ − φ , using this observation and the conserved quantity N =∣ c ∣ + ∣ c ∣ + ∣ c ∣ to reduce from three degrees of freedom to two, and, finally, converting the systemback to complex coordinates.In the new variables, the Hamiltonian is independent of φ and contains no square roots becauseof the quadratic and quartic dependence of Hamiltonian (3.3) on c and ¯ c : H ( z , ¯z ; (cid:15), N ) = H ( z , ¯z ) + (cid:15)H ( z , ¯z ; (cid:15), N ) + H ( z , ¯z ) ; (3.6)where H = − ∆ ∣ z ∣ + ∆ ∣ z ∣ ; (cid:15)H = (cid:15) (∣ z ∣ + ∣ z ∣ ) + AN ( ( z ¯ z + ¯ z z ) − ( z − z ) − ( ¯ z − ¯ z ) ) ; H = A [ (∣ z ∣ + ∣ z ∣ ) − ( z ¯ z + ¯ z z ) + (∣ z ∣ + ∣ z ∣ ) (( z − z ) + ( ¯ z − ¯ z ) − ( z ¯ z + ¯ z z ))] . (3.7)Fixed points of this system correspond to RFP’s of system (3.3) and periodic orbits in the reducedsystem (3.6) correspond to RPOs in the full system. Remark 3.2.
Since c appears in Hamiltonian (3.3) only in even powers, the square root introducedin the change of variables (3.5) disappears from the Hamiltonian (3.6) . While many of the sourcescited in Sections 1.2 and 1.3 perform some sort of reduction, few take advantage of this feature. Inorder to apply this observation to the DST system (1.8) , the linear part of the equations must firstbe diagonalized, which results in a Hamiltonian more closely resembling equation (3.3) . .3 Reduction on the even subspace B even Reduction (3.5) is not valid for representing the two NNM solutions (3.4b), nor for solutions on B even more generally, because the angle φ is ill-defined. The dynamics on B even can be simplifiedby setting c = c and c into canonical polar coordinates c j = √ J j e iφ j and using the relation (1.4), yielding H = Ω J + Ω J − A ( J + J J + J )− A √ J J ( J + J ) cos ( φ − φ ) − AJ J cos ( φ − φ ) The only dependence on the angles comes in terms of the phase difference φ − φ , so we make thecanonical change of variables φ = ( θ − θ ) , φ = ( θ + θ ) , J = ρ − ρ , J = ρ + ρ . The angle θ is cyclic, so its conjugate ρ is conserved. We note that ρ = N / θ and ρ , yielding a Hamiltonian H = ρ − AN √ N − ρ cos θ − A ( N − ρ ) cos θ − AN + ∆ N + N Ω ; (3.8)and evolution equations ˙ θ = + AN ρ cos θ √ N − ρ + Aρ cos θ ; (3.9a)˙ ρ = − A ( ( N − ρ ) cos θ + N √ N − ρ ) sin θ. (3.9b) The nonlinear bound states become fixed points in the reduced systems corresponding to (3.6)and (3.8). We describe the fixed points of the two reductions separately. (3.6)
Since the leading order terms of H ( z , ¯z ) are quadratic, z = U − of system (3.3). The other fixed points arereal, z j = ¯ z j ≡ ζ j , and satisfy (− AN − ∆ + (cid:15) ) ζ + AN ζ + Aζ − Aζ ζ − Aζ ζ − Aζ = AN ζ + (− AN + ∆ + (cid:15) ) ζ − Aζ − Aζ ζ − Aζ ζ + Aζ = . (4.2)13n order to count its roots, we use Mathematica to eliminate ζ , which yields an odd parity ninth-degree polynomial for ζ alone, b ζ + b ζ + b ζ + b ζ + b ζ = , (4.3)with coefficients depending on the parameters ( A, ∆ , (cid:15) ) and on the conserved quantity N . This canbe factored as ζ times a quartic polynomial in Z = ζ . We again use Mathematica to compute thediscriminant of the Z polynomial. We find numerically the values of N for which the discriminantvanishes, which indicates where bifurcations occur in the solutions of equation (4.3). For thedefault values of ( A, ∆ , (cid:15) ) , bifurcations occur at N a1 ≈ .
246 (creating the RFPs U ++0 and U +00 described in Remark 1.2) and N a2 ≈ .
667 (the RFPs U − ++ and U − +0 ). Figure 3 shows numericallycalculated fixed points ( ζ , ζ ) of equation (4.2), at three values of N confirming the existence oftwo saddle-node bifurcations. Note that (− ζ , − ζ ) is the same relative fixed point as ( ζ , ζ ) . - - - - z z H a L - - - - z z H b L - - z z H c L Figure 3: The red and blue curves are, respectively, the zero level-sets of the two equations insystem (4.2), so intersections solve the system. (a) N = . < N a1 , (b) N a1 < N = . < N a2 , (c) N = . > N a2 . Since all intersections are inside the circles z + z = N (dashed), they correspondto physical solutions. Since all intersections lie near the circle, the PDE solution is dominated bythe even eigenfunctions, at least for small N . The origin in this figure corresponds to the solution U − and the boundary circle to the invariant subspace B even . B even Setting the right side of the ˙ ρ equation (3.9b) to zero implies that either sin θ = θ =
0. Therefore θ = θ = π . Insertingcos θ = ± A ρ + A ∆ ρ + ( ∆ − A N ) ρ − A ∆ N ρ − ∆ N = . (4.4)At small values of N , this system has two fixed points, one with θ = U , and the other with θ = π , corresponding to U − + − . Both are stable to perturbationswithin B even . A saddle-node bifurcation occurs when the discriminant of equation (4.4), consideredas a function of ρ , vanishes, which occurs for exactly one real, positive value of N , N SN , even = √ ( + ⋅ / + ⋅ / ) ∆32 A ≈ . A ≈ . . (4.5)14or N > N SN , even , there are two additional solutions with θ = π . Numerical examples are shown inFigure 4 for N = . < N SN , even , for N = . < N SN , even , and for N = > N SN , even . (a) (b) (c) Figure 4: The phase plane of the even subsystem (3.9) with (a) N = .
1, (b) N = . N = U , while those near the top edge are close to second nonlinear excited state U − + − .For N ≪ ( c , c ) coordinates U = √ N ( − O ( N ) O ( N ) ) and U − + − = √ N ( O ( N ) − O ( N )) . However, for large N , the behaviors of the two NNM solutions are quite different. The ρ coordinateof solution U approaches ρ =
0, meaning that ρ − ρ →
0. Since θ =
0, the two solutions arein phase. Using this in ansatz (3.1) to construct the approximate PDE solution implies that themodes U and U interfere destructively near z = ± L but constructively near z =
0, so that for large N , the solution concentrates in the middle well.By contrast, the solution U − + − stays near the upper boundary of the region ρ ≈ N /
2, so that thethe shape of the PDE standing wave is not greatly altered for large N , so this mode is labeled U − + − ;again, see Remark 1.2. Similarly for the two modes that appear in the saddle-node bifurcation,the one labeled U +++ stays near the top edge, with a profile that resembles U ( x ) , while the onelabeled U +0+ has destructive interference near z = z = ± L .Thus far we have described nine nonlinear bound states: three NNM’s and six more that comeinto existence in three saddle node bifurcations. This is summarized in Figure 5. U The stability of the U solution, the trivial solution (4.1) in the reduced dynamics, can be foundby examining the quadratic part of the reduced Hamiltonian. Separating the ODE i ˙ z j = ∂H∂ ¯ z j | z | / N (+++)(++0)(+0+)(0+0)(+00)(-++)(-+0)(-+-)(-0+) (a) N | z | / N (+++)(++0)(+0+)(0+0)(+00)(-++)(-+0)(-+-)(-0+) Figure 5: (Color online) The nine branches of nonlinear bound states, showing (a) the fraction of N given by ∣ z ∣ and (b) the fraction of N given by ∣ z ∣ . The stability of each branch is denotedby line style. Thick blue: four imaginary eigenvalues. Thin red: two real, two imaginary. Dashedblack: Krein quartet. Dash-dot green: four real eigenvalues.into real and imaginary parts z j = x j + iy j gives linearized evolution equations ddt v = M v where v = ⎛⎜⎜⎜⎝ x x y y ⎞⎟⎟⎟⎠ and M = ⎛⎜⎜⎜⎝ n − ∆ + (cid:15) n n n + ∆ + (cid:15)n + ∆ − (cid:15) − n − n n − ∆ − (cid:15) ⎞⎟⎟⎟⎠ , (5.1)with n = AN . The trivial solution may change stability when it has multiplicity-two eigenvalueson the imaginary axis, in which case the characteristic polynomial P ( λ ; n, (cid:15), ∆ ) = λ + ( n + + (cid:15) ) λ + ( n ∆ − n (cid:15) + ∆ − (cid:15) + (cid:15) − n (cid:15) ) has double roots. Because P contains only even powers of λ , this is equivalent to the conditionthat P (√ q ; n, (cid:15), ∆ ) have double roots, which is that its discriminant vanishes. The discriminant ofthis polynomial is D ( n ; ∆ , (cid:15) ) = n + (cid:15)n + ( (cid:15) − ∆ ) n + ∆ (cid:15) . (5.2)Since ∆ ≫ (cid:15) , when n = O ( (cid:15) ) , the two leading terms are ∆ ( (cid:15) − n ) , and the leading order solutionis n = (cid:15) + O ( (cid:15) ) , i.e. N HH1 = n A = (cid:15) A + O ( (cid:15) ) (5.3)which is equivalent to the condition found in [9]. Using ∆ ≫ (cid:15) , when n = O ( ∆ ) (large) one findsthe largest positive root of the discriminant is N HH2 = ∆ − (cid:15) A + O ( (cid:15) ) , so that this solution is unstable if N HH1 ⪅ N ⪅ N HH2 , as is shown numerically in [6] and [9]. Directlysolving equation (5.2) using Remark 1.1, we find N HH1 ≈ . N HH2 ≈ . , N HH1 ≈ . N HH2 ≈ . , in reasonable agreement with the values above. This branch is the bottom curve of both branchesof Figure 5. HH bifurcations are visible at the two computed values. (3.6) To find the stability of the other fixed points z ∗ = ( ζ , ζ ) , we linearize, letting z j = ζ j + ( x j + iy j ) and ¯ z j = ζ j + ( x j − iy j ) . The linearized equations are ddt v = ( M + AM ( z ∗ )) v where M and v aregiven in equation (5.1) and M = ⎛⎜⎜⎜⎝ ζ − ζ ζ + ζ − ζ − ζ ζ − ζ − ζ − ζ ζ − ζ ζ − ζ ζ + ζ − ζ + ζ ζ + ζ ζ + ζ ζ + ζ ζ + ζ ζ + ζ ζ + ζ ζ − ζ ⎞⎟⎟⎟⎠ . The stability is calculated numerically and shown in Figure 5. The branches U ++0 and U − ++ alwayshave two real and two purely imaginary eigenvalues. When created in saddle node bifurcations, thebranches U − +0 and U +00 both have four imaginary eigenvalues and are stable. While branch U +00 isstable for all N , the branch U − +0 loses and regains stability in a pair of HH bifurcations. B even The reduced system (3.6) is ill-defined on B even , and system (3.8) describes only the motion within B even and cannot describe the stability of perturbations in directions complementary to B even . Todetermine the stability of the solutions on B even requires a reduction similar to equation (3.5),except one that eliminates z or z instead of z . The matrix describing the linearization about thesolutions on B even is reducible to two smaller systems: one describing the motion of perturbationswithin B even and one for perturbations orthogonal to B even . It is clear from Figure 4 that theeigenvalues corresponding to in-plane motion for modes U − + − , U , and U +0+ are imaginary whilethose for U +++ are real. The remaining two eigenvalues, corresponding to perturbations out of theplane B even , come from the linearized equations of motion for z = r + is in a neighborhood of thefixed points, ( ˙ r ˙ s ) = ⎛⎜⎝ − ∆ − (cid:15) + A ( N + ρ ) + A ( N − ρ )√ ρ cos θ √ N − ρ ∆ + (cid:15) + A ( N − ρ ) − A ( N − ρ )√ ρ cos θ √ N − ρ ⎞⎟⎠ ( rs ) . If the product of the off-diagonal terms is positive, then the associated fixed point has real eigen-values, and if it is negative, these eigenvalues are imaginary. This is found to be negative for U and U − + − and positive for U +0+ and U +++ . See [6] for a more thorough discussion.17 Normal form for N = O ( (cid:15) ) At the HH bifurcation for N = N HH1 , the number of Lyapunov families of periodic orbits attachedto the RFP U − drops from two to zero, and we would like to know their fate. Then, at the secondHH bifurcation at N HH2 = O ( ∆ ) , the number of Lyapunov families jumps back up to two. Theinterval between these two bifurcations is marked by the dashed black curve segment in Figure 5.Later we will address the question of how and whether the branches that exist for N > N HH2 relateto those that exist for N < N HH1 and will discover via numerical experiments that they are entirelyseparate.We turn first to the HH bifurcation at N = N HH1 . We apply a result due to Chow and Kim [28]which allows us to find the periodic orbits of Hamiltonian, and which requires knowledge only ofthe truncated normal form of equation (3.6) H N , trunc ( z , ¯z ; (cid:15) ) = H ( z , ¯z ) + H N2 ( z , ¯z ; (cid:15), N ) + H N4 ( z , ¯z ) . (6.1)We proceed in three steps. First, we state the the theorem to be applied. Second, we calculate theterms in the truncated normal form (6.1). Third, we we apply the theorem to the normal form,thereby enumerating the periodic orbits. This leads to further insights into the bifurcation thatoccurs at N = N HH1 from equation (5.3), which we explore via further normal-form analysis.Chow and Kim point out that, according the Moser-Weinstein reduction theorem [30], thereexists for sufficiently small (cid:15) , a C ∞ normal form E ( z , ¯z ; (cid:15) ) that depends smoothly on (cid:15) for 0 < ∣ z ∣ ≪ H N ( z , ¯z ; (cid:15) ) up to any finite order in (cid:15) . Hence, they conclude, it is sufficient tolook for periodic orbits in the truncated normal form system H N , trunc . Theorem 1.
Consider a Hamiltonian of the form (3.6) whose leading order term H is in semisim-ple -1:1 resonance and which has truncated normal form H N , trunc in equation (6.1) . Let z =( z , z , ¯ z , ¯ z ) and define the matrices L and B by H ( z , ¯z ) = ⟨ z , L z ⟩ and H N2 ( z , ¯z ) = ⟨ z , B z ⟩ Let z ∗ be a critical point of H N , trunc , subject to the constraint H = h , i.e. assume that H ( z ∗ ) = h and ∇ ( H N , trunc − ηH )∣ z = z ∗ = , (6.2) then system (3.6) has a π / ∆ η periodic orbit on the level set H = h z ( t ) = e ηJLt z ∗ , where the Lagrange multiplier η is given by η = + ⟨ L z , (cid:15)B z + ∇ H N4 ⟩∣ z ∣ RRRRRRRRRRR z = z ∗ (6.3) for sufficiently small values of (cid:15) , and ∣ z ∗ ∣ . In [9], we carry out this reduction using von Zeipel averaging, and the periodic orbits are givenas fixed points of the corresponding averaged equations.18 .1 Calculating the truncated normal form
The leading order Hamiltonian H given by equation (3.7) is in a suitable form to apply Theorem 1,so we can skip the first step described in Section 2.2.4 and move right on to construct the adjointoperator ad H . In the ( z , ¯z ) coordinates, the Poisson bracket is { F, G } = i ∑ j ∈{ , } ( F ¯ z j G z j − F z j G ¯ z j ) . so the adjoint operator of H isad H G = i ∆ ( G z − G z − G ¯ z + G ¯ z ) . When applied to a monomial z α ¯z β ∈ P j ,ad H z α ¯z β = i ∆ ( α − α − β + β ) ⋅ z α ¯z β . The operator ad ( j ) H acts diagonally on the basis of monomials, and so is symmetric andker ad ( j ) H T = ker ad ( j ) H = { z α ¯z β ∈ P j ∣ α − α − β + β = } . Thus the leading order normal form is simply the projection of equation (3.6) onto the span of themonomials whose exponents are positive integer solutions of ( − − ) ⎛⎜⎜⎜⎝ α α β β ⎞⎟⎟⎟⎠ = ( j ) , with general solution α + β = j α + β = j . This has positive integer-valued solutions only for j ∈ Z + . The resonant monomials can be enu-merated by specifying α and α to be integers drawn from the set { , . . . , j / } . The resonantquadratic and quartic monomials are tabulated in Table 1. α / α z ¯ z ∣ z ∣ ∣ z ∣ z z (a) Degree Two α / α z ¯ z ∣ z ∣ ¯ z ¯ z ∣ z ∣ ∣ z ∣ ¯ z ¯ z ∣ z ∣ ∣ z ∣ ∣ z ∣ z z ∣ z ∣ ∣ z ∣ z z z z (b) Degree Four Table 1: Resonant monomials of degree two and fourThus, we can read off the resonant terms that make up the normal form of the Hamiltonian (3.6) H N2 ( z , ¯z ; (cid:15), N ) = (cid:15) (∣ z ∣ + ∣ z ∣ ) + AN ( z z + ¯ z ¯ z ) ; (6.4a) H N4 ( z , ¯z ) = A [ ∣ z ∣ − ∣ z ∣ ∣ z ∣ + ∣ z ∣ − (∣ z ∣ + ∣ z ∣ ) ( z z + ¯ z ¯ z )] . (6.4b)19 .2 Finding the periodic orbits We find it convenient to eliminate the Lagrange multiplier from system (6.2), and to present theresults in canonical polar coordinates z j = √ J j e iθ j .The constrained critical points satisfy the equations √ J J ( (cid:15) − A ( J + J )) + A ( N ( J + J ) − J − J J − J ) cos Θ =
0; (6.5a) √ J J ( N − J − J ) sin Θ = , (6.5b)where Θ = θ + θ . The physically meaningful solutions lie in a solid torus, with triangular cross-section defined by the simplex 0 < J ; 0 < J ; J + J < N ; Θ ∈ S . (6.6)On the boundaries of the triangular cross-section one or both of the angles θ , θ is undefined, asis Θ. Therefore some solutions of system (6.5) do not correspond to periodic orbits of system (3.6).Since equation (6.5b) is factored, it is convenient to solve this first, and to use its factors toclassify the periodic orbits. The equations are equivariant to the transformation that interchanges J and J , so the diagram containing the solution curves of this system is symmetric across theplane J = J . The factor ( N − J − J ) lies on the boundary of the open domain of definition (6.6),and solutions in which this term is zero do not correspond to periodic orbits of the system (3.2). We find that setting J = J = (a) ( J , J ) = ( , ) , (b) ( J , J ) = ( N, ) , and (c) ( J , J ) = ( , N ) (6.7)for any value of Θ. We note that Θ is ill-defined for these solutions, but that these solutionsare exactly the NNMs. Solution (a) corresponds to U − using (4.1), while solutions (b) and (c)correspond to the NNM solutions U and U − + − of equation (3.4b), with which they agree to theorder of the truncation. J + J < N Two continuous branches of solutions are defined by setting sin Θ = = nπ , making cos Θ = ± S with Θ = π and S π on which Θ = π mod 2 π . It is useful to look at the intersections of these twobranches with the symmetry axis J = J = J , 0 < J < N . On this axis, equation (6.5a) becomes J ( (cid:15) − AJ + A ( N − J ) cos Θ ) = S , this has solution J = (cid:15) + AN A which is in the interval ( , N ) when N > N even = (cid:15) A . J (a) J J (b) J J (c) Figure 6: The locations in ( J , J , Θ ) space of periodic orbits of system H N , trunc . Solutions areshown for N j ∈ { . , . , . } , chosen such that N < N even < N < N HH1 < N . The trivialsolutions from Section 6.2.1 are at the corners and the solutions S π in blue and on S in dashedred from Section 6.2.2.On the family S π , equation (6.8) has solution J = J = (cid:15) − AN A which is positive when N > (cid:15) / A , which we recognize as the right-hand side of N HH1 given by theleading-order term of condition (5.3) for the instability of the trivial solution. The birth of thissolution would appear to be a consequence of the Hamiltonian Hopf bifurcation which is discussedfurther in Section 7.The results of these calculations are shown in Figure 6. The triangle of physical solutions (6.6)is shaded gray, with the three NNMs (6.7) at the corners: U − at the origin, U to its right, and U − + − above it. The branch S π is shown as a solid line, and the branch S dashed. All the RPOson the branch S are unphysical for N < N even , and lie outside the triangle in Figure 6(a). For N > N even , as in images (b) and (c), some of the solutions on S are physical, having crossed intothe triangle. For N < N HH1 , as in images (a) and (b), the branch S π consists of two pieces, the firstinterpolating between U − and U , and the second between U − and U − + − . At N = N HH1 , thetwo pieces merge with each other and detach from U − +0 , interpolating directly between U and U − + − . U − The goal in this section is to understand the two HH bifurcations of U − in more detail by focusingon a neighborhood of the bifurcation points. After describing the normal form for a Hamiltoniannear a non-semisimple -1:1 resonance, we will further normalize the normal-form Hamiltoniandefined by equations (6.1) and (6.4) where it becomes unstable at N = (cid:15) / A ≈ N HH1 . In Section 9,we numerically check the bifurcation type at N HH2 in terms of this analysis.21 .1 The general normal form
It is well known [29] that a two-degree-of-freedom quadratic Hamiltonian system with a deficientpair of imaginary eigenvalues ± i Ω can be put in the form H ( ξ ) = Ω ( ξ η − ξ η ) + σ ( ξ + ξ ) , (7.1)with σ = ±
1. This leads to a linear evolution equation ddt ξ = B ξ where B = Ω S − σ J , S = ⎛⎜⎜⎜⎝ − − ⎞⎟⎟⎟⎠ , and J = ⎛⎜⎜⎜⎝ ⎞⎟⎟⎟⎠ . (7.2)Constructing the normal form for the perturbed system H ( ξ ) = H ( ξ ) + δH ( ξ ) + H ( ξ ) (7.3)requires projecting H ( ξ ) and H ( ξ ) onto the null spaces, respectively, of the operators ad ( ) H T0 andad ( ) H T0 . Unlike in Section 6, however, the operator ad H does not act diagonally on the monomialbasis, so determining its adjoint null-space, and constructing the projection operator (2.4) requiresexplicit construction of the operators ad ( ) H and ad ( ) H in terms of the monomial bases of P ( ) ( R ) and P ( ) ( R ) . These are 10 ×
10 and 35 ×
35 matrices, so computer algebra is very helpful.Defining the quadratic combinationsΓ = ξ η − ξ η , Γ = ξ + ξ , and Γ = η + η H ( ξ ) = ΩΓ + σ Γ and resonant terms of the formker ( ad ( ) H T0 ) = span { Γ , Γ } and ker ( ad ( ) H T0 ) = span { Γ , Γ Γ , Γ } . Projecting terms from the Hamiltonian (7.3) onto these bases yields the normal form: H N = H ( Γ , Γ ) + δH N2 ( Γ , Γ ) + H N4 ( Γ , Γ ) (7.4)where H N2 = κ Γ + κ Γ , and H N4 = κ , Γ + κ , Γ Γ + κ , Γ . After making a scaling (which is symplectic with multiplier µ ) with µ ≪ ξ j → µ ξ j ; η j → µη j ; δ → µ s (with s = ± H → H / µ , the Hamiltonian becomes H = ΩΓ + µH pert + O ( µ ) where H pert = σ Γ + s κ Γ + κ , Γ . u ( t ) = u ( τ ) + µ u ( τ ) + O ( µ ) where τ = ( ω + µω + O ( µ )) t and u j ( τ + π ) = u j ( τ ) . the sequence of equations is O ( ) ∶ ω d ξ dτ − Ω S ξ = O ( µ ) ∶ ω d ξ dτ − Ω S ξ = − ω d ξ dτ + J ∇ H pert . (7.6)The sequence continues for higher orders in µ , but the essential information, i.e. the topology ofthe periodic orbits is determined by the expansion to this order.The solution to equation (7.5) is ξ = ⎛⎜⎜⎜⎝ cos τ sin τ − sin τ cos τ τ sin τ − sin τ cos τ ⎞⎟⎟⎟⎠ ⎛⎜⎜⎜⎝ Ξ Ξ Y Y ⎞⎟⎟⎟⎠ and ω = Ω . The conditions to avoid secular terms in equation (7.6) are Ξ = − ω Y , Ξ = ω Y , and σω + κ , r = s κ . (7.7)The bifurcation can thus take two forms depending on the sign of σκ , . Assume that the originis stable for δ < δ >
0. Define the hyperbolic HH bifurcation to occur when σκ , <
0. In this case, the two Lyapunov families of periodic orbits lie on separate branches of ahyperbola for δ <
0; for δ = ( r, ω ) = ( , ) ; and for δ > r =
0. See Figure 7(a). If σκ , >
0, then the two Lyapunov families lie on a single, continuous,approximately elliptical branch for δ <
0. This shrinks to a point at δ = δ > elliptical HH bifurcation . r ω δ >0 δ <0 δ =0 (a) r ω (b) Figure 7: The branches of time-harmonic solutions in a neighborhood of the origin. (a) HyperbolicHH (b) Elliptical HH. 23 .2 The HH bifurcation at N HH1
In this section, we apply the above analysis to the normal form Hamiltonian (6.1) with N = (cid:15) / A + δ with δ ≪ (cid:15) using real coordinates x j + iy j = z j . At δ =
0, the linearized evolution equations under Hamilto-nian H N , trunc are ddt ⎛⎜⎜⎜⎝ x x y y ⎞⎟⎟⎟⎠ = ⎛⎜⎜⎜⎝ (cid:15) − ∆ − (cid:15) − (cid:15) ∆ + (cid:15) ∆ − (cid:15) − (cid:15) − (cid:15) − ∆ − (cid:15) ⎞⎟⎟⎟⎠ ⎛⎜⎜⎜⎝ x x y y ⎞⎟⎟⎟⎠ ≡ A x This matrix has multiplicity-two eigenvalues ± i ∆ and is non-diagonalizable. To put it in theform (7.1), an algorithm of Burgoyne and Cushman [31] is used to construct a symplectic matrix P —one satisfying P T J P = J —such that in the new coordinates ξ = ⎛⎜⎜⎜⎝ ξ ξ η η ⎞⎟⎟⎟⎠ = P − x , the leading order linear part is given by B = P − A P in equation (7.2). The sign σ is determinedby the algorithm. The Burgoyne algorithm gives P = ⎛⎜⎜⎜⎜⎜⎝ √ (cid:15) √ (cid:15) √ (cid:15) −√ (cid:15) − √ (cid:15) √ (cid:15) √ (cid:15) √ (cid:15) ⎞⎟⎟⎟⎟⎟⎠ and the normal form algorithm gives σ = , κ = , κ = − Aδ, κ , = A, κ , = , and κ , = A(cid:15) Since σκ , = A(cid:15) >
0, the HH bifurcation is hyperbolic.The solid curve in Figure 8 is a numerical analog to Figure 7 and contains three plots generatedfrom the same data as Figure 6. On the x -axis is r = ρ + ρ and on the y -axis we plot ( η − ) fromequation (6.3), which in these coordinates is − ∆ ( ρ − ρ ) (− A √ ρ ρ cos Θ + A ( ρ + ρ ) + (cid:15) ) ( ρ + ρ ) . This figure clearly shows the pinch-off associated with the hyperbolic HH bifurcation.We have analyzed the bifurcation at N HH1 by applying a second normal form reduction to thetruncated normal form Hamiltonian derived in the previous section. We do not have this availableat the other HH bifurcation that occurs at N = N HH2 . In Section 9.2 we determine the type of thisbifurcation numerically. 24 +J ω × -5 -4-3-2-101234 (a) J +J ω × -5 -6-4-20246 (b) J +J ω × -5 -8-6-4-202468 (c) Figure 8: Graphs of the solution amplitude vs. the frequency correction, with the family S solidand the family S π dashed, reproducing the bifurcation picture in Figure 7. The saddle-node bifurcations resulting in the pairs of RFPs U ++0 / U +00 and U − +0 / U − ++ , depicted inFigures 3 and 5, both show the creation of a pair of fixed points, one with two pairs of imaginaryeigenvalues, and one with one pair of imaginary and one pair of real eigenvalues. This bifurcationis also known at the Hamiltonian 0 iω bifurcation [32, 33]. In this section we consider a normalform for such behavior. The system undergoes the bifurcation as the parameter µ crosses zero: H = ( q + p ) + α ( p + µq − q ) + βIq + H higher ( q , I ) + R ∞ ( q , q , p , p ) (8.1)where I = ( q + p )/ α >
0, but β may have either sign. This normal form can be formally deduced in thesame manner that we determined the resonant terms in Sections 6.1 and 7, where H = ( q + p + p ) and the remaining terms lie in ker ( ad H T0 ) .The term H higher represents a formal power series expansion and contains only resonance terms.The term R ∞ is a beyond-all-orders small term. Because H higher is independent of the angle θ inthe ( q , p ) plane the action I is conserved in any truncation of the system to finite degree. Thus,the act of computing a normal form for this bifurcation introduces a new conserved quantity, thusrendering the normal-form equations integrable when the underlying dynamics may not be. Thiswill be important in the numerical calculations of Section 9.4.3.In this preliminary analysis, we consider the truncated system where the terms H higher ( q , I ) and R ∞ ( q , q , p , p ) are ignored. The equations of motion are˙ q = ( + βq ) p ; ˙ q = αp ; (8.2a)˙ p = −( + βq ) q ; ˙ p = − α ( µ − q ) − βI. (8.2b)When µ > q µ ± = ( q , q , p , p ) ≈ ( , ±√ µ, , ) , (8.3)with q µ − stable and q µ + unstable. Each has a pair of imaginary eigenvalues λ = ± iω ≈ ± i comingfrom the ( q , p ) directions. The stable fixed point has an additional pair of imaginary eigenvalues25 = ± iω = ± O (√ µ ) i . Therefore both fixed points have a family of short -periodic orbits comingfrom ω and the stable fixed point q µ − has an additional family of long -periodic orbits coming from ω when ω − ∉ Z . A third family of periodic orbits that is not Lyapunov plays an important rolein this bifurcation. We describe the global character of the three families.The subspace defined by ( q , p ) = ( , ) is invariant. Off of this subspace, we introduce canonicalaction-angle coordinates in the ( q , p ) degree of freedom q = √ I cos θ, p = √ I cos θ. We find that I and θ satisfy ˙ I = θ = −( + βq ) , which has solution I = ν Q θ ( t ) − θ ( ) = − ∫ t ( + βq ( t ′ )) dt. (8.4)Let µ = ν s with ν > s = ± , q = νx, and γ = βQ α + s. Then x ( t ) solves ¨ x + να ( γ − x ) = . (8.5)When γ > s = µ >
0) and for any b satisfying b < γ this system has a solution x ( t ) = c + ( b − c ) sn ( ωt, k ) , (8.6)where a = − b + √ γ − b c = − b − √ γ − b k = √ b − ca − c ; ω = √ ( a − c ) α ν . (8.7)The above calculation shows that equations (8.4) and (8.6) are the general solution to sys-tem (8.2). Among these orbits are three families of periodic orbits, which we describe next. These orbits correspond to fixed points of the ODE (8.5), i.e. x ( t ) = X = γ or X − β α Q = s. (8.8)This is reminiscent of equation (7.7), showing that this bifurcation comes in two types. If β < elliptical iω bifurcation . There are no periodic orbits when s = − s =
1, which is depicted as the heavy solid curve inFigure 9(a). Note that s = µ >
0, the condition for the existence of fixed points, and thatthese are just degenerate periodic orbits along this curve. The curve of periodic orbits shrinks toa point as ν →
0. 26he condition β > hyperbolic iω bifurcation . Regardless of whether s = ± s = − s =
1. Only in the case s = µ =
0, the two branches merge into a pair of crossed lines, and re-connect. These areshown as the heavy solid curves in Figure 9(b) and (c). -0.02 0 0.02 q -0.08-0.06-0.04-0.0200.020.04 q (c) q p (d) Figure 9: System (8.1) has three families of periodic orbits—short, long, and mixed, representedby heavy solid, dashed, and light red curves, respectively—and two types—elliptic and hyperbolic.In the elliptic case, the branches only exist for µ >
0, shown in (a) for ( α, β, µ ) = ( , − , − ) .In the hyperbolic case, the short and mixed branches exist regardless of the sign of µ , with theshort-periodic branches undergoing a reconnection bifurcation at µ = µ n . The long periodic branch only exists for µ > (b) and (c) for ( α, β, µ ) = ( , , − − ) and ( α, β, µ ) = ( , , − ) , respectively. (d) The long-period dynamics due to the saddle-node bifurcation.
The Lyapunov family of long-period orbits exists only when µ >
0. It lies in the invariant plane ( q , p ) = ( , ) and consists of the bounded solutions to equation (8.5) with Q =
0; see Figure 9(d).Representing each member of a family of periodic orbits by the values of ( q , q ) when p = ( q , q ) = ( , √ µ ) to the furthest point on thehomoclinic loop at ( , − √ µ ) . Each periodic orbit is represented by two points on this line onopposite sides of the elliptic fixed point at ( q , q ) = ( , −√ µ ) ; these are the marked points anddashed lines in Figure 9(a) and (c). The period increases to infinity toward the endpoints of thesegment. The two families of Lyapunov periodic orbits have trivial dynamics in one or the other degree offreedom. In the third family, both degrees of freedom evolve in time.The component x ( t ) given by equation (8.6) has period T = K ( k ) ω where K ( k ) is a complete elliptic integral of the first kind. It oscillates between x ( ) = c and x ( T / ) = b . 27he orbit is periodic if ( q , p ) oscillates an integer number of times during one period of x ( t ) i.e. if ∆ θ is integer multiple of 2 π over one period of x ( t ) . Integrating equation (8.4), we find∆ θ = − T − νβ ∫ T x ( t ′ ) dt ′ = − √ ( a − c ) α ν ( K ( k ) + νβ [ aK ( k ) − ( a − c ) E ( k )]) . Thus the condition that an orbit be periodic is that √ ( a − c ) α ν ( K ( k ) + νβ [ aK ( k ) − ( a − c ) E ( k )]) = πn (8.9)for some n ∈ Z . This condition depends on the initial condition when ˙ q =
0, i.e. on b or c , and onthe value Q . These orbits are shown in Figure 9 for α = β = ±
2, and µ = ± − .The topology of these branches mirrors that of the short-periodic orbits. The upper-limitingenvelope satisfies equation (8.8) with X = b = √ γ in equation (8.6). The lower envelope is relatedto b by the equation for c in equation (8.7) and is given by c = − √ γ . Each of the branches ofmixed periodic orbits is a π -times-integer-valued level-set curve of the form (8.9).The most interesting question to ask is how the branches change as µ changes. The most obviousthing to note is that for β < µ < µ >
0. The behavior for they hyperbolic case β > µ < µ > µ where the linearization about the fixedpoint is resonant, i.e. where the linearization of system (8.2) about the fixed point q µ − has a longperiod that is n times its short period, i.e. when1 α √ ν = n − βν . This occurs at ν = √ µ = α n + β − √ α n + n α ββ = α n − β α n + O ( n − ) . To analyze the bifurcation, we perform a perturbation expansion. We let ν = α n − β α n + C α n , with C to be determined. We further assume that Q = ν / ˜ Q and assume the solution (8.6) issmall amplitude with b = −√ γ + O ( ν ) , in particular choosing b = −√ γ + ν ˜ B . Inserting these assumptions into the resonance relation (8.9), we find that to eliminate terms of O ( n − ) requires 4 ˜ B − βα ˜ Q = C. This quadratic form equation mimics equation (8.8) for the short-periodic orbits. The case β < C >
0, where as the case β > C =
0. This confirms the elliptical and hyperbolicconfigurations of the branches in Figure 9. 28 emark 8.1.
For the orbit to be periodic, it is sufficient that the ratio of the two frequencies berational. However at the resonant bifurcation the Lyapunov center theorem fails to hold, and thishas important consequences for the periodic orbits, as shown by the numerics below.
The effect of the term H higher ( q , I ) in the Hamiltonian (8.1) is minor, at least for small initialconditions. The dynamics in ( q , p ) remains a circular rotation with an I -dependent frequency, andthe ( q , p ) dynamics remain a Newtonian evolution with a potential that depends parametricallyon I . The effect of the neglected term R ∞ ( q , q , p , p ) is more significant as it destroys thesequalitative properties and introduces nontrivial interactions between the two degrees of freedom.In particular, the structure of the normal form ignoring R ∞ is such that the phase difference betweenthe two oscillators is arbitrary. The simplest effect of R ∞ is to select a phase difference. Its othereffects will be described in detail in Section 9.4 where periodic orbits are calculated numerically. Our analysis has proceeded as a sequence of simplifications: from PDE to ODE and then to normal-form ODE. Our numerical investigations will proceed in the opposite direction. We have alreadyvisualized the bifurcation diagram determined by our analysis of the normal form. We comparethese with branches of numerically-computed periodic orbits of the “full” two degree-of-freedomODE system (3.6). N HH1
Numerical continuation of periodic orbits of system (3.6)In this section we compare numerically computed periodic orbits of the Hamiltonian system (3.6)with the approximate periodic orbits obtained from the normal form equations, satisfying sys-tem (6.5). The ODE orbits are computed using an algorithm by Viswanath that combines aspectsof the Poincar´e-Lindstedt algorithm and Newton iteration [34]. Because the periodic orbits lieon one-parameter families, we use pseudo-arclength continuation to generate them. Due to thesystem’s symmetries, all periodic orbits can be chosen such that z ( ) and z ( ) are real.Figure 10 shows the period of the numerical solution as a function of the initial conditionmagnitude, analogous to Figure 8. These solutions confirm the predictions based on the normalizedsystem (6.1), but lack the symmetry present in the solutions of the truncated normal form equations. Numerical calculation of near-periodic solutions to NLS/GP
Exact time-quasiperiodic orbits of system (1.1) on the whole real line have been shown not toexist [35]. Numerically, Marzuola find such solutions with a two-well potential slowly lose energyto radiation [12]. Nonetheless the decay of energy can be very small, and we are able find solutionsthat return very close to their initial conditions and match closely the periodic solutions to thereduced system above.We take initial conditions of the form (3.1), with c j ( ) ≪
1, run the PDE solver and then projectthe numerical solutions onto the eigenfunctions, giving approximate projections c j ( t ) . Cancelling29 ( x (0) + x (0) ) / P e r i od (a) ( x (0) + x (0) ) / P e r i od (b) AB C ( x (0) + x (0) ) / P e r i od (c) Figure 10: The period of solutions plotted as a function of the maximum magnitude of the nu-merically computed periodic orbit of system (3.3) over one period, with the same parameter valuesas previous figures. Lines represent solutions of the Hamiltonian system (3.6), while black circlesdenote near-periodic orbits of NLS/GP.out the phase of c ( t ) gives slowly evolving envelopes z ( t ) and z ( t ) . Our procedure to find near-periodic orbits works entirely with these time series and not with the full solution, and thus relieson the fact that the numerical solution is well-described by the solution ansatz.We note that the periodic orbits to the reduced ODE system are all symmetric across both thereal and imaginary z j -axes, and that the solution components z and z are always either in-phaseor 180 ○ out of phase. Thus at instants t when z ( t ) ∈ R , so is z ( t ) , and, similarly, both arepurely imaginary at the same instants. We use this observation to construct a shooting procedure.Starting with real-valued initial conditions, integrate the PDE until a time τ such that z ( τ ) ispurely imaginary. Then use a root finding algorithm to adjust the initial conditions such that z ( τ ) is purely imaginary to some small tolerance. The computed solution is then one quarterof a periodic orbit. Depending on the shape of the periodic orbits, it is sometimes numericallypreferable to start with purely imaginary initial data and integrate until the solutions are real, orelse to run the simulation until z is imaginary and adjust until z is also imaginary. Some periodicorbits discussed in later sections have fewer symmetries, so we calculate half-periods instead ofquarter-periods.The program to implement this scheme numerically solves the NLS/GP system using a codewritten by T. Dohnal. It uses operator splitting to separate the stiff linear terms from the non-stiffnonlinear term, and perfectly matched layers to absorb outgoing radiation at the boundaries [36, 37].While the PDE system of course loses energy to radiation escaping the computational domain,we find at these amplitudes that the numerical solution ( z , z ) returns to within 10 − of the initialcondition. Using this initial condition for a long numerical integration, we find that the solutionremains close to periodic for times of the order of thousands. The black circles in Figure 10 indicateapproximate periodic orbits obtained in this way, showing good qualitative agreement.Figure 11 shows numerical solutions of NLS/GP corresponding to the solutions marked A-C in Figure 10(b). These three solutions are chosen because they have near complete cancellationleading to dark spots in the PDE solution field. In cases A and B , these dark spots occur in thewaveguides on the edge whereas in C , the dark spot occurs in the middle waveguide.30
10 0 10 x -0.0500.05 t=0 (a) -10 0 10 x -0.0500.05 t=T/4 -10 0 10 x -0.0500.05 t=T/2 -10 0 10 x -0.0500.05 t=3T/4 -0.1 0 0.1 Real(z j ) -0.1-0.0500.050.1 I m ag ( z j ) (a) -10 0 10 x -0.0500.05 t=0 (b) -10 0 10 x -0.0500.05 t=T/4 -10 0 10 x -0.0500.05 t=T/2 -10 0 10 x -0.0500.05 t=3T/4 -0.1 0 0.1 Real(z j ) -0.0500.05 I m ag ( z j ) (b) -10 0 10 x -0.0500.05 t=0 (c) -10 0 10 x -0.0500.05 t=T/4 -10 0 10 x -0.0500.05 t=T/2 -10 0 10 x -0.0500.05 t=3T/4 -0.05 0 0.05 Real(z j ) -0.0500.05 I m ag ( z j ) (c) Figure 11: Left: The illumination patterns of the three periodic orbits marked ‘A’, ‘B’, and ‘C’in Figure 10 over four complete periods. Center: The real and imaginary parts of the solutionnormalized such that the projection onto U is always real and non-negative. Right: The variables z (solid, blue) and z (dashed, red) calculated from the numerical solution. N HH2
We turn our attention to the behavior near the bifurcations at N HH2 and N HH2 . Consideringfirst the dynamics of the reduced ODE system with N − N HH2 positive, but small, we find abranch of periodic orbits that shrinks to a point as N ↘ N HH2 ; an elliptical HH bifurcation. Ournumerical experiments show that the periodic orbits along these branches reach their maximumamplitudes when z and z are purely imaginary. We therefore set the phase of our orbits suchthat ( z ( ) , z ( )) = i ( y , y ) , and use y and y in the subsequent plots.Figure 12 shows the new branch of Lyapunov periodic orbits that arises as N increases acrossthe bifurcation value N HH2 . We show both the analog of Figure 7b, and also y vs. y . Figure 13is the equivalent for the PDE. Note we use different values of N and N in the two figures, sincethe bifurcation values are not exactly the same.In the above computations, the total amplitude is close to the bifurcation value, and the values of ∣ z ( t )∣ and ∣ z ( t )∣ remain small in comparison with N or N , meaning that all periodic orbits alongthese branches are essentially small perturbations of the dipole mode U − . Over the entire branchesdepicted in these bifurcation diagrams, z and z vary nearly harmonically as their evolution is onlyweakly nonlinear. B even The two RFPs U ( N ) and U − + − ( N ) are each stable for all values of N and thus should eachpossess two families of Lyapunov periodic orbits. In fact, these have already been calculated. Eachhas one family within the even invariant subspace B even itself, as seen in Figure 4. Each also has31 p y + y P e r i od (a)N=0.43N=0.45N=0.5 -0.2 0 0.2 y -0.3-0.2-0.100.10.20.3 y (b) Figure 12: Branches of periodic orbits of the reduced system that arise in the HH bifurcation at N HH2 . Note that because of symmetries, the point on the curve at ( y , y ) and (− y , − y ) representthe same periodic orbit. p y + y P e r i od (a) N=0.45N=0.5N=0.55 -0.2 0 0.2 y -0.3-0.2-0.100.10.20.3 y (b) Figure 13: Branches of periodic orbits of the NLS/GP system that arise in the HH bifurcation at N HH2 . 32 family of Lyapunov periodic orbits that extends out of this subspace, and these families havebeen computed in Figure 10. The points on these curve at the boundary between the shaded andunshaded regions are precisely the NNM solutions and the branches extending out from them arethe other family of Lyapunov periodic orbits.Of the other two RFPs on the B even discussed in Section 4.2, RFP U +++ has no imaginaryeigenvalues and thus no Lyapunov periodic orbits, while RFP U +0+ has two imaginary eigenvaluesand thus one family of Lyapunov periodic orbits, which encircle this point within B even as shownin Figure 4. The remaining four RFPs arise in a pair of saddle-node bifurcations as described in Section 4.1.The normal form analysis near these bifurcations—see Section 8—partially describes the branchesof RPOs. The eigenvalue count in Figure 5 shows that the branches U +00 and U ++0 should beaccomapanied by, respectively, two branches and one branch of Lyapunov periodic orbits. Similarly,the branches U − +0 ( N ) and U − ++ ( N ) , have, respectively two branches and one branch of Lyapunovperiodic orbits near the saddle-node bifurcation at N a2 . In addition the branch U − +0 ( N ) firstdestabilizes and then restabilizes in a pair of HH bifurcations.The normal form argument of Section 8 misses an important feature of the dynamics which wedescribe below. U +00 and U ++0 We first consider the Lyapunov branches of periodic orbits through U +00 and U ++0 , which appearin a saddle-node bifurcation at N = N a ≈ . U +00 . Figure 14(b) shows a magnified view of the small blackrectangle at the lower left. It shows that the long-period and mixed periodic branches do notsimply intersect, and cannot be considered as separate entities. Each apparent intersection is infact a pseudo-intersection: the numerically generated branch follows the path of the long-periodbranch before taking a sharp turn and following the mixed-periodic branch.Figure 14(a) shows three such long/mixed periodic branches, each drawn in a different colorand each a simple closed curve. These branches do intersect the short-periodic branches in a seriesof marked points. Since each point on these curve denotes the initial condition for a periodic orbitof an ODE with unique solutions, the two branches must pass through the same periodic orbit(or fixed point). Since the orbit on the long/mixed periodic branch has a longer period, this orbitmust consist of an integer number m of copies of the short-period orbit. We call such orbits m -foldorbits.The mixed-periodic branches were shown in Section 8.3 to arise in bifurcations when the longperiod is an integer multiple of the short period. The resonant values of N at which this occurs forthis system are provided in Table 2. How the branches of fast/periodic branches change at thesebifurcations is explored in the computations summarized in the rest of Figure 14.33 .56 0.57 0.58 0.59 0.6 0.61 0.62 x (0) -0.18-0.16-0.14-0.12-0.1-0.08 x ( ) (a) N=0.2496 x (0) -0.171-0.17-0.169-0.168-0.167-0.166-0.165-0.164 x ( ) (b) N=0.2496 x (0) -0.25-0.2-0.15-0.1-0.05 x ( ) (c) N=0.258 x (0) -0.4-0.3-0.2-0.10 x ( ) (d) N=0.299 x (0) -0.4-0.3-0.2-0.10 x ( ) (e) N=0.3004 x (0) -0.4-0.3-0.2-0.100.1 x ( ) (f) N=0.31 Figure 14: (Color online) Branches of Lyapunov periodic orbits related to the points U +00 at variousindicated values of N . Coloring of branches is consistent between the subfigures. The dashed curveshows the short-period orbits. (a) A value of N < N < N . (b) Magnification of the region outlinedby a black rectangle, showing avoided intersections. (c) N < N < N , a kink and two 3-fold orbitshave appeared on the long-period branch. (d) With N slightly less than N an isola has appearedcontaining two 2-fold orbits. (e) At N even closer to N the isola approaches the long-period branchat the point U +00 . (f) For N > N , the point U +0+ has swapped branches with a 2-fold orbit. Key: ◻ : U +00 , ○ : U ++0 , ∗ : 2-fold orbits, ✡ : 3-fold, ☆ : 4-fold, ▽ : 5-fold, △ : 6-fold, + : 7-fold, and × :8-fold. 34ifurcation SN . . . N k . . . N k for which the long and short periods of the mode U +00 are in k ∶ N between N and N . The long-periodbranch through this point, marked in red, intersects the short-period branch both at U +00 and, atthe marked point, a 4-fold orbit. Similarly the next (yellow) curve contains a 4-fold orbit, two5-fold orbits, and a 6-fold orbit, and the next (green) curve has a 6-fold, two 7-folds, and one 8-foldorbit. There are (presumably) an infinite number of such branches not depictedAs N is further increased, these closed branches deform, and additional branches appear at theresonant values N k . The behavior is different depending on the parity of k . No branches disappearas N is increased, but we drop some from the figures for clarity of exposition. The behavior at N is typical of resonances when k is odd. The straight portion of long/mixed period branch developsa kink, and, with it, two additional intersections with the short-period branch; see Figure 14(c)which correspond to 3-fold orbits.Typical behavior for even values of k is shown next for k =
2. At some critical value of N justbelow N , an isola of periodic orbits bifurcates into existence. This branch intersects the short-period branch in two 2-fold orbits; Figure 14(d). At N = N , this isola touches the branch U +00 ;Figure 14(e). Finally for N > N , the point U +00 and the 2-fold orbit exchange branches. This isseen in the computation of the bifurcations at the k ∶ k .The bifurcations that occur in a neighborhood of N ≈ N k where the two frequencies are resonant,is long studied, e.g. by Schmidt or Duistermaat [38, 39]. However, we do not know of a study of asequence of such bifurcations as seen here. U − +0 and U − ++ At N = N a ≈ . U − +0 and U − ++ arise in a saddle-node bifurcation. The branch U − +0 has both short- and long-period Lyapunov families of periodic orbits while the branch U − ++ has only the short-period family. Figure 15 shows some of the branches of periodic orbits that existfor N = . N = .
67, and N = . iω bifurcation.At resonant values of N , displayed in Table 3, these branches collide with and reconnect with theLyapunov family of long-period orbits that emerges from U − +0 .Bifurcation SN . . . N k . . . N k for which the long and short periods of the mode U +00 are in k :1 resonance.The two 1:1 resonances are the HH bifurcations.Figure 15(a), for which N < N a , features only the short period and mixed-periodic families,each of which, locally, is a hyperbola. Figure 15(b) shows a computation with N < N < N . Thelong-period branch is, as in the elliptic case, broken up by a sequence of pseudo-intersections withmixed-period branches. Figure 15(c) shows a computation with N < N < N , demonstrating howthe branches re-connect following a resonant bifurcation.While these branches of periodic orbits look locally like hyperbolas, they extend far outside the35 .2 0.25 0.3 0.35 0.4 x (0) x ( ) (a) N=0.664 x (0) x ( ) (b) N=0.67 x (0) x ( ) (c) N=0.673 Figure 15: (Color online) Branches of Lyapunov periodic orbits through the points U − +0 and U − ++ .As N increases through the resonant values, the vertical branches, each of which includes a pairof k :1 periodic orbits, collide at the point U − +0 and reconnect into a pair of horizontally-orientedhyperbolas, with no k :1 periodic orbits. (a) N < N a . (b) N a < N < N < N . (c) N < N < N .Markers as in Figure 14. Since this bifurcation involves re-connection, not all curves are identifiedby color between panels. 36egion of the points U − +0 and U − ++ , and exist on both sides of the saddle-node bifurcation. Thesebranches extend far outside a region where the local normal form expansion of Section 8 applies.The branch U − +0 also undergoes a pair of HH bifurcations at larger values of N . These followthe pattern seen in Sections 7, 9.1, and 9.2. First, the periods of the short- and long-periodbranches approach a common value, and the angle between them at U − +0 decreases to zero at theHH bifurcation value N ≈ . N ≈ . The phenomena seen in this section has been largely described in two separate streams of research,but the pieces have not been assembled in this manner before, as far as we can tell. Schmidt [38]considers the bifurcations of the short and long periodic branches that occur in a neighborhood of k : l bifurcations where k > l ≥
1. While in Section 8.3, we considered only mixed periodic orbitswhere the long and short period orbits were in n :1 resonance, but the argument there works with n replaced by any rational kl >
1. Schmidt shows that, in fact, those periodic orbits exist in thefull system for k / l > l = k ≥
4, he provesa similar result, but with the appearance of a gap in the branch of long period orbits as seen inour numerics. For l = k = k =
3, the results are a bit more delicate but can be put incorrespondence with our numerical results. In particular, he explains the development of the kinkin the long period branch in the 3:1 resonance bifurcation shown in Figure 14(c). Duistermaat laterrecast and extended Schmidt’s results in a more geometric setting [39].The above results take a localized view in a neighborhood of the individual k :1 resonances. Gel-freich and Lerman [40] take a more global view of the Hamiltonian 0 iω bifurcation of system (8.1)for fixed parameter. In a re-scaled system, the plane ( q , p ) = ( , ) becomes a slow manifold thatpersists for small nonzero µ and is foliated by long periodic orbits, except for a sequence of annulargaps. These gaps occur precisely at the pseudo-intersections between the long periodic and mixedperiodic branches shown in Figures 14 and 15. They illustrate their result with an example thatis equivalent to Figure 14. The results of this paper are confined to a small neighborhood of theslow manifold and do not address the global topology of the branches nor the difference betweenthe elliptic and hyperbolic cases.
10 Conclusion
In the context of of a finite-dimensional reduced system, we have used normal forms to show howbranches of relative periodic orbits behave as the total amplitude in the system is raised. Wehave supplemented this with numerical continuation studies tracing out families of relative periodicorbits that lie beyond the reach of the local normal forms. Our goal has been to illustrate some ofthe behaviors found in this system beyond simple standing waves or numerically-generated Poincar´esections. An unexpected find has been the complex arrangement of periodic orbits near the twosaddle-node bifurcations. Further work will be necessary to more fully specify and to combineprevious results into a rigorous global description of the bifurcations of these branches.The NLS equation with a multi-well potential is an excellent model system for understandingin detail the nonlinear dynamics that occur in nonlinear dispersive equations. Previous studies of37he double-well potential showed rigorously that periodic orbits of a reduced system are shadowedfor long time by nearly-periodic solutions to NLS [11–13]. A motivation for this paper has beento enumerate some of the analogous orbits in the triple-well problem, setting the stage for similarrigorous work in this more complicated case. Sigal has shown that true quasiperiodic orbits toNLS/GP cannot exist [35], so the long-time breakdown of these nearly quasiperiodic solutions isalso an important question.While this paper has avoided talk of chaotic dynamics, they are of course ubiquitous, and wereexplored numerically in our previous work [9]. One place where we should be able to say somethingrigorous is in a neighborhood of the HH bifurcations. The normal form (7.4) is integrable toall orders, so standard Melnikov methods will not show any splitting of separatrices leading tohomoclinic chaos; however Gaivao and Gelfreich have showed that an exponentially small splittingdoes indeed occur [41] in systems with a nonsemisimple -1:1 resonance.
Acknowledgments
The author received support from NSF DMS–0807284. He is grateful for conversations with Casayn-dra Basarab, Panayotis Kevrekidis, Richard Kollar, Jeremy Marzuola, Dmitry Pelinovsky, ArndScheel, and Eli Shlizerman, and to Stephen Shipman for a reading of the final manuscript. Theauthor gratefully acknowledges the support and hospitality provided by the IMA during his visitwhich took place from September to December 2016.
References [1] A. Newell and J. Moloney,
Nonlinear Optics . Advanced Book Program, Westview Press, 2003.[2] X. Chen and J. Holmer, “Focusing quantum many-body dynamics: The rigorous derivation ofthe 1D focusing cubic nonlinear Schr¨odinger equation,”
Arch. Ration. Mech. Anal. , vol. 221,pp. 631–676, 2016.[3] L. Erd˝os, B. Schlein, and H.-T. Yau, “Derivation of the cubic non-linear Schr¨odinger equationfrom quantum dynamics of many-body systems,”
Invent. Math. , vol. 167, no. 3, pp. 515–614,2007.[4] L. Pitaevskii and S. Stringari,
Bose-Einstein condensation . No. 116 in Int. Ser. Monogr. onPhys., Oxford University Press, 2003.[5] E. W. Kirr, P. G. Kevrekidis, E. Shlizerman, and M. I. Weinstein, “Symmetry-breaking bifur-cation in nonlinear Schr¨odinger/Gross-Pitaevskii equations,”
SIAM J. Math. Anal. , vol. 40,pp. 566–604, 2008.[6] T. Kapitula, P. G. Kevrekidis, and Z. Chen, “Three is a crowd: Solitary waves in photore-fractive media with three potential wells,”
SIAM J. Appl. Dyn. Syst. , vol. 5, pp. 598–633,2006.[7] E. W. Kirr, P. G. Kevrekidis, and D. E. Pelinovsky, “Symmetry-Breaking Bifurcation in theNonlinear Schr¨odinger Equation with Symmetric Potentials,”
Commun. Math. Phys. , vol. 308,pp. 795–844, Oct. 2011. 388] E. W. Kirr, “Long time dynamics and coherent states in nonlinear wave equations,” arXivmath.AP , May 2016.[9] R. H. Goodman, “Hamiltonian Hopf bifurcations and dynamics of NLS/GP standing-wavemodes,”
J. Phys. A: Math. Theor. , vol. 44, p. 425101, 2011.[10] A. Sacchetti, “Nonlinear Schr¨odinger equations with multiple-well potential,”
Phys. D , vol. 241,pp. 1815–1824, 2012.[11] R. H. Goodman, J. L. Marzuola, and M. I. Weinstein, “Self-trapping and Josephson tunnelingsolutions to the nonlinear Schr¨odinger / Gross-Pitaevskii equation,”
Discrete Contin. Dyn.Syst. , vol. 35, pp. 225–246, 2015.[12] J. L. Marzuola and M. I. Weinstein, “Long time dynamics near the symmetry breaking bifurca-tion for nonlinear Schr¨odinger/Gross-Pitaevskii equations,”
DCDS-A , vol. 28, pp. 1505–1554,2010.[13] D. Pelinovsky and T. Phan, “Normal form for the symmetry-breaking bifurcation in the non-linear Schr¨odinger equation,”
J. Diff. Eq. , vol. 253, pp. 2796–2824, 2012.[14] R. Fukuizumi and A. Sacchetti, “Bifurcation and stability for nonlinear Schr¨odinger equationswith double well potential in the semiclassical limit,”
J. Stat. Phys. , vol. 145, pp. 1546–1594,2011.[15] M. Albiez, R. Gati, J. F¨olling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, “Directobservation of tunneling and nonlinear self-trapping in a single Bosonic Josephson junction,”
Phys. Rev. Lett. , vol. 95, p. 010402, 2005.[16] J. Yang, “A normal form for Hamiltonian–Hopf bifurcations in nonlinear Schr¨odinger equationswith general external potentials,”
SIAM J. Appl. Math. , vol. 76, pp. 598–617, 2016.[17] P. G. Kevrekidis, D. E. Pelinovsky, and A. Saxena, “When linear stability does not excludenonlinear instability,”
Phys. Rev. Lett. , vol. 114, p. 214101, 2015.[18] J. C. Eilbeck, P. S. Lomdahl, and A. C. Scott, “The discrete self-trapping equation,”
Phys. D ,vol. 16, pp. 318–338, 1985.[19] J. Carr and J. C. Eilbeck, “Stability of stationary solutions of the discrete self-trapping equa-tion,”
Phys. Lett. A , vol. 109, pp. 201–204, 1985.[20] H. Susanto, “Few-lattice-site systems of discrete self-trapping equations,” in
The DiscreteNonlinear Schr¨odinger Equation (P. G. Kevrekidis, ed.), Berlin, Heidelberg: Springer, 2009.[21] P. Baˇnack´y and A. Zajac, “Theory of particle transfer dynamics in solvated molecular com-plexes: analytic solutions of the discrete time-dependent nonlinear Schr¨odinger equation. I.Conservative system,”
Chem. Phys. , vol. 123, pp. 267–276, 1988.[22] V. M. Kenkre and D. K. Campbell, “Self-trapping on a dimer: Time-dependent solutions of adiscrete nonlinear Schr¨odinger equation,”
Phys. Rev., B Condens. Matter , vol. 34, pp. 4959–4961, 1986. 3923] B. Liu, L.-B. Fu, S.-P. Yang, and J. Liu, “Josephson oscillation and transition to self-trappingfor Bose-Einstein condensates in a triple-well trap,”
Phys. Rev. A , vol. 75, p. 033601, 2007.[24] S. Zhang and F. Wang, “Interference effects between three coupled Bose–Einstein conden-sates,”
Phys. Lett. A , vol. 279, pp. 231–238, 2001.[25] M. Johansson, “Hamiltonian Hopf bifurcations in the discrete nonlinear Schr¨odinger trimer,”
J. Phys. A , vol. 37, pp. 2201–2222, 2004.[26] P. Panayotaros, “Instabilities of breathers in a finite NLS lattice,”
Phys. D , vol. 241, pp. 847–856, 2012.[27] C. H. Basarab,
Hamiltonian Bifurcations in Schr¨odinger Trimers . PhD thesis, New JerseyInstitute of Technology, 2016.[28] S.-N. Chow and Y.-I. Kim, “Bifurcation of periodic orbits for non-positive definite Hamiltoniansystems,”
Appl. Anal. , vol. 31, pp. 163–199, 1988.[29] K. Meyer, G. Hall, and D. Offin,
Introduction to Hamiltonian Dynamical Systems and theN-Body Problem , vol. 90 of
Applied Mathematical Sciences . Springer, 2nd ed., 2010.[30] J. Moser, “Periodic orbits near an equilibrium and a theorem by Alan Weinstein,”
Commun.Pure Appl. Math. , vol. 29, pp. 727–747, 1976.[31] N. Burgoyne and R. Cushman, “Normal forms for real linear Hamiltonian systems with purelyimaginary eigenvalues,”
Celest. Mech. Dyn. Astr. , vol. 8, pp. 435–443, 1974.[32] H. W. Broer, S.-N. Chow, Y.-I. Kim, and G. Vegter, “A normally elliptic Hamiltonian bifur-cation,”
Z. Angew. Math. Phys. , vol. 44, pp. 389–432, 1993.[33] V. Gelfreich and L. M. Lerman, “Separatrix splitting at a Hamiltonian 0 iω bifurcation,” Regul. Chaotic Dyn. , vol. 19, pp. 635–655, 2014.[34] D. Viswanath, “The Lindstedt-Poincar´e technique as an algorithm for computing periodicorbits,”
SIAM Rev. , vol. 43, pp. 478–495, 2001.[35] I. M. Sigal, “Non-linear wave and Schr¨odinger equations I. Instability of Periodic andQuasiperiodic Solutions ,”
Commun. Math. Phys. , 1993.[36] T. Dohnal and T. Hagstrom, “Perfectly matched layers in photonics computations: 1D and2D nonlinear coupled mode equations,”
J. Comput. Phys. , vol. 223, pp. 690–710, 2007.[37] C. A. Kennedy and M. H. Carpenter, “Additive Runge–Kutta schemes for convection–diffusion–reaction equations,”
Appl. Numer. Math. , vol. 44, pp. 139–181, 2003.[38] D. S. Schmidt, “Periodic solutions near a resonant equilibrium of a Hamiltonian system,”
Celestial Mech. , vol. 9, pp. 81–103, 1974.[39] J. J. Duistermaat, “Bifurcations of periodic solutions near equilibrium points of Hamiltoniansystems,” in
Bifurcation Theory and Applications , pp. 57–105, Springer Berlin Heidelberg,1984. 4040] V. Gelfreich and L. M. Lerman, “Long-periodic orbits and invariant tori in a singularly per-turbed Hamiltonian system,”
Phys. D , vol. 176, pp. 125–146, 2003.[41] J. P. Gaivao and V. Gelfreich, “Splitting of separatrices for the Hamiltonian-Hopf bifurcationwith the Swift-Hohenberg equation as an example,”