You are an expert at summarizing long articles. Proceed to summarize the following text: 
the old - new interest in the problem of soliton - soliton intertaction and soliton molecules has been increasingly accumulating particularly over the past few years . this is mainly motivated by the application of optical solitons as data carriers in optical fibers @xcite and the realization of matter - wave solitons in bose - einstein condensates @xcite . one major problem limiting the high - bit rate data transfer in optical fibers is the soliton - soliton interaction . on the one hand , soliton - soliton interaction is considered as a problem since it may destroy information coded by solitons sequences . on the other hand , it is part of the problem s solution , since the interaction between solitons leads to the formation of stable soliton molecules which can be used as data carriers with larger `` alphabet '' @xcite . the interaction force between solitons was first studied by karpman and solovev using perturbation analysis @xcite , gordon who used the exact two solitons solution @xcite , and anderson and lisak who employed a variational approach @xcite . it was shown that the force of interaction decays exponentially with the separation between the solitons and depends on the phase difference between them such that in - phase solitons attract and out - of - phase solitons repel . this feature was demonstrated experimentally in matter - wave solitons of attractive bose - einstein condensates @xcite where a variational approach accounted for this repulsion and showed that , in spite of the attractive interatomic interaction , the phase difference between neighboring solitons indeed causes their repulsion @xcite .    for shorter separations between the solitons , malomed @xcite used a perturbation approach to show that stationary solutions in the form of bound states of two solitons are possible . however , detailed numerical analysis showed that such bound states are unstable @xcite . stable bound states were then discovered by akhmediev @xcite and a mechanism of creating robust three - dimensional soliton molecules was suggested by crasovan _ recently , soliton molecules were realized experimentally by stratmann _ _ in dispersion - managed optical fibers @xcite and their phase structure was also measured @xcite . perurbative analysis was used to account theoretically for the binding mechanism and the molecule s main features @xcite . quantization of the binding energy was also predicted numerically by komarov _ _ @xcite . in refs.@xcite , a hamiltonian is constructed to describe the interaction dynamics of solitons . the mechanism by which the relative phase between the solitons leads to their force of interaction , and hence the binding mechanism , is understood only qualitatively as follows . for in - phase ( out - of - phase ) solitons , constructive ( destructive ) interference takes place in the overlap region resulting in enhancement ( reduction ) in the intensity . as a result , the attractive intensity - dependent nonlinear interaction causes the solitons to attract ( repel ) @xcite . a more quantitative description is given in refs . @xcite .    in view of its above - mentioned importance from the applications and fundamental physics point of views , we address here the problems of soliton - soliton interaction and soliton molecule formation using the exact two solitons solution . this approach has been long pioneered by gordon @xcite where he used the exact two solitons solution of the homogeneous nonlinear schr@xmath0dinger equation to derive a formula for the force of interaction between two solitons , namely @xmath1 where @xmath2 is the solitons separation and @xmath3 is their phase difference . this formula was derived in the limit of large solitons separation and for small difference in the center - of - mass speeds and intensities , which limits its validity to slow collisions . with appropriately constructed hamiltonian , wu _ et al . _ have derived , essentially , a similar formula that gives the force between two identical solitons and reliefs the condition on slow collisions @xcite .    here , we present a more comprehensive treatment where we derive the force between two solitons for arbitrary solitons intensities , center - of - mass speeds , and separation . we also generalize gordon s formula to inhomogeneous cases corresponding to matter - wave bright solitons in attractive bose - einstein condensates with time - dependent parabolic potentials @xcite and to optical solitons in graded - index waveguide amplifiers @xcite . many interesting situations can thus be investigated . this includes the various soliton - soliton collision regimes with arbitrary relative speeds , intensities , and phases . most importantly , soliton - soliton interaction at short solitons separations will now be accounted for more quantitatively than before . specifically , soliton molecule formation is clearly shown to arise from the time - dependence of the relative phase which plays the role of the restoring force . in this case , the force between the two solitons is shown to be composed of a part oscillating between attractive and repulsive , which arises from the relative phase , and an attractive part that arises from the nonlinear interaction . the time - dependence of the relative phase results in a natural oscillation of the molecule s bond length around an equilibrium value . the various features of the soliton molecule , including its equilibrium _ bond length _ , _ spring constant _ , _ frequency _ and _ amplitude _ of oscillation , and _ effective mass _ , will be derived in terms of the fundamental parameters of the solitons , namely their intensities and the nonlinear interaction strength . the two solitons solution is derived here using the inverse scattering method @xcite . although the two solitons solution of the homogeneous nonlinear schr@xmath0dinger equation is readily known @xcite , here we not only generalize this solution to inhomogeneous cases , but also present it in a new form that facilitates its analysis . the solution will be given in terms of the four fundamental parameters of each soliton , namely the initial amplitude , center - of - mass position and speed , and phase . the main features of the solution will be shown clearly such as the contribution of the nonlinear interaction to the actual separation and phase difference between solitons where it turns that the separation between the two solitons increases with logarithm of the difference between the amplitudes of the two solitons . furthermore , the general statement that _ a state of two equal solitons with zero relative speed and finite separation does not exist as a stationary state for the homogeneous nonlinear schr@xmath0dinger equation _ , will be transparently and rigorously proved . stability of soliton molecules is an important issue since , in real systems , perturbations caused by various sources such as losses , raman scattering , higher order dispersion , and scattering from local impurities , tend to destroy the molecules . to investigate the stability of the soliton molecules described by our formalism , we have considered three situations . first , we studied the reflection of the molecule from a hard wall and a softer one . while for the hard wall the molecule preserves its molecular structure after reflection , it generally breaks up for the softer ones due to energy losses at the interface . secondly , the scattering of the molecule by a potential barrier was also investigated . we show that the molecular structure is maintained only for some specific heights of the barrier . this suggests a quantization in the binding energy as predicted by komarov _ the oscillation period of the reflected molecule is noticed to be smaller than for the incident one . in addition , the outcome of scattering depends on the phase of the molecule s oscillation at the interface of the barrier . for instance , a dramatic change in the scattering outcome takes place if the _ coalescence _ point of the molecule lies exactly at the interface . in such a case , the otherwise totally transmitting molecule will now split into reflecting and tunneling solitons . thirdly , we have considered the collision between a single soliton with a stationary soliton molecule . the effects of different initial speeds , amplitudes , and phases of the scatterer soliton were studied . it turns out that for slower collisions , it is easier for the scatterer soliton to break up the soliton molecule , while for fast collisions the scatterer soliton expels and then replaces one of the solitons in the molecule . the phase of the scatterer soliton plays also a crucial rule in preserving or breaking the bond of the molecule , which can be used as _ key _ tool to code or uncode data in the molecule . the rest of the paper is organized as follows . in section  [ exact_sec ] , we use the inverse - scattering method to derive the two solitons solution of the inhomogeneous nonlinear schr@xmath0dinger equation and present the solution in the above - mentioned appealing form . the main features of the solution will be discussed in subsection  [ mainsec ] . the center - of - mass positions and relative phases will be derived in subsections  [ com_sec ] and [ phases_sec ] , respectively . the force between solitons will be derived in section  [ force_sec ] where gordon s formula will be extracted as a special case in subsection  [ gordon_sec ] and our more general formula will be derived in subsection  [ our_sec ] . in subsection [ num_sec ] , we compare our formula with the numerical calculation . in section  [ molsec ] , we show the possibility of forming soliton molecules , derive their main features in subsection  [ formation_sec ] , and investigate their stability in subsection  [ stability_sec ] . we end in section  [ conc_sec ] with a summary of results and conclusions . the details of the derivation of the two solitons solution and the center - of - mass positions are relegated to appendices a and b , respectively . matter - wave solitons of trapped bose - einstein condensates and optical solitons in optical fibers can be both described by the dimensionless gross - pitaevskii equation @xmath4 , time to @xmath5 , and the wavefunction @xmath6 to @xmath7 , where @xmath8 and @xmath9 are the characteristic frequencies of the quasi one - dimensional ( @xmath10 ) trapping potential in the axial and radial directions , respectively . in these units , the strength of the interatomic interaction will be given by the ratio @xmath11 , where @xmath12 is the @xmath13-wave scattering length . for the case of optical solitons , the function @xmath6 represents the beam envelope , @xmath14 is the propagation distance , @xmath15 is the radial direction , and the intensity - dependent term represents the kerr nonlinearity . in this case , scaling is in terms of the characteristic parameters of the fiber , as for instance in ref . the specific form of the prefactors of the inhomogeneous and nonlinear terms guarantees the integrability of this equation @xcite . for the special case of @xmath16 , the homogeneous case is retrieved . other interesting special cases have also been considered @xcite .    as outlined in appendix  [ appa ] , we use the darboux transformation method to derive the two solitons solution of this gross - pitaevskii equation , which can be put in the form @xmath17\nonumber\\ & + & \sqrt{\frac{n_2 \alpha_{22}(t)}{2}}e^{i ( \phi_{01}+\phi_{02}+\phi_2(x , t)+ \tan^{-1}({\alpha_2\over\alpha_1 } ) + \tan^{-1}({\alpha_4\over\alpha_3}))}\nonumber\\ & \times&{\rm{sech}}\left[\alpha_{22}(t)(x - x_{\rm cm2}(t ) ) + \frac{1}{2}\log\left(\frac{\alpha_1 ^ 2+\alpha_2 ^ 2}{\alpha_3 ^ 2 + \alpha_4 ^ 2}\right)\right]\label{exactsol},\end{aligned}\ ] ] where +   + @xmath18 , @xmath19 , +   + @xmath20 , +   + @xmath21 , +   + @xmath22 , +   + @xmath23 , +   + @xmath24 , +   + @xmath25 , +   + @xmath26 , +   + @xmath27 , +   + @xmath28 , +   + @xmath29 , +   + @xmath30 , +   + @xmath31 , +   + @xmath32 , +   + @xmath33 , +   + @xmath34 , +   + @xmath35 , +   + @xmath36 . +   + the solution is put in this suggestive form to facilitate its analysis . the first sech part corresponds to the exact single soliton solution with center - of - mass position @xmath37 , width @xmath38 , phase @xmath39 , and normalization @xmath40 . hence , @xmath41 and @xmath42 correspond to the initial center - of - mass position and speed , respectively . the second sech term contains the same features in addition to a shift in both the center - of - mass position and phase . it should be noted , however , that @xmath43 , @xmath44 , @xmath45 , and @xmath46 correspond to the center - of - mass position and speed , normalization , and phase of the single noninteracting solitons . due to the interaction between solitons , these four characteristic quantities may not correspond exactly to the values of the same physical quantities as they did for the single soliton solution . for instance , @xmath37 will not correspond to the center - of - mass of one of the solitons . instead , the soliton may be shifted from that position due to the interaction with the other soliton . in the following , we present a detailed analysis of the locations and phases of the two solitons . inspection shows that there are two main regimes for the two solitons solution , namely the regime of resolved solitons and the regime of overlapping solitons . in the former case , the center - of - mass concept is well defined and analysis of the relative dynamics becomes feasible . the solitons are considered resolved as long as the two main peaks are not overlapping , which means that partial overlap may occur in this regime . the analysis in this section assumes the resolved solitons regime .    in the resolved solitons regime , the argument of the second sech term of eq . ( [ exactsol ] ) , namely @xmath47 $ ] , simplifies to a function with three roots . the fact that the sech function is peaked at the roots of its argument , leads to that the second sech term corresponds to three `` solitons '' . this is shown in fig  [ fig1 ] . we denote these solitons as the `` left '' , `` central '' , and `` right '' solitons with peak locations at @xmath48 , @xmath49 , and @xmath50 , respectively . we notice that the central soliton is located at the position of the soliton of the first sech term , namely near @xmath37 . further inspection shows that the two solitons at this location are out of phase and interfere destructively such that they do not appear in the total profile . therefore , the two solitons that our solution of eq . ( [ exactsol ] ) describes are in fact the left and right solitons arising from the second sech term . this is different from what one aught to conclude from the form of the exact solution , namely that the first sech term corresponds to one soliton and the second sech term corresponds to the other soliton . in general , the center - of - mass locations of the left and right solitons , @xmath48 and @xmath50 , do not match @xmath37 and @xmath51 . typically , @xmath48 will be shifted to the left of @xmath37 , @xmath50 will be shifted to the right of @xmath51 , while @xmath52 remains near @xmath37 . the amount of shift depends mainly on the normalization difference @xmath53 and the relative speed @xmath54 , as will be shown in the next subsection . an interesting general and exact result is that , for the homogeneous case @xmath16 , a state of two equal solitons , @xmath55 , with zero relative speed , @xmath56 , and finite separation , does not exist as an exact solution of the gross - pitaevskii equation . this can be proven by substituting @xmath55 and @xmath56 in @xmath57 to find that the right and left solitons migrate to @xmath58 and @xmath59 , respectively , while the center - of - mass of the central soliton matches exactly @xmath37 . furthermore , we show in section [ phases_sec ] that the phase difference between this central soliton and that of the first sech term equals , in this case , @xmath60 ; guaranteeing their destructive interference . thus , the three solitons disappear in such a special case and @xmath6 becomes the trivial solution .    in view of the above , a need arises to derive formulae for the center - of - mass positions , @xmath48 and @xmath50 in terms of the solitons parameters , which will then be used to derive the force between the two solitons .      in this section , we derive formulae for the three roots of @xmath57 which correspond to the locations of the left , central , and right solitons , @xmath48 , @xmath52 , and @xmath50 , respectively . to facilitate the derivation , we define @xmath61}$ ] , @xmath62}$ ] , @xmath63 , @xmath64 , @xmath65 , and @xmath66 . the equation @xmath67 is a third - order polynomial in both @xmath68 and @xmath69 . in principle , this equation can be solved algebraically for @xmath68 or @xmath69 . however , extracting @xmath15 from the resulting three roots will not be possible analytically for general @xmath70 . alternatively , and as can be seen from fig . [ fig1 ] , we can exploit the simple linear behavior of @xmath57 near its roots . assuming , without loss of generality , @xmath71 , and noting that in the resolved solitons regime the solitons separation @xmath72 is large , we argue in appendix  [ appb1 ] that @xmath73 for all @xmath15 , @xmath74 for @xmath75 , @xmath76 for @xmath77 , and @xmath78 for @xmath79 . based on this , the center - of - mass position @xmath50 can be derived from a taylor expansion of @xmath57 in powers of large @xmath68 and @xmath69 . expanding @xmath57 in powers of large @xmath68 only and leaving @xmath69 arbitrary , accounts for @xmath52 and @xmath48 simultaneously . keeping terms up to first order of @xmath80 , we show in appendix  [ appb1 ] that the position of the right soliton is given by @xmath81 where @xmath82 and the position of the the left soliton is given by @xmath83 where @xmath84 , @xmath85 , and @xmath86 and @xmath87 are given in appendix [ appb1 ] . the second and third terms on the right hand side of eqs . ( [ xreq2 ] ) and ( [ xleq2 ] ) account for the shift in the center - of - mass position with respect to the single soliton ones . the third terms are much smaller than the second ones since they decay exponentially with the solitons distance @xmath88 . in the limit @xmath89 and @xmath90 , both @xmath91 and @xmath92 take the form @xmath93}$ ] . thus , it is obvious that for @xmath94 , @xmath95 and @xmath96 . this agrees with our earlier result that two equal solitons with zero relative speed and finite separation , do not exist as an exact solution .      in this section , we calculate the phases of the left , central , and right solitons in reference to the phase of the soliton of the first sech part in the two solitons solution . for simplicity , the special case of @xmath97 and @xmath98 will be assumed . the phase difference between the left , central , and right solitons on one hand and the soliton of the first sech term on the other hand is generally given by @xmath99 which by observing that for all @xmath15 @xmath100 reduces to @xmath101 this expression gives the phases of the right soliton @xmath102 , the central soliton @xmath103 , and the left soliton @xmath104 , for @xmath75 , @xmath52 , and @xmath48 , respectively . to calculate these phases we express the parameters @xmath105 in terms of @xmath68 and @xmath69 , as follows @xmath106 @xmath107 @xmath108 @xmath109 in general , @xmath73 for all @xmath15 , but @xmath74 only for @xmath110 , and @xmath111 for @xmath112 , as shown in appendix  [ appb1 ] . +   + * i. phase of the right soliton @xmath102 : * +   + in this case @xmath113 which leads to @xmath114 and thus @xmath115 , @xmath116 , @xmath117 , @xmath118 . therefore +   +   + @xmath119 +   +   + and @xmath120 , which gives @xmath121 where @xmath122 is a function that depends on the ratio @xmath123 for nonzero @xmath124 and @xmath125 . +   + * ii . phase of the central soliton @xmath103 : * +   + in this case @xmath126 which gives @xmath127 , and hence @xmath128 , and @xmath129 . finally , we get @xmath130 the last result shows that , for @xmath131 , i.e. , two equal solitons with zero relative speed , the central soliton and the soliton of the first sech term of the two solitons solution are out of phase and therefore interfere destructively . +   +   + * iii . phase of the left soliton @xmath104 : * +   + in this case @xmath132 which results in @xmath78 and @xmath133 , @xmath134 , @xmath135 , @xmath136 . therefore , @xmath128 and @xmath137 , which gives @xmath138 the phase difference between the left and right solitons is thus given by @xmath139 noting that @xmath140 for @xmath141 and small but finite @xmath125 , and @xmath142 for @xmath143 and small but finite @xmath124 , we finally conclude that the phase difference between the two solitons is given by @xmath144 the above results are verified in fig . [ fig2 ] , where we plot @xmath145 , @xmath102 , @xmath104 , and @xmath103 versus @xmath15 . the agreement between our estimated values and the exact curve is evident .      in this section , we use the results of the previous two subsections to derive the force between the left and right solitons . the force is proportional to the acceleration of the solitons separation @xmath146\,e^{-\frac{1}{4 } g_0 { n_s } { x_d}(t ) e^{\gamma ( t)}}\label{deltaeq},\end{aligned}\ ] ] where @xmath147 and the coefficients @xmath148 are given in appendix  [ appb2 ] . the first two terms on the right hand side of eq . ( [ deltaeq ] ) are the dominant ones since they correspond to the noninteracting solitons separation , @xmath149 , and their logarithmic shifts @xmath150 arising from the interaction between solitons . the time - dependence of @xmath151 originates from @xmath152 , @xmath153 , and @xmath149 . the acceleration can thus be derived @xmath154     \ , e^{-\frac{1}{4 } g_0 n_s { \delta}(t ) e^{\gamma ( t)}+3 \gamma ( t)}\label{force2},\ ] ] where the coefficients @xmath155 are given in appendix  [ appb2 ] . in the last equation , we have used eq . ( [ deltaeq ] ) with only the first two terms of its right hand side to substitute for @xmath149 in terms of @xmath2 in the exponential factor . the first term on the right hand side of eq . ( [ force2 ] ) corresponds to the force due to the external potential which vanishes for the homogeneous case . the rest of the terms correspond to the force of interaction between the solitons . the interaction force depends , as expected , on the phase difference of the two solitons and decays exponentially with their separation . it should be noted that this equation is a generalization of gordon s formula @xcite in two aspects . first , it is derived for a time - dependent inhomogeneous medium . secondly we have , essentially , no restriction on the difference between the two solitons amplitudes and speeds ; apart from some extreme cases which were mentioned in appendix [ appb1 ] and will be discussed further below .      for the homogeneous case , @xmath16 , and in the limits @xmath89 and @xmath90 , the acceleration formula , eq .  ( [ force2 ] ) , simplifies considerably . an apparent inconsistency occurs when switching the order of these two limits , namely @xmath156 while @xmath157 which differ by an overall minus sign . the conflict is resolved by invoking eq . ( [ phasediff ] ) where it is shown that in the first case : @xmath158 , while in the second case : @xmath159 . therefore , the two approaches agree on the following result @xmath160 this is essentially gordon s result , eq . ( [ gordoneq ] ) , since in his derivation gordon took @xmath161 and soliton amplitude @xmath162 , @xmath163 , as can be seen in eqs.(1 ) and ( 6 ) of ref . @xcite . substituting @xmath161 and @xmath164 in the last equation , it becomes identical to eq . ( [ gordoneq ] ) . it should be mentioned here that eq . ( [ gordoneq ] ) was also derived in ref . @xcite using a perturbation analysis based on the inverse scattering method , and in ref . @xcite using a variational calculation .      for nonzero @xmath152 and in the limits @xmath89 and @xmath90 , the acceleration formula takes the form @xmath165 this is a generalization to gordon s formula for the inhomogeneous case as modeled by eq . ( [ gp ] ) . depending on the specific form of @xmath152 , the two force terms , namely the external ( first term ) and the interaction ( second term ) , can be repulsive , attractive , or oscillatory . in addition , the phase difference @xmath3 also depends on @xmath152 . it is established in the homogeneous case , as will also be shown in section [ molsec ] , that the time - dependence of the phase difference is responsible for binding the two solitons in the soliton molecule . here , in addition to the possibility of forming soliton molecules , the dependence of @xmath3 on @xmath152 allows for controlling the parameters of the molecule such as its equilibrium bond length , period , and spring constant . this and possibly other interesting phenomena will be left for future investigation .      to obtain an estimate of the accuracy of the general acceleration formula , eq . ( [ force2 ] ) , we calculate numerically the acceleration from the exact two solitons solution , eq . ( [ exactsol ] ) , and compare the two results . the distance between the two solitons of the function @xmath166 is determined using a numerical algorithm that employs our formulae for @xmath48 and @xmath50 given by eqs . ( [ xreq2 ] ) and ( [ xleq2 ] ) to calculate seed values . the distance is then differentiated numerically twice at @xmath167 . in fig . [ fig3 ] , we compare the two results . good agreement is obtained for @xmath168 . the analytical solution diverges at @xmath169 . the value at which divergence takes place is set by the specific choice of parameters in fig . as pointed out in section [ com_sec ] , the divergence occurs due to the merging of the central and left solitons . this artifact divergency can be remedied by associating the location of the local maximum of @xmath57 to @xmath48 once this maximum has reached the @xmath15-axis from above . restricting our study to the region where agreement is obtained , we interestingly notice that the acceleration is oscillating between positive and negative values . this means that the force between the solitons is oscillating between attractive and repulsive . the possibility of having attractive forces for finite @xmath124 is particularly interesting ; for two solitons with nonzero relative positive speed , i.e. the solitons are initially diverging from each other , the force between them is attractive . this suggests that , if the force remains attractive for sufficient time , the two solitons will slow down and eventually converge at some point . if true , this should occur at small distances since the force decays exponentially with distance and when the two solitons are allowed to diverge even for a short while , the force might be weakened such that the two solitons can not return back . to be able to judge on such a possibility , we need to know what happens to the acceleration at later times . to that end , we calculate numerically the acceleration in terms of @xmath124 and @xmath14 . the result is plotted in fig . [ fig4 ] , where it is clear that the acceleration indeed decays with time for all @xmath124 . this leads to that any nontrivial effect of the oscillating force is most likely to take place at short solitons separations . this is what we find in the next section where the possibility of forming stable soliton molecules is pointed out . we have shown in the previous section that , as a result of the solitons time - dependent relative phase , the force of interaction between solitons is oscillating between repulsive and attractive . since the force decays exponentially with the solitons separations , this oscillation will have a tangible effect only when the two solitons are close to each other . in this section , we investigate the force of interaction between solitons for short solitons separation . in such a special case , eq .  ( [ deltaeq ] ) takes a simple form that accounts for the solitons separation in terms of their relative phase . using this formula , we show that the solitons will be bound to oscillate around some equilibrium distance where the phase plays the role of the restoring force . comparison with exact numerical calculations shows that this formula is accurate for almost the full range of the solitons separation , except at the coalescence point ( if any ) . in subsection  [ formation_sec ] , we discuss the main features of the resulting soliton molecules , and in subsection  [ stability_sec ] , we investigate numerically their stability in different scattering regimes .    to focus on the role of relative phase , we simplify the analysis by restricting our treatment to the homogeneous case , @xmath16 , and zero relative speed , @xmath170 . we also set @xmath171 so that any separation between the solitons to be as small as possible which in this case arises only from the logarithmic shifts ( @xmath172 in eq . ( [ deltaeq ] ) ) . in this case , the solitons separation @xmath151 , given by eq . ( [ deltaeq ] ) , simplifies in the limit @xmath173 to @xmath174-\frac{4}{g_0 n_s}\log \left(\frac{g_0     n_d^2}{n_s}\right)\label{delmol},\ ] ] and the acceleration is given by @xmath175 it is noted that for @xmath176 or @xmath177 , gordon s formula is retrieved , but here with an explicit time - dependence of the phase , @xmath178 . this acceleration formula deviates considerably from gordon s formula for @xmath179 . specifically , for @xmath180 , @xmath151 diverges to @xmath59 at @xmath181 , which indicates that the two solitons coalesce . this is confirmed below by examining the exact solution at this condition . we note here that an approximate expression for the solitons separation was also derived in refs . in addition , our predicted molecule s oscillation frequency ( see eq . ( [ freqeq ] ) below ) agrees with these references .    to verify this feature , we calculate numerically the distance between the two solitons directly from the exact solution , eq . ( [ exactsol ] ) , for different values of @xmath182 . for @xmath183 , the density plot in fig . [ fig5]a , shows a soliton molecule of two clearly resolved solitons with a separation oscillating around some nonzero equilibrium distance . approaching the solitons coalescence point with @xmath184 , the density plot in fig . [ fig6]a , shows the two solitons approaching each other more than the previous case . furthermore , this figure shows a slight bounce back by one of the solitons in the region of collision . approaching further the coalescence condition with @xmath185 , we indeed observe in fig . [ fig7]a that the two solitons merge almost completely . for a more quantitative comparison , we calculate numerically the center - of - mass trajectories of the two solitons . we show the trajectory curves in the density plots of figs . [ fig5]a-[fig7]a . in figs . [ fig5]b-[fig7]b , we plot the solitons separation obtained from formula ( [ delmol ] ) and the numerical trajectories obtained from the exact solution . it is clear from these figures that this formula agrees well with the exact soliton separation except near the collision region . in fig . [ fig5]b , the two solitons remain away from each other during the collision , and therefore good agreement is obtained with the exact result even in the collision region . in fig . [ fig6]b the two solitons approach each other further such that formula ( [ delmol ] ) does not account for the above - mentioned slight bounce of one of the solitons . in fig . [ fig7]b , agreement with the exact solution in the collision region is qualitative . we found that at the condition @xmath180 and for @xmath173 , the analytic curve overlaps with that of the exact solution ; apart from the horizontal segments where formula ( [ delmol ] ) diverges to @xmath59 . further insight is obtained by plotting the density profile of the soliton molecule at some specific times , as shown in figs . [ fig5]c-[fig7]c . in fig . [ fig5]c , we observe that the initial amplitude imbalance is never removed during the dynamics . instead , it becomes maximum when the two solitons are closest to each other . in addition , we notice that the oscillation amplitude of the larger soliton around its equilibrium position is larger . figure  [ fig6]c shows clearly the soliton bounce which takes place in the time interval @xmath186 to @xmath187 . in these subfigures we plot two vertical dashed lines that indicate the position of the solitons at the closest approach . it is clear that after first closest approach at @xmath186 , the right soliton bounces back with a maximum displacement at @xmath188 . in fig . [ fig7]c , it is shown that , although the two solitons coalesce , two small symmetric @xmath189 appear . a detailed examination of these wings shows that they are the remnants of the two solitons after they coalesce and they both bounce back in the collision region similar to the case of fig . [ fig6 ] . it is also instructive to show the dynamics of the phase profile during the molecule s oscillation . this is shown with the contour plots of fig . [ fig8 ] which correspond to the molecules of figs . [ fig5]-[fig7 ] . in fig . [ fig8]a , which corresponds to fig . [ fig5 ] , the two solitons start initially in phase . by time the phase of the right soliton , which is the one with higher intensity amplitude and larger oscillation displacement , starts to exceed that of the left soliton . at the point of closest approach , the phase difference is exactly @xmath60 . after that point , the two solitons diverge again , the phase difference starts to decrease , and the cycle is repeated . similar behavior is seen in fig . however , in fig . [ fig8]c , where the two solitons coalesce for a considerable amount of time , the phase difference during the coalescence time is zero . it is thus not completely understood why , in this case , the two solitons still repel each other and eventually split . having established the existence of the soliton molecule from the exact two solitons solution and derived a formula that describes its bond length , here we use this formula to examine more closely the propertie of the soliton molecule and its mechanism of binding . it is clear from eq . ( [ delppmol ] ) that the sinusoidal time - dependence of the solitons relative phase leads to a force of interaction that oscillates between attractive and repulsive and hence allowing for soliton molecule formation . further details of the mechanism of binding will be uncovered by expressing the acceleration , @xmath190 in terms of @xmath2 by substituting for @xmath191 from eq . ( [ delmol ] ) into eq . ( [ delppmol ] ) , to get @xmath192 this shows that the interaction force between the two solitons is the resultant of an attractive part and a repulsive part . the equilibrium bond length , defined by @xmath193 , is given by @xmath194 in consistence with our previous result , the equilibrium bond length diverges as @xmath195 . solving the last equation for @xmath196 and then substituting in eq . ( [ delpp2mol ] ) , @xmath197 simplifies to @xmath198.\ ] ] for small amplitude oscillations , @xmath199 , the last equation gives @xmath200 the restoring force ( @xmath201 ) originates from the phase - dependent terms , @xmath202 . this appealing form of the force of interaction shows that the force between the solitons is of hooke s law type with a spring constant @xmath203 where @xmath204 is the _ bare mass _ of the molecule . expressed in terms of @xmath205 and @xmath125 , the spring constant takes the form @xmath206 which shows that @xmath207 for @xmath143 ; corresponding to a soliton molecule of infinite bond length . furthermore , @xmath208 diverges for @xmath180 , which signifies soliton coalescence , as we have pointed out in the previous subsection . since the frequency of the soliton molecule is given by @xmath209 and the spring constant is given in eq . ( [ keq ] ) , the effective mass @xmath210 will be given by @xmath211 which again diverges at the soliton coalescence condition , @xmath180 . having determined the main properties of the soliton molecule , we can now return back to eq . ( [ delmol ] ) to express @xmath151 as @xmath212 where @xmath213 is the initial value of @xmath151 , which is given by the solitons parameters through @xmath214 the amplitude of the oscillation @xmath215 is thus given by @xmath216 which gives an elastic potential energy @xmath217 it should be noted here that this is equal to the mechanical energy since the initial speed vanishes , @xmath218 . the fact that the potential energy diverges at the coalescence condition @xmath180 is a gain an artifact of the calculation , but it at least indicates that the bond is tighter than cases where @xmath176 or @xmath177 .      here , we investigate the stability of the soliton molecule against break up in the following three collision regimes : @xmath219 reflection by a hard wall , @xmath220 crossing a finite potential barrier , and @xmath221 collision with a single soliton . to that end we solve the gross - pitaevskii equation , eq . ( [ gp ] ) , numerically . as an initial state , we use , for cases @xmath219 and @xmath220 the two solitons solution , eq . ( [ exactsol ] ) , which represents the soliton molecule . for case @xmath221 , we use the superposition of the exact single soliton , eq . ( [ psi1_exact ] ) , with the two solitons solution .    before starting the discussion of results , we point out that in figs . [ fig9]-[fig14 ] , we present the results of this section using spacio - temporal density plots . since the solitons are too thin compared to the spacial range that we consider , a density plot with full spacial and time ranges will not show a clear solitons peak density or center - of - mass path , as fig . [ fig9]c shows . to solve this problem , we restricted the density plotting to a finite range of @xmath166 , namely between 0.025 and 0.15 corresponding to the upper part of the solitons peaks . this results in an easier tracking of both the solitons peak density and center - of - mass path , as shown in fig . [ fig9]a , b , d and the rest of subsequent figures .    for reflection from a hard wall we solve eq . ( [ gp ] ) with a potential step of the form @xmath222 where @xmath223 and @xmath224 are the hight and location of the potential wall , respectively . the result of reflection from this hard wall , with @xmath225 , is shown in fig . the soliton molecule preserves its molecular structure but with different characteristics . the solitons in the reflected molecule do not coalesce as in the incident molecule . in other words , the equilibrium bond length becomes larger . the density plot shows that initially the two solitons are of comparable intensities . after reflection , the brighter color of the left soliton and darker color of the right soliton indicate that the left soliton acquires higher intensity on the expense of the right soliton . we also notice that the left soliton performs two reflections from the potential interface . after the first reflection , it collides with the right soliton and then collides with the potential interface for the second time . the picture becomes different when the hight of the wall is reduced to @xmath226 , as shown in fig . the soliton molecule breaks up after reflection . this is due to loss of energy at the interface of the potential . part of the soliton molecule transmits as a nonsolitonic pulse that broadens and decays in intensity by time . by plotting @xmath166 in fig . [ fig9]c with its full range , we can see the nonsolitonic part as the left- and right - going two red ejections corresponding to the transmitted and reflected nonsolitonic pulses , respectively . in fig . [ fig9]d , we combine fig . [ fig9]b and fig . [ fig9]c to show the locations of the nonsolitonic ejections with respect to the solitons centers . in the case of reflection from a hard wall , the nonsolitonic ejections are essentially not present which results in the stability of the molecular structure .    for reflection from a potential barrier we solve eq . ( [ gp ] ) with the potential @xmath227 where @xmath223 , @xmath228 , and @xmath224 , are the height , width , and location of the right side of the barrier , respectively . in fig . [ fig10 ] , we show the many different possibilities that result when the height of the barrier is changed . the free evolution case with @xmath229 is shown as a reference plot . the full reflection case is shown for @xmath225 , which is similar to the previous case of reflection from a hard wall . reducing the height of the barrier to @xmath230 , we notice that the soliton molecule breaks up after reflection . as pointed above , this is due to the nonsolitonic ejections taking place at the interfaces of the potential . reducing the height of the barrier to @xmath231 , a sign of solitons recombining appears in the form of a soliton molecule of a short lifetime . at @xmath232 a stable molecule is remarkably formed with a considerably shorter period than for the incident molecule . we have confirmed numerically that this molecule remains stable for much longer time provided that the soliton molecule remains sufficiently far from the boundaries of the spacial grid . this unique structure remains for some small domain around @xmath232 , but is lost for @xmath233 , where the molecule breaks for long evolution times . decreasing the height of the barrier to @xmath234 , the molecule breaks at the interface and splits into a reflected and transmitted solitons . for @xmath226 , the molecule breaks at the interface , but both solitons transmit through the barrier . for @xmath235 , the transmitted solitons show a sign of recombining again but with shorter period than for the free evolution case and larger than in the case of @xmath232 .    motivated by the fact that at the coalescence point the intensity of the molecule is considerably higher than at other instants , we expect to find different scattering dynamics when the soliton molecule meets the interface of the potential at different phases of its periodic oscillation . in fig . [ fig11 ] , this is investigated by fixing the height of the potential barrier at @xmath232 and changing the initial launching position of the molecule . starting at @xmath236 , the molecule breaks after reflection . a shortly - lived molecule is obtained at @xmath237 , and a stable molecule is found for @xmath238 , which corresponds to the @xmath232 case of fig . [ fig10 ] . at @xmath239 , the soliton splits at the interface into transmitted and reflected solitons . the coalescence point is , in this case , located at the interface . transmission takes place due to the high intensity of the soliton at the coalescence point . at @xmath240 , the two solitons still split as in the previous subfigure but with a weak transmitted soliton intensity less than 0.025 and hence will not be shown in our density plots which are restricted to the intensities between 0.025 and 0.15 , as pointed out previously . for @xmath241 , the coalescence point takes place before the molecule reaches the interface and both solitons reflect but the molecule breaks up . for @xmath242 the two reflected solitons start to recombine forming a stable molecule at @xmath243 . at @xmath244 the reflected molecule starts to break up since the second coalescence point becomes close to the interface . thus , the conclusion from this figure is that the soliton molecule is more vulnerable to break up when it meets the interface at the coalescence point . equivalently , soliton molecules with larger equilibrium bond length , such that coalescence does not occur , will be more stable against breakup post reflection from barriers .    finally , we present in figs . [ fig12]-[fig14 ] the results of scattering of a soliton molecule by a single soliton described by eq . ( [ newsol1 ] ) with normalization @xmath245 , center - of - mass position and speed @xmath246 and @xmath247 , respectively . the effects of the phase , speed , and amplitude of the injected soliton are investigated separately . in fig . [ fig12]a , a soliton initially at @xmath248 is launched towards a stationary soliton molecule near @xmath249 . at the impact , the molecule brakes up , its right soliton is ejected in the direction of the positive @xmath15-axis , and the left soliton combines with the scatterer soliton to form a new stationary molecule shifted by about a bond length to the left . we point out here that for such an outcome to occur , it is essential that the amplitude of the scatterer soliton is nearly equal to that of the right soliton of the molecule . otherwise , a different outcome , as that of fig .  [ fig14 ] , will be obtained . in fig . [ fig12]b , the same numerical experiment is repeated but with adding a @xmath60 to the phase of the scatterer soliton . clearly , this phase addition prevents the formation of a new molecule resulting in three solitons diverging from each other . from applications point of view , the phase of @xmath60 could be used as a `` key '' to `` unlock '' the molecule for the purpose of extracting stored data . in fig . [ fig13 ] , we show the effect of the initial speed of the injected soliton . in contrary to one s first judgment , the molecule preserves its structure for fast collisions , as in fig . [ fig13]a , and breaks up for slower collisions , as in fig . [ fig13]b . in fig . [ fig14 ] , the injected soliton has an amplitude that is approximately two times larger than any of the two solitons of the molecule . the injected soliton penetrates the molecule leaving it almost unchanged apart from a center - of - mass shift to the left . we have used the inverse scattering method to derive the two solitons solution of a nonlinear schr@xmath0diner equation with a parabolic potential and cubic nonlinearity with time - dependent coefficients , as given by eq . ( [ gp ] ) . the solution was then simplified and put in a suggestive form in terms of the fundamental parameters of the two solitons , namely their amplitudes , center - of - mass positions and speeds , and their phases . in this form , two different regimes of the solution , namely the resolved solitons and overlapping solitons , were distinguished and the main features such as the solitons separation and relative phase were extracted . from the expression for the solitons separation we find that for the homogeneous case and zero solitons relative speed , the solitons separation diverges logarithmically with the solitons amplitude difference such that , for equal solitons , the trivial solution is obtained . the force of interaction was then derived , essentially , for arbitrary solitons parameters . this resulted in generalizing gordon s formula @xcite to _ i _ ) the generalized inhomogeneous case considered here , _ ii _ ) arbitrary solitons relative speed and amplitudes , and _ iii _ ) short solitons separations ( compared to their width ) . with this formula , the possibility of forming soliton molecules emerged naturally , where the force at short distance was shown to be composed of an attractive part , resulting from the nonlinearity , and another part that oscillates between repulsive and attractive resulting from the time - dependent relative phase . the main features of the soliton molecule , including its equilibrium bond length and bond spring constant , frequency and amplitude of oscillation , effective mass , and its elastic potential energy , where then calculated in terms of the solitons parameters . it turns out that the amplitudes difference @xmath63 plays an important role in determining these quantities . furthermore , we show that at the condition @xmath250 , the two solitons coalesce while away from this condition , the solitons approach each other but remain resolved . at this condition , the molecule s effective mass and spring constant have maximum values . in our expressions eqs . ( [ keq ] ) and ( [ masseq ] ) diverge because these formulae were derived assuming the solitons remain resolved .    to have a sense of its stability we investigated numerically _ i _ ) the collision of the soliton molecule with a hard wall and softer one , _ ii _ ) scattering by a potential barrier , and _ iii _ ) collision with a single molecule . the first case showed that while the molecular structure is preserved after reflection from a hard wall , it breaks when reflecting from a softer one . reflection from a finite barrier showed that the molecular structure is preserved only for specific heights of the barrier . for an incident molecule with the coalescence condition satisfied , the molecule will be most vulnerable to break up when the coalescence point takes place at the interface of the barrier . this is simply understood by the fact that the intensity of the soliton molecule is maximum when the two solitons coalesce such that tunneling becomes possible . stability of the molecule was also investigated by scattering the molecule with a single soliton . it turned that slower collisions tend to break up the molecule more easily than faster ones . in addition , the outcome of the collision depends on the phase of the incoming soliton such that a scatterer soliton which is in - phase with the molecule will typically preserve its molecular structure , but for an out - of - phase soliton , the molecule breaks up . the two solitons solution presented here and the analysis that shows how to extract the solitons separation and relative phase may constitute the basis for a more accurate and detailed investigation of the origin of the soliton - soliton force , especially for short separations and at coalescence . the results of this paper will be hopefully of relevance to possible future applications of soliton molecules as data carriers or memories . the author acknowledges helpful discussions with vladimir n. serkin . the lax pair associated with the gross - pitaevskii equation , eq . ( [ gp ] ) , is obtained using our lax pair search method @xcite and reads @xmath251 @xmath252 where +   + @xmath253 , @xmath254 , @xmath255 , @xmath256 , +   +   + @xmath257 , +   +   + @xmath258 , @xmath259 , and @xmath260 and @xmath261 are arbitrary constants . here , @xmath6 is the solution of eq . ( [ gp ] ) and @xmath262 is the _ auxiliary _ field .    the _ compatibility condition _ @xmath263 of the linear system , eqs . ( [ phix ] ) and ( [ phit ] ) , is equivalent to the gross - pitaevskii equation , eq . ( [ gp ] ) , and its complex conjugate . for a known _ seed _ solution , @xmath264 , of eq . ( [ gp ] ) the linear system will have the solution @xmath265 . the darboux transformation is defined as @xmath266={\bf\phi}\cdot{\lambda}-{\bf\sigma}\,{\bf\phi}$ ] , where @xmath266 $ ] is the transformed field and @xmath267 . requiring the linear system to be _ covariant _ under the darboux transformation imposes the transformation @xmath268={\bf p}+{\bf j}\cdot{\bf\sigma}-{\bf\sigma}\cdot{\bf j}$ ] , where @xmath268 $ ] is the transformed @xmath269 . this gives the new solution in terms of the seed solution as @xmath270    using the trivial solution @xmath271 as a seed , the darboux transformation generates the well - known sech - shaped single soliton solution @xcite @xmath272 where @xmath273 @xmath274 @xmath275 @xmath276 , @xmath277 , @xmath278 , and the constant @xmath279 corresponds to an arbitrary overall phase . this solution corresponds to a soliton density profile that is localized at @xmath37 , moving with center - of - mass speed @xmath280 , and normalized to @xmath40 .    using this single soliton solution as a seed , the darboux transformation generates a two solitons solution . the solution of the linear system , eqs . ( [ phix ] ) and ( [ phit ] ) , can in this case be derived and simplified to the following form @xmath281 @xmath282 @xmath283 @xmath284 @xmath285 where , @xmath286 @xmath287 @xmath288 @xmath289    finally , the two solitons solution is obtained by substituting for @xmath290 and @xmath291 in eq . ( [ newsol1 ] ) , which upon substituting @xmath292 and further simplification can then be put in the form of eq . ( [ exactsol ] ) . since @xmath68 and @xmath69 are functions of @xmath15 , we start by examining their values near the roots of @xmath57 . this task can be simplified by rewriting @xmath68 as @xmath293 . for @xmath294 , we get @xmath295 , and for @xmath296 , we have @xmath297 . for @xmath113 the first term in the exponent of @xmath68 is larger than the second one , provided that @xmath40 is not much larger than @xmath298 , as remarked at the end of this section , hence @xmath73 . for @xmath132 the magnitude of the first term in the exponent of @xmath68 is smaller than the magnitude of the second one , which again leads to @xmath73 . for the region @xmath299 , the condition @xmath73 is always satisfied . in conclusion , @xmath73 for all @xmath15 , apart from situations with extreme values of @xmath300 . the situation is simpler for @xmath69 : @xmath74 for @xmath113 , @xmath301 for @xmath302 , and @xmath303 for @xmath132 . thus , to find @xmath50 , we expand @xmath57 in powers of large @xmath68 and @xmath69 and to find @xmath48 and @xmath52 , we expand @xmath57 in powers of large @xmath68 . expanding @xmath57 in powers of @xmath68 and @xmath69 around @xmath58 up to first order in @xmath80 and @xmath304 , we find @xmath305 where @xmath306 the root of this equation gives the center - of - mass position of the right soliton . the first two log terms equal @xmath307 . the third log term is constant and diverges for @xmath94 . the last term is small since it is proportional to @xmath80 , but it is needed because it contains the phase - dependent contributions . due to the combination of @xmath308 and @xmath68 , finding an algebraic expression for the root of this equation will not be possible . instead , we ignore at first the @xmath80 term to obtain the dominant contribution which will then be used to find the @xmath80 contribution . this gives @xmath309 substituting back in eq . ( [ q1 ] ) and solving for @xmath15 , we finally get @xmath310 where @xmath311 .    for the left and central solitons , we expand @xmath57 in powers of @xmath68 ; leaving @xmath69 arbitrary . in this manner , we account for the two roots , @xmath48 and @xmath52 , simultaneously . expanding @xmath57 in powers of large @xmath68 , we get @xmath312 where @xmath313\nonumber\\&/&{\left(g_0 ^ 2 e^{2 \gamma_0 }     ( { n_d}+{n_s } { y } ) ^2 + 4 ( { y } + 1)^2 { \eta_d}^2\right)}.\end{aligned}\ ] ] similar to the above procedure for @xmath50 , we solve first eq .  ( [ q2 ] ) without the @xmath80 term , which gives @xmath314 - 8     { \eta_d}^2}{2 g_0 ^ 2 { n_s}^2 e^{2 \gamma_0}+8 { \eta_d}^2},\ ] ] where @xmath315 in the limit @xmath89 and @xmath90 , the solution @xmath87 approaches 0 , which corresponds to @xmath316 . the solution @xmath317 approaches 1 , which corresponds to @xmath318 . thus , the solutions @xmath87 and @xmath317 correspond the left and central solitons , respectively . since we will be interested only in the left and right solitons , we take for the rest of this section the @xmath87 . to find the contribution of the @xmath80 term , we substitute for @xmath87 in the @xmath80 term of eq . ( [ q2 ] ) , and then solve for @xmath69 to get a corrected expression for @xmath87 @xmath319     -8 { \eta_d}^2}{2 g_0 ^ 2 { n_s}^2 e^{2     \gamma_0}+8 { \eta_d}^2},\ ] ] where @xmath320 , @xmath321 , @xmath322 , @xmath323 , and @xmath324 . expanding for large @xmath325 and then solving for @xmath15 , we finally get @xmath326    final remarks about the validity of the above derivation are in order . the condition @xmath73 will be met in the region @xmath113 only for @xmath327 . note that the numerator of the right hand side of this inequality is less than the denominator by at least @xmath328 . for @xmath132 , the condition @xmath73 will be met only for @xmath329 . here , the numerator of the right hand side of this inequality is larger than the denominator by at least @xmath328 . a more quantitative estimate for the ratio @xmath300 can be obtained using the above results for @xmath330 and @xmath331 . furthermore , the above - derived formula for @xmath48 is limited to values of the parameters for which the quantities @xmath332 and @xmath333 are positive . at @xmath334 , the two roots @xmath48 and @xmath52 coincide . in fig . [ fig1 ] , this corresponds to the maximum of the @xmath57-curve lying on the @xmath15-axis . for @xmath335 , the maximum of the curve is below the @xmath15-axis and the two roots become nonreal .        noting that : + @xmath341 + @xmath342 + @xmath343 + @xmath344 + @xmath345 + @xmath346 , +   + the acceleration , @xmath190 , can be calculated to take the form of eq . ( [ force2 ] ) with coefficients @xmath347 @xmath348 @xmath349 @xmath350      a. hasegawa and y. kodama , _ solitons in optical communications _ ( oxford university press , new york , 1995 ) ; l.f . mollenauer and j.p . gordon , _ solitons in optical fibers _ ( acadamic press , boston , 2006 ) ; g.p . agrawal , _ nonlinear fiber optics _ ( academic press , san diego , 2001 ) , 3rd ed . ; n.n . akhmediev and a. ankiexicz , _ solitons : nonlinear pulses and beams _ ( chapman and hall , london , 1997 ) . l. khaykovich , f. schreck , g. ferrari , t. bourde , j. cubizolles , l.d . carr , y. castin , c. salomon , science * 296 * , 1290 ( 2002 ) . m. stratmann , t. pagel , and f. mitschke , phys . * 95 * , 143902 ( 2005 ) .                                              ) , @xmath47 $ ] . thick ( green ) curve : the soliton intensity @xmath351 . light ( red ) curve : the soliton intensity with only the second sech term of eq . ( [ exactsol ] ) . the parameters used are : @xmath16 , @xmath352 , @xmath353 , @xmath354 , @xmath355 , @xmath356 , @xmath357 , @xmath358 , and @xmath167.,width=566 ]    ) . dotted ( blue ) curve corresponds to the phase of the right soliton , @xmath102 , calculated from eq . ( [ phaser ] ) . thick dashed curve corresponds to the phase of the central soliton , @xmath103 , calculated from eq . ( [ phasec ] ) . dashed - dotted curve corresponds to the phase of the left soliton , @xmath104 , calculated from eq . ( [ phasel ] ) . the solid curves correspond to the density profiles as in fig . the parameters used are : @xmath16 , @xmath352 , @xmath353 , @xmath354 , @xmath355 , @xmath359 , @xmath360 , @xmath357 , @xmath358 , and @xmath167.,width=566 ]    , versus @xmath124 . thick ( blue ) curve is calculated from eq . ( [ force2 ] ) . light ( red ) curve is calculated numerically from the exact two solitons solution , eq . ( [ exactsol ] ) . the parameters used are : @xmath16 , @xmath352 , @xmath353 , @xmath354 , @xmath355 , @xmath359 , @xmath357 , @xmath358 , and @xmath167.,width=566 ]      ) . * ( b ) * thick ( green ) curve is solitons separation calculated from the exact solution ( [ exactsol ] ) . light ( red ) curve is the solitons separation calculated from formula ( [ delmol ] ) . * ( c ) * density profile at some specific times . the parameters used are : @xmath16 , @xmath361 , @xmath362 , @xmath363 , @xmath356 , @xmath364 , @xmath358.,width=566 ]        -[fig7 ] . subfigures * ( a ) * , * ( b ) * , and * ( c ) * , correspond to fig . [ fig5 ] , fig . [ fig6 ] , and fig . [ fig7 ] , respectively . the blue ( lower ) and green ( upper ) curves correspond to the center - of - mass trajectories of the left and right solitons , respectively.,width=566 ]    ) . the dashed vertical line represents the interface of the potential step . the parameters used are : @xmath364 , @xmath366 , @xmath367 , @xmath368 , @xmath369 , @xmath370 , @xmath371 . * ( a ) * @xmath225 . * ( b),(c),(d ) * @xmath226 . ,
the problem of soliton - soliton force is revisited . from the exact two solitons solution of a nonautonomous gross - pitaevskii equation 
 , we derive a generalized formula for the mutual force between two solitons . 
 the force is given for arbitrary solitons amplitude difference , relative speed , phase , and separation . 
 the latter allows for the investigation of soliton molecules formation , dynamics , and stability . 
 we reveal the role of the time - dependent relative phase between the solitons in binding them in a soliton molecule . 
 we derive its equilibrium bond length , spring constant , frequency , effective mass , and binding energy of the molecule . 
 we investigate the molecule s stability against perturbations such as reflection from surfaces , scattering by barriers , and collisions with other solitons .