Effective dynamics for N -solitons of the Gross-Pitaevskii equation
EEFFECTIVE DYNAMICS FOR N -SOLITONS OF THEGROSS-PITAEVSKII EQUATION TREVOR POTTER
Abstract.
We consider several solitons moving in a slowly varying external field. Weshow that the effective dynamics obtained by restricting the full Hamiltonian to the finitedimensional manifold of N -solitons (constructed when no external field is present) providesa remarkably good approximation to the actual soliton dynamics. That is quantified asan error of size h where h is the parameter describing the slowly varying nature of thepotential. This also indicates that previous mathematical results of Holmer-Zworski [8]for one soliton are optimal. For potentials with unstable equilibria the Ehrenrest time,log(1 /h ) /h , appears to be the natural limiting time for these effective dynamics. Wealso show that the results of Holmer-Perelman-Zworski [7] for two mKdV solitons applynumerically to a larger number of interacting solitons. We illustrate the results by applyingthe method with the external potentials used in Bose-Einstein soliton train experimentsof Strecker et al [14]. Introduction
In many situations a wave moving in a slowly varying field, that is, a field describedby a potential whose derivatives are much smaller than the oscillations/width of the wave,can be described using classical dynamics. This is the basis of the semiclassical/short waveapproximation, perhaps best known in the case of the linear Schr¨odinger equation,(1.1) ih∂ t u = − h ∂ x u + V ( x ) u , where V is an infinitely differentiable potential. A typical result concerns a propagation ofa coherent state u ( x,
0) = exp (cid:18) ih (cid:0) ( x − x ) ξ + i ( x − x ) / (cid:1)(cid:19) , maximally concentrated near the point ( x , ξ ) in the position-momentum space. In thatcase,(1.2) u ( x, t ) = a ( x, t ) exp (cid:18) ih ϕ ( x, t ) (cid:19) + O ( h ) , < t < T ( h ) , where Im ∂ x ϕ >
0, Im ϕ ≥
0, andIm ϕ ( x, t ) = 0 ⇒ x = x ( t ) , ∂ x ϕ ( x, t ) = ξ ( t ) , where ( x ( t ) , ξ ( t )) satisfy Newton’s equations:(1.3) x (cid:48) ( t ) = ξ ( t ) , ξ (cid:48) ( t ) = − V (cid:48) ( x ) , x (0) = x , ξ (0) = ξ . a r X i v : . [ m a t h . A P ] O c t TREVOR POTTER (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) Figure 1.
A side-by-side comparison of the effective dynamics versus theexact solution of (1.6) for 4 solitons with the potential W ( x ) = − e cos x .The plot on the left shows the absolute value of the solutions up to time t = 1. The plot on the right shows the real part of the solutions at times t = 0 . t = 1. Compared to the solutions in Figures 4 and 5, much lessdiscrepancy between the two solutions is visible.The time of the validity of (1.2), T ( h ), depends on the properties of the flow (1.3), and ingeneral it is limited by the Ehrenfest time ,(1.4) T ( h ) ∼ log (cid:18) h (cid:19) , see [1] for a recent discussion on the case of one dimension.The approximation (1.2) means that the solution is concentrated for logarithmically longtimes on classical trajectories. The phase ϕ and the amplitude a can be described veryprecisely and a can be refined to give an asymptotic expansion – see [6] for an earlymathematical treatment and [13] for more recent developments and references.In this paper we consider the Gross-Pitaevski equation, which is the cubic non-linearSchr¨odinger equation with a potential:(1.5) ih∂ t u = − h ∂ x u − u | u | + V ( x ) u . It provides a mean field approximation for the evolution of Bose-Einstein condensate in anexternal field given by the potential V ( x ) – see the monograph [12] and references giventhere. Questions about propagation of localized states are also natural in the setting of(1.5) and have been much studied. One direction is described in a recent monograph [2].In this note we present a numerical study of multiple soliton propagation for (1.5) andshow that it can be described very accurately using a natural effective dynamics – seeFigure 1. That effective dynamics is based on mathematical results of Holmer-Zworski[8] and Holmer-Perelman-Zworski [7] and we refer to those papers for pointers to earliermathematical works on that subject. FFECTIVE DYNAMICS FOR N -SOLITONS 3 Following the convention of earlier papers – see Fr¨ohlich et al [4] – we rescale equation(1.5) so that the parameter h is in the potential which is now slowly varying:(1.6) i∂ t u = − ∂ x u − u | u | + V ( x ) u , V ( x ) = W ( hx ) . For V ≡ N -soliton solutions: u ( x, t ) = q N ( x, a + tv, v, θ + t µ + v ) , µ ) ,a, v ∈ R N , θ ∈ ( R / π Z ) N , µ ∈ R N + , (1.7)where the construction of q N = q N ( x, a, v, θ, µ ) will be recalled in § V (cid:54)≡ u ( x,
0) = q N ( x, a, v, θ, µ ) , the exact dynamics (1.7) is replaced by u ( x, t ) = q N ( x, a ( t ) , v ( t ) , θ ( t ) , µ ( t )) + O ( h ) , (1.8)where the precise meaning of the error, and its optimality, will be described below. Theparameters of the multisoliton approximation solve the system of ordinary differential equa-tions: ˙ v j = − µ − j ( ∂ a j V N + v j ∂ θ j V N ) , ˙ a j = v j + µ − j ∂ v j V N , ˙ µ j = ∂ θ j V N , ˙ θ j = v j / µ j / µ − j v j ∂ v j V N − ∂ µ j V N , (1.9)and where V N ( a, v, θ, µ ) def = 12 (cid:90) R V ( x ) | q N ( x, a, v, θ, µ ) | dx . Although somewhat complicated looking, the equations (1.9) have a natural interpretationin terms of Hamiltonian systems: they are the Hamilton-Jacobi equations for the fullHamiltonian of (1.6) restricted to the symplectic 4 N -dimensional manifold of N -solitons –see § V ≡ < t < C log(1 /h ) /h .In other words, the equations (1.9) give the minimal exact effective dynamics valid up tothe Ehrenfest time log(1 /h ) /h , where h is the parameter controlling the small variation ofthe potential, see (3.1). In this work, we show that the approximation errors O ( h ) and theEhrenfest time bound are sharp. See [11] for a survey of soliton dynamics under integrablesystems that have been perturbed.One motivation for this study is the experimental and theoretical investigation of solitontrains in Bose-Einstein condensates [14]. We show that the effective dynamics described in § § N -soliton solutionsfor V ≡ §
3, the Hamiltonian structure of the equation and the derivation of theeffective equations of motion. In § TREVOR POTTER
Figure 2.
Four solitons with alternating phases bunching up and spreadingout in the potential V ( x ) = ( x/ . The full solution is plotted on the rightwith a bird’s eye view. The figures on the left are snapshots of that solution.Due to their alternating phases, the solitons repel and never pass througheach other.of solutions to (1.6) and draw some quantitative conclusions. Specific potentials similarto those in [14] are then discussed in §
5. We investigate effective dynamics for the mKdVequation in §
6. Finally, in § N -solitons for cubic NLS When V ≡
0, we recover the nonlinear cubic one dimensional Schr¨odinger equation,which has N -soliton solutions with explicit formulas that we now recall – see [3] for adetailed presentation of this completely integrable equation.We will construct functions q N ( x ) that depend on 4 N parameters: positions, velocities,phases, and masses:(2.1) q N ( x ) = q N ( x, a, v, θ, µ ) , a, v, θ ∈ R N , µ ∈ (0 , ∞ ) N . Put(2.2) λ j = v j + iµ j , γ j ( x ) = e iλ j x e i ( θ j − v j a j ) e µ j a j , and define matrices M ( x ) ∈ R N × N , M ( x ) ∈ R ( N +1) × ( N +1) , FFECTIVE DYNAMICS FOR N -SOLITONS 5 (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239)
500 x (cid:239)
100 exp(cos(x)) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239)
100 sech(5x) + x (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) Figure 3.
A gallery of potentials used for the numerical experiments. Sincethe solitons in the experiments have width approximately 1 /
10, the inter-esting potentials should have size approximately 100. This is suggested bythe rescaling (3.9). The potentials vary on a scale comparable to 1, hencethe effective h is approximately 1 /
10. The exception is the upper right plotwhere we intentionally chose a potential which will exhibit some failures ofeffective dynamics. In the analysis of errors, for instance in Figure 6, onlyrelative sizes of h matter.by(2.3) M jk ( x ) = 1 + γ j ( x )¯ γ k ( x ) λ j − ¯ λ k , M = M ( x ) γ(cid:126) where(2.4) γ = [ γ , · · · , γ N ] T , (cid:126) , · · · , ∈ R N . Finally,(2.5) q N ( x ) def = det M ( x )det M ( x ) . Remarkably, this gives a solution to (1.5) with V ≡
0, the N -soliton solution:(2.6) u ( x, t ) = q N ( x, a + tv, v, θ + t µ + v ) , µ ) . As one can see from the formula, some restrictions on the parameters apply, see [3].
TREVOR POTTER Effective Dynamics Equations
We consider potentials defined on R that are slowly varying in the sense that(3.1) V ( x ) = W ( hx )where W ( x ) is C in x , and | ∂ kx W ( x ) | ≤ C (1 + | x | ) N , k ≤ . where C and N are independent of h . This means that h is the parameter controlling theslow variation of V .To obtain an effective dynamics for the evolution we use the Hamiltonian structure ofthe equation. In the physics literature an approach using Lagrangians is more common –see for instance Goodman-Holmes-Weinstein [5] and Strecker et al [14]. In the mathematicstreatments [4],[8],[7] the Hamiltonian approach was found easier to use, which we followhere.The basic claim is that an approximate evolution of q N is obtained by restricting theHamiltonian flow generated by the Gross-Pitaevskii equation to the manifold of N -solitonsdescribed in §
2. The Hamiltonian associated with the Gross-Pitaevskii equation is(3.2) H V ( u ) def = 14 (cid:90) ( | ∂ x u | − | u | ) dx + 12 (cid:90) V | u | with respect to the symplectic form(3.3) ω ( u, v ) = Im (cid:90) u ¯ v . The manifold of solitons, M N , is 4 N -dimensional and equipped with the restricted sym-plectic form given by the sum of forms for single solitons:(3.4) ω M def = ω | M = N (cid:88) j =1 ( µ j dv j ∧ da j + v j dµ j ∧ da j + dθ j ∧ dµ j ) .H V restricted to M N is(3.5) H N def = H V | M N ( a, v, θ, µ ) = N (cid:88) j =1 (cid:18) µ j v j − µ j (cid:19) + V N ( a, v, θ, µ ) , (3.6) where V N ( a, v, θ, µ ) def = 12 (cid:90) R V ( x ) | q N ( x, a, v, θ, µ ) | dx . The effective dynamics is given by the flow of the Hamilton vector field of H N on themanifold M N . That vector field, Ξ H N , is defined using the symplectic form (3.4):(3.7) dH N = ω M ( · , Ξ H N ) . FFECTIVE DYNAMICS FOR N -SOLITONS 7 A computation based on this gives the following ordinary differential equation for theparameters a, v, θ and µ , called the effective dynamics:˙ v j = − µ − j ∂ a j H N − µ − j v j ∂ θ j H N = − µ − j ( ∂ a j V N + v j ∂ θ j V N ) , ˙ a j = µ − j ∂ v j H N = v j + µ − j ∂ v j V N , ˙ µ j = ∂ θ j H N = ∂ θ j V N , ˙ θ j − = µ − j v j ∂ v j H N − ∂ µ j H N = v j / µ j / µ − j v j ∂ v J V N − ∂ µ j V N . (3.8)We remark that one can scale the Gross-Pitaevskii equation (1.6) in the following way:Consider any function u ( x, t ), scaling parameter α , let ˜ x = αx , ˜ t = α t , and define the newfunction(3.9) ˜ u (˜ x, ˜ t ) def = 1 α u ( x, t ) . Then if u ( x, t ) satisfies (1.6) with the potential V ( x ), ˜ u (˜ x, ˜ t ) also satisfies (1.6) with thenew potential(3.10) ˜ V (˜ x ) = 1 α V (cid:18) ˜ xα (cid:19) . This means that if we deal with a soliton of width comparable with α , the potentials forwhich interesting dynamics should appear should have size approximately α − and theslowly varying factor replaced by h/α .The effective dynamics equations (3.8) scale similarly: if ( a ( t ) , v ( t ) , θ ( t ) , µ ( t )) satisfies(3.8) and we define ˜ x, ˜ t, and ˜ V as above, then (cid:16) ˜ a (˜ t ) , ˜ v (˜ t ) , ˜ θ (˜ t ) , ˜ µ (˜ t ) (cid:17) also satisfies (3.8) with(3.11) ˜ a (˜ t ) = αa ( t ) , ˜ v (˜ t ) = v ( t ) α , ˜ θ (˜ t ) = θ ( t ) , ˜ µ (˜ t ) = µ ( t ) α . The scalings (3.9) and (3.11) are related in the following way: if u ( x, t ) = q N ( x, a ( t ) , v ( t ) , θ ( t ) , µ ( t )) , then ˜ u (˜ x, ˜ t ) = q N (˜ x, a (˜ t ) , v (˜ t ) , θ (˜ t ) , µ (˜ t )) /α = q N ( x, ˜ a (˜ t ) , ˜ v (˜ t ) , ˜ θ (˜ t ) , ˜ µ (˜ t )) . Comparison of effective and exact dynamics
For given values a , v , θ , µ in R N , we consider the solution u ( · , t ) of (1.6) with ini-tial data q N ( · , a , v , θ , µ ) and the solutions a ( t ) , v ( t ) , θ ( t ) , µ ( t ) of the effective dynamicsequations (3.8) with initial values a , v , θ , µ . In the following discussions we will refer to u ( · , t ) as the exact solution and q N ( · , a ( t ) , v ( t ) , θ ( t ) , µ ( t )) as the effective dynamics.Holmer and Zworski [8] proved that in the case N = 1,(4.1) (cid:107) u ( · , t ) − q N ( · , a ( t ) , v ( t ) , θ ( t ) , µ ( t )) (cid:107) H = Ch − δ , for t < δ log(1 /h ) Ch ,
TREVOR POTTER
Figure 4.
A side-by-side comparison of the absolute value of the exactsolution of (1.6) versus the effective dynamics (3.8) for 1,2,3, and 4 solitonswith potential W ( x ) = −
100 sech (5 x )+10 x . The sharpness of the sech (5 x )term creates clearly visible discrepancy between the two solutions.where δ ∈ (0 , /
2) can be chosen, and where C depends only on the potential and initialvelocity of the soliton, but not on δ . The H norm measures the size of the function andits spatial derivative in L : (cid:107) v (cid:107) H def = (cid:107) v (cid:107) L + (cid:107) ∂ x v (cid:107) L , (cid:107) v (cid:107) L def = (cid:90) R | v ( x ) | dx . This norm measures the energy of the solution.The limiting time log(1 /h ) /h is the Ehrenfest time discussed int § N >
1. This is suggested by [7], whichproves the analagous theorem for the modified Korteweg-de Vries (mKdV) equations forcase N = 1 ,
2, see § N = 2, multiple FFECTIVE DYNAMICS FOR N -SOLITONS 9 soliton interactions are not theoretically understood. All these considerations provided astrong motivation for this numerical study.We note that for N >
1, it was conjectured in [7] that the error bounds (4.1) will holdnot only in the H norm, but for the H N norm, which measures the the size of a functionand its first N derivatives, where N is the number of solitons. This has been proven in themKdV case with N = 2 [7]. For our numerical experiments, we consider only the H norm.We present numerical simulations to show that the result (4.1) holds for N > § § § O ( h − δ ) error estimate in (4.1) holds as h → /h ) /h timescale, or Ehrenfest timescale, in § A numerical case study.
We consider initial data q N ( · , ¯ a N , ¯ v N , ¯ θ N , ¯ µ N ), where ¯ a N =( a , . . . , a N ) , N = 1 , , , v N , ¯ θ N , and ¯ µ N are similarly defined with( a , a , a , a ) = ( − , − . , , v , v , v , v ) = ( − , , , θ , θ , θ , θ ) = ( π/ , , − , − µ , µ , µ , µ ) = (17 , , , q N ( · , ¯ a N , ¯ v N , ¯ θ N , ¯ µ N )is close to 0 outside of ( − π, π ), our numerical domain (see § V ( x ) = − e cos x , see Figure 3. The factor −
100 is chosen to create a deep enough well so that the solutionsremain in ( − π, π ), but the potential is otherwise chosen arbitrarily.We compute the exact solution and the effective dynamics solution for N = 1 , , , t = 1, which is chosen to allow for multiple soliton interactions. We plot thesolution for 4 solitons in Figure 1, where we observe very little difference between the exactand effective dynamics solutions. An equally small amount of discrepancy between thesolutions was observed for N = 1 , , V ( x ) = −
100 sech (5 x ) + 10 x , see Figure 3. This potential is chosen to be outside of the slowly varying regime for whichthe effective dynamics give good approximations. This is due to the − (5 x ) term,which creates a sharp dip roughly the width of the solitons we are studying; thus we donot expect the exact solution to maintain its soliton structure very well. This causes the (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239)
505 xt = .35 (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) Figure 5.
The plot on the left is a different view of the exact solution with3 solitons shown in Figure 4. The plot on the right compares the real partsof the N = 3 exact solution with the effective dynamics solution at times t = 0 .
35 and t = 0 . x term ensures the solutions remain on the interval ( − π, π ).We compute the solutions for N = 1 , , , t = 0 .
7, which is again chosen toallow for multiple soliton interactions. Figure 4 displays these solutions and demonstratesthat the effective dynamics captures the true motion, regardless of the number of solitonsand regardless of multiple soliton interactions. In the experiments presented in this paper,we only consider N ≤
4, but we have observed good agreement between the exact solutionand the effective dynamics for N ≤
7. We did not investigate futher due to increasingcomputational time needed to solve (3.8).We note that for N ≥ N = 3, at approximately t ≈ . Quantitative study of the error as h → . We investigate the O ( h − δ ) errorbetween the exact solution and the effective dynamics on a fixed time interval. The estimatethat gives rise to (4.1) is(4.3) (cid:107) u ( · , t ) − q N ( · , a ( t ) , v ( t ) , θ ( t ) , µ ( t )) (cid:107) H ≤ Ch e Cht . If t = δ log(1 /h ) / ( Ch ), then the RHS reduces to Ch − δ . When dealing with fixed length oftime or even time of size O (1 /h ) the RHS is O ( h ) and that form of error will be shown tobe optimal. FFECTIVE DYNAMICS FOR N -SOLITONS 11 (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239)
10 log10(h) l og10 ( H e rr o r) a t T = . Convergence of solution of effective dynamics equations to exact solution as h tends to zero N=1N=2N=3N=4
Figure 6.
A log-log plot of the H error, relative to the H norm of theinitial data, between the exact solution of 1.6 and the N -soliton evolvingaccording to the effective dynamics equations 3.8, as function of h . Here,the potential is V ( x ) = W ( hx ), where W ( x ) = −
100 sech (5 x ) + 10 x . Forsmaller values of h , the slope of the lines approaches 2, in agreement withthe theoretical upper bound on the error in (1.8).We reconsider our second potential from § h : V ( x ) = W ( hx ) , W ( x ) = −
100 sech (5 x ) + 10 x , and explore the H error, relative to the H norm of the initial data, between the exactsolution and effective dynamics as h →
0. We expect that as h becomes small enough, theequation will enter the slowly varying regime and display O ( h ) error. Indeed, the log-logplot in Figure 6 demonstrates the error is bounded by C N h as h →
0, where the constant C N varies only slightly between different values of N .We fit the data from Figure 6 to a line using the 6 smallest values of h .N 1 2 3 4Slope 1.86 1.76 1.91 1.86 C N -1.07 -1.46 -1.27 -1.37Thus, we conclude that the error is approximately O ( h ).4.3. Ehrenfest time.
We now investigate the length of time for which the effective dynam-ics approximation is accurate. In (4.3), we recalled that the error is bounded by Ch e Cht and hence the approximation breaks down at the Ehrenfest time, T ( h ) ∼ log(1 /h ) /h . We have already verified for small h and fixed time this error behaves as O ( h ). Thus wefocus on observing exponential growth in the error as a function of time, and verifying thatit is of the form O ( e Cht ).For this we must choose a potential and initial data to exhibit exponential instability.We are motivated by Newton’s equations for V ( a ) = − a / a = v , ˙ v = a , v = v cosh t + a sinh t , a = a cosh t + v sinh t . r e l a t i v e H e rr o r timeerror between exact solution and effective dynamics h=.3h=.4h=.5h=.6h=.7h=.8h=.9h=1.0h=1.1 0.2 0.4 0.6 0.8 1 1.201020304050607080 B hlinear dependence of B on h Figure 7.
The plot on the left shows H error, relative to the H normof the initial data, between the exact solution of (1.6) and the effectivedynamics for a single soliton sliding down the concave potential V ( x ) = W ( hx ) , where W ( x ) = 1000 x . The error is plotted as a function of timeand for several values of the parameter h . On the right sight, B is plotted asa function of the parameter h , when the errors from the plot on the left arefitted to a curve of the form A ( e Bt + C ). We expect B to depend linearly on h .In this case, we have exponential instability of classical dynamics. This suggests choosingpotentials with a non-degenerate maximum and working near the unstable equilibriumpoints.Hence we will investigate solutions to (1.6) with potential V ( x ) = W ( hx ) , where W ( x ) = − x . Figure 7 below demonstrates exponential divergence between the exact solution to (1.6)and the effective dynamics for several values of h and a single soliton initial condition q ( x, . , , , A ( e Bt + C ), for the timeperiod when the soliton’s position was between x = . /h and x = 1 . /h . This range wasobserved to be a region where exponential increase dominated the error and before thesoliton approached the numerical boundary.In Figure 7 we observe a linear dependence of B on h , in agreement with (4.3). Thisindicates that for certain potentials the Ehrenfest time C log(1 /h ) /h is the appropriatebound for the length time we expect the effective dynamics to give a good approximationto (1.6). We note that in our experiments with other potentials we often observe a linearincrease in error which would correspond to a timescale of C/h instead of the Ehrenfesttime C log(1 /h ) /h. FFECTIVE DYNAMICS FOR N -SOLITONS 13 Figure 8.
Surface plots of the error between the exact solution of (1.6) andthe effective dynamics for a single soliton sliding down the concave portionof the potential V ( x ) = W ( hx ) , where W ( x ) = 1000 x , as in Figure 7.We have plotted the absolute value of the difference between the spatialderivatives between the two solutions.5. Application to Bose-Einstein Condensates
Strecker, Partridge, Truscott and Hulet [14] discovered that Bose-Einstein condensatesform stable soliton trains while confined to one-dimensional motion. When set into motionin a suitably chosen optical trap, a Bose-Einstein condensate forms multiple soliton forma-tions which exist for multiple oscillatory cycles without being destroyed by dispersion ordiffraction. We can observe this same behavior numerically, using the effective dynamicsequations. We choose a potential of the type described in [14] V ( x ) = (cid:16) x (cid:17) and set N = 4, which was the most frequent case in their experiment. Strecker et al inferredthat the repulsive behavior of the solitons indicated alternating phases. Their argumentwas based on considering a certain reduced Lagrangian.In our numerical experiment we put ¯ θ = (0 , π, , π ) and then set the four solitons inmotion with the same velocity near the center of the potential. Similar to [14], we observebunching and spreading of the soliton train for several oscillations. See Figure 2.6. Effective dynamics for the mKdV equation
The mKdV equation(6.1) ∂ t u = − ∂ x ( ∂ x u + 2 u ) , Figure 9.
The left plot shows a side-by-side comparison of the exact solu-tion of the mKdV equation (6.2) and the effective dynamics solution for 3solitons with potential b ( x ) = 300 cos x . No discrepancy between the twosolutions is visible. The right plot displays the exact solution from a differentangle.like the nonlinear Schr¨odinger equation, has soliton solutions and a Hamiltonian structure.Holmer, Perelman, and Zworski [7] derived effective dynamics equations for the mKdVequation with a slowly varying potential(6.2) ∂ t u = − ∂ x ( ∂ x u − b ( x, t ) u + 2 u ) , b ( x, t ) = b ( hx, ht )and proved a result analogous to the (4.1) for N = 1 ,
2: the H N error between the solutionof (6.2) and its associated effective dynamics with N -soliton initial data is bounded by(6.3) Ch e Cht , for t < Ch log 1 h , Similarly to § H error is O ( h ) as h → Numerical methods
We now describe the numerical methods we employ to compute the Gross-Pitaevskii PDE(1.6) and the ODE (3.8) arising from the effective dynamics. When comparing a solutionof (1.6) with a solution of (3.8), we refine our numerical solutions until the error betweensucessive refinements of solutions to the same equation is several orders of magnitudesmaller than the error between solutions of the two equations.
FFECTIVE DYNAMICS FOR N -SOLITONS 15 (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239)
10 log10(h) l og10 ( H e rr o r) a t T = . Convergence of solution of effective dynamics equations to exact solution as h tends to zero N=1N=2N=3N=4
Figure 10.
A log-log plot of the H error, relative to the H norm of theinitial data, between the exact solution to the mKdV equation (6.2) and the N -soliton evolving according to the effective dynamics, as a function of h .For smaller values of h , the slope of the lines approaches 2, in agreement withthe theoretical upper bound on the error. The theoretical upper bound hasonly been proven for N = 1 ,
2, but this figure gives evidence that it holds forall N .Numerically solving the ODE arising from the effective dynamics (3.8) necessitates com-puting q N ( x, a, v, θ, µ ) and its derivatives with respect to the parameters a, v, θ ,and µ effi-ciently. For this we note that an equivalent definition of q N in (2.5)(7.1) q N ( x ) = − (cid:126) M − γ , where (cid:126) , M , and γ are as in (2.5). Since iM is Hermitian, M − γ can be efficiently computedusing the Cholesky factorization. Differentiating q N numerically, for larger values of N , istoo costly. Instead we used (7.1) to obtain explicit formulas for the derivatives of q N andagain used Cholesky factorizations to efficiently compute them. With this we compute theintegrand in (3.6), and then numerically integrate it using the trapezoidal method. Oncewe can efficiently compute the RHS of the effective dynamics equations (3.8), the standardfourth-order Runge Kutta method was found to be suitable to solve to the ODE.In order to solve the PDE (1.6) we used a Fourier spectral method to study the evolutionon the numerical domain ( − π, π ). This requires our solution u ( x, t ) to be periodic in space,so we choose initial conditions such that u ( x, t ) decays to zero, to machine precision, beforethe endpoints − π and π . Arbitrary initial data can be handled by either extending thenumerical domain or rescaling the equation (see (3.9)).One difficulty arises in that a non-trivial potential W ( hx ) cannot be periodic for all h .However, the potential W ( hx ) need not be periodic on ( − π, π ) so long as the product W ( hx ) u ( x, t ) is periodic on ( − π, π ), which is achieved if u ( x, t ) decays fast enough atthe endpoints − π and π . If W ( x ) is periodic on ( − π, π ) and one wishes to consider a solution u ( x, t ) that doesn’t decay before the endpoints − π and π , the rescaling (3.9) maybe employed with α = h :(7.2) ˜ x = hx , ˜ t = h t , ˜ u (˜ x, ˜ t ) def = 1 h u ( x, t ) . Then if ˜ u (˜ x, ˜ t ) satisfies (1.6) with periodic potential ˜ V (˜ x ) = W (˜ x ) /h , u ( x, t ) also satisfies(1.6) with potential V ( x ) = W ( hx ).This rescaling also makes it clear that as h →
0, a soliton solution becomes sharperrelative to the potential. This requires higher resolution in order to apply our numericalmethod to solve the PDE (1.6), while the effective dynamics equations (3.8) are unaf-fected. Indeed, our numerical experiments confirmed that increased computational effortwas needed to resolve the PDE as h →
0, but not the effective dynamics ODE.We now describe the method to solve a general solution u ( x, t ) of (1.6) on a periodicdomain with periodic initial data and potential V ( x ). The Fourier modes ˆ u k ( t ) of a solution u ( x, t ) to (1.6) evolve according to(7.3) ∂ t ˆ u k = − i k ˆ u k + i ( (cid:99) uv ) k , v = | u | − V Discretizing space and replacing the Fourier Transform with the Discrete Fourier Transformgives rise to a finite dimensional system of ODE, which we now represent in the generalform(7.4) u t = L u + N ( u )where L is a stiff linear transformation corresponding to the first term of (7.3) (representedby a diagonal matrix in our case) and N is a non-linear operator from the second termof (7.3). To solve (7.4) we compared the fourth order implicit-explicit (IMEX) methodARK4(3)6L[2]SA proposed by Kennedy and Carpenter [10] with the exponential time dif-ferencing (ETD) method ETDRK4 proposed by Kassam and Trefethen [9]. The IMEXscheme update formula is(7.5) u n +1 = u n + ∆ t ( b ( k + l ) + · · · + b s ( k s + l s )) , where ∆ t is the time step and k i and l i are chosen such that(7.6) k i = L ( u n + ∆ t ( A i k + · · · + A is k s + ˆ A i l + · · · + ˆ A is l s ))(7.7) l i = N ( u n + ∆ t ( A i k + · · · + A is k s + ˆ A i l + · · · + ˆ A is l s )) . Here A, ˆ A are s × s lower triangular matrices with ˆ A having zeros along its diagonal.. Thisallows us to solve for the k i and l i one stage at a time, only inverting the diagonal linearoperators I − ∆ tA ii L . The implicit treatment of the L term mitigates the stiffness arisingfrom the k factor in (7.3), while the l i can be computed explicitly, so non-linear equationsinvolving N need not be solved.The ETD method, on the other hand, uses an exact formula for obtaining the next step u n +1 from u n based on solving the linear portion exactly:(7.8) u n +1 = e L ∆ t u n + e L ∆ t (cid:90) ∆ t e −L τ N ( u ( t n + τ ) , t n + τ ) dτ FFECTIVE DYNAMICS FOR N -SOLITONS 17 (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) (cid:239) l og10 ( e rr o r) a t T = Convergence of numerical solution of PDE as the timestep tends to zero 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 (cid:239) (cid:239) (cid:239) (cid:239) l og10 ( e rr o r) a t T = Convergence of numerical solution of PDE as function of computational time ETDRK4ARK4(3)6L[2]SA ETDRK4ARK4(3)6L[2]SA
Figure 11.
Log-log plots of the convergence of the fourth order schemesETDRK4 and ARK4(3)6L[2]SA as a function of the timesteps and compu-tational time, respectively. ARK4(3)6L[2]SA is slightly more efficient pertimestep, but ETDRK4 is significantly more computationally efficient.The integral in (7.8) can then be numerically approximated using matrix exponents of L and evaluations of N . Thus as with the IMEX method, we do not need to solve non-linearequations, and computations involving L (namely computing e ∆ t L ) are efficient because L is diagonal. Stiffness is mitigated by solving the linear portion of (7.4) exactly. Wefound that the ETDRK4 scheme computed a solution of a desired accuracy nearly twiceas fast as the ARK4(3)6L[2]SA scheme. While the ARK4(3)6L[2]SA scheme had a slightlysmaller error rate per step, more computations per step made it significantly less efficient.Neither method demonstrated any instability in the range of step sizes required for oursolutions. Below, we plot the convergence of the two schemes as the timestep goes to zero.To obtain the results in the figures, we used the same potential and initial data as in § W ( x ) = −
100 sech (5 x ) + 10 x . Acknowledgements.
The author was supported through the National Science Founda-tion through grant DMS-0955078 and by the Director, Office of Science, Computationaland Technology Research, US Department of Energy under Contract DE-AC02-05CH11231.The author would like to thank Jon Wilkening and Maciej Zworski for helpful discussionsand comments.
References
1. S. De Bi`evre and D. Robert,
Semiclassical propagation on | log (cid:126) | time scales , Int. Math. Res. Not. (2003), 667–696.2. R. Carles, Semi-classical analysis for nonlinear schr¨odinger equations , World Scientific Publishing Co.Pte. Ltd., Hackensack, NJ, 2008.3. L.D. Faddeev and L.A. Takhtajan,
Hamiltonian methods in the theory of solitons , Springer, Berlin,1987.4. J. Fr¨ohlich, S. Gustafson, B.L.G. Jonsson, and I.M. Sigal,
Solitary wave dynamics in an externalpotential , Comm. Math. Physics (2004), 613–642.5. R.H. Goodman, P.J. Holmes, and M.I. Weinstein,
Strong nls soliton-defect interactions , Physica D (2004), 215–248.6. G. Hagedorn,
Semiclassical quantum mechanics , Ann. Phys. (1981), 5870.
7. J. Holmer, G. Perelman, and M. Zworski,
Soliton interaction with slowly varying potentials , preprint. arXiv:0912.5122 (2008).8. J. Holmer and M. Zworski,
Soliton interaction with slowly varying potentials , IMRN (2008), Art.ID runn026, 36 pp.9. A. Kassam and L. N. Trefethen,
Fourth-order time-stepping for stiff pdes , Siam J. Sci. Comput. (2005), 1214–1233.10. C. A. Kennedy and M. H. Carpenter, Additive runge-kutta schemes for convection-diffusion-reactioneqautions , Applied Numerical Mathematics (2003), 139–181.11. Y.S. Kivshar and B.A. Malomed, Dynamics of solitons in nearly integrable systems , Rev. Mod. Phys. (1989), 763–915.12. L.P. Pitaevskii and S. Stringari, Bose-Einstein condensation , Oxford: Clarendon Press, 2003.13. D. Robert,
On the herman-kluk semiclassical approximation , arXiv 0908.0847, 2009.14. K.E. Strecker, G.B. Partridge, A.G. Truscott, and R.G. Hulet,
Formation and propagation of matter-wave soliton trains , Nature (2002), 150–153.
Department of Mathematics, University of California, Berkeley, CA 94720, USA
E-mail address ::