{'fraction_non_alphanumeric': 0.05088434624497217, 'fraction_numerical': 0.013995261446911676, 'mean_word_length': 4.1707977207977205, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Soliton Interactions in Perturbed Nonlinear Schroedinger Equations

We use multiscale perturbation theory in conjunction with the inverse scattering transform to study the interaction of a number of solitons of the cubic nonlinear Schroedinger equation under the influence of a small correction to the nonlinear potential. We assume that the solitons are all moving with the same velocity at the initial instant; this maximizes the effect each soliton has on the others as a consequence of the perturbation. Over the long time scales that we consider, the amplitudes of the solitons remain fixed, while their center of mass coordinates obey Newton's equations with a force law for which we present an integral formula. For the interaction of two solitons with a quintic perturbation term we present more details since symmetries -- one related to the form of the perturbation and one related to the small number of particles involved -- allow the problem to be reduced to a one-dimensional one with a single parameter, an effective mass. The main results include calculations of the binding energy and oscillation frequency of nearby solitons in the stable case when the perturbation is an attractive correction to the potential and of the asymptotic"ejection"velocity in the unstable case. Numerical experiments illustrate the accuracy of the perturbative calculations and indicate their range of validity.


Introduction
This paper is concerned with the asymptotic behavior of solutions of the initial-value problem for the perturbed nonlinear Schrödinger equation (NLSE) subject to the initial condition ψ(x, 0) = ψ 0 (x) for certain initial fields ψ 0 (x), in the limit when the perturbation term p[ψ, ψ * ] becomes formally small. The unperturbed problem, when p[ψ, ψ * ] ≡ 0 in (1), is well-known to be solvable [1] by an inverse scattering transform, one consequence of which is the existence of finite energy soliton solutions that are dynamically stable and robust with respect to collisions. The unperturbed NLSE arises in two different physical situations in modern optics [2]. Firstly, for high-speed telecommunications, (1) describes the propagation of light wavepackets along an optical fiber. In this interpretation, t is the spatial coordinate along the fiber and x is the retarded time variable for the signal; accordingly the solitons of the unperturbed problem (and usually also the solitary waves of the perturbed problem, when they exist) are called temporal solitons. The suggestion by   [3,4] that temporal solitons, being immune to dispersion, might serve as bits in a high-speed data stream has since generated a large body of work, much of which is comprehensively reviewed in [5,6]. Secondly, for photonic switching devices, (1) describes the stationary envelope of monochromatic light waves in a planar waveguide under the paraxial approximation. Here x and t are both spatial variables, with t being the propagation direction and x being the transverse direction; accordingly the solitons of the unperturbed problem are called spatial solitons. In both of these applications, the cubic term |ψ| 2 ψ in (1) models the Kerr effect in which the refractive index of the material depends linearly on the local intensity of light. For weakly nonlinear effects, when intensities are not too large, this effect is dominant in isotropic materials like glass. This fact, along with the integrability afforded by neglecting p[ψ, ψ * ], makes the unperturbed NLSE one of the most important models in modern optics.
Of course, real materials can have a complicated dependence of refractive index on intensity, for which the Kerr effect is only an idealization. Modeling such phenomena requires introducing corrections to the coefficient |ψ| 2 in the cubic term of the NLSE. The perturbative term p[ψ, ψ * ] might also include corrections related to higher-order dispersion, the Raman effect, self-steepening of pulses, etc. In this paper we will consider only the influence of higher-order nonlinearity on solitons of the unperturbed NLSE. For spatial solitons in photorefractive media, such a perturbation can be the main factor influencing propagation. In particular, we take the perturbation in (1) in the form of a quintic term p[ψ, ψ * ] = σǫ 2 |ψ| 4 ψ , where ǫ > 0 is a small parameter and σ = ±1.
In view of the possibility of using solitons as bits in optical fibers or dynamically controllable switches in planar waveguides, it is of some interest to determine the effect of such a perturbation on the solitons of the unperturbed problem. If one considers an initial condition ψ 0 (x) that is a "snapshot" of a simple soliton solution of the unperturbed problem, then there are many approaches available to study the perturbed evolution. Because the unperturbed soliton is stationary in some Galilean frame, the main effect of p[ψ, ψ * ] will be an adiabatic adjustment of the soliton's amplitude and phase parameters. This fact, together with the simplicity of the form of the soliton solution, means that direct perturbative methods can be used to study their slow evolution. In particular, variational methods and multiscale methods applied directly to (1) often give valid results. These perturbative methods are dynamical in origin and capture effects on finite but long scales. Other methods can be used to answer infinite time questions concerning the persistence of solitary waves. In fact in the presence of quite general perturbations solitary waves continue to exist for arbitrary ǫ [7,8] and these can be expressed in closed form in some cases [2].
The presence of more than one soliton complicates the analysis. If the solitons are isolated then the field may be approximated as a sum of solitons plus a small error term, and the adiabatic coupling among the solitons may be calculated by several methods. Note that if the solitons are moving with respect to each other then they will always be in isolation except possibly for a short time. An early analysis of this kind was carried out by Gordon [9], who studied the exact two-soliton solution of the unperturbed NLSE for equal velocities. When the solitons are well-separated, there is an effective force between them (even in the unperturbed NLSE) that varies sinusoidally with their phase difference. This phase difference grows linearly in t if the solitons differ in amplitude. The force is therefore zero on average [10] and one expects periodic motion. This is a physical explanation of the mathematical fact that the intensity |ψ| 2 of the exact two-soliton solution for equal velocities is a periodic function of t. An extension of this argument to perturbed problems was given by Ankiewicz [11], who obtained a simple description of soliton interactions with the use of complex averaged potentials. Again, the essential assumption is that the solitons are well-separated in x, so that the field may be approximated as a sum of solitons. If the solitons are close to each other, nonlinear interference effects cause the field to adopt a form very different from the linear superposition of individual solitons, and therefore a different approach is needed. Often, one turns to numerics to study the interactions of solitons in various media (see, for example, [12,13,14]) without the restriction of the solitons being isolated.
In the scattering transform domain, where the dynamics of the unperturbed NLSE are trivial, a state in which two solitons are close to each other in x has the same spectrum as a state in which they are far apart. This suggests that for studying the influence of perturbations on multisoliton bound states (that is, several solitons traveling with the same velocity, represented by a collection of eigenvalues of the Zakharov-Shabat equations with the same real part) it is best to carry out the analysis in the transform domain using soliton perturbation theory [7,15,8]. With p[ψ, ψ * ] = 0, the evolution of the scattering data is no longer trivial, and thus the scope of possible dynamics in near-integrable systems like (1) is much greater than in the unperturbed NLSE, including effects like repulsion, attraction, and energy exchange among bound or colliding solitons. Other techniques that have been used to study these effects include the judicious use of conserved quantities [2], variational methods [16,17], "equivalent particle" approaches [18,19], and of course, numerics.
In this paper, we use soliton perturbation theory to study perturbations of the nonlinear potential in (1), for initial conditions ψ 0 (x) that are snapshots of multisoliton bound states of the unperturbed NLSE. With respect to treating the solitons in isolation, this is a worst-case scenario since in the unperturbed NLSE a tightly-bound state of solitons will remain so for all time. Nonetheless, it is a scenario of some interest, in particular for the quintic perturbation (2). If σ = +1, then it is known that the solution remains bound, and this case has been studied using conservation laws [20]. If σ = −1, then the bound state becomes destabilized. Recently it was shown [21] by simulations of (1) that the instability causes the bound state to divide into isolated solitons that are ejected from the origin with nonzero relative velocities. On the time scales over which this splitting occurs, the solitons do not appear to exchange energy. In mathematical terms, each eigenvalue in the bound state ensemble, originally confined to the imaginary axis (zero velocity), appears to slowly "grow" a real part while its imaginary part remains fixed. Once the solitons escape, they no longer interact and the velocities no longer change. The wave guidance properties of Y-junctions engineered from such splittings of spatial solitons have also been analyzed [22].
By considering the relative velocities to be small, we will find an integral formula that expresses the asymptotic velocity difference between a pair of initially co-propagating solitons destabilized by the quintic perturbation (2) with σ = −1. Along the way, we will write down a coupled system of differential equations that describes the interaction of any number of solitons under more general perturbations over long time scales. These equations are just Newton's equations for a system of interacting particles in one space dimension; the particle coordinates have the interpretation of the soliton centers of mass. The force is translationally invariant, conserves the total momentum, and is also proportional to σ, so the forces giving rise to attraction and repulsion are related just by a change of sign. For the interaction of two solitons, the problem may be reduced to a single degree of freedom, the relative separation of the solitons. The force law scales simply with the (fixed) amplitudes, which have the interpretation of masses. The result is a one-parameter family of problems indexed by a normalized effective mass. If the separation is small in the attractive case σ = +1, the force is nearly linear and the frequency of motion becomes a function of the normalized effective mass. We calculate this frequency, a quantity that is connected with the vibrations of solitons that are infinitely close, a limit opposite to the well-separated case.
Our paper begins in §2 with a review of the theory of the scattering transform for the Zakharov-Shabat eigenvalue problem and of the inverse theory that holds in the reflectionless case. We also recall the derivation of the exact equations of motion in the transform domain corresponding to the perturbed NLSE (1). Then, in §3 we consider perturbations of the form p[ψ, ψ * ] = ǫ 2 W (|ψ| 2 )ψ and apply multiscale perturbation theory to find asymptotic solutions of the equations of motion in the transform domain. The approximations are uniformly valid as ǫ ↓ 0 on expanding time intervals of length ǫ −1 , and are given in terms of solutions of Newton's equations for particles interacting in one dimension under a force law that has several universal features. In §4 we focus on the quintic perturbation (2) and study the interaction of two solitons. We reduce the problem to the motion of a single particle and then explicitly perform the averaging required to remove secular terms from the asymptotic expansion. This leaves the force law in the form of a 1D integral that we study numerically. We use it to compute the "ejection" velocity observed by Artigas et. al. [21] in the unstable case and the harmonic frequency of tightly-bound solitons in the stable case. Finally, we compare the results of perturbation theory with direct simulations of (1). The Appendix contains the more cumbersome formulae that nonetheless are among our main analytical results.
Regarding notation, we will use stars for complex conjugation, and matrices will be written with bold letters, except for the Pauli matrices

Exact Inverse Scattering Theory For The Perturbed NLSE
Here, we review the known inverse scattering theory for the Zakharov-Shabat eigenvalue problem to fix our notation. In general, we wish to consider (1) where p[ψ, ψ * ] is a polynomial in ψ, ψ * , and their x-derivatives. The field ψ is taken to be in the Schwartz space as a function of x.

Scattering data.
We will work with the scattering transform of ψ, a map that associates to the complex field ψ at each fixed time a set of "scattering data" from which ψ can be reconstructed by inverting the map. As is well-known, the advantage of this is that the time evolution of the scattering data corresponding to the time evolution of ψ is trivial when p ≡ 0. Consequently, when |p| ≪ 1, this proves to be a useful setting for perturbation theory. Fix t, and assume the complex function ψ(x, t) to be given. For λ ∈ R denote by M ± (x, t, λ) the 2 × 2 matrix solutions of the linear differential equation satisfying the boundary conditions M ± (x, t, λ) exp(iλσ 3 x) → I as x → ±∞ Since L is traceless, these boundary conditions guarantee that these matrices are unimodular for all x. For each λ there can only be two linearly independent column vector solutions of (4); therefore there is a matrix S(t, λ), λ ∈ R, the scattering matrix, such that The first column of M − (x, t, λ) and the second column of M + (x, t, λ) turn out to be boundary values of analytic functions for ℑ(λ) > 0, while the second column of M − (x, t, λ) and the first column of M − (x, t, λ) are the boundary values of analytic functions for ℑ(λ) < 0. Adjoining the second column of M + (x, t, λ) on the right of the first column of (5) and taking determinants gives which is therefore the boundary value of a function analytic for is the boundary value of a function analytic for ℑ(λ) < 0. Fix λ ∈ R. Then, from (4), M ± * = σ 2 M ± σ 2 , and thus S * = σ 2 Sσ 2 , so that S 22 = S * 11 and S 21 = −S * 12 . In particular, this means that as an analytic function for ℑ(λ) < 0, S 22 (t, λ) = S 11 (t, λ * ) * . Also, for λ ∈ R the fact that det(S) = 1 implies the normalization condition |S 11 | 2 + |S 12 | 2 = 1.
The analytic function S 11 (t, λ) with ℑ(λ) > 0 may have zeros λ 1 (t), . . . , λ N (t). The determinant formula (6) then shows that there exist complex numbers γ 1 (t), . . . , γ N (t) such that The conjugation symmetry of M ± (x, t, λ) for λ ∈ R, when extended to the complex plane, implies that at the complex conjugate points Hilbert transform theory can be used in conjunction with the normalization condition to express S 11 (t, λ) for ℑ(λ) > 0 in terms of its zeros and the values of S 12 (t, λ) on the real axis [23]: The so-called "trace formulae" that equate certain functionals of the potential ψ to functionals of the scattering data will be useful below. In particular, we will use the formula: This functional (not to be confused with the perturbation p[ψ, ψ * ]) has the interpretation of the total momentum of the wave function ψ(x, t). For the unperturbed problem, as well as in the presence of many physically important perturbations, the total momentum is a constant of motion.

Reconstruction of the potential in the reflectionless case.
The miracle of inverse scattering theory is that for each fixed t, the potential ψ(x, t) can be recovered from its scattering data, namely the reflection coefficient S 12 (t, λ) for λ ∈ R, the eigenvalues {λ k (t)} with ℑ(λ k ) > 0, and the proportionality constants {γ k (t)}. The reconstruction is particularly simple if S 12 (t, λ) ≡ 0 as a function of λ for some t, since it then follows from (8) that which extends to ℑ(λ) < 0 as a meromorphic function. Similarly one sees that S 22 (t, λ) = 1/S 11 (t, λ) and that S 21 (t, λ) ≡ 0. Since S(t, λ) is diagonal in this case, the solution matrices M ± (x, t, λ) can be expressed in terms of a common solution matrix U(x, t, λ) by setting The columns of U(x, t, λ) = (U 1 (x, t, λ), U 2 (x, t, λ)) necessarily satisfy the relations for all k = 1, . . . , N . It follows that U(x, t, λ) takes the simple form that is, a polynomial in λ times an exponential, where the matrix coefficients U (p) (x, t) are determined uniquely from (12). This means that (12) can be viewed as a linear algebraic system of 4N equations in 4N unknowns, the matrix elements of U (p) (x, t). Moreover, it can be shown that U constructed in this way is satisfies ∂ x U = LU if and only if the potential function in L is This formula reconstructs ψ(x, t) from the discrete scattering data {λ k (t)} and {γ k (t)} in the "reflectionless" case when S 12 (t, λ) ≡ 0. This treatment of multisoliton potentials via the matrix U follows Krichever [25], Manin [26], and Date [27]. See [28] for a relevant application.

Dynamics of the scattering data.
We now recall how the data evolve in t when ψ satisfies (1). The motivating observation [8] is that (1) can be cast in matrix form: where the matrix L is the one appearing in the linear scattering problem (4), and where Using the fact that M ± satisfies (4), multiply (15) on the right by M ± and find This equation is solved for (i∂ t − B)M ± by variation of parameters. Introducing a new unknown J ± (x, t, λ) defined through the relation (i∂ t − B)M ± = M ± J ± , one finds that J ± satisfies ∂ x J ± = −M ±−1 PM ± . We now integrate to find J ± explicitly, taking into account the boundary conditions satisfied by M ± as x → ±∞ and the fact that in both limits B → λ 2 σ 3 . With the use of these explicit formulae for J ± the equations (i∂ t − B)M ± = M ± J ± become equations of motion for the matrices M ± : As written, (18) does not make sense for ℑ(λ) = 0. But for λ ∈ R, the columns M − 1 and M + 2 are the boundary values of functions analytic for ℑ(λ) > 0, and we will also need equations for them that hold for ℑ(λ) > 0. To this end, we introduce the matrix M(x, t, λ) := (M − 1 , M + 2 ), and as before define the new unknown J(x, t, λ) = (J 1 , J 2 ) through the relation (i∂ t − B)M = MJ, and then integrate: As before, these expressions are used in (i∂ t − B)M = MJ to yield the equation of motion for M, valid for ℑ(λ) > 0 except at {λ k }, where M fails to be invertible. Each singularity is, however, removable, since det M = S 11 and hence (writing We make the natural assumption that the (isolated) zeros λ = λ k (t) of the denominator S 11 (t, λ) are simple [23]. But then the numerator of each entry is analytic at λ = λ k (t) and is easily seen to vanish there, thus cancelling the singularity. Hence, the evolution equation for M makes sense as λ → λ k (t). We accordingly introduce the notation The equations of motion for M ± and M determine the evolution of the scattering data. Using S := M +−1 M − , for real λ one finds Substituting from (18) yields Finally, since S does not depend on x, it may be brought inside the integrals. With the use of its definition the integrals are combined, giving the equation of motion: Note that since P is off-diagonal, the equation for S 11 (t, λ) only involves quantities analytic for ℑ(λ) > 0. Likewise, the equation for S 22 (t, λ) only involves quantities analytic for ℑ(λ) < 0. The equation of motion for the reflection coefficient S 12 (t, λ) is contained in that for S: The integrand here is p[ψ, ψ * ]M + 22 M − 22 − p[ψ, ψ * ] * M + 12 M − 12 , evaluated at x ′ , t, and λ, which generally only makes sense for λ ∈ R, as required. Now, the expression defining the zeros λ k (t) of S 11 (t, λ) is S 11 (t, λ k (t)) = 0. Differentiating with respect to t gives Using the equation of motion for S, one therefore finds The integrand here is p[ψ, , evaluated at x ′ , t, and λ = λ k (t). As remarked above, this makes sense with ℑ(λ k (t)) > 0. It remains to find an equation for {γ k (t)}. Differentiating the defining relation M + 2 (x, t, λ k (t)) = γ k (t)M − 1 (x, t, λ k (t)) with respect to t and using the evolution equation for M taken in the limit λ → λ k (t), yields the equation of motion with λ = λ k (t). The equations (25), (27), and (28) describe the evolution of the scattering data, but are coupled to the equations for M and M ± . This coupling disappears for P ≡ 0: for k = 1, . . . , N , as was first observed by Zakharov and Shabat [1]. From this simple system, it is possible to introduce the coupling perturbatively, leading to closed systems order by order.

Perturbation Theory with Nearly Bound Solitons
We now suppose that p[ψ, ψ * ] = ǫ 2 W (|ψ| 2 )ψ for some real-valued function W , taking ǫ > 0 to be a small parameter, and seek a perturbative solution of the equations of motion for the scattering data. We want a description of the solution up to an O(ǫ 2 ) error, containing important physical information, and valid uniformly over time scales of length O(ǫ −1 ). The initial data we consider is Proposition 1 The solution of the initial-value problem of (25), (27), and (28) with initial conditions (30), is given asymptotically for small ǫ by S 12 (t, λ) = O(ǫ 2 ) and where x k (T ), v k (T ), and θ k (t) are certain functions to be specified below. They satisfy x k (0) = x 0 k , v k (0) = 0 and θ k (0) = θ 0 k . This approximation is uniformly valid for times t = O(ǫ −1 ).
We develop the expansion using the multiscale formalism. Introducing the slow time variable T = ǫt, and assuming all quantities to depend functionally on both t and T , we replace the time derivatives in (25), (27), and (28) according to the chain rule: ∂ t → ∂ t + ǫ∂ T . Observe that for the initial conditions (30), there is no enforced magnitude for ℜ(λ k ) or S 12 (λ). We may thus select the scaling of these quantities to achieve a dominant balance. We choose to scale S 12 (λ) as ǫ 2 and ℜ(λ k ) as ǫ. Thus, setting λ k = ǫa k + ib k and γ k = exp(∆ k + iξ k ), we assume the expansions: Substituting into the equations of motion and collecting powers of ǫ, we examine the resulting equations order by order. First, from the leading-order terms in (27) we find for k = 1, . . . , N that so that these quantities do not depend on the fast time t. Similarly, looking at (28) we see that The description we desire will follow upon determining the T dependence of these leading-order quantities. The O(ǫ) contribution in the equation for b k , the imaginary part of (27), is At O(ǫ), equation (28) gives Again, avoid secular growth of ∆ and then take ∆ (1) k = 0. If we now define: then (38) takes the simple form: An equation for v k (T ) is found at O(ǫ 2 ) in the real part of (27). We find where G k (t, T ) is the leading term, divided by ǫ 2 , of the right hand side of (27). In more detail, from (8) and the leading-order behavior of {λ k }, we first see that To find the leading-order behavior of M ± , recall that S 12 = O(ǫ 2 ) so that we can use the "reflectionless" construction of M ± and ψ in terms of U, which in turn is constructed from Now, it is clear from (12) that all of the x and t dependence in U and ψ enters through the products where angled brackets denote averaging over t with T fixed. The force functions F k depend parametrically on the masses m k . Equations (40) and (45) imply Newton's equations for a system of interacting particles of mass m k and coordinate x k : It is easy to see that F k (x 1 + dx, x 2 + dx, . . . , x N + dx) = F k (x 1 , x 2 , . . . , x N ) so that the forces only depend on the relative coordinates. There is also a symmetry for (46) coming from the conservation of momentum that holds exactly (and thus to all orders of expansion) in (1) with p[ψ, ψ * ] = ǫ 2 W (|ψ| 2 )ψ. This symmetry follows from the trace formula (9) and shows that the total force on the system is zero: The dynamical system (46) describes the evolution of the scattering data. Since the reflection coefficient vanishes to second order on the time scales of interest, solutions of (46) can be used to build, at each fixed t, the N -soliton potential as in §2.2. This allows a direct comparison between numerics for (1) and the predictions of (46).

Two Particles
Consider the case N = 2. The aforementioned symmetries imply that the system takes the form for some function F . The relevant quantity is then the relative distance y(T ) := x 2 (T ) − x 1 (T ), which has the simple-looking equation of motioñ where the effective mass is defined bym :

Writing down the force function.
We begin our study of the force functions by simplifying the integrand in (43) to isolate terms that are exact x-derivatives and do not contribute. In this context, consider the Then, these quantities again satisfy a linear system of equations Using the quadratic forms associated with U, f as defined by (44) is seen to be an exact x-derivative: where the last equality follows from the fact that the determinant of any solution of (4) is independent of x because L is traceless. For N = 2, we use the relations (12) and the parameters λ 1 = im 1 , λ 2 = im 2 , γ 1 = exp(−2m 1 x 1 (T ) + i(θ 0 1 − 2m 2 1 t)), and γ 2 = exp(−2m 2 x 2 (T ) + i(θ 0 2 − 2m 2 2 t)) to find U 12 = e iλx (λψ/(2i) + ϕ) and where ψ is the well-known two-soliton "breather" solution, and using Since W (·) ∈ R, only ℜ(f (x, t, im k )) is needed to find ℑ(G k (t, T )). From (52) one finds f = −∂ x (λ 2 |ψ| 2 /4 + λℑ(ψϕ * ) + |ϕ| 2 ), and therefore ℜ(f (x, t, im k )) = m 2 k ∂ x |ψ| 2 /4 − ∂ x |ϕ| 2 . Using this in the formula (43) for ℑ(G k (t, T )), one finds that the first term is an exact derivative of a rapidly decreasing function and hence integrates away. In terms of the two quantities |ψ| 2 and |ϕ| 2 obtained directly from (53) we finally obtain In particular, it follows that −2m 1 ℑ(G 1 (t, T ))−2m 2 ℑ(G 2 (t, T )) = 0 so that the total instantaneous (that is, before averaging over t) force vanishes. Specializing further to the quintic perturbation (2) by taking W (ρ) := σρ 2 and writing we have found the following explicit formula for h: where the individual terms h k are given in the Appendix. They depend on a dummy integration variable z that differs from x by a simple translation. Note that, by the periodicity with respect to the "fast" function θ = θ 2 (t) − θ 1 (t), averaging over t is equivalent to averaging over θ.

Scale invariance.
From (56) and the explicit formulae for the terms h k in the Appendix, note the important symmetry: for all nonzero ξ ∈ R. Setting y = ξq and S = ξ −3 T , the equation of motion takes the form: For arbitrary masses m 1 and m 2 , we may then set ξ = (m 1 m 2 ) −1/2 . Becausem is homogeneous of degree one in m 1 and m 2 , it is convenient to use the normalized masses Here, M 1 and M 2 satisfy M 1 M 2 = 1 and may therefore be expressed in terms of the normalized effective massM by solving 2(M −1 2 + M −1 1 ) −1 =M subject to this constraint to find: assuming without loss of generality that M 2 > M 1 . From now on, we will work exclusively with the normalized masses, in which case the force depends only on S andM .

General features of the force function.
Unfortunately, (65) cannot be evaluated in closed form because the integrand generally involves both exp(ζ 1 ) and exp(ζ 2 ). Even if M 2 /M 1 ∈ Q so that the integrand becomes a rational function of, say, exp(ζ 1 ), the denominator is irreducible already for the simplest resonance, M 2 = 2M 1 . In spite of these difficulties, certain elementary features of the force law can be extracted: • F (q;M ) is proportional to the constant σ = ±1, as is clear from (65).
• F (q;M ) is an odd function of q, since the integrand satisfies H(−q, z) = −H(q, −z) and moreover this symmetry holds term by term in the formula for H.
• F (q;M ) decays to zero for large q. This follows from the fact that the denominator of each term H mn in the integral is bounded and the corresponding numerator vanishes for large q whenever g mn = 0. The result then follows from a dominated-convergence argument.
• F (q;M ) only vanishes exactly for q = 0. Thus it is strictly of one sign for q > 0.
• The normalized effective massM enters the dynamics both as a mass parameter multiplying the acceleration q ′′ (S) and as a parameter in F (q;M ) itself.
The force F (q;M ), as computed from the integral formula (65), is plotted in Figure 1 for several different values of the normalized effective massM . A formula for the spring constant k(M ) can be found by simply differentiating with respect to q in (65) and setting q = 0, however it seems less useful to present than a plot, shown in Figure 3, of the (numerically) evaluated formula. In Figure 4 we plot the corresponding frequency (on the time scale S), the latter being a directly observable quantity. It is noteworthy here that the dynamics of solitons can be described by a linear theory even though their amplitudes are not at all small. The parameter linearizing the theory is the distance between the solitons, rather than the soliton amplitude. We also remark that the limit in which this linear behavior holds is that of infinitessimally-separated solitons, a limit in which methods assuming the solitons to be well-separated are invalid.

Repulsive case. Asymptotic velocity.
For σ = −1, the force and displacement q have the same sign, resulting in q always becoming large.
Solitons that are near each other at t = 0 are ejected from the origin as observed by Artigas et. al. [21]. This effect is captured accurately by our theory, as shown in Figure 5. The work done by the force in moving the particle from q = q 0 to q = ∞ determines the asymptotic velocity of an initially stationary particle upon ejection. Taking q 0 = 0 corresponds to the ultimate velocity of a  Figure 2: A two-soliton bound state affected by an attractive perturbation. Here, ǫ = 0.0387, m 1 = 0.6, and m 2 = 1. Left: an approximation to |ψ| 2 found by solving Newton's equations and then constructing the field using reflectionless inverse theory. Right: the corresponding numerics for (1). The bound state has too much energy for the harmonic approximation to hold, and the period of motion, about 120 time units, is longer than the harmonic period.
stationary particle that is slightly perturbed from (unstable) equilibrium at the origin. With zero initial velocity, one equates the asymptotic kinetic energy with the work done: to find a formula for the asymptotic velocity difference: (68) Figure 6 shows the asymptotic velocity difference q ′ (∞) for q 0 = 0 found from (68) as a function of the normalized effective massM . To apply the graph in Figure 6 to problems with unnormalized masses, it is useful to unravel the changes of variables made so far. Given m 1 and m 2 , the scaling parameter is ξ = (m 1 m 2 ) −1/2 and the effective mass ism = 2·(m −1 1 +m −1 2 ) −1 . Then, the normalized effective mass used in Figure 6 isM = ξm. Next, from the graph one finds the asymptotic velocity q ′ (∞). The true velocity in the original coordinates is then dy/dt = ǫξ −2 q ′ (∞). For example, the parameters used in Figure 5 imply a normalized effective mass ofM ≈ 0.9. From Figure 6 one finds q ′ (∞) ≈ 5.0, and thus dy/dt ≈ 0.15. This value agrees well with the pictures in Figure 5.
In the attractive case, the integral (67) also has physical meaning as the binding energy of the two-soliton state. A relative velocity in excess of q ′ (∞), the escape velocity, will "ionize" the state.

Discussion
Multiscale asymptotics shows that under certain conditions the behavior of a multisoliton initial condition in a perturbed NLSE reduces to Newton's equations for a system of interacting particles, one particle per soliton. The theory applies over time scales of length O(ǫ −1 ) for perturbations of size ǫ 2 , when the initial velocities of the solitons mutually differ by an O(ǫ) amount. Our calculations make very concrete the often-cited analogy between solitons and particles. We want to emphasize that the limit considered here is one in which the relative velocities of the solitons are small but the solitons may be strongly nonlinearly superimposed, precisely the limit in which methods exploiting large distances between solitons fail.
For a quintic perturbation of the NLSE and an initial condition composed of two solitons, the resulting dynamical system can be analyzed. When the perturbation is attractive (σ = +1), the system describes a nonlinear oscillator with all solutions q(S) being periodic. If the energy associated with q(S) is small (that is, if q(0) and q ′ (0) are both small), then the periodic motion is nearly harmonic, and formulae for the associated spring constant and frequency of motion can be found; in this limit the model for the soliton interaction linearizes even though the soliton amplitudes are not at all small. The latter are determined by the masses m 1 and m 2 and are not related to the coordinate q(S). For larger energies, the spring "softens" and the frequency decreases with increasing energy. The pictures in Figure 2 show oscillations in the nonlinear regime, where the frequency of motion is smaller than the linear frequency. Of course even in the nonlinear regime, the dynamics still obey the simple modelM q ′′ = F (q;M ). Although the periodic motion is predicted and observed over long time scales of size O(ǫ −1 ), it is not likely to persist for all time, due to the influence of higher-order resonant coupling effects.
On the other hand, when the quintic perturbation is repulsive (σ = −1), the nodal point at the origin in the phase plane gets replaced with an unstable saddle point. All orbits apart from the fixed point itself represent the nonlinear development of the instability. Because the force vanishes fast enough for large q, the velocity q ′ (S) ultimately saturates as the two-soliton state becomes "ionized". From the force function F (q;M ) this "ejection" velocity may be calculated, giving excellent agreement with direct simulations of the perturbed NLSE. This analysis explains the observations reported in [21]. The symmetry-breaking that determines which soliton ends up on the right and which on the left can be traced to the location of the initial phase point in relation to the separatrix connected to the saddle. Unlike in the attractive (σ = +1) case, the approximation obtained from multiscale asymptotics for the repulsive (σ = −1) case is expected to be uniformly valid for all time, since as the solitons separate, further effects due to resonant coupling diminish.
Given the formula for the force F (q;M ), it is possible to compute the harmonic frequency and ejection velocity, more explicitly than we have done here. For example, the formulae would be expected to simplify in the limitsM ↓ 0 (corresponding to two solitons differing very much in amplitude) andM ↑ 1 (corresponding to two solitons with nearly the same amplitude). The calculation of the ejection velocity is challenging because it may require uniform approximation of F (q;M ) for all q in the limit of interest; pointwise asymptotics for fixed q are not enough to approximate the work integral without further information.

Acknowledgements
PDM acknowledges the support of NSF grant DMS 9304580 while a member of the School of Mathematics at the Institute for Advanced Study. During the preparation of this paper, JAB and NNA were affiliated with the Australian Photonics Cooperative Research Centre.

Appendix: Formulae for the Two-Particle Force Function Integrand
Here, we record the details of the formulae for the two-particle force function needed to calculate or approximate for special values ofM the force and related quantities to any desired accuracy.