Self-trapping and Josephson tunneling solutions to the nonlinear Schr\"odinger / Gross-Pitaevskii Equation

We study the long-time behavior of solutions to the nonlinear Schr\"odinger / Gross-Pitaevskii equation (NLS/GP) with a symmetric double-well potential, continuing work of the 2nd and 3rd authors. NLS/GP governs nearly-monochromatic guided optical beams in weakly coupled waveguides with both linear and nonlinear (Kerr) refractive indices and zero absorption. The optical power ($L^2$ norm) is conserved with propagation distance. At low optical power, the beam energy executes beating oscillations between the two waveguides. There is an optical power threshold above which the set of guided mode solutions splits into two families of solutions. One type of solution corresponds to an optical beam which is concentrated in either waveguide, but not both. Solutions in the second family undergo tunneling oscillations between the two waveguides. NLS/GP can also model the behavior of Bose-Einstein condensates. A finite dimensional reduction (system of ODEs) well-approximates the PDE dynamics on long time scales. In particular, we derive this reduction, find a class of exact solutions and prove the very long-time shadowing of these solutions by applying the approach of the 2nd and 3rd authors.

1. Introduction. We study the long-time behavior of solutions to the nonlinear Schrödinger / Gross-Pitaevskii equation (NLS/GP) with a symmetric double-well potential in R 1 . Equations of NLS/GP-type, in one or more spatial dimensions, arise as models of many physical systems, notably (a) nearly-monochromatic guided optical beams in weakly coupled waveguides with both linear and nonlinear (Kerr) refractive indices and no absorption [2,21], and (b) the time-evolution of a Bose-Einstein condensate (BECs) confined by a magnetically-induced linear potential [4,5,6,7,24].
Specifically, we consider NLS/GP in one space dimension where V (x) denotes a double well potential. In the context of optics, u(x, t) denotes the complex-valued slowly varying envelope of the electric field for a nearly monochromatic stationary beam propagating in the "t" direction. The potential V (x) is obtained from the transverse refractive index profile of the two coupled waveguides and g < 0 is proportional to the Kerr nonlinear coefficient. In the macroscopic quantum setting, u(x, t) denotes a macroscopic (mean-field) wave-function, V , the magnetic trap and g, which is proportional to the microscopic two-body scattering length, can be positive or negative, depending on the underlying atomic species.
A key quantity conserved by solutions of (1) is the squared L 2 norm, corresponding in the two systems to (a) the optical power, conserved with propagation distance, and (b) the particle number, conserved under the time evolution. The aim of this paper is to describe how long term dynamics of specific solutions of equation (1) vary with N , extending the result of [17].

Background and prior results.
Linear theory of double wells. We take the double-well, V (x), to be bimodal, i.e. the sum of translates of a unimodal potential V 0 (x), We shall assume that the basic well, V 0 (x), is spatially localized and even with respect to x = 0: and supports exactly one discrete eigenpair, (Ω , ϕ (x)) which solves −∂ 2 x + V 0 (x) ϕ = Ω ϕ . Without loss of generality, we assume ϕ L 2 = 1.
For large and positive, detailed information for the double-well Schrödinger operator H = −∂ 2 x + V . can be deduced from the properties of the basic single well, V 0 (x) [12]. In particular, there is a well-separation distance, 0 > 0, such that if > 0 (weak-coupling), then the linear operator H has two simple eigenvalues, Ω 0 = Ω 0 ( ) and Ω 1 = Ω 1 ( ), with Ω 0 ( ) < W < W 1 ( ) and Ω 1 ( ) − Ω 0 ( ) = O e −c0 , for some c 0 > 0. By the symmetry of V , the corresponding eigenfunctions ψ j = ψ j (x; ), satisfying H ψ j = Ω j ψ j , j = 0, 1 may be taken to be L 2 -normalized and possess, respectively, even and odd spatial symmetry: For sufficiently large , these eigenfunctions are bimodal and satisfy Nonlinear standing waves. There has been a great deal of interest in the existence, stability, and bifurcations of stationary solutions to NLS/GP (1) of standing wave type: For fixed value N > 0, U = U (x; N ) and ω = ω(N ) denotes a solution of the nonlinear eigenvalue problem defined jointly by equation (2) and by Questions of both mathematical and physical interest are: 1. Given N > 0, classify the solutions of (4). 2. How does the set of solutions vary as N is increased?
In [15], it is shown that there exist solution branches which are continuations of the linear solutions as N → 0, i.e. such that for j = 0, 1, U j (x; N ) ∼ √ N ψ j (x) and ω → Ω j , and that these branches possess the same symmetries as their linear counterparts. In the case of focusing nonlinearity (g < 0), they find that the "ground state solution," the continuation of the linear ground state ψ 0 (x), undergoes a symmetry-breaking bifurcation as N is increased. Fix > 1 sufficiently large. Then, there exists a threshold power / particle number, N cr ( ), such that for 0 < N < N cr ( ), the only solutions to equation (4) are U 0 (x; N ) and U 1 (x; N ), but for N > N cr , there exists another pair of solutions U ± (x; N ), concentrated in the left or in the right well but not symmetrically in both. As N is increased further, the solution becomes concentrated more strongly in one well or the other. Moreover, for N > N cr , stability is transferred from the ground state U 0 (x; N ) to the asymmetric states U ± (x; N ).
The analysis of [15] is based on a Lyapunov-Schmidt reduction. For N > 0 small, solutions are decomposed into their projections onto the span of {ψ 0 , ψ 1 } and their projection onto the orthogonal complement. For large, one may view this system as a two-dimensional system of nonlinear algebraic equations, which is weakly coupled to an infinite dimensional system. The two-dimensional truncated system has a bifurcation diagram of the type described above and it can be shown that neglected (infinite-dimensional) corrections are small for large. In particular, if N cr ( ) denotes the approximate symmetry breaking threshold obtained from the 2-dimensional reduction, then we have that for some κ > 0, see [15], Equation (1.1). Remark 1. Throughout this paper we shall assume g < 0, the case of focusing nonlinearity. For the defocusing nonlinearity g > 0, an analogous result holds with one difference. It is the mode U 1 (x, N ) that undergoes a symmetry-breaking bifurcation as N is increased; see for example [26].
The standing-wave result has been generalized in several ways. Kirr, Kevrekidis, and Pelinovsky [14] perform a global bifurcation analysis for the class of symmetric double well potentials with a non-degenerate maximum. Yang, in [27,28], studies the detection and classification of symmetry-breaking bifurcations in NLS equations with more general nonlinearities. In [13], Kapitula et al. study the richer family of bifurcations for NLS with triple well potentials.
Time-dependent dynamics. To study the dynamics of solutions to NLS/GP near the symmetry breaking bifurcation it is natural to express the solution to the initial value problem as where c 0 (t) and c 1 (t) are time-dependent complex-valued amplitudes and ρ(x, t), the projection of u(x, t) onto span ⊥ {ψ 0 , ψ 1 }. The resulting system consists of ODEs for c 0 (t) and c 1 (t) coupled to a PDE solved by ρ(x, t) and is mathematically equivalent to NLS/GP.
In [17], the latter two authors pursue the following strategy. First, they derive a finite-dimensional model by neglecting the terms involving ρ(x, t) in the evolution equations for c 0 (t) and c 1 (t). This is reduced by symmetry to an ODE resembling the Duffing equation, which is studied by phase-plane methods. Second, they show that solutions to the full PDE system shadow the ODE solutions on very long time-scales.
The standing wave solutions U j (x; N ) and U ± (x; N ) correspond to fixed points of this reduced system. The stability of the standing waves corresponds to that of these fixed points, as shown in [15]. Since the reduced system is conservative and two-dimensional, stable fixed points are surrounded by families of nested periodic orbits. The result proven in [17] is that the periodic orbits sufficiently close to these stable fixed points are shadowed by solutions to (1) over long but finite times. The result of the present paper is to extend the shadowing theorem to other periodic solutions of the ODE that are not confined to small neighborhoods of stable fixed points. This is explained further using the phase portrait in Section 1.2 below.
In recent related work, Pelinovsky and Phan [23] use a different reduction ansatz that leads to the standard Duffing oscillator. They control shadowing of a wider class of orbits than in [17] using only energy-type estimates and Gronwall's inequality in large data and arbitrary regimes provided there exist two distinct eigenvalues for H . On the other hand, in the decomposition (5), the representation of the initial conditions is more straightforward. In the present work, we furthermore extend the results of [17] to include orbits outside the separatrix, expanding the class of orbits we can shadow to be much closer to that of [23].

Qualitative discussion of results.
The finite-dimensional model. In [17], by neglecting coupling to ρ(t), equation (1) is viewed as a perturbation of the following two degree of freedom Hamiltonian ODE system for the evolution of (c 0 (t), c 1 (t)). Under the change of variables [16,25] c 0 (t) = A(t)e iϑ(t) , the evolution of the overall phase ϑ(t) decouples from the evolution of the other three quantities, due to the phase invariance of the underlying physical system, to give for We can assume A ≥ 0 with A = 0 only on the invariant circle α 2 + β 2 = N . As a shorthand for these coordinates, we define a four-dimensional vector and its three-dimensional truncation A direct consequence of (7a) is the conserved quantity corresponding to the L 2 invariance of (1). The constraint (9) can be used to further reduce (7a) to the 2-dimensional system d dt where H DW (α, β) = Ω 10 2 β 2 + σ 2 α 2 + α 4 + α 2 β 2 and J = 0 1 −1 0 for Note that (11) is a family of Hamiltonian systems parametrized by N . The Hamiltonian H DW differs from the standard Duffing Hamiltonian, by the presence of the mixed term α 2 β 2 . The phase plane of system (10) is displayed in Figure 1. It is topologically equivalent to that of the Duffing system (13) for either sign of σ. For σ < 0, the phase plane is foliated by concentric closed curves enclosing the fixed point at the origin; see Figure 1(a). For σ > 0, the origin becomes unstable and two stable fixed points appear, one on either side of the origin; Figure 1(b). In this regime, the dynamics feature two types of periodic orbits. There are two families of smaller periodic orbits, encircling exactly one of the stable fixed points and a second family of periodic orbits encircling all three fixed points.   [17] and the pink (lighter) shaded regions, the domains of validity in this paper.
Control of large time behavior for solutions to (1) requires precise estimates of the period and amplitude of the periodic orbits in the reduced dynamical system. These estimates are straightforward for separable Hamiltonians of the form such as the Duffing oscillator. Although H DW is not separable, we nevertheless obtain closed-form periodic orbits in Section 3. These are used to obtain appropriate estimates on the periods and amplitudes required for the shadowing analysis in Section 4.
Interpretation of the ODE solutions. Our results describe orbits which have an interpretation in quantum and electromagnetic contexts. In the electromagnetic context, the orbits we study represent nearly monochromatic beams within neighboring wave-guides exchanging energy.
For quantum settings, by examining the form of the eigenfunctions (3), the solution ansatz (5), and the reduction (6), notice that when α > 0 and β = 0, the basis functions ψ 0 and ψ 1 interfere constructively in the left well, centered at x = −L and destructively in the right well, centered at x = L, so that most of the optical power is concentrated in the left well. When α < 0, this is reversed and most of the power is concentrated in the right well. Thus, for the two families of periodic orbits encircling only one of the two fixed points in figure 1(b), most energy stays in one or the other of the two potential wells. Solutions of this type are called self-trapped in the Bose-Einstein condensate literature [1].
For the two families of periodic orbits in which α(t) changes sign, both for N < N cr and for solutions outside the separatrix for N > N cr , the location and magnitudes of the maximum and minimum of the reconstructed solution alternate with a fixed period. In the subcritical case N < N cr of Figure 1(a), the phase difference between c 0 and c 1 changes slowly due to the closeness of the frequencies Ω 0 and Ω 1 . For small-amplitude periodic orbits, this is a manifestation of the common beating phenomenon, while for solutions further from the origin, we refer to this as nonlinear beating. Periodic orbits encircling all three fixed points in Figure 1(b) are known as Josephson tunneling solutions in the BEC literature, where the word tunneling refers to the fact that these orbits must cross a local maximum of the potential energy to travel between the two wells. Note that the nonlinear beating solutions speed up as they cross α = 0, while the Josephson tunneling solutions slow down. At large amplitude, there is no practical difference between Josephson tunneling and nonlinear beating solutions. Both self-trapped and Josephson tunneling solutions have been directly observed in Bose-Einstein condensates by Albiez et al. [1].
It is instructive, also, to construct approximate periodic solutions using equations (5) and (6) Figure 2 shows two nonlinear beating solutions in the subcritical regime N < N cr . The first is a small periodic orbit around to the origin, in which case the solution stays very close to the symmetric mode ψ 0 . In the second, the values of A and |α(t) + iβ(t)| are close, and the solution migrates almost completely from the right well to the left well each period.
Bottom row: absolute value of reconstructed field, c 0 ψ 0 + c 1 ψ 1 , in (5). Subfigure (a) shows a solution near the stable fixed point, and (b) shows a larger periodic orbit. Figure 3 shows four orbits in the supercritical regime N > N cr . The first two are self-trapped, and the last two are Josephson tunneling solutions. Subfigure (a) shows a solution close to the right fixed point, the field makes small oscillations about a steady asymmetric profile. Subfigures (b) and (c) show solution trajectories near the separatrix. Both these trajectories spend long times close to the hyperbolic fixed point α = β = 0, where the field is nearly symmetric, with short bursts to asymmetric states. In (b), inside the separatrix, all these bursts move toward one asymmetric state, but in (c), the field alternates between the two. Finally in subfigure (d), the solution makes larger swings between the two asymmetric states without pausing near the symmetric state. The earlier result of [17] essentially shows the existence of the solutions of type (a) in both the subcritical and supercritical regimes, while the present result allows for patterns similar to those seen in the rest of the figures. Direct numerical simulations of NLS/GP can be seen, for comparison, in [17]. We now give a non-technical statement of the main result generalizing Theorem 5.1 of [17], which applies only to periodic orbits lying sufficiently close to stable equilibria of (7a). A precise statement appears in Sec. 2, Theorem 2.1. Main Result: Let 0 < |N − N cr | be sufficiently small. For periodic solutions to the ODE system (10) of sufficiently small amplitude for N < N cr or N > N cr , as long as the periodic solutions are sufficiently bounded away from the separatrix in the case N > N cr , there are corresponding solutions to the NLS equation (1) which shadow these orbits of the finite dimensional reduction on very long, but finite, time scales.
The remainder of this paper is organized as follows: Section 2 contains a precise formulation of the main result, Theorem 2.1. Section 3 presents exact formulae for the periodic solutions of (10), from which we derive period and amplitude estimates in Section 4. These bounds are used to show that the periodic orbits or the finite dimensional dynamical system satisfy the assumptions of Theorem 2.1. The section concludes with a discussion of how to apply these estimates together with the analysis of [17] to prove the main theorem. Section 5 contains concluding remarks and a discussion of future directions. The proof in this paper takes advantage of the exact solutions, which are convenient, but are not an essential feature for such a shadowing result to hold. Appendix A details the use of Lie transforms to construct a normal form for system (10). This is a possible first step in completing the proof in the absence of exact solutions.
2. The main theorem. We begin by describing the regions of parameter space and phase space in which the results hold. Our results apply to the regime Recall from [17] that when 1, we have Hence, we will define an asymptotic parameter σ such that with 7/9 < γ < 1 as specified in Theorem 2.1 and take initial N such that The 7/9 here is not sharp, but is a remnant of the asymptotic techniques in [17].
The main result of this paper is the following theorem. Note that the key difference between Theorem 5.1 of [17] and our main result, Theorem 2.1 (below), is that we now omit Assumption 3 of [17]. We will use the statement and part of the proof of [17], Theorem 5.1 in Section 4.3.
Remark 2. The full system we solve in Theorem 2.1 can be found in equations (2.5) − (2.6) and (2.11) − (2.17), with all error terms written out in Appendix A of [17]. Up to a phase shift, the phase term θ(t) in equation (16) evolves under an equation similar to (7b) for ϑ(t) in equation (6), though it differs by a small interaction term stemming from error terms in the full PDE expansion. This leads to O(1) differences on the time scales considered. This is discussed in greater detail in appendix A and Section 3.1 of [17] and revisited in Section 4.3 below. The variables δ 0 and δ 1 are sufficiently small, but non-zero constants chosen such that the bootstrapping arguments of [17] hold. In particular, they arise from the fact that the error terms must remain of lower order than the dominant dynamics over the time scale we study.
The proof that there exist periodic orbits of the finite dimensional model that are shadowed by solutions to NLS is then provided by the following Lemma: Lemma 2.2. There exists σ 0 > 0 and δ > 0 (chosen as in [17]), such that for 0 < σ < σ 0 , the system (10) has periodic solutions (α,β)(t) which through equation (9) determine a periodic function A(t) such that: The periods of these orbits satisfy the bound (14) required in Theorem 2.1. Furthermore, the fundamental solution matrix M (t) for the linearized system (see systems (3.37), (3.39) of [17]) about these periodic orbits satisfies the bound (15) from Theorem 2.1.
We next embark on the proof of this theorem, beginning with constructing periodic orbits and bounding their periods.
3. Explicit construction of periodic orbits. To solve the ODE system (10), we first obtain a simpler form of the equations by the rescaling: In terms of the new variables p, q, the Hamiltonian (up to the addition of a constant) is given by The potential energy V (q) has double-well structure in the supercritical case ζ > 0 (N > N cr ). Note that the kinetic energy depends both on the position and momentum. When ζ > 0, the system has three fixed points, a saddle at q = 0, and centers at q = ± ζ 2 . Fromq = ∂H ∂p = (1 + q 2 )p, we obtain p =q/(1 + q 2 ). By conservation of Hamiltonian, a solution with initial condition (q(0), p(0)) = (q 0 , 0) satisfies Solving forq 2 yieldṡ We now proceed with the exact integration of (20).
Josephson tunneling and nonlinear beating solutions. Let us first analyze orbits exterior to the separatrix for N > N cr . By (19), this implies that q 0 − √ ζ > 0. Define a by a 2 ≡ q 2 0 − ζ. Equation (20) may be rewritten aṡ

Separation of variables and integration yields
where we have defined , and u 2 = q 2 1 + q 2 . This integral can be found in a table [11, integral 3.152.4] or further simplified by a trigonometric substitution Q = B sin θ. In either case, we find where F (δ, k) is the incomplete elliptic integral of the first kind given by Some calculation yields the formula To obtain q as a function of τ requires an identity that inverts the equation θ = F (δ, k). This inverse defines the Jacobi elliptic function: sn(θ, k) = sin δ.
An introduction to these quantities in the language of modern dynamical systems is [19], while a comprehensive handbook is by Byrd and Friedman [3]. The functions cn and dn are defined in a similar manner. From here, there remains only the algebra to invert all the above changes of variables, which makes use of the identity sn 2 (δ, k) + cn 2 (δ, k) = 1, and eventually arrives at q(τ ) = q 0 cn (ωτ, k) 1 + q 2 0 sn 2 (ωτ, k) .
The orbit q(t) has period given in terms of K(k), the complete elliptic integral of the first kind This is a decreasing function of q 0 for q 0 > √ ζ and fixed ζ. For the case N < N cr , we note that ζ < 0, and hence the analysis follows in an identical fashion.
Self-trapping solutions. The periodic orbits encircling only one fixed point in the supercritical case ζ > 0 may be found in the same manner. These orbits are of the form .
Then, once again we have .

SELF-TRAPPING AND JOSEPHSON TUNNELING SOLUTIONS TO NLS 237
4. Bounds on the monodromy matrix and the period for Josephson tunneling orbits.

4.1.
Verification of Hypothesis (H1) from Theorem 2.1 and amplitude bounds from Lemma 2.2. We now prove the period estimate (14) from the expression (22). Each orbit of (7a) has an initial condition of the form The exact solution α(t) and β(t) to system (7a) is given by equations (18). Then, A(t) is computed from (9). The values of q 0 and ζ in terms of σ are read off from these expressions.
To satisfy (23), we take for some C(σ) bounded from above but which we do not attempt to optimize here.
The complete elliptic integral has the asymptotic behavior K(k) → ∞ logarithmically as k 2 → 1, where k 2 − 1 in effect encodes the distance of the exact solution from the separatrix. Hence, we must choose q 0 sufficiently far from the separatrix such that K(k) is uniformly bounded with a computable constant dependent upon k 0 when |k| < k 0 < 1. Indeed, we observe via a Taylor expansion of k in the parameter ζ/2q 2 0 in equation (21) that More specifically, given the scaling √ ζ q 0 σ, equation (21) yields the approximation which for q 0 as in (24) and σ sufficiently small implies It follows that At this stage we note that since K is a decreasing function as k → 0, orbits that lie sufficiently close to the separatrix violate the period bound (14). Also, the choice of 7/8 is not sharp, but suffices to bound the term 3(1 + σ)/4 arising as the leading order of the Taylor expansion of k using (24).
To verify the amplitude bounds from Lemma 2.2 for q 0 satisfying estimate (24), we directly evaluate (18) and observe |α|, |β| σ 1 2 (25) which, using the conservation equation (9) (7) about an arbitrary periodic solution, with period T * . Separating the linearized evolution component into the coupled (α, β, A) system (7a) and the independent θ evolution equation yields Floquet's theorem says that there exists a matrix-valued T * -periodic function P (t) and a constant matrix B * , such that the fundamental solution matrix of subsystem (27a) with initial condition M (0) = I is M (t) = P (t)e B * t with P (0) = P (T * ) = I.
In [17], Section 3.3, the authors construct M * = e B * t at a stable equilibrium point and then restrict their analysis to nearby periodic orbits, for which they could prove bounds on M (t) perturbatively. Here, we explicitly construct the matrix M at any periodic orbit using the exact solution derived in Section 3. The linearly independent column vectors generating the matrix M are found by differentiating equation (7a) with respect to the canonical system parameters ζ, q 0 and t. These represent translation in the mutually orthogonal directions through energy space (translation perpendicular to the phase plane through energy space), within the energy plane onto a nearby orbit (translation in the normal direction to an orbit) and along a heteroclinic orbit respectively (translation in the tangential direction to an orbit).
To proceed, we evaluate the periodic solution χ(t) and then its derivatives ∂ t χ, ∂ q0 χ, and ∂ ζ χ using the exact solution. Each of these vectors solves equation (27b). A matrix whose columns are given by these three vectors, or, in fact by any three independent linear combinations of these solutions, defines a fundamental solution matrix. Of these three vectors, the first is T -periodic in time, while the second and third exhibit linearly-growing oscillations.
We first compute

SELF-TRAPPING AND JOSEPHSON TUNNELING SOLUTIONS TO NLS 239
At t = 0, only the second component is nonzero, so we can define the renormalized vector The vectors ∂ q0 χ and ∂ ζ χ can be computed similarly. As the expressions for them in general is rather long, we display here only their values at t = 0: Using from Section 4.1 that A(0) = O σ γ 2 , these may be renormalized to where C is a positive O(1) constant. We can then define the fundamental solution matrixM (t) = χ ren To prove (15) for the fundamental solution operator M (t), we prove first that the growth over one full period is bounded. Then, we prove that the variation within a single period is bounded. Proof. This lemma is essentially Proposition 3.3 in [17], but we recall the proof here for completeness. It is clear thatχ = (α(t),β(t),Ȧ(t)) is a solution to (27a), giving at least one Floquet multiplier λ 1 = 1. Similarly, differentiation with respect to the period gives a similar result, implying λ 1 = λ 2 = 1.
Note, using a similar approach to above, there exist at least two linearly independent vectors by using differentiation in t and T .
The above lemma shows that the n-times iterated monodromy operator satisfies the bound for some positive constant C. To show that the solution to (1) stays close to the finite dimensional orbit χ * on the time scale of shadowing, we require that the norm of M (t) for 0 < t < T does not grow too large to prevent us from applying the bootstrapping arguments in the proof of the main theorem from [17]. Hence, we require slightly better control on the growth of M than we have currently shown. Indeed, we require the following lemma bounding the operator within a the evolution of a single period.
Proof. Note, it suffices to prove this theorem forM in (28). To prove this, we may simply calculate the entries inM (t) and find the maximum magnitude of each as a function of time over all values of 0 ≤ t < T . To begin, we recall the functional form the exact solution for α(t): .
There are similar expressions for β(t) and A(t). We will make use of the fact that cn(z, k), sn(z, k), and their z-derivatives are bounded functions for all k and z, and that derivatives with respect to k exhibit linearly-growing oscillations for fixed k and increasing z. These bounds take advantage of the extensive literature on Jacobi elliptic functions [3,22]. Since we have chosen q 0 and ω such that k, K and K = dK dk are uniformly bounded, we can find similar bounds for β and A but must be mindful of the dependence on the asymptotic parameter σ in constructing the fundamental solution matrix. We observe that the fundamental solution matrix can be represented in terms of the column vectors defined in (29). Note, under such a choicẽ For t > 0, the implicit differentiation in the columns of M (t) introduces terms that are purely periodic as well as terms that have linear growth as observed in the orbits in citeMW, Proposition 3.3. This calculation relies upon boundedness of cn, sn and their derivatives in k.

SELF-TRAPPING AND JOSEPHSON TUNNELING SOLUTIONS TO NLS 241
and have bounds maximum bounds given by The terms that do grow linearly in t are of the form: ∂k ∂q 0 ωt, and Ω 3/2 Noticing that on the time scale of interest and keeping track of the renormalization constants, we observe the fundamental solution bound is given by, for instance, the largest of the time-dependent terms. Hence, in addition to (32). This is the desired bound for M (t) given 0 ≤ t ≤ T .

4.3.
Brief recap of the perturbative methods from [17]. Lemma 4.2 implies that the evolution within one period is also bounded, hence the Duhamel evolution operator, M (t)M −1 (s) for 0 < t − s < T , used in the perturbation theoretic arguments from [17] has the same bound. Coupling this implicit Duhamel evolution bound, the period bound from (23), and the amplitude bounds from (25) and (26) shows that assumptions (14), (15) from Theorem 2.1 and the amplitude bounds in Lemma 2.2 all hold for the exact solution. Hence, the dynamical solutions satisfy the assumptions of Theorem 5.1 in [17] for proving shadowing of the dynamical system for orbits outside the separatrix in (10) and large orbits can be shadowed for long times in (1)! Following the analysis in [17], we now write the system (27) with the orbit χ * given by an exact periodic orbit of the finite dimensional truncation: where the η A , η α and η β will account for the error terms in gluing the finite dimensional model equations into the full infinite dimensional model. This corresponds to a solution of the form: u(x, t) = e iθ(t) (Ã(t) + η A (t)) ψ 0 + [(α(t) + η α (t)) + i(β(t) + η β (t))]ψ 1 + R(x, t) with initial conditions u 0 (x) = e iθ(0) Ã (0) ψ 0 (x) + [α(0) + iβ(0)]ψ 1 (x) . Centered about χ * , equation (1) becomes the system: Here, F FD and F b represent the finite dimensional and infinite dimensional projections of the terms dependent only upon χ(t) in the ansatz. Also, G FD and G θ represent the interaction terms in the finite dimensional between R and χ, as given by the terms Error A , Error α , Error β , Error θ from Appendix A of [17]. The term F R is the interaction term in the infinite dimensional evolution, as it is in Appendix A of [17]. In particular, the evolution equations for η are similar to those for the finite dimensional linearization (δα, δβ, δA) in (27a) but now with an additional forcing term given by G FD . See Appendix A of [17] for full details. The estimates are derived in [17]. We decompose R =R + w, whereR the leading order part of R, driven by the periodic solution χ * (t) as defined in (17), and w is a correction. IntroducingM (t), a fundamental solution matrix for the system of ODEs with time-periodic coefficients ∂ t η = D χ F (χ * (t)) η.
allows system (34) to be written as a system of integral equations for η(t), w(x, t), and θ(t): We view a solution, ( η, w), of this system of integral equations as fixed point of a mapping M: ( η, w) = M( η, w).
Given the bounds proven above in the proof of Theorem 2.1, M is a contraction map in a particular Banach space in both x and t designed to optimally measure the dynamics of both η and w. To that end, we recall the Strichartz space where 2 p + 1 q = 1 2 and T * is given in Theorem 2.1. Following [17], we define the space We define B σ (I) ⊂ X(I) such that ( η, R) ∈ B σ (I) if and only if where δ 1 > 0 must be chosen in the course of the analysis.
Then the following proposition makes the desired contraction mapping precise. 2. There exists κ < 1 such that given ( η j , w j ) ∈ B σ (I) for j = 1, 2, Thus, there exists a unique solution ( η, w) in B σ (I).
The main result then follows by applying the asymptotic analysis from Section 5 of [17] to prove for example bounds of the form 5. Discussion, further problems, and a remark. We consider NLS/GP with a symmetric double well potential having sufficiently separated wells. We have shown that there are periodic solutions of NLS/GP which, on a long time scale, shadow periodic orbits of the finite dimensional ODE model both inside and outside the separatrix. An obvious question is whether similar behavior can be shown in systems with multiple potential wells. It is shown in [9,10] that periodic orbits and a class of quasi-periodic orbits arise due to a Hamiltonian Hopf bifurcation. Work is underway to show that these periodic orbits are shadowed in the full PDE model. Remark 3. The extension of the proof of [17] to the more general class of periodic orbits considered here was facilitated by a closed-form solution of the reduced system. It is instructive to consider how this was used: 1. Bounds on the solution's period are needed in order to obtain bounds on a Duhamel operator in the proof of [17], Page 21, Equation (4.4). The necessary bounds follow from the explicit solution. In the case where there is no explicit solution, for trajectories near the separatrix the dominant contribution to the period can be computed by an asymptotic analysis near the saddle fixed point at the origin. 2. In Lemma 4.1, the symmetries of the exact solution are used to show that the monodromy matrix has an eigenvalue λ = 1 of algebraic multiplicity three and geometric multiplicity two. This information is obtainable from the symmetry of the equation and does not require the exact solution. 3. The exact solution simplifies the estimate of the bounds of the various terms in the fundamental solution operator, equations (31)-(33). In the absence of exact formulas, these bounds require terms beyond first-order in the expansion of M (t). An alternate way to determine them is to first put the Hamiltonian system (10) into normal form up to some order, which agrees with the Hamiltonian of the Duffing oscillator to leading order. This is demonstrated in the Appendix using the method of Lie Transforms. Truncation in the normal form transformation introduces a secondary source of error that must be controlled carefully in any estimates.
Appendix A. Outline of approach using Hamiltonian normal forms. A key ingredient in the proof of shadowing is prove bounds on the fundamental solution matrix operator for a reduced system. While these can be studied via the explicit solutions of the reduced system arising for double wells, more generally explicit solutions are unavailable and therefore a different approach must be taken. For example, in the case where V (x) is a triple well potential, the orbits can only be calculated via perturbative methods [9,10]. In this appendix we give a brief sketch of the general normal form method and the results of its implementation for (10), in the regime where σ is a small parameter; see [18] for details. Consider a general Hamiltonian of the form H(Q, P ) = H 0 (Q, P ) + σf (Q, P ; σ), where σ 1, H 0 is quadratic (so that the associated differential equations are linear and thus integrable) and f has a finite or convergent power series. A normal form for the Hamiltonian (36) is an equivalent but simpler Hamiltonian obtainable by a suitable canonical change of variables (Q, P ) = Φ(q, p) = (q, p) + . . . under which new Hamiltonian K(q, p) = H(Φ(q, p)) = H 0 (q, p) + σg(q, p; σ). Roughly, the function g(q, p; σ) should contain as few terms as possible at each order in perturbation theory. It is well-known [18, §10.4] that for a Hamiltonian like (11) with leading-order part H 0 = Ω10 2 P 2 , the normal form Hamiltonian is K(q, p) = Ω 10 2 p 2 + σg(q; σ).
That is, the normal form equation is separable, with p-dependent kinetic and qdependent potential energy. Using perturbative technique based on the Lie transform, we can calculate the change of coordinates and the transformed Hamiltonian to arbitrary order using The three leading terms of K(q, p) are equivalent to the Duffing Hamiltonian (13), which has well-known solutions involving Jacobi elliptic functions [22]. These solutions then form the leading order parts of Poincaré-Lindstedt approximations to true periodic orbits of K(q, p), which can be computed to arbitrary order using computer algebra.
While the asymptotic series defining a normal form change of variables does not generally converge, there exist, in some cases, useful error estimates for finite truncations of the series, for example Giorgilli and Galgani [8], who prove an estimate of Nekhoroshev type [20]. By truncating the series to order r, they show that the error due to the normal form approximation can be controlled to within O(σ) for times of O(σ −r ). Further, as σ, the order r can be chosen in order to control the error over exponentially long time scales in 1/σ. Using the normal form approximation introduces approximations into the fundamental solution operator of the finite-dimensional system, which must be controlled in order to verify the bounds (14) and (15) which are discussed in Section 4. The time scales on which these estimates are proven must be reconciled with the time scales on which the normal form approximation is valid.