Long time dynamics near the symmetry breaking bifurcation for nonlinear Schr\"odinger/Gross-Pitaevskii Equations

We consider a class nonlinear Schr\"odinger / Gross-Pitaevskii equations (NLS/GP) with a focusing (attractive) nonlinear potential and symmetric double well linear potential. NLS/GP plays a central role in the modeling of nonlinear optical and mean-field quantum many-body phenomena. It is known that there is a critical $L^2$ norm (optical power / particle number) at which there is a symmetry breaking bifurcation of the ground state. We study the rich dynamical behavior near the symmetry breaking point. The source of this behavior in the full Hamiltonian PDE is related to the dynamics of a finite-dimensional Hamiltonian reduction. We derive this reduction, analyze a part of its phase space and prove a {\it shadowing theorem} on the persistence of solutions, with oscillating mass-transport between wells, on very long, but finite, time scales within the full NLS/GP. The infinite time dynamics for NLS/GP are expected to depart, from the finite dimensional reduction, due to resonant coupling of discrete and continuum / radiation modes.


Introduction and Outline
The cubic nonlinear Schrödinger / Gross-Pitaevskii (NLS / GP) equation i∂ t u = (−∆ + V (x))u − |u| 2 u (1.1) plays a central role in the mathematical description of nonlinear optical and quantum many-body phenomena. In the context of nonlinear optics, u denotes the slowly varying envelope of a nearly monochromatic electromagnetic field propagating in a wave-guide, t the distance along the wave-guide and x ∈ R 2 the dimensions transverse to the wave guide [26,8]. At low intensity, light is guided to a higher refractive index region, corresponding to a potential well V (x). The Kerr nonlinear effect gives rise to an increase in the refractive index in regions of higher intensity, |u| 2 , and therefore a "deeper" effective potential well V (x) − |u| 2 . In the context of quantum many-body physics, NLS/GP emerges in the mean-field limit of many weakly interacting identical quantum particles obeying Bose statistics, as the number of particles tends to infinity [28,10,30] and other recent works. The potential V (x) governs the confining trap for the quantum particles. 1 In this paper we focus on a class of symmetric double-well potentials. A model to keep in mind is the two parameter family of symmetric double well potentials: which converges as σ → 0 + to a double-delta well at ±L. In figure 1 the case L = 3, σ = 1 is shown. 1 NLS / GP falls into a larger class of models: allowing for more general nonlinear terms, e.g. more general functions of |u| 2 or K nonlocal; see, for example, [21]. For g < 0 (focusing case) our analysis yields very similar results to those obtained above. For g > 0 (defocusing) our methods can be adapted to prove a shadowing result for trajectories of the finite dimensional reduction which arises.  Figure 1. Double-well potential, V (x; 3, 1). Superposed are the two eigenstates, even ground (ψ 0 ) and odd excited (ψ 1 ) states, of the Schrödinger operator, −∂ 2 x + V (x; 3, 1) Double-well potentials are of particular interest in optics as models of coupled parallel wave guides channeling light, which interacts through through evanescent tails. In the quantum context, these are natural and simple systems in which to study quantum tunneling. As discussed below, the combined effects of a confining double-well potential with focusing cubic nonlinearity lead to the phenomenon of spontaneous symmetry breaking of the ground state at sufficiently high optical power or particle number. See, for example, [40,4,19] for experimental studies of symmetry breaking. Our goal in this paper is to investigate the phase space dynamics of NLS/GP near the symmetry breaking point.
The equation NLS/GP, (1.2), is a Hamiltonian system and expressible in the form: (1.4) i∂ t u = δH δu * . H denotes the conserved Hamiltonian energy functional: The conserved squared L 2 norm (particle number / optical power) is denoted: We are interested in the dynamics near special classes of nonlinear bound states of NLS/GP. Nonlinear bound states are solutions of the form where Ψ Ω is spatially localized: Consider first the linear case of the linear eigenvalue problem: (1.7) (−∆ + V (x))Ψ Ω = Ω Ψ Ω , Ψ Ω ∈ H 1 (R).
In this case, there is a least energy ground state, ψ 0 with corresponding simple eigenvalue Ω 0 [29]. If the separation between wells is sufficiently large, then the ground state eigenfunction is a bimodal positive symmetric state, which reflects the discrete symmetry of the potential [12,34,17]; see figure 1. In addition, there is an anti-symmetric (odd) state, ψ 1 with energy Ω 1 , such that Ω 0 < Ω 1 < 0.
In the attractive / focusing nonlinear case, (1.1), the character of solutions, and the solution set varies with the solution norm. Indeed, if we consider the set of solutions of (1.6) on the level set (1.8) |Ψ Ω | 2 = N , we find that for large enough well-separation, there is a symmetry breaking threshold N cr [21]; see figure2: (1) If N < N cr there is a unique positive, symmetric and bimodal state. (2) For N > N cr , (modulo phase) there are three positive localized states: a symmetric state (which exists for all N > 0) and two are asymmetric states, biased respectively to the right and left wells. (3) As N increases beyond N cr this symmetry broken state becomes increasingly concentrated in one of the wells [7,21]. (4) The symmetric (bimodal) state is dynamically stable for N < N cr and unstable for N > N cr . For N > N cr the asymmetric states are stable.
That symmetry breaking occurs at sufficiently large values of N was studied variationally in [7] for the nonlinear Hartree equation. Their method can be adapted to a large class of equations, for which a ground state can be realized as a minimizer of a Hamiltonian, H, subject to fixed N . The above detailed portrait of the symmetry breaking transition and the exchange of stability among branches was studied in detail in [21]. In particular, for large well-separation, L,  Our goal is to explore the detailed general dynamics near the symmetry breaking transition. Toward a formulation of precise results, we first introduce a class of double-well potentials in dimension one 2 . Following [12], start with a single rapidly decaying potential well centered at 0, V 0 (x), for which the Schrödiner operator H 0 = −∂ 2 x + V 0 has exactly one (simple) eigenvalue ω.  Figure 3. A numerical plot of the symmetry breaking in the soliton curve. Specifically, we plot here N versus Ω for a double-(Gaussian) well potential with separation parameter L = 3 and σ = 1. Note the stable curve is invariant under reflection about the y axis, hence the slightly thicker curve represents the resulting degeneracy for the asymmetric states. potential and define the Schrödinger operator (1.11) H L = −∂ 2 x + V L . There exists L 0 > 0 such that for L > L 0 , H L has a pair of simple eigenvalues, Ω 0 = Ω 0 (L) and Ω 1 = Ω 1 (L) and corresponding eigenfunctions ψ 0 (even) and ψ 1 (odd): Hψ j = Ω j ψ j , j = 0, 1; ψ j ∈ L 2 Ω 0 < Ω 1 < 0.
As noted above, the symmetry breaking threshold, N cr (L) (equation (1.9)) is exponentially small for large well-separation. Therefore, to study the dynamics in a neighborhood of the symmetry breaking point, it is natural use coordinates associated with the linear operator H 0 = −∂ 2 x + V . Throughout this result, we will assume that V is such that ψ 0 and ψ 1 are the only discrete eigenfunctions of H 0 and in addition that V is sufficiently smooth and decaying as defined in [39] and discussed in Appendix B as to guarantee dispersive estimates to do perturbation theory. We note that these assumptions will be satisfied by double delta wells, for which we have done our numerical computations. These ideas will be explored further in a forthcoming note [9].
Thus, we expand the solution as follows: u(x, t) = c 0 (t)ψ 0 (x) + c 1 (t)ψ 1 (x) + R(t, x), ψ j , R(·, t) = 0, j = 0, 1. (1.12) By orthogonality, we have Referring to this decomposition, we now give an overview of the paper: (1) In section 2 we express the NLS / GP as an equivalent dynamical system governing c 0 (t), c 1 (t) and R(x, t). This Hamiltonian system, equivalent to NLS / GP, has the form of two equations governing discrete nonlinear "oscillators", coupled to an equation for a wave field, R(x, t).
(2) In section 3 we study the finite dimensional reduction governing c 0 and c 1 , obtained by dropping all terms which couple the oscillator and field variables. This reduction is a finite dimensional Hamiltonian system with conserved Hamiltonian: Remark 1.1. Note also that the reduction above is what we get if we make the Ansatz (1.12) into the Lagrangian of NLS/GP, that is restrict L N LS−GP to the bound state manifold, and obtain the equations of motion as in [15].
Use of these two invariants facilitates an analysis of the finite dimensional phase space; trajectories with prescribed values of H and N lie on a two-dimensional surface. We then discuss some of the rich dynamics of this finite dimensional reduction and our goal is to prove their persistence, for non-trivial time scales , for with the full PDE, NLS/GP.
For the dynamics on level set N ∼ N F D cr ≈ N cr , we establish the following behavior: (a) N < N F D cr : There is an elliptic fixed point, corresponding to the stable symmetric state of NLS/GP for N < N cr . A neighorhood of this fixed point is foliated by stable time-periodic solutions. (b) N > N F D cr : The fixed point for N < N F D cr persists, but transitions from being a stable elliptic point to an unstable saddle. This corresponds to the unstable symmetric-bimodal state for N > N cr . At N = N F D cr , two new elliptic equilibria bifurcate and, a neighborhood of each is foliated by stable time-periodic solutions. These stable time-periodic oscillations correspond to stable oscillations around the stable asymmetric standing waves for N > N cr . (c) N > N F D cr : There are periodic solutions, outside a separatrix, which encircle both new equilibria. In the physical configuration space, these correspond to soliton transport from one well to the other and back, continuing periodically. Furthermore, one can quantify the "energy barrier" that must be exceeded to dislodge a "soliton" from localization about one of the wells. Energy thresholds, such as that described above play an important role in transport of energy in inhomogeneous and discrete systems. Natural directions to pursue beyond this work are transport in systems with many wells and the Peireles-Nabbaro barrier for motion of localized coherent structures discrete lattice systems.
(3) In section 5 we prove that certain phase space structures for the finite dimensional dynamical systems, persist in character for the full NLS/GP system, on long, but finite time scales. Stated nontechnically, Theorem: For any sufficiently small amplitude periodic solution about equilibrium an equilibrium state of the finite dimensional reduction (N above or below the bifurcation threshold), there is a solution of the PDE (NLS/GP), whose projection into the finite dimensional phase space, shadows this finite dimension orbit on very long time scales.
A precise statement is given in Theorem 1; see also Figure 4. The time scales on which these results hold enable us to see nearly periodic oscillations for the PDE on long time scales, i.e. through many, many oscillations, but not on an infinite time scale. Indeed, we do not expect the persistence of such oscillations on infinite time scales due to nonlinear coupling of bound to to radiation modes for the full system, for example, [36,11]. We conjecture that in the infinite time limit, the soliton executes a (radiation) damped oscillation to some stable nonlinear bound state, corresponding to the damped oscillatory decay to a stable equilibrium of the finite dimensional reduction. Evidence is presented in section 6, where numerical simulations are discussed.
We remark that our current theorem does not apply to perturbations of general "large" periodic orbits of the finite dimensional reduction. More detailed information on the Floquet theory of the linearized equations about general periodic orbits is still needed. This is currently being investigated. (1) The spaces H s (R n ), L p (R n ), W k,p are the standardly defined Sobolev integration spaces.
(2) We have the L 2 inner product: f, g = R n fḡ.

Formulation of NLS/GP as a coupled finite-infinite dimensional system
In this section we derive an equivalent formulation of NLS/GP, appropriate for studying the exchange of energy between the bound and radiative parts of the solution. We substitute the decomposition (1.12) into NLS / GP and, to the resulting equation, apply the projection operators P 0 , P 1 and P c to obtain the coupled system Here, We have that a ijkl = 0 ⇐⇒ i + j + k + l = 0 mod2. For simplicity we take Initial conditions for (2.1) are: c j (0) = ψ j , u(·, 0) , j = 0, 1; R(·, 0) = P c R(·, 0).
The system (2.1) may be viewed as an infinite dimensional Hamiltonian system, comprised of two coupled subsystems: one finite dimensional, governing bound state degrees of freedom described by c 0 and c 1 , and a second, infinite dimensional, governing a dispersive wave field, R(·, t) = P c R(·, t). In the next section, we focus on the finite dimensional truncation of (2.1) obtained by setting R = 0. After obtaining a detailed description of the phase space of this truncation in section 3, we then turn toward proving the long-time persistence of structures within the finite dimensional system, within the full infinite dimensional problem.
Before embarking on this path, we conclude this section with an alternative coordinate description of (2.1). These coordinates prove very useful in understanding the bifurcations within the finite dynamical subsystem, as N is varied, and its participation within the infinite dimensional dynamics.
2.1. Alternative coordinates. We shall require the following decomposition, proved in Appendix A.
Substitution of this ansatz into NLS/GP, we see where F is determined as in (2.4). Hence, A = Im(π 0 F ), α =θβ + Ω 1 β + Im(π 1 (F )), Note, there will be linear terms in R andR contained in F , which we must handle very carefully in our eventual iteration argument.
Proof. These follow directly from the computations in Appendix A. 9 For convenience, we rewrite the system (2.11)-(2.15) in compact form as: For simplicity, we take F b = P c F b and F R = P c F R in the sequel.

Phase space of the finite dimensional Hamiltonian truncation of NLS/GP
In this section we study the finite-dimensional system obtained by setting the dispersive part of the solution, R(x, t), equal to zero in (2.1). We denote the solution of the resulting system by (ρ 0 (t), ρ 1 (t)) ∈ C 2 : Symmetry breaking for a system of this type arising from a general class of defocusing nonlinearities was considered recently in [31].
This is a two degree of freedom Hamiltonian system, with time-translation and phase invariances inherited from NLS/GP. The associated time-conserved Hamiltonian and L 2 (optical power or particle number) functionals are: In terms of H, system (3.1) can be expressed in Hamiltonian form In terms of the alternative coordinates of section (2.1), hence the system (3.1) takes the form: Recall that, for simplicity, we have set a ijkl = 1 for all i, j, k, l = 0, 1. Note that θ completely decouples from the A, α, β equations, meaning we have where now the θ equation is decoupled to givė One can verify changing coordinates in (3.2-3.3), or directly from (3.6), the time-conserved quantities: We obtain a closed system for (α, β) using that N = A 2 + α 2 + β 2 is conserved. Then, we may reduce the system toα In addition, the system has the conserved quantity Now, let us define the matrices B eq andB eq to be those related to linearization about the equilibrium solution (A eq , α eq , β eq , θ eq ) and the decoupled reduced system (A eq , α eq , β eq ) respectively. Similarly, let M (t) andM be the resulting monodromy matrices for nearby time dependent periodic orbits (A(t), α(t), β(t), θ(t)) and the decoupled reduced system (A(t), α(t), β(t)) respectively.
In the following section, we actually discuss the relevant bounds on the operators e Beqt and eB t .
3.1. Bifurcation of equilibria for (3.7) and Symmetry Breaking in NLS/GP. Recall that N = N [A, α, β] is a constant of the motion for the system (3.6), corresponding to the physical quantities optical power or particle number. Thus it is natural to explore the nature of the phase space restricted to the level sets of N . We are interested in time-periodic states of frequency Ω, corresponding in the physical space to solutions of NLS/GP of the form u(x, t) = e −iΩt U . Thus, we transform the system to a rotating frame by setting (3.12) θ(t) = Θ(t) − Ωt and obtainα The states we seek are equilibria in this rotating frame. Thus, we have whose solutions we consider on the level set It is easy to observe an equilibrium corresponding to the Symmetric states: Via (2.10) we identify this with the symmetric ground state of NLS/GP; Due to its correspondence with the symmetric state, we refer to this equilibrium of the finite dimensional reduction, the symmetric equilibrium.
A second, bifurcating family can be found explicitly as follows. Define Remark 3.1. Had we set the nonlinearity coefficient g = −1 and the "interaction weights", a ijkl equal to one, we would have which is easily observed by comparison to a single well potential symmetric state for L sufficiently large (see [21]). 3 (3.20) again to eliminate Ω from (3.18). Thus we have, since β = 0

Using (3.20) we first eliminate Ω from (3.17) and obtain Ω
Solving for A and α we have the following equilibria, corresponding to symmetry broken states, which bifurcate for at N = N F D cr : Symmetry broken states: Via (2.10) we identify this with the asymmetric ground state of NLS/GP; Due to its correspondence with the asymmetric (symmetry-broken) states of NLS/GP, we refer to this equilibria of the finite dimensional reduction as asymmetric equilibria.

3.2.
Stability of equilibria; finite dimensional analysis. We consider the stability of the various solution branches obtained in the previous section. We rewrite the system (3.13), using the last equation to eliminate −Ω +Θ from the equations for α and β. Thus we have Note that in these coordinates the equations for α, β and A decouple from the equation for Θ. The finite dimensional system has a phase portrait, equivalent to (3.1), see Figure 5. In particular, we observe elliptic and hyperbolic equilibria, and periodic orbits. We now embark on detailed linear stability analysis of these states.
Linearization about an arbitrary solution gives the linearized perturbation equation Since the evolution of α, β and A, decouple from that for Θ, we consider the behavior of the reduced system Then,B is the 3 by 3 block in matrix of equation (3.28). In obtaining (3.30), we used that N cr = Ω 10 /2.
Linearized dynamics about the symmetric equilibrium state: For the symmetric equilibrium, as displayed in (3.22) we have Hence, For the reduced system, we have whose eigenvalues and implied linear stability character is as follows: Fig. 5. Thus, the symmetric state transitions from stable to unstable as N increases beyond N F D cr . Furthermore,B − can be diagonalized and the linear evolution is given by In other words, we have the bound The full linearized dynamics are governed by the matrix where B has the same eigenvalues asB, with λ = 0 now a generalized eigenvalue of multiplicity two, implying linear growth of δθ(t). Explicitly, we have Hence, it is clear Note here the leading order behavior is then a system which oscillates at the correct period and grows linearly only in the phase term. On a component basis, we may say Alternatively, we can use the spectrum of B − to define the matrix Then, we have Hence, we see Linearized dynamics about asymmetric equilibrium states, N > N F D cr : For either asymmetric equilibrium, displayed in (3.22) we have In this case, we substitute (3.25) into the expression forB in (3.30) and obtain: The eigenvalues ofB are: Therefore, the bifurcating asymmetric states are stable elliptic points. Also,B + is diagonalizable, resulting in In a manner analogous to the case of symmetric bound states, we have Hence, it follows Once again, we see the dominant behavior be oscillations in α, β, A with a linear growth accompanied by a factor of order N 1 2 in θ. As before, we can use the spectrum of B + to define the matrix Then, we have Hence, we see Remark 3.2. The analysis of e B±t provides partial understanding of the inherent instability in the phase of the finite dimensional periodic orbits. Specifically, an orbit (A 1 , α 1 , β 1 ) will oscillate with period T 1 while a different orbit (A 2 , α 2 , β 2 ) will oscillate with period T 2 . As a result, if two orbits begin quite close in phase, they will naturally oscillate out of phase with one another, which from the analysis above gives the linear shift in the δθ component. However, the oscillations will be purely described by the (A, α, β) system, meaning we must decouple the phase from the equation to get stability.
We have can the following Proposition 3.1. The system (3.11) has periodic orbit solutions.
Proof. The proof follows from standard level set techniques. Namely, we have the convex geometric curve H, which we slice at a particular value of N . Since the α, β coordinates are constrained to live on that level set, the only possible orbits are those that orbit with a period determined by N . Note, the equilibrium solution is the minimum of H with respect to N .
First of all, we have , 0) are local minima when N < N F D cr . Hence, for a fixed N , we select a bounded, closed curve for (α, β). By the Poincare-Bendixson Theorem, the possible curves are asymptotic to a limit cycle. However, the orbit is a closed curve, so it is necessarily periodic. Linearizing about the equilibrium solutions of (3.6) for N − N F D cr > 0, one may define where h 1 , h 2 small perturbations. Plugging this ansatz into (3.6) gives Hence, the period of oscillation near the equilibrium point is of the form This leads us to the following: Take (ρ eq 0 , ρ eq 1 ) be the corresponding equilibrium solution for (3.1). For any periodic solution ρ 0 (t), ρ 1 (t) of (3.1) such that we have (ρ 0 (t + T ), ρ 1 (t + T )) = (ρ 0 (t), ρ 1 (t)), (3.50) where Remark 3.4. For small perturbations of the equilibrium point (for either (N −N F D cr ) < 0 or (N −N F D cr ) > 0, it is for precisely the period T ± (N − N F D cr ) on which we must control the coupling to the continuous spectrum for the full solution to (1.1) in order to prove these finite dimensional structures are observable over many oscillations, see Section 5. In order to generalize our result to any periodic solution predicted by the finite dimensional dynamics, we must understand fully the period of each full oscillation. This will be discussed further in Section 7.

An Ansatz for the Coupled System
Starting with the finite dimensional system given in (3.7), we consider a periodic orbit, (α(t),β(t)), near equilibrium point and construct a periodic orbit of the extended system (3.6): Below we shall specify how near the equilibrium σ * need be.

See Appendix A for the explicit expressions in the system.
We denote byR the leading order part of R, driven by the periodic solution σ * (t): The correction toR is given by w, which satisfies: Introduce,M (t), a fundamental solution matrix for the system of ODEs with time-periodic coefficients: . We shall study the following system of integral equations for η(t), θ(t), w(x, t): We view a solution, ( η, R), of this system of integral equations as fixed point of a mapping, M: In the next section we formulate and solve this fixed point problem in a function space, which yields the existence of solutions to NLS/GP which "shadow" the periodic orbit, σ * (t), on the time scale of many periods. 21

Main Results
Recall that the period of the orbits described in Proposition 3.1 satisfies In this section we shall construct solutions to the full PDEs, which "shadow" the finite dimensional orbits for many periods.
Let τ > 0 denote a number to be chosen sufficiently small. And define the region of parameter space, about the symmetry breaking point, in which we conduct our study by For now 0 < γ < 1, but we shall place further constraints will be placed on γ. In terms of τ , we now describe the periodic solutions discussed in Propositions 3.1 and 3.2: Proposition 5.1. There exists τ 0 > 0, such that for 0 < τ < τ 0 , the system (3.7) has periodic solutions, which are small perturbations of the equilibria: where δ > 0 is to be chosen below.
The period of oscillations of the periodic solutions of Proposition 3.1 is We shall seek to contruct solutions on a time interval of the form
Take initial data for NLS/GP of the form where θ(0) ∈ R is chosen arbitrarily.
Then there exists a solution u(x, t) of (1.1) of the form whereR as in (4.2).

Proof of theorem.
In the sequel, we will often use the contracted notation where T * is given by (5.5), (5.1).
We now show that the map defined in section 4 is a contraction.
Thus, there exists a unique solution ( η, w) in B τ (I).

From Proposition 5.1 and the finite dimensional conservation laws
where |ǫ 0 | ∼ τ 1+δ 2 . Based on the expansion in the ansatz, we have . The term of largest order in this expansion isÃ The bounds on the remaining terms will follow similarly, so we look at In particular, we will show there exists δ 0 > δ, such that Hence, the leading order constant terms from the derivative are N F D cr − Ω 0 . We write By selecting γ > 2 3 , we ensure that all terms resulting from integration by parts are bounded by τ 1+δ1 for all t ∈ I = [0, τ − 1+γ 2 −ǫ ].
Using dispersive estimates, we have additional bounds forR given by the following for any Strichartz pair (p, q).
Proof. To begin, we first note that by the L 2 boundedness of wave operators discussed in Appendix B, we need only consider H R L 2 x . As in Lemma 5.3, we look at the the worst term, namely where we have used the Strichartz estimate (B.10) with dual Strichartz norm L 4 3 (I; L 1 ). The Strichartz estimate follows similarly.
Remark 5.2. Though such estimates do not arise in this work, let us also point out the following simple estimate 19) which is proved in [32] and may be applicable when trying to prove long or infinite time results on similar problems to the one studied here.
Note, by standard Sobolev embeddings and the contraction assumption, we have Since we have proper bounds onR, we must now bound From Appendix B, we have for any Strichartz pair (p, q) that where (p,q) is a dual Strichartz pair. In one dimension, it is useful to note we may takep = 4 3 + µ 1 and q = 1 + µ 2 for µ j > 0 small, j = 1, 2. In other words, we can be as close to the endpoint estimate ofp = 4 3 andq = 1 as we need to be. By reducing the system and including the lowest order terms from the expansions above, we can reduce to controlling a model problem of the form where χ ∈ S, and For our model problem, we now take (η, w) ∈ X, (η, w) X ≤ τ 1 2 +δ1 . In addition, take (η j , w j ) ∈ X, (η j , w j ) X ≤ τ 1 2 +δ1 for j = 1, 2. 27 It follows provided we choose In addition, it is clear from the analysis there that Next, for T erm η 2 we get Similarly, we have For T erm η 3 , we have Since T erm η 3 is independent of (η, w), it follows easily that T erm η1 3 − T erm η2 3 L ∞ = 0, meaning this term does not factor into the contraction.
The bound on T erm η 4 is of the form It follows immediately that For T erm η 5 , we have , which follows from previous constraints on γ. Once again, we have as well For T erm η 6 , we have using Proposition 5.3 T erm η For T erm η 7 , we have Furthermore, For T erm η 8 , we have 29 by previous constraints on γ. Again, it follows that We now study the map on the dispersive part, w. For T erm w 1 ,using the Strichartz estimates (B.10), we have by previous constraints on γ. Hence, For T erm w 2 , we have using bounds onR in H 1 (Lemma 5.4) which follows using γ > 7 9 + 14 9 ǫ. Hence, it follows that For T erm w 3 , using Lemma 5.4 once again we have which follows easily from previous constraints on γ. Hence, it follows directly that For T erm w 4 , we have which follows from previous constraints on γ. Furthermore, For T erm w 5 , we have which follows easily by previous constraints on γ. Moreover, For T erm w 6 , we have which follows easily from previous constraints on δ 1 . Hence, it follows directly that For T erm w 7 , we have from Lemma 5.4 T erm w which follows from previous constraints on γ. It follows directly that The bounds for T erm w 8 follows in a similar manner. For T erm w 9 and T erm w 10 , we have and T erm w Both of these terms are independent of (η, w), meaning the contraction mapping follows easily.
Once we have solved for ( σ, R), it is clear that the resulting function θ ∈ C(I) by construction. 6. Numerical Simulations 6.1. Polar Coordinates. As it will simplify the process of building initial conditions for numerically solving (1.1) that display the behaviors we study above, let us discuss here an alternative set of coordinates for (3.1). Namely, we set ρ 0 = r 0 e iθ0 and ρ 1 = r 1 e iθ1 . This leads to the following system of ODE's: where ∆θ = θ 1 − θ 0 . Given the system above, we can say that the bifurcation of stability occurs at As we are interested in the behavior quite near the bifurcation point, we define new parameters ǫ 0 , ǫ 1 and n such that (1 + cos(2∆θ)). (6.9) It should be noted, using the conservation laws we can write the system for ǫ 1 and ∆θ independently (1 + cos(2∆θ)) (6.10) and hence analyze phase plane diagrams, see Figures 8 and 9. In particular, the behavior of ∆θ in these regions leads to interesting oscillatory behavior as shown in the phase diagrams featured in Figures 8 and  9. For n > 0, a simple calculation shows that the equilibrium solutions occur for ǫ 1 = n 2 , ∆θ = kπ for k ∈ Z.
For n > 0, we have trapped orbits near the values ∆θ = kπ for k ∈ Z. This is a manifestation of the orbital stability of these mixed states. However, oscillations between wells can be generated by a large enough phase shift to leave the region where such trapped orbits occur. These oscillations are then large, in particular ǫ 1 must reach some ǫ max before decreasing.
For n < 0, we see the oscillations still attain a maximum at ∆θ = kπ, however their amplitude approaches 0 with the initial ǫ 1 (0). This is a manifestation of the stability of the symmetric state in this regime. One may ask if such finite dimensional Hamiltonian dynamics would appear in the infinite dimension dynamics of the PDE. For this, see Figures 10,11 and 12 for numerical evidence of their existence for long times. To begin, we locate the symmetry breaking point for a particular system. To do this, we use spectral renormalization with an asymmetric initial function to find the symmetry breaking point, see Fig. 3.
In the remainder of this section, we run simulations with wells of the form ). (6.11) Similar results hold and the same numerical analysis tools are applicable for potentials with more regularity. Knowing the bifurcation point as discussed in [17], we can numerically integrate using a finite element method similar to that in [14], where scattering of soliton solutions across single delta function potentials was analyzed (see [2] for analysis of finite element methods for nonlinear Schrödinger equations without potential). It is quite simple to adapt the method presented there to allow for a double-well potential (for delta functions or smoother potentials). The initial data is generated by finding the lowest entries of the spectrum of the discretized representation of H = −∆ + V in the Galerkin approximation, which is an operation embedded in many numerical software programs. For simplicity, we use the eig function from Matlab. Then, we may numerically solve the PDE system (1.1) with initial data corresponding to that necessary for the three types of oscillation described in Section 3. Note, one could also use the solitons from the spectral renormalization code (see [6]) as initial data quite easily, however these represent true nonlinear structures and we wish to observe structures derived from the finite dimensional dynamics, which we only expect to persist on finite time scales due to the nonlinear structure. The orbital stability of the 33 nonlinear objects is an interesting question in its own right and was explored in [21]. The equilibrium point of our dynamical system is in fact the finite dimensional part of a soliton solution, so there is no question the orbital stability of the soliton and the long time existence of oscillations near an equilibrium point are related. The phase plane diagrams for the finite dimensional dynamics are plotted using the MATLAB software program pplane7 [3].  Figure 11. At top, a numerical plot of the phase plane diagram for n = N − N F D cr > 0 with specific points chosen along a closed orbit which shows localization of the mass on one side of the well. Below, numerical plots of the absolute value of the solution to Equation (1.1) at various times with initial data such that n = N − N F D cr > 0 and ∆θ(0) = 0. The plots correspond to points a, b, c, d respectively from the specified orbit.
Finally, by taking a system such that |Ω 0 − Ω 1 | is comparably large (or well-separation distance, L, comparably close but sufficiently large to guarantee the hypothesized discrete spectrum), we can observe for large perturbations coupling to the continuous spectrum and hence decay of a fully oscillatory solution to a ground state in a finite time. In general, the mass dispersion is rather rapid, hence after a prescribed number of time steps determined by the computational domain, we cut off the solution near the origin and continue solving with the cut-off initial data. See Figures 13, 14, 15 and 16 for various computed time evolutions of the (radiation) damped oscillations observed in such a system.

Conclusion and discussion
The fact that we may observe oscillations between potential wells in systems such as (1.1) is not a new result, however we have been able to give a representation of the phenomenon in terms of a classical oscillatory system. In future work on double well potentials, the authors hope to prove long time stability 35 for oscillations far from the equilibrium point and optimize the time of existence proof by having better control of the damping caused by coupling to the continuous spectrum. In addition, these pseudo-bound states represent possible solutions resulting from the problem of scattering of solitons across double well potential wells, which in the high velocity limit have been recently studied in [1].
Finally, the authors would like to point out this question of oscillation and the resulting dynamical systems becomes more challenging and interesting as the number of wells is increased, see [20]. In particular,    as one increases the number of wells to ∞, the phase shift required to see oscillation from one well to the next might give insight into the celebrated Peireles-Nabbaro barrier for discrete nonlinear Schrödinger systems, see [27], [22], [18].