Wan-Yu Tsai

Peng-Cheng Gao

Barbara Daffos

Pierre-Louis Taberna

Carlos R Perez

Yury Gogotsi

Frederic Favier

Patrice Simon
email: simon@chimie.ups-tlse.fr

Ordered mesoporous silicon carbide-derived carbon for high-power supercapacitors

Keywords:

Micro/mesoporous carbon was prepared by chlorination of ordered mesoporous silicon carbide derived from magnesio-thermal reduction of templated carbon-silica precursors. These materials were then used as active materials for electrochemical capacitors and characterized in 1.5 M NEt 4 BF 4 /AN. The electrodes showed outstanding rate capability (90% of capacity retention at 1 V/s and time constant of 1 s) with high specific areal capacitance (0.5 F/cm 2 of electrode), that makes such hierarchical porous carbons promising for high power and energy density supercapacitors.

Introduction

Electrochemical double layer capacitors (EDLCs), also known as supercapacitors, store charge by simple electrostatic attraction between electrolyte ions and charges at the electrode surface. Due to a different nature of energy storage, supercapacitors hold higher power density (15 kW/kg) and much better cyclability (> 10 6 ) as compared to batteries [1,2].

Carbon-based materials with high specific surface area (SSA) have been widely studied as supercapacitor electrodes due to their chemical stability and high conductivity. While activated carbon is the material of choice for commercial supercapacitors, carbide derived carbons (CDCs) [3], onion-like carbons (OLCs) [4,5], carbon nanotubes (CNTs) [6], and graphene [7] have also been investigated. Among those, microporous CDCs are highly attractive thanks to their tunable structure leading to high capacitance and high power capabilities [8][9][10]. CDCs are prepared by chlorination of metal carbides; depending on the choice of precursor material and synthesis conditions, the resulting microstructure can be nicely controlled with a narrow pore size distribution, while retaining the original carbide shape and volume (a conformal transformation) [3]. We recently developed an original route for the synthesis of ordered mesoporous carbides by a magnesiothermal reaction on silica-carbon mesoporous composites [11]. The present work investigates the electrochemical performance of such ordered micro/mesoporous silicon carbide-CDC (SiC-CDC) materials in supercapacitor applications.

Experimental

Synthesis of SiC-CDC

Mesoporous carbon-silica composites were prepared with slight modifications of the synthetic route described by Liu et al. [12]. A triblock copolymer P123 (Aldrich) was chosen as soft template. In a typical procedure, under vigorous stirring, 1 g of 0.2 M HCl was added to 9.6 g of 20 wt.% P123 alcoholic solution. After stirring at 40 °C for 4 h, a clear solution was obtained. Subsequently, 2.5 g of 20 wt.% resol alcoholic solution and 2.08 g tetraethoxysilane (TEOs, 98% Alfa Aesar) were added. The mixture was continuously stirred for 2 h and painted afterwards onto a glass plate. The resulting film was kept at ambient temperature overnight to evaporate residual ethanol before being polymerized at 100 °C for 24 h. The as-made film was scraped and ground into fine powder before being carbonized in a tube furnace at 800 °C for 3 h under N 2 flow. C/Si molar ratio was measured at 1.04:1 by residual silica content after burning SiC precursors at 550 °C for 2 h.

The mesoporous silicon carbide was synthesized by using Mg as reducing agent at moderate temperatures as described in our recent work [11]. After sealing a carbon/SiO 2 composite mixed with Mg powder (10% Mg molar excess over SiO 2 ) under Ar in a 316 L stainless steel tube, the magnesio-thermal reduction was performed at 800 °C for 24 h at a heating rate of 1 °C/min. SiC product was purified using 1 M HCl and 1 M HNO 3 aqueous solutions in sequence.

CDC powders were produced by chlorination of the SiC precursor, as previously described [3]. The precursor was placed in a horizontal tube furnace and heated to 800 °C in Ar, followed by Cl 2 flow (10-15 cm 3 /min) for 3 h. After flushing residual Cl 2 with Ar and lowering the temperature to 600 °C, the sample was annealed for 2 h using H 2 .

Material characterization

Transmission electron microscopy (TEM) images were obtained using a JEOL 1200 EX2 TEM operating at 100 kV. Porosity characteristics were calculated from argon sorption isotherms measured at 87 K with Micromeritics ASAP 2020 porosimeter. The specific surface area was estimated by using both Brunauer-Emmett-Teller (BET) and density functional theory (DFT) method while the pore volumes and the pore size distributions were calculated from adsorption isotherms by using quenched solid density functional theory (QSDFT) method [13].

Electrochemical characterization

SiC-CDC electrode films were prepared by mixing 95 wt.% of SiC-CDC powder with 5 wt.% of polytetrafluoroethylene (PTFE; 60 wt.% dispersion in water) using ethanol as the medium for mixing the slurry.

Electrode film thickness was about 210 μm and the SiC-CDC loading was 7.3 mg/cm 2 , close to that of commercial electrodes (100-200 μm or about 10 mg/cm 2 ) [14], leading to an apparent electrode density of 0.35 g/cm 3 . Electrochemical characterization of SiC-CDC electrodes was conducted in 1.5 M tetraethylammonium tetrafluoroborate (NEt 4 BF 4 ) (Acros Organics, CAS #429-06-1) in acetonitrile (AN) (Acros Organics, CAS #75-05-8, H 2 O b 10 ppm) in a two-electrode Swagelok cell [15] assembled in a glove box under argon atmosphere, at room temperature. The electrode films were laminated onto treated aluminum current collectors [8], and two layers of 25 μm-thick porous cellulose were used as separator.

Cyclic voltammetry (CV) tests at various scan rates were performed between 0 and 2.3 V. Electrochemical impedance spectroscopy (EIS) measurements were performed at the rest potential, using a sinusoidal signal of ±5 mV from 50 kHz to 10 mHz. Electrochemical tests were performed using a multichannel VMP3 potentiostat/galvanostat (Bio-Logic, France).

Results and discussion

The TEM image, gas adsorption/desorption isotherm and the pore size distribution obtained from QSDFT (Fig. 1) show that the mesoporous structure of the silicon carbide was retained after chlorination, with a significant fraction (~30%) of the pores having sizes greater than 2 nm. Additionally, the SiC-CDC contains narrowly distributed micropores as a peak at 1 nm is observed in Fig. 1b (inset). The specific surface areas of SiC-CDC calculated from BET method and DFT method are 1619 m 2 /g and 1335 m 2 /g, respectively.

Fig. 2a shows cyclic voltammograms of the SiC-CDC at various scan rates from 20 to 300 mV/s in 1.5 M NEt 4 BF 4 in AN electrolyte. All the CVs exhibit a rectangular shape, which is characteristic of a pure capacitive behavior. When the potential scan rate increased, the rectangular profile was nicely preserved up to 300 mV/s, although a slight distortion appeared when the potential was switched back. The perfect symmetry in the distortion evidences an ohmic origin, which is linked to the ohmic drop in the bulk electrolyte. The specific capacitance of the prepared porous carbon is 80 F/g. This value is certainly not among the highest reported in the literature for porous carbons in organic electrolyte [16]; however, taking into account the carbon weight loading of 7.3 mg/cm 2 , it corresponds to a specific capacitance of more than 0.5 F/cm 2 per electrode, which is far beyond most values reported in the literature [4,6,17,18] at the same scan rate. The specific capacitance per cm 2 of the electrode is representative to the energy that could be stored in a real device [14].

The EIS Nyquist plot of the cell (Fig. 2b) shows a typical capacitive behavior of a porous blocking electrode [19] as well as an extremely low resistance: an equivalent series resistance (ESR) of 0.25 Ω cm 2 and an ionic resistance of the electrolyte inside the porous carbon electrode of 0.35 Ω cm 2 . Both ESR and ionic resistance obtained are smaller than those found in the literature for other carbon-based materials (activated carbon [19], other CDCs [20], CNT [4], carbon onion [4], graphene [7]) in the same electrolyte under similar testing conditions. While micropores provide a good charge storage capability, a high power is anticipated from the low intrinsic resistance of SiC-CDC and high ion accessibility through the mesopores.

Fig. 3a shows CVs at high scan rates from 20 to 3000 mV/s and Fig. 3b shows the gravimetric capacitance (derived from the slope of the Q-V curve during cell discharge) and the areal capacitance (obtained by dividing the electrode capacitance by the geometric area of the cell) as a function of the potential scan rate. The capacitive CV shape was preserved up to 1 V/s, with a gravimetric capacitance of 75 F/g, corresponding to > 90% of capacitance retention; this already evidences the high power performance of the material.

The frequency behavior of the cell was analyzed using a complex capacitance model [19], based on the modeling of the capacitance in real part C′(ω) and imaginary part C″(ω), both functions of the frequency. Fig. 4a shows the change of the real part of the capacitance C′(ω) with frequency. The graph shows a transition from purely resistive behavior at high frequency (phase angle close to 0) to purely capacitive behavior at low frequency (phase angle of -87°). The whole capacitance of the electrode is reached at around 0.2 Hz since a plateau at high frequency is observed; the capacitance value for a phase angle of -45°obtained at 1 Hz is 0.25 F/cm 2 . This corresponds to half of the total capacitance of the electrode [19]. The change of C′ vs. frequency first evidences the high power capability of the prepared porous carbon, since the transition frequency is high, at 1 Hz. Higher power can be reached using thin films in microdevices for instance, or exohedral carbons with accessible but limited surface area. However, it is achieved at the expense of the energy density. As such, the areal capacitance (F/cm 2 ) of such devices or electrodes is low, from μF/cm 2 to few mF/cm 2 [17,18,21,22]. In contrast, thanks to the hierarchical porous structure of the SiC-CDC prepared in this work, an attractive areal capacitance of 0.5 F/cm 2 was obtained, demonstrating a higher energy density.

Fig. 4b shows the imaginary part of the capacitance, C″(ω), vs. frequency. The relaxation time τ 0 in this system is 1 s, which can be calculated from τ 0 = 1/f 0 , where f 0 is the frequency corresponding to the maximum of the curve of C″ vs. f. These micro/mesoporous carbons thus show high power performance with high specific capacitance. These results suggest that the magnesio-thermal reduction of SiC is an efficient way for producing carbon for high power applications when high energy density is also required.

Conclusions

Evaluation of the ordered mesoporous SiC-CDC as an electrode material for supercapacitors was done in 1.5 M NEt 4 BF 4 /AN electrolyte, giving respectable gravimetric capacitance (80 F/g), and high capacity retention up to 1 V/s, a high rate compared to other high-SSA porous carbons. EIS study showed a low intrinsic resistance of CDC, high ion accessibility to the micro/mesoporous hierarchical structure, and thus a remarkably low relaxation time of the system, evidencing that the use of such hierarchical porous carbons for high power supercapacitors is very promising.

Fig. 1 .Fig. 2 .

 12 Fig. 1. TEM image of SiC-CDC (a) and gas adsorption/desorption isotherms of SiC-CDC (SSA ~1619 m 2 /g); inset shows the pore size distribution of the SiC-CDC powder derived from QSDFT (b).

Fig. 3 .Fig. 4 .

 34 Fig. 3. CVs at scan rates up to 3 V/s (a) and capacitance vs. scan rate (b).

Acknowledgments P.C.G. would like to thank the Chinese Scientific Council for fellowship. W.-Y.T. was supported by the European Research Council (ERC, Advanced Grant, ERC-2011-AdG, Project 291543 -IONACES). P.S. acknowledges the funding from the Chair of Excellence of the EADS foundation "Embedded multi-functional materials". C.P. and Y.G. are thankful to the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Science which supported this research under Award Number DESC0001912. The US-France collaboration was enabled through support from the Partnership University Fund (PUF).


Radoniaina Andriatsimanefitra

Thomas Genet

Laurent Guillo

Jean-Francois Lalande

David Pichardie

Valérie Viet

Triem Tong

Projet Kharon du Labex CominLabs Découvrir, comprendre et reconnaître des malware Android par suivi de flux d'information

L'avènement des téléphones et tablettes ces dernières années a favorisé le développement d'un nouveau modèle économique qui repose sur la livraison d'un téléphone nu sur lequel l'utilisateur peut installer des applications via des magasins d'applications ou marchés. Ces marchés sont aussi bien officiels car approuvés par les entreprises développant les systèmes d'exploitations de ces téléphones comme Google Play ou AppStore ou non officiels. Dans ce qui suit, nous nous consacrons aux téléphones et tablettes utilisant le système d'exploitation Android sur lesquelles un utilisateur peut installer des applications via Google Play. Sur ce marché, nous trouvons beaucoup d'applications (plus de un million d'applications) téléchargées par beaucoup d'utilisateurs. Une même application peut être téléchargée plus de 100 millions de fois. Les développeurs de ces applications sont d'origines très variées puisque pour la modique somme de 25 USD n'importe qui peut créer un compte développeur pour ce marché et distribuer sa propre application. Il est aussi possible d'installer des applications sans passer par Google Play. Android est ainsi devenu rapidement une cible de choix pour les développeurs malveillants et le nombre de malware n'a cessé d'augmenter pour atteindre 10 millions d'applications malveillantes en 2012-2013 selon Kaspersky [START_REF]Kaspersky web site[END_REF].

Ces malware essaient principalement de soutirer de l'argent à un utilisateur. Par exemple, le malware Simplelocker [START_REF]Android community web site[END_REF] chiffre sur le téléphone les données de l'utilisateur comme ses photos, ses documents les rendant ainsi inutilisables et demande ensuite une rançon pour déchiffrer ces mêmes données. D'autres comme DroidKungFu1[3] installent des applications à l'insu de l'utilisateur. D'autres encore envoient des sms à des numéros sur-taxés à l'insu, là encore, de l'utilisateur. Ces malware sont souvent repackagés dans différentes applications pour atteindre un plus grand nombres d'utilisateurs et se disséminer plus aisément dans l'offre disponible.

Pour débarrasser le marché de ces applications malveillantes, il existe deux grandes approches : l'analyse statique du bytecode des applications et l'analyse dynamique des exécutions de ces applications.

Le but du projet Kharon est de développer une plateforme d'analyse dynamique utilisant le suivi de flux d'information. Cette analyse dynamique surveillera le comportement d'un malware et étudiera comment ce malware se dissémine dans le système d'exploitation. Nous avons montré dans [START_REF] Andriatsimandefitra | Capturing android malware behaviour using system flow graph[END_REF] que cette approche permet de créer des signatures comportementales de malware, c'est à dire des descriptions précises du comportement de ces malware. Cette approche se base sur l'hypothèse forte que le malware s'exécute effectivement durant l'analyse dynamique. Cette hypothèse n'est pas acceptable en pratique car les expériences que nous avons menées montrent que la plupart des malware ne s'exécutent pas automatiquement au lancement de l'application. Ces malware attendent des évènements particuliers comme l'envoi d'une commande par un serveur distant, une date précise, une série d'actions de l'utilisateurs. Dans le projet Kharon, nous proposons utiliser l'analyse statique pour découvrir quels sont ces évènements déclencheurs et développer une méthode de fuzzing utilisant ces données.

Ce projet regroupe des membres de l'équipe Celtique et Cidre de l'IRISA. L'équipe Celtique est spécialisée en analyse statique de bytecode java [START_REF] Hubert | Sawja : Static Analysis Workshop for Java[END_REF] et l'équipe Cidre a proposé les signatures comportementales utilisées dans ce projet. Les membres à temps plein de ce projet sont actuellement un ingénieur expert et un étudiant en master 2. Ce projet a débuté en janvier 2015.

Dans le cadre de la conférence RESSI, nous proposons de faire une présentation des objectifs de Kharon et de les illustrer au travers d'une démonstration.

Références


G Araujo Pimentel
email: <guilherme.araujopimentel@umons.ac.be

A Rapaport
email: rapaport@supagro.inra.fr

A Vande Wouwer

Nonlinear Model Predictive Control of a Wastewater Treatment Process Fitted with a Submerged Membrane Bioreactor

Keywords: process control, optimization, mathematical modeling, membrane bioreactor

Submerged membrane bioreactors are increasingly applied for wastewater treament but requires a tight control of the membrane fouling so as to ensure safe and efficient operation. The objective of this paper is to design a nonlinear model predictive control to minimize the irreversible resistance while keeping the trans-membrane pressure, which is a good indicator of membrane fouling, at an acceptable level. To this end, the manipulated variables are the permeate flow and the air scouring flow, which allows the material layer formed on the membrane (in short the "cake") to be detached. The NMPC structure is tested in simulation considering a detailed simulator as the reference process, and a reduced-order model as the predictor. The results show that the process can be regulated until the irreversible resistance takes the main role in the fouling resistance. When this state is reached, a chemical cleaning is required, or a larger trans-membrane pressure has to be accommodated.

INTRODUCTION

Submerged membrane bioreactors (sMBRs) are a relatively new technology for wastewater treatment. Their main advantages are related to the outstanding water quality, the solid effluent concentration, the separation between solid retention time (SRT) and hydraulic retention time (HRT) and the smaller process footprint [START_REF] Cicek | A review of membrane bioreactors and their potential application in the treatment of agriculture wastewater[END_REF][START_REF] Bella | An integrated model for physical-biological wastewater organic removal in a submerged membrane bioreactor: Model development and parameter estimation[END_REF][START_REF] Sarioglu | Model evaluation of simultaneous nitrification and denitrification in membrane bioreactor operated without an anoxic reactor[END_REF]. However, one of the main drawbacks of this technology is the membrane fouling that hampers the efficiency of the process [START_REF] Le-Clach | Fouling in membrane bioreactors used in wastewater treatment[END_REF][START_REF] Zarragoita-González | Modeling of submerged membrane bioreactor: Conceptual study about link between activated sludge biokinetics, aeration and fouling process[END_REF].

Nowadays, the operation of an sMBR plant is controlled by splitting it in two distinct parts: (i) the filtration efficiency control; (ii) the biological degradation process. In this study, a simple model integrating both physical and biological parts is implemented. The macroscopic model structure, characterized by a modest size, allows the application of a nonlinear model predictive control (NMPC) that has been used in industrial process control with an impressive success rate in the last decades (Kouvaritakis and Cannon, 2001). In the present case, NMPC can be used to act on the irreversible fouling resistance of the membrane while maintaining the trans-membrane pressure (T M P ) at an acceptable level using the combination of permeate flow and air scouring flow (or air cross-flow). This paper is organized as follows. In Section 2, a detailed model of an sMBR is presented along with the reducedorder dynamic model. A nonlinear model predictive controller is designed in Section 3 while Section 4 discusses the numerical results. Finally, Section 5 draws some conclusions and outlines possible future investigation.

sMBR MODELS

A schematic view of an sMBR is presented in

Integrated Models

Many single process models have been developed for sM-BRs. On the one hand, the biological aspects are usually based on the activated sludge model (ASM) framework [START_REF] Henze | Activated sludge models ASM1[END_REF], while dedicated models have been proposed to describe the membrane fouling mechanism. However, these processes are strongly linked and an integrated model is needed [START_REF] Mannina | An integrated model for biological and physical process simulation in membrane bioreactor (MBRs)[END_REF]. [START_REF] Naessens | Critical review of membrane bioreactor models -part 2: Hydrodynamic and integrated models[END_REF] point out the importance of developing an integrated model and the risk of using these separate models for optimization of the entire system, since optimal conditions for one process may not at all be optimal for another process (e.g. optimizing the air cross flow can modify biomass properties resulting in a faster/slower degradation of the substrate).

Integrated sMBR models have been proposed and applied for the first time by [START_REF] Lee | Modeling of submerged membrane bioreactor process for wastewater treatment[END_REF] and [START_REF] Wintgens | Modeling of membrane bioreactor system for municipal wastewater treatment[END_REF]. The former presented an ASM1-SMP model (activated sludge model No.1 with soluble microbial products, which are free in mixed liquor) and the latter an ASM3-EPS model (activated sludge model No.3 with extracellular polymeric substances that are biopolymers attached on the microorganisms). Several other models were proposed afterwards, as for instance Li andWang (2006), Zarragoita-González et al. (2008) and [START_REF] Bella | An integrated model for physical-biological wastewater organic removal in a submerged membrane bioreactor: Model development and parameter estimation[END_REF].

Usually, these models have a large number of parameters, as evidenced by the model proposed by [START_REF] Mannina | An integrated model for biological and physical process simulation in membrane bioreactor (MBRs)[END_REF] that includes 45 parameters. In this same study, a sensitivity analysis showed that only around 30% of the parameters are crucial to process dynamics emulation. Instead of considering a detailed model, it is therefore more appropriate to develop a simplified integrated model for control purposes. This will be dealt with in the next subsection.

The biological sub-model proposed by [START_REF] Mannina | An integrated model for biological and physical process simulation in membrane bioreactor (MBRs)[END_REF] is a modified version of the well-known ASM1 [START_REF] Henze | Activated sludge models ASM1[END_REF], which takes the influence of SMPs into account. In short, the model uses Monod-type kinetics for the degradation of the different substrates and mass balance equations for the different substrate modeling biomass growth, biomass decay, ammonia and carbon removal processes. The filtration sub-model is an improved version of the model proposed by [START_REF] Bella | An integrated model for physical-biological wastewater organic removal in a submerged membrane bioreactor: Model development and parameter estimation[END_REF]. This model is a resistance in series model with R total modeled as follows:

R total = R m + R p + R sc + R dc (1) R irr R rev
where R m is the intrinsic membrane resistance, R irr is the irreversible resistance, which can only be removed by chemical cleaning while R rev denotes the reversible resistance that is affected by the air cross-flow.

The specific pore-blocking resistance R p is proportional to the amount of permeate produced and is computed through the specific pore fouling resistance (r p ) and in terms of the filtrate volume.

R p = r p ∑ Jt f (2) J [m 3 /m 2 d]
is the filtration flux, t f is the duration of the filtration period. The stable sludge cake resistance R sc and the dynamic sludge cake resistance R dc are presented in equations (3a) and (3b).

R sc = r sc • M sc R dc = r dc • M dc (3a) (3b)
These expressions include M sc [kg/m 2 ] that denotes the stable sludge cake mass and M dc [kg/m 2 ] the dynamic sludge cake mass onto membrane surface.

         dM sc dt = 24C ss J 2 24J + C d d p G - β b (1 -α)GM 2 sc γ b V f t + M sc dM dc dt = -β b (1 -α)GM 2 dc γ b V f t + M dc (4a) (4b)
where C ss [kg/m 3 ] is the sludge concentration, C d is the coefficient of the drag and lift force, d p [m] is the particle size, β b is the erosion rate coefficient of the dynamic sludge film, α is the stickiness of the biomass particle, γ b [kg/m 3 s] is the compression coefficient for the dynamic sludge film,

V f = Jt f [m 3 /m 2 ] is water production within a filtration period of an operation cycle, t [d] is time and G [s -1 ] is the shear intensity G = √ ρ s • g • Q a µ S ( 5
)
where ρ S and µ S are the density and viscosity, respectively, of the sludge mixture, Q a is the air cross-flow and g is the gravitational constant.

Simplified Integrated Model

In a previous work [START_REF] Pimentel | Design, analysis and validation of a simple dynamic model of a submerged membrane bioreactor[END_REF], the simplified model ( 6) is presented and validated by comparison with a more complex biological and filtration model.

                 dβ dt =γβ dS dt =- 1 Y µ(S)X+ Q in V (S in -S) dX dt = ( µ(S) - Qw V ) X+ Q in V X in - Qout V X+ J air V µ air (m)m dm dt =QoutX -J air µ air (m)m (6a) (6b) (6c) (6d) with µ(S) = µ S,max S K S + S , µ air (m) = β m K air + m
The membrane filtration is modeled by equations ( 6a) and (6d). These equations represent the sludge cake dynamic on the membrane surface. The long-term sludge cake evolution β is a key factor of the model and represents the ease (or difficulty) of detaching the sludge cake from the membrane by air scouring, influenced by the air cross-flow J air . Considering a process with constant permeate flow (Q out ≈ constant), the capacity of J air to detach the sludge cake decreases with time, due to the drag force on the solid particles, which results in a decay behavior of β (difficulty of detachment), and thus γ [d -1 ] is a negative constant. It must be noted that if the process has constant trans-membrane pressure (T M P ≈ constant), the permeate flow decreases with time. Hence, β increases, and γ is positive, therefore the efficiency of J air increases as a consequence of the loss of the drag force of the membrane to the particle deposition. However, in this case, the model should include a term limiting the growth of β to a practical maximum. Depending on the process setup (amplitude of permeate, waste and air flow, bulk concentration, process relaxation and filtration cycles), the value of γ changes. Based on the observation of experimental data, it can be concluded that this phenomenon has a slow behavior (Merlo et al., 2000), thus γ has a small value. The state vector includes the long-term sludge cake mass evolution β [m -1 ] with β(0) = β 0 > 0 and the cake mass m [g] with m(0) = m 0 > 0. Both states are positive and bounded if γ < 0; Q out ≤ Q out ≤ Q out is the permeate flow and J air ≤ J air ≤ J air is the air cross-flow.

The sludge cake mass dynamic (6d) is divided in two terms. The first term represents the attachment of suspended particles on the membrane surface, ruled by the effluent flow rate Q out and suspended solids bulk concentration X. The second term represents the sludge cake detachment proportional to air cross-flow, governed by J air and also by the sludge cake mass m and K air [g], the half-saturation coefficient.

The biological activity is described using a simple chemostat reactor [START_REF] Smith | The theory of the chemostat[END_REF] extended to an sMBR, involving one biomass growing on a limiting substrate. It is important to note that this simple biological model structure can be easily extended to more than one biological reaction (see, for instance, [START_REF] Dochain | Dynamical modeling and estimation in wastewater treatment processes[END_REF]).

Equation (6b) represents the consumption of the substrate by the free biomass, ruled by a Monod law µ(S) = µ S,max S K S +S and the transportation of incoming and outgoing substrate through the tank. Note that the substrate is not affected by the membrane, keeping in mind that only solid matter are retained. Equation (6c) shows that there is an interaction between the suspended solid matter and sludge cake mass buildup. The first part of the equation represents the growth of the free biomass that consumes the substrate. Material transportation involves the cake attachment -QoutX V and detachment and the instantaneous "conversion" in suspended solids due to the air cross-flow + Jair V µ air (m)m. The free particle matter is transformed in sludge cake mass and vice-versa depending on the process input values. The waste flow is represented by Q w and the influent is defined as

Q in = Q w + Q out .
The biological time scale is governed by the consumption rate of substrate and, consequently, the growth of biomass. This rate is represented by a Monod law equation (µ(S)) and is normally measured in days.

Usually, the sludge cake mass cannot be measured due to the lack of reliable sensors [START_REF] Dalmau | Development of a decision tree for the integrated operation of nutrient removal MBRs based on simulation studies and expert knowledge[END_REF]. Thus, the trans-membrane pressure is used to indirectly control the sludge cake mass build-up, according to the following equation:

T M P = Q out ηR total . ( 7
)
Q out is the permeate pump flow selected by the operator (or controller) and η [mbar

• d] is the apparent viscosity. The total fouling resistance, R total,s [m -1 ], is modeled by R total,s = R m + R cake + δ R , ( 8
)
where R m [m -1 ] is the intrinsic resistance (assumed constant) and δ R is used to represent the total resistance disturbance, resulting from pore-blocking, biofilm, concentra-tion polarization and scaling. R cake is the cake resistance and, as previously reported [START_REF] Lee | Modeling of submerged membrane bioreactor process for wastewater treatment[END_REF][START_REF] Khan | Prediction of membrane fouling in MBR systems using empirically estimated specific cake resistance[END_REF], it is the most important term in the fouling buildup. Since R m is proportionally small, the total resistance can thus be modeled as the cake resistance:

R cake = ρ m + m 0 A , ( 9
)
where ρ [m • g -1 ] is the specific cake resistance, m 0 [g] is the initial cake mass, m [g] is the cake mass and A is the membrane area.

This simplified model allows the use of advanced control theory such as NMPC and has also been validated with real data from a pilot plant with data set at 30 days long, reaching an accuracy correlation coefficient of R 2 = 0.95 [START_REF] Pimentel | Validation of a simple fouling model for submerged membrane bioreactor[END_REF]. These results encourage us to cast the model in an NMPC optimization framework, but before, the controllability study of the model is achieved.

Controllability Study of the Simplified sMBR Model

The study of controllability is important in order to understand in which way the process behavior can be influenced through certain inputs. According to [START_REF] Slotine | Applied nonlinear control[END_REF], one of the most used mathematical tools for nonlinear systems is the Lie Brackets, briefly presented in the following.

Assuming Q out and J air are the manipulated variables and rewriting system (6) in the following form ẋ = f (x) + g(x)u (10) Definition 2.1. Let f and g be two vector fields on R n . The Lie Brackets of f and g is a third vector field defined by

[f , g] = ▽g f -▽f g (11)
It can be defined as (ad 1 f ) ≡ [f, g] where "ad" means "adjoint". Theorem 2.1. The system defined by:

ẋ = f (x) + m ∑ i=1 g i (x)u i (12)
is locally accessible at about x 0 if the accessibility distribution C spans an n space, where n is the rank of x and C is defined as follows:

C = [ g 1 , g 2 , . . . , gm, [g i , g j ], . . . , [ad k g i , g j ], . . . , [f, g i ], . . . , [ad k f , g i ], . . .

]

The system is therefore controllable.

Applying the theorem 2.1 and considering Q out and J air as process actuators the system is controllable.

NMPC-sMBR PROCESS CONTROL

The main motivation behind the development of a simple dynamic model is the potential of applying advanced model-based control as NMPC. The advantage of this control technique is its ability to handle model nonlinearity and various types of constraints on the actuators and state variables [START_REF] Santos | Nonlinear model predictive control of fedbatch cultures of micro-organisms exhibiting overflow metabolism[END_REF].

NMPC uses a model to predict the trajectory of the system on a prediction horizon and computes an optimal control Fig. 2. NMPC-Scheme sequence on a control horizon [START_REF] Allgöwer | Nonlinear model predictive control: from theory to application[END_REF]. The first important element is therefore a nonlinear model in the form:

ẋ = f (x(t), u(t)), x(0) = x 0 , ( 13
)
together with constrains in the form u(t) ∈ U, ∀t ≥ 0, x(t) ∈ X, ∀t ≥ 0 where x(t) ∈ R n and u(t) ∈ R m are the vector of states and inputs, respectively. The sets U and X are compact and can be represented by

U := u ∈ R m |u min ≤ u ≤ u max , and X := x ∈ R n |x min ≤ x ≤
x max with the constant vector u min , u max and x min , x max .

The NMPC control moves are usually the results of a finite horizon open-loop optimal control problem, which is solved at every sampling instant. In generic notation, the NMPC problem can be expressed as:

min ϕu(•) J N M P C (x(t), ϕ u (•)) s.t. φx (τ ) = f (ϕ x (τ ), ϕ u (τ )), ϕ x (t) = x(t), ϕ u (τ ) = ϕ u (t + T c ), ∀τ ∈ [t + T c , t + T p ], (14a) (14b)
with the cost function

J N M P C (x(t), ϕ u (•)) := t+Tp ∫ t F (ϕ x (τ ), ϕ u (τ ))dτ ;
where

ϕ u (τ ) ∈ U, ∀τ ∈ [t, t + T c ], ϕ x (τ ) ∈ X, ∀τ ∈ [t, t + T p ],
T p and T c are the prediction and control horizon with T c ≤ T p . ϕ x (•) denotes the new value of the state x(•) computed by the closed loop equation φx using the new input value ϕ u found by the optimization problem for each instant over the moving finite horizon T c (see figure 2).

The cost function, equation ( 15), is chosen based on the process desired performance and the first choice for the cost function is often the quadratic function. Positive weighting matrices (Ω 1 and Ω 2 ) can also be included in the cost function:

Table 1. Input concentrations of the sMBR and physical model parameters

Compound Units Value S I [gCOD m -3 ] 30.00 S S [gCOD m -3 ] 0.77 X I [gCOD m -3 ] 3554.43 X S [gCOD m -3 ] 59.35 X BH [gCOD m -3 ] 3572.44 X BA [gCOD m -3 ] 311.33 X P [gCOD m -3 ] 2373.11 S O [g m -3 ] 2.19 S N O [gN m -3 ] 11.54 S N H [gN m -3 ] 0.33 S N D [gN m -3 ] 0.63 X N D [gN m -3 ] 4.40 V [m 3 ] 0.19 A [m 2 ] 0.93 rsc [m kg -1 ] 1.0e15 µ [P a s] 0.001 Rm [m -1 ] 1e11 γ b [kg/(m 3 s)] 2.5e -5 β b [-] 3.5e -4 α [-] 0.5 dp [m] 1e -6 C d [-] 1 K 1 [-] 4e -6 ρs [kg/m 3 ] 1000 g [m/s 2 ] 9.81 rp [m -2] 1.4e14 F (x, u) = (x-x ref ) T Ω 1 (x-x ref )+(u-u ref ) T Ω 2 (u-u ref ). (15)
where x ref and u ref are the desired reference of a state and an input, respectively.

SIMULATION RESULTS

Both models, the descriptive model from [START_REF] Mannina | An integrated model for biological and physical process simulation in membrane bioreactor (MBRs)[END_REF] and the simplified model, equation ( 6), are implemented in the Matlab environment. The biological process inputs and the descriptive model parameters can be found in Table 1.

First, the parameters of the simplified model are identified based on some prior experiments (see [START_REF] Pimentel | Design, analysis and validation of a simple dynamic model of a submerged membrane bioreactor[END_REF] for the details of the parameter estimation procedure, which is based on a time-scale separation). The identified parameters are presented in Table 2. Following the model calibration, the NMPC methodology is applied using the simplified model ( 6) as a predictor. The cost function is defined as: where the effluent multiplied by the square of the time is minimized in order to significantly reduce the irreversible resistance and at the same time maintain the transmembrane pressure at desired setpoint T M P * and Ω i are the scaling factors. Note that the irreversible resistance, represented by the equations ( 2) and (3a), depend on the permeate flow and time that should be minimized to enlarge the periods between chemical cleaning. The following constraints are added: (i) Q out ≥ 0 for the physical range of the permeate pump and (ii) J air ≥ 0 for the cross-flow range. The methodology is applied using the Matlab code presented by [START_REF] Grüne | Nonlinear Model Predictive Control: Theory and Algorithms[END_REF].

F (x, u) = ( Ω 1 Q out A t 2 ) 2 + (Ω 2 (T M P -T M P * )) 2 (16)
The results presented in Figure 3 are obtained assuming that all the state variables are measured. The NMPC uses a sampling time of one day, while a prediction horizon T p = 3 days and a control interval of T c = 1 day.

The first plot in Figure 3, represents T M P computed with the total resistance, equation (1). Note that, in blue, it is possible to see the influence of the reversible resistance and, in green, the influence of the irreversible resistance on the T M P value. The setpoint T M P * is represented by a red line, which is set to 100 mbar. To maintain the desired setpoint, the controller increases the air cross-flow (Q a /J air ) and, at the same time, decreases the permeate flow (Q out ). These input values are presented in the last two plots of Figure 3. The decay of sludge cake mass, observed in the second plot, shows that the sludge cake resistance (R rev ) is much more important than the irreversible resistance at the beginning of the process, matching the observations of [START_REF] Mannina | An integrated model for biological and physical process simulation in membrane bioreactor (MBRs)[END_REF]. It is important to highlight that the controller actions maintain the desired setpoint even though some nonlinearities are not modeled in the simplified model. In order to test the robustness of the control scheme, a simulation study is achieved by considering an additive white noise with a standard deviation of 15% on X in , so as to emulate the variation in the daily inflow concentrations, represented in Figure 4. The results show that the NMPC is robust enough to maintaining in the desired operating range until a chemical cleaning is required. Additionally, to emulate errors in the actuators, it is added 5% random noise in the input Q out . Even though this input affects all the Fig. 4. The NMPC acting with 15% random noise added to X in .

states and the irreversible and reversible layers, the process maintains its closed-loop performance.

The behavior of the controller when the irreversible resistance becomes more important than the sludge cake is shown in Figure 5. This occurs around day 30 where the sludge cake mass is extremely small, see second plot in Figure 5. To keep the T M P at the selected setpoint the controller acts vigorously by changing the air cross-flow, as seen in the third plot. This action does not results in a decay of the T M P value, and around day 34, the controller significantly changes the permeate flow to maintain the setpoint. This action is explained recalling the equation (2), that links permeate flow to the irreversible resistant, when the fractions of the R sc and R dc are extremely small. When this state is reached, a chemical cleaning is required, or a larger trans-membrane pressure has to be selected. Note that the chemical cleaning, when applied, deeply changes the properties of the membrane filtration, resulting in the need for a new identification of the model parameters before the NMPC framework is applied to the process.

The simulation in Figure 3 shows the importance of airscour to maintain the reversible sludge layer at an admissible level until the irreversible layer becomes more important. Comparing the process without and with noise, it is apparent that more air-cross flow is needed in the second case. Figure 5 shows the inefficiency of air-cross flow to detach the irreversible cake, resulting into a waste of energy used by the air-blowers that are the key factor in the energy consumption in a sMBR plant. Furthermore, the process performance is controlled using air-scouring and permeate flows. These two inputs profiles are good indicators about the process health. Analyzing these profiles one can understand the process. In real applications, alarms could be set when these inputs reach a certain values or if its derivatives have large value, predicting the maintenance of the membrane process.

CONCLUSION AND FUTURE WORK

A NMPC has been applied to a submerged membrane bioreator process intended to stabilize the trans-membrane pressure at a desired value. The irreversible and sludge Fig. 5. Blue is the sludge cake influence on T M P and green is the irreversible resistance on T M P . cake layer resistances have been used in order to optimize the process. The results show that after a certain amount of time, the process cannot be stabilized anymore due to the irreversible resistance. The process input values of Q out and J air lead us to believe that the setpoint should be changed or the chemical cleaning procedure should be carried out. Hence, it is possible to control the process and at the same time monitor the input variations that try to predict the next cleaning procedure needed. Future research will focus on the recirculating (Q r ) and waste (Q w ) flows that will be incorporated in the NMPC framework for a better control of biological degradation in the process while modifying the hydraulic and soluble retention time factors.

  Figure (1). The following manipulated variables are usually available: the permeate flow Q out [m 3 /d], the waste flow Q w [m 3 /d], the air cross-flow J air [m 3 /m 2 d], as well as the input concentrations of biomass X in and substrate S in and the trans-membrane pressure T M P [mbar].

Fig. 1 .

 1 Fig. 1. Schematic representation of a submerged membrane bioreactor (sMBR).

Fig. 3 .

 3 Fig. 3. The NMPC acting on the descriptive process.

Table 2 .

 2 Parameter Calibration

	Parameters	Values
	K air [g]	9.3367 × 10 6
	ρ [m/g]	2.0229 × 10 10
	m 0 [g]	4.4869
	γ [day]	0.0194
	Y [-]	1.8558
	µ S,max	0.0781


Fanny Delebecque

AND Stefan Le Coz

Rada M Weish Äupl

MULTI-SPEED SOLITARY WAVES OF NONLINEAR SCHR ÖDINGER SYSTEMS: THEORETICAL AND NUMERICAL ANALYSIS

Keywords: . 2010 Mathematics Subject Classification. 35Q55(35C08, 35Q51, 37K40) solitons, solitary waves, nonlinear Schrödinger systems

We consider a system of coupled nonlinear Schrödinger equations in one space dimension. First, we prove the existence of multi-speed solitary waves, i.e solutions to the system with each component behaving at large times as a solitary wave. Then, we investigate numerically the interaction of two solitary waves supported each on one component. Among the possible outcomes, we find elastic and inelastic interactions, collision with mass extraction and reflexion.

Introduction

We consider the following nonlinear Schrödinger system:

i∂ t u 1 + ∂ xx u 1 + µ 1 |u 1 | 2 u 1 + β|u 2 | 2 u 1 = 0, i∂ t u 2 + ∂ xx u 2 + µ 2 |u 2 | 2 u 2 + β|u 1 | 2 u 2 = 0, ( NLS
)
where for j = 1, 2 we have u j : R × R → C, µ j > 0, and β ∈ R \ {0}.

When µ 1 = µ 2 = β, system (NLS), also called Manakov system has been introduced by Manakov (see [START_REF] Manakov | On the theory of two-dimensional stationary self-focusing of electromagnetic waves[END_REF] for example) as an asymptotic model for the propagation of electric fields in waveguides. In this particular case, it is to be noticed that the usual roles of x and t are inverted to study the evolution of the electrical field along the propagation axis.

It has also been used later on to model the evolution of light in optical fiber links. One of the main limiting effects of transmission in optical fiber links is due to the polarization mode dispersion (PMD). It can be explained by the birefringence effect, i.e the fact that the electric field is a vector field and that the refraction index of the medium depends on the polarization state (see e.g [START_REF] Ablowitz | Discrete and continuous nonlinear Schrödinger systems[END_REF][START_REF]Nonlinear fiber optics[END_REF]). The evolution of two polarized modes of an electrical field in a birefringent optical fiber link can indeed be modeled by (NLS) in the case where µ 1 = µ 2 and β measures the strength of the cross phase modulation which depends of the fiber (see [START_REF] Manakov | On the theory of two-dimensional stationary self-focusing of electromagnetic waves[END_REF]). Randomly varying birefringence is studied adding random coefficients in both nonlinearity and coupling terms of (NLS) (see for example [START_REF] Garnier | Effective pulse dynamics in optical fibers with polarization mode dispersion[END_REF])

In higher dimensions, systems of nonlinear coupled Schrödinger equations appears in various physical situations such as the modeling of the interaction of two Bose-Einstein condensates in different spin states.

Systems of type (NLS) have also been studied from the mathematical point of view. When µ 1 = µ 2 = β, in dimension 1, the system (NLS) has the particularity to be completely integrable. Hence explicit calculations of solutions are possible and one can exhibit a variety of "truly" nonlinear solutions like solitons, multi-solitons or breathers (see e.g. the book [START_REF] Ablowitz | Discrete and continuous nonlinear Schrödinger systems[END_REF]). The integrability property is however not robust, and the slightest change in the parameters µ 1 , µ 2 and β destroys it. Many works (see, among many others, [START_REF] Ambrosetti | Standing waves of some coupled nonlinear Schrödinger equations[END_REF][START_REF] Bartsch | Note on ground states of nonlinear Schrödinger systems[END_REF][START_REF] De Figueiredo | Solitary waves for some nonlinear Schrödinger systems[END_REF][START_REF] Sirakov | Least energy solitary waves for a system of nonlinear Schrödinger equations in R n[END_REF]) have been devoted to the study of the stationary version of (NLS)

∂ xx φ 1 + µ 1 |φ 1 | 2 φ 1 + β|φ 2 | 2 φ 1 = ω 1 φ 1 , ∂ xx φ 2 + µ 2 |φ 2 | 2 φ 2 + β|φ 1 | 2 φ 2 = ω 2 φ 2 , (1)
that one obtains when looking for standing waves solutions (u 1 , u 2 )(t, x) = (e iω1t φ 1 (x), e iω2t φ 2 (x)).

When standing waves exist, it is natural to study their stability and again many works have been devoted to this problem (see, again among many others, [START_REF] Hajaiej | Orbital stability of standing waves of some ℓ-coupled nonlinear Schrödinger equations[END_REF][START_REF] Maia | Orbital stability property for coupled nonlinear Schrödinger equations[END_REF][START_REF] Montefusco | Energy convexity estimates for non-degenerate ground states of nonlinear 1D Schrödinger systems[END_REF][START_REF] Ohta | Stability of solitary waves for coupled nonlinear Schrödinger equations[END_REF]). Existence and stability of standing waves are often proved using variational techniques. The analysis of (1) through variational techniques is very subtle and the introduction of new ideas is necessary to understand the full picture (see [START_REF] Nguyen | Existence and stability of a two-parameter family of solitary waves for a 2-coupled nonlinear Schrödinger system[END_REF]). Note that (NLS) is Galilean-invariant. Hence a Galilean transform modifies a standing wave into a solitary wave traveling at some non-zero speed. Our goal in this paper is to provide a new point of view on the study of this system. We aim at understanding better the behavior in large times of solutions starting at initial time as two scalar solitary waves carried by the two different components. We will use a mixture of theoretical and numerical tools, a combination seldom seen when dealing with this kind of problems.

First, we propose to push further a study initiated in [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF] on the multi-speed solitary waves of system (NLS) and followed up for a different nonlinearity in [START_REF] Wang | Multi-speed solitary wave solutions for a coherently coupled nonlinear schrödinger system[END_REF]. A multi-speed solitary wave is a solution of (NLS) which behaves at large time as two solitary waves. Here and as in [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF], we restrict ourselves to the case where the composing solitary waves are each carried on only one component of the system. In other words, taken independently, each component behaves as a scalar solitary wave at large time. Our first aim is to remove the high speed assumption under which the main result in [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF] was proved. We therefore consider the system (NLS) in dimension 1 and benefit from the fact that scalar solitary waves are in that case orbitally stable (see e.g [START_REF] Cazenave | Orbital stability of standing waves for some nonlinear Schrödinger equations[END_REF]).

Our next aim is to investigate further the properties of multi-speed solitary waves solutions when they are crossing at positive time. Our theoretical result (Theorem 1) indeed only guarantees existence of multi-speed solitary wave solutions to (NLS) when no interaction can occur at large time between the composing waves. But what happens when two solitary waves carried by different components collide? Due to the possible complexity of the phenomenon and the lack of appropriate theoretical tools to study it, we proceed the following numerical experiment. We take as initial data solitons on each components, both away from 0 but facing each other for the direction of propagation. Among the possible outcomes, we find elastic and inelastic interactions, interaction with mass exctraction, and reflexion.

1.1. The theoretical result. Before stating our main theoretical result, let us give a few preliminaries.

Let Q ω ∈ H 1 (R) be the unique positive radial ground state solution to

-∂ xx Q ω + ωQ ω -|Q ω | 2 Q ω = 0, Q ω > 0, Q ω ∈ H 1 rad (R). (2)
From simple calculations we note that the following scaling occurs:

Q ω (x) := √ ω Q 1 ( √ ωx). (3)
For j = 1, 2, consider ω j > 0, γ j ∈ R, x j , v j ∈ R and define

R j (t, x) = e i(ωj t-1 4 |vj | 2 t+ 1 2 vj •x+γj) 1 µ j Q ωj (x -v j t -x j ). (4)
The function R j is a solitary wave solution to

i∂ t u + ∂ xx u + µ j |u| 2 u = 0. (5)
In this paper, we want to investigate the existence of solutions to (NLS) where each component behaves like a solitary wave R j solution to the scalar equation [START_REF] Bao | Ground states and dynamics of multicomponent Bose-Einstein condensates[END_REF]. Our main theoretical result is the following. Theorem 1. Let µ 1 , µ 2 > 0 and β ∈ R \ {0}. For j = 1, 2, take v j , x j , γ j ∈ R, ω j > 0 and consider the ground state profile Q ωj solution to (2) and the soliton R j defined in (4). Then, there exist C > 0, T 0 > 0 and u 1 , u 2 solution to (NLS) on the time interval [T 0 , +∞) such that for all t ∈ [T 0 , +∞), we have

u 1 (t) u 2 (t) - R 1 (t) R 2 (t) H 1 ×H 1 Ce -√ ω * v * t ,
where [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF]Theorem 1], the main differences are the following. Our result is valid for any speeds, whereas the one in [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF] required a high speed assumption. We restrict ourselves to dimension 1 to have stable solitons (in [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF], any dimension was allowed). The overall proof strategy is similar, but in our case we need to perform several technical refinements which include in particular working with localized momenta and modulated waves.

ω * = 1 2304 min{ω 1 , ω 2 } and v * = |v 1 -v 2 |. Remark 1. Compare to
In addition, we are introducing the technical artefact consisting into introducing arbitrary constants in the definition [START_REF] Garnier | Effective pulse dynamics in optical fibers with polarization mode dispersion[END_REF] of the global action. This is a new feature for this type of analysis, which is quite surprising as usually such a flexibility is not allowed by the algebra of the problem.

The scheme of the proof is inspired by the one developed for the study of multisolitons in scalar nonlinear Schrödinger equations in [START_REF] Côte | High-speed excited multi-solitons in nonlinear Schrödinger equations[END_REF][START_REF] Côte | Construction of multi-soliton solutions for the L 2supercritical gKdV and NLS equations[END_REF][START_REF] Martel | Multi solitary waves for nonlinear Schrödinger equations[END_REF][START_REF] Merle | Construction of solutions with exactly k blow-up points for the Schrödinger equation with critical nonlinearity[END_REF] (see also [START_REF] Bellazzini | Multi-solitary waves for the nonlinear Klein-Gordon equation[END_REF] for a similar approach applied to Klein-Gordon equations). It consists in solving (NLS) backward in time, taking as final data a couple of solitary waves R 1 (T n ), R 2 (T n ) , for an increasing sequence of times T n → +∞. Thus we get a sequence u n 1 , u n 2 of solutions to (NLS) on a time interval (-∞,

T n ] such that u n 1 (T n ), u n 2 (T n ) = R 1 (T n ), R 2 (T n ) .
We then have to prove the existence of a time T 0 , independent of n such that for n large enough,

u n 1 , u n 2 is close to R 1 , R 2 on [T 0 , T n ].
The key tools at hand to prove Theorem 1 are

• uniform in n estimates ∀t ∈ [T 0 , T n ], u n 1 (t) u n 2 (t) - R 1 (t) R 2 (t) H 1 ×H 1 Ce -√ ω * v * t ,
• a compactness argument that gives the existence of

u 0 1 , u 0 2 ∈ H 1 (R) × H 1 (R) such that u n 1 , u n 2 converges strongly in H s (R) (s ∈ [0, 1)) towards u 0 1 , u 0 2 . Remark 2.
The method used to obtain Theorem 1 is a powerful tool to obtain sharp existence results for multi-solitons composed of ground states. Another approach relying on a fixed point argument has been developed for nonlinear Schrödinger equations in [START_REF] Le Coz | Fast-moving finite and infinite trains of solitons for nonlinear Schrödinger equations[END_REF][START_REF] Coz | Infinite soliton and kink-soliton trains for nonlinear Schrödinger equations[END_REF]. This approach is very flexible and allows to prove existence of solutions more complicated than the ones in Theorem 1 like infinite trains of solitons or multi-kinks. The main drawback is that it always requires a large speed assumption.

We also have tested numerically if the multi-speed solitary wave configuration was stable provided the starting waves are well-ordered and well-separated. In other words, we took the interaction to be small at the origin and the composing waves going away from each other. With such a well-prepared initial configuration, we remain close to a similar configuration in large time. This suggests that the multi-speed solitary waves are stable (no matter the coupling parameter). Note that this is expected due to the fact that each wave taken individually is stable. We have however no theoretical mean to verify this conjecture. Similar difficulties arise in the analysis of the stability for multi-solitons in nonlinear scalar Schrödinger equations (see e.g. [START_REF] Martel | Stability in H 1 of the sum of K solitary waves for some nonlinear Schrödinger equations[END_REF]).

1.2. The numerical experiments. We solve the system (NLS) in one dimension by adapting the time-splitting spectral method described in [START_REF] Bao | Numerical solution of the Gross-Pitaevskii equation for Bose-Einstein condensates[END_REF]. This method is unconditionally stable, time reversible, of spectral-order accuracy in space and second-order accuracy in time, and it conserves the discrete total mass [START_REF] Bao | Ground states and dynamics of multicomponent Bose-Einstein condensates[END_REF]. One can refer to [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] for other possible schemes and their properties.

We will also compute the real valued ground state (minimizer of the energy on fixed L 2 mass constraints) of the system (1) using a normalized gradient flow approach. This will be used to make the comparison between the outcome of the interaction between two solitary waves and a solitary wave with profile (φ 1 , φ 2 ).

As already mentioned, the experiment consists in taking as initial data the initial data of two solitary waves facing each other, each on one component.

We considered four cases, the first one being the integrable case, where we expect the solitons after the interaction to move with the same velocity and amplitude. Apart when µ 1 = µ 2 = β, the system is not integrable, hence we do not expect pure elastic interaction between solitons. However, there are still regimes where we expect the outcome of interaction between solitons to be also a multi-speeds solitary wave, different from the input at two level: first, there are modifications in the speeds and amplitudes of the composing solitons. Second, there is a loss of a bit of energy, mass and momentum into a small dispersive remainder. In certain cases, we have been able to identify the profiles of the outcome of the interactions as ground states of the stationary system [START_REF] Ablowitz | Discrete and continuous nonlinear Schrödinger systems[END_REF].

The rest of this paper is organized as follows. In Section 2, we prove Theorem 1 assuming uniform estimates. In Section 3, we prove the uniform estimates. The numerical methods are described in Section 4, and the numerical experiments are presented in Section 5. Appendix A contains the proof of a modulation result.

Existence of multi-speed solitary waves

This short section is devoted to the proof of Theorem 1, assuming uniform estimates proved in the next section.

In this section and in the next one, we assume that µ 1 , µ 2 > 0 and β ∈ \{0} are fixed constants and that we are given for j = 1, 2 soliton parameters ω j , v j , x j , γ j ∈ R. Denote by Q ωj and R j the corresponding profile and soliton.

We make the assumption that

0 < v 1 = -v 2 . (6)
Since (NLS) is Galilean invariant, this assumption can be done without loss of generality. This will simplify calculations later on. Note that it follows from classical arguments (see [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF]) that the Cauchy problem for (NLS) is globally well-posed in the energy space

H 1 (R) × H 1 (R) and also in L 2 (R) × L 2 (R). In particular, for any initial data (u 0 1 , u 0 2 ) ∈ H 1 (R) × H 1 (R) there exists a unique global solution (u 1 , u 2 ) of (NLS) in C(R, H 1 (R) × H 1 (R)) ∩ C 1 (R, H -1 (R) × H -1 (R)).
Let (T n ) be an increasing sequence of times such that T n → +∞ as n → +∞. Let (u n 1 , u n 2 ) be the sequence of solutions to (NLS) defined by solving (NLS) backward on (-∞,

T n ] with final data (u n 1 , u n 2 )(T n ) = R 1 , R 2 (T n ).
The proof of Theorem 1 then relies on the following two ingredients.

First, we have uniform estimates on the distance between the sequence u n 1 , u n 2 and the multi-speed solitary wave profile R 1 , R 2 .

Proposition 3 (Uniform estimates). There exist T 0 > 0, n 0 ∈ N such that, for all n n 0 and for all t ∈ [T 0 , T n ] we have

u n 1 u n 2 (t) - R 1 R 2 (t) H 1 ×H 1 e -√ ω * v * t ,
where, as in Theorem 1,

ω * = 1 2304 min{ω 1 , ω 2 } and v * = |v 1 -v 2 |.
The proof of Proposition 3 is rather involved and we postpone it to Section 3. The next ingredient is a compactness result on the initial data u n 1 , u n 2 (T 0 ). This result was already present in this form in [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF] and we recall it without proof.

Proposition 4 (Compactness). There exists u

0 1 , u 0 2 ∈ H 1 (R) × H 1 (R) such that, up to a subsequence, u n 1 , u n 2 (T 0 ) converges strongly towards u 0 1 , u 0 2 in H s (R) × H s (R) for all s ∈ [0, 1).
With these two ingredients in hand, we can now conclude the proof of Theorem 1.

Proof of Theorem 1. Let (u 1 , u 2 ) be the solution on R of the Cauchy problem (NLS) with initial data

(u 0 1 , u 0 2 ) at t = T 0 . By H 1 (R) × H 1 (R) boundedness and local well-posedness of the Cauchy problem in H s (R) × H s (R) for all s ∈ [0, 1), we have weak convergence in H 1 (R) × H 1 (R) of (u n 1 , u n 2 )(t) towards (u 1 , u 2 )(t) for any t ∈ R. Combined with the uniform estimates of Proposition 3, this implies for all t ∈ [T 0 , +∞) that u 1 u 2 (t) - R 1 R 2 (t) H 1 ×H 1 lim inf n→+∞ u n 1 u n 2 (t) - R 1 R 2 (t) H 1 ×H 1 Ce -√ ω * v * t .
This concludes the proof of Theorem 1.

Uniform estimates

This section is devoted to the proof of Proposition 3. In all this section, T n and u n 1 , u n 2 are given as in the beginning of Section 2. 3.1. The bootstrap argument. We first reduce the proof of Proposition 3 to the proof of the following bootstrap result.

Proposition 5 (Bootstrap argument).

There exist T 0 > 0 and n 0 ∈ N such that for all n n 0 and for any t 0 ∈ [T 0 , T n ] the following property is satisfied. If for all t ∈ [t 0 , T n ] we have

u n 1 u n 2 (t) - R 1 R 2 (t) H 1 ×H 1 e -√ ω * v * t , (7)
then for all t ∈ [t 0 , T n ] we have

u n 1 u n 2 (t) - R 1 R 2 (t) H 1 ×H 1 1 2 e -√ ω * v * t . (8)
The proof of Proposition 5 will occupy us for most of the rest of this section. We divided it into several steps. We first perform a geometrical decomposition of the sequence (u n 1 , u n 2 ) onto the manifold of multi-speed solitary waves in order to obtain orthogonality conditions. We then introduce an action-like functional, which turns out to be coercive due to our orthogonality conditions. This functional is not a conserved quantity, but since it is made with localized conservations laws it is almost conserved. Using that property and a control on the geometrical modulation parameters, we are able to conclude the proof of Proposition 5.

Before going on with the details of the proof of Proposition 5, let us show how it implies Proposition 3.

Proof of Proposition 3. Since we have (u n 1 , u n 2 )(T n ) = (R n 1 , R n 2 )(T n
) at the final time T n , by continuity there exists a minimal time t 0 such that for all t ∈ [t 0 , T n ] we have

u n 1 u n 2 (t) - R 1 R 2 (t) H 1 ×H 1 e -√ ω * v * t . (9)
We prove that t 0 = T 0 by contradiction. Assume that t 0 > T 0 . By Proposition 5, for all t ∈ [t 0 , T n ] we have

u n 1 u n 2 (t) - R 1 R 2 (t) H 1 ×H 1 1 2 e -√ ω * v * t .
Therefore by continuity there exists t 00 < t 0 such that on [t 00 , T n ] estimate ( 9) is satisfied. This however contradicts the minimality of t 0 . Hence t 0 = T 0 and this concludes the proof.

For the rest of Section 3, T 0 > 0 and n 0 ∈ N will be large enough fixed numbers, and we assume the existence of t 0 T 0 such that for all t ∈ [t 0 , T n ] the bootstrap assumption ( 7) is verified, i.e. we have

u n 1 u n 2 (t) - R 1 R 2 (t) H 1 ×H 1 e -√ ω * v * t . (10)
Our final goal is now to prove that in fact ( 8) holds for all t ∈ [t 0 , T n ].

Modulation.

Let us start with a decomposition lemma for our sequence of approximated multi-speed solitary waves.

Lemma 6 (Modulation).

There exist C > 0 and C 1 functions

ωj : [t 0 , T n ] → (0, +∞), xj : [t 0 , T n ] → R, γj : [t 0 , T n ] → R, j = 1, 2,
such that if for j = 1, 2 we denote by Rj the modulated wave

Rj (t, x) = e i( 1 2 vj •x+γj(t)) 1 √ µ j Q ωj(t) (x -xj (t)), (11)
then for all t ∈ [t 0 , T n ] the functions defined by

ε 1 ε 2 (t) = u n 1 u n 2 (t) -R1 R2 (t)
satisfy for j = 1, 2 and for all t ∈ [t 0 , T n ] the orthogonality conditions

ε j (t), Rj (t) 2 = ε j (t), i Rj (t) 2 = ε j (t), ∂ x Rj (t) 2 = 0. ( 12
)
Moreover, for all t ∈ [t 0 , T n ], we have

2 j=1   |∂ t ωj (t)| 2 + |∂ t xj (t) -v j | 2 + ∂ t γj (t) + v 2 j 4 -ωj (t) 2   C ε 1 ε 2 (t) 2 H 1 ×H 1 + Ce -3 √ ω * v * t . ( 13
)
Remark 7. It is to be noticed that estimate [START_REF] Côte | Construction of multi-soliton solutions for the L 2supercritical gKdV and NLS equations[END_REF] clearly implies that, for T 0 large enough, for j = 1, 2, and for all t ∈ [t 0 , T n ], we have

xj (t) v * 2 √ 2 t > 2L and ωj (t) 1152ω * .
Moreover, a better estimate can be obtained for ωj and will be stated later on in Lemma 12.

This type of modulation result is classical in the literature dealing with solitary waves of nonlinear dispersive equations (see e.g. the fundamental paper of Weinstein [START_REF] Weinstein | Modulational stability of ground states of nonlinear Schrödinger equations[END_REF] for an early version or [START_REF] Martel | Stability in H 1 of the sum of K solitary waves for some nonlinear Schrödinger equations[END_REF] for a recent approach). Its proof consists essentially in the application of the implicit function theorem combined with the use of the evolution equation to find equation ( 13) for the evolution of the modulation parameters. We refer to the appendix for the details of the proof.

Energy estimates and coercivity.

In this subsection, we analyze the different quantities that are conserved or almost-conserved in our coupled-vectorial problem. Remember that in the case of the scalar equation ( 5) the energy, mass and momentum, defined as follows, are conserved along the flow of ( 5):

E(u, µ j ) := 1 2 ∂ x u 2 L 2 - µ j 4 u 4 L 4 , M (u) := 1 2 u 2 L 2 , P (u) := 1 2 Im R u∂ x udx.
The solution Q ω we chose of equation ( 2) is known to be the unique positive radial ground state of the action S := E(•, 1)+ωM . Consequently, each soliton R j defined by ( 4) is a critical point of the scalar functional S j defined by

S j := E(•, µ j ) + ω j + v 2 j 4 M + v j P. ( 14
)
Coercivity properties of linearizations of S j -like functionals are the key tool of the analysis of multi-solitons interaction (see for example [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF][START_REF] Martel | Multi solitary waves for nonlinear Schrödinger equations[END_REF]).

In the vectorial case we are interested in here, the coupled system (NLS) admits its own conservation laws. In particular, the mass of each component is preserved, as in the scalar case. However, the coupling does not preserve conservation of the scalar energy and momentum for each component and we only have conservation of the total energy (made of individual energies plus a coupling term) and total momentum (sum of the scalar momenta). More precisely, total energy, total momentum, and scalar masses of the whole system (defined as follows) are conserved quantities for the flow of system (NLS):

E u 1 u 2 := E(u 1 , µ 1 ) + E(u 2 , µ 2 ) - β 2 R |u 1 | 2 |u 2 | 2 dx, P u 1 u 2 := P (u 1 ) + P (u 2 ), M j u 1 u 2 := M (u j ), j = 1, 2.
In order to use the conservation of the momentum as in the scalar case, we here need to localize the momentum of each soliton, as was done in [START_REF] Côte | High-speed excited multi-solitons in nonlinear Schrödinger equations[END_REF][START_REF] Côte | Construction of multi-soliton solutions for the L 2supercritical gKdV and NLS equations[END_REF] for the scalar mass and momentum. Note that this was not needed for the analysis in [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF]. Let us define the cut-off functions

χ 1 L (x) = χ x L , χ 2 L = 1 -χ 1 L ,
where L > 0 is arbitrary but fixed and χ is a C 3 function such that

0 χ 1 on R, χ(x) = 0 for x -1, χ(x) = 1 for x > 1, χ ′ 0 on R,
and satisfies, for some positive constant C and for all x ∈ R the estimates

(χ ′ (x)) 2 Cχ(x), (χ ′′ (x)) 2 Cχ ′ (x).
Localized momenta P j loc are defined by:

P j loc u 1 u 2 = 1 2 Im R u 1 ∂ x u 1 + u 2 ∂ x u 2 χ j L dx, j = 1, 2.
Remark that P = P 1 loc + P 2 loc . Note that since we are assuming (6) the momenta above defined are localized around each composing solitary wave of the profile. The advantage of having made assumption ( 6) is that the cut-off does not depend on time. This will simplify our next calculations.

In the sequel, we are interested in the following global action:

S u 1 u 2 = E u 1 u 2 + j=1,2 ωj (t) + v 2 j 4 M j u 1 u 2 + j=1,2 v j P j loc u 1 u 2 + C 1 (β, v 1 )M 2 u 1 u 2 + C 2 (β, v 2 )M 1 u 1 u 2 , (15)
where C j (β, v j ) are positive constants depending only on β and v j and whose exact values will be decided later on. Note that the action implicitly depends on t via ωj . It is to be noted that, in this work, we have the freedom to add these two coupled-mass terms that do not appear in the usual definition of S-like functionals. This is a key point in our analysis.

Let us now state in the following lemma several estimates related to the localization of R1 and R2 and which will be of great use in the sequel. Lemma 8. For j = 1, 2, if T 0 is large enough, then for all t ∈ [t 0 , T n ] and for all x ∈ R we have

| Rj (t, x)| + |∂ x Rj (t, x)| χ 3-j L (x) C(1 + |v j |)e -3 √ ω * v * t e -√ ω * |x| , (16)
k=1,2

| Rk (t, x)| + |∂ x Rk (t, x)| C(1 + |v 1 | + |v 2 |) 2 e -3 √ ω * v * t e -√ ω * |x| . ( 17
)
Lemma 8 follows from the support properties of the cut-off function and the exponential localization of the solitons profiles. Indeed, recall that in fact, the profile

Q = Q 1 is explicitly known Q(x) = 2 sech(x),
and it follows that Q and its derivatives are exponentially decaying, i.e. for any η < 1 we have

|Q xx | + |Q x | + |Q| C η e -η|x| .
Proof of Lemma 8. We prove only [START_REF] Hajaiej | Orbital stability of standing waves of some ℓ-coupled nonlinear Schrödinger equations[END_REF], the proof of (17) following from similar (simpler) arguments. For simplicity in notation, assume j = 1, the case j = 2 being perfectly symmetric. Due to the exponential decay of the soliton profiles, we have

| R1 (t, x)| Ce -3 4 √ ω1|x-x1|
The cut-off function χ 2 L is supported on (-∞, L], and since for T 0 large enough x1 > 2L, for x ∈ (-∞, L] we have (we recall here Remark 7)

|x -x1 | |x|, and |x -x1 | 1 2 |x 1 |.
As a consequence, we have

| R1 (t, x)|χ 2 L (x) Ce -1 4 √ ω1|x1| e -1 4 √ ω1|x| .
In addition, as noticed in Remark

S u 1 u 2 = S R1 R2 + H ε 1 ε 2 + O(e -3 √ ω * v * t ), (18)
H ε 1 ε 2 = H free ε 1 ε 2 + H coupled ε 1 ε 2 , with H free ε 1 ε 2 = j=1,2 1 2 ∂ x ε j 2 L 2 + 1 2 ωj + v 2 j 4 ε j 2 L 2 -µ j R |ε j | 2 | Rj | 2 dx - µ j 2 Re R ε 2 j R2 j dx + 1 2 v j Im R ε j ∂ x ε j χ j L dx
and

H coupled ε 1 ε 2 = C 1 (β, v 1 ) ε 2 2 L 2 + C 2 (β, v 2 ) ε 1 2 L 2 - β 2 R |ε 1 | 2 | R2 | 2 + |ε 2 | 2 | R1 | 2 dx + 1 2 v 1 Im R ε 2 ∂ x ε 2 χ 1 L dx + 1 2 v 2 Im R ε 1 ∂ x ε 1 χ 2 L dx.
Proof of Lemma 9. First note that, according to the definition (15) of S,

S u 1 u 2 = j=1,2 Sj (u j ) + v j • P j loc u 1 u 2 -P (u i ) + C 1 (β, v 1 )M (u 1 ) + C 2 (β, v 2 )M (u 2 ) - β 2 R |u 1 | 2 |u 2 | 2 dx, ( 19
)
where Sj denote the same functional as S j (see the definition ( 14)) with ωj instead of ω j . For j = 1, 2, let us expand u j (t) = Rj (t) + ε j (t) in the expression of Sj . A simple computation leads to:

Sj ( Rj + ε j ) = Sj ( Rj ) + S′ j ( Rj )ε j + S′′ j ( Rj )ε j , ε j + O( ε j 3 H 1 )
Now, as Rj is a critical point of the scalar functional Sj , we have S′ j ( Rj ) = 0, and thus Sj ( Rj

+ ε j ) = Sj ( Rj ) + S′′ j ( Rj )ε j , ε j + O( ε j 3 H 1 ), (20)
where

S′′ j ( Rj )ε j , ε j = 1 2 ∂ x ε j 2 L 2 + 1 2 ωj + v 2 j 4 ε j 2 L 2 + 1 2 v j Im R ε j ∂ x ε j dx - µ j 2 Re R ε 2 j R2 j dx -µ j R |ε j | 2 | Rj | 2 dx. ( 21
)
Let us now develop the remaining terms in [START_REF] Coz | Infinite soliton and kink-soliton trains for nonlinear Schrödinger equations[END_REF]. As far as the momentum part is concerned, we write the expansion for j = 1 for simplicity:

P 1 loc R1 + ε 1 R2 + ε 2 -P ( R1 + ε 1 ) = -Im R R1 ∂ x R1 χ 2 L dx + Im R R2 ∂ x R2 χ 1 L dx -Im R R1 ∂ x ε 1 + ε 1 ∂ x R1 χ 2 L dx + Im R R2 ∂ x ε 2 + ε 2 ∂ x R2 χ 1 L dx -Im R ε 1 ∂ x ε 1 χ 2 L dx + Im R ε 2 ∂ x ε 2 χ 1 L dx. ( 22
)
Concerning the β coupling part in [START_REF] Coz | Infinite soliton and kink-soliton trains for nonlinear Schrödinger equations[END_REF], we get:

R | R1 + ε 1 | 2 | R2 + ε 2 | 2 dx = R | R1 | 2 | R2 | 2 dx + 2Re R | R1 | 2 R2 ε 2 + | R2 | 2 R1 ε 1 dx + R |ε 1 | 2 | R2 | 2 + |ε 2 | 2 | R1 | 2 + 4Re(ε 1 R1 )Re(ε 2 R2 ) dx + 2 R |ε 2 | 2 Re(ε 1 R1 ) + |ε 1 | 2 Re(ε 2 R2 ) dx + R |ε 1 | 2 |ε 2 | 2 dx. (23)
Finally, the extra-masses terms in [START_REF] Coz | Infinite soliton and kink-soliton trains for nonlinear Schrödinger equations[END_REF] expand into

M ( Rj +ε j ) = 1 2 Rj +ε j 2 L 2 = M ( Rj )+ Rj , ε j 2 +M (ε j ) = M ( Rj )+M (ε j ), ( 24
)
where we have used the orthogonality conditions ( 12) to obtain the last equality. In ( 22) and ( 23), all terms containing a product of solitons or cut-off functions with different indices/exponents are of order O(e -3 √ ω * v * t ) by Lemma 8. The terms containing a degree 3 or higher term in (ε 1 , ε 2 ) are also of order O(e -3 √ ω * v * t ) by the bootstrap assumption [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF]. Therefore, gathering ( 19)-( 20)-( 21)-( 22)-( 23)-( 24) together gives:

S R1 + ε 1 R2 + ε 2 = S R1 R2 + H free ε 1 ε 2 + H coupled ε 1 ε 2 + O(e -3 √ ω * v * t ).
This concludes the proof.

Lemma 10 (Coercivity of H).

There exists λ > 0 such that, for any t 0 ∈ [T 0 , T n ] and for all t ∈ [t 0 , T n ]:

H free ε 1 ε 2 2λ ε 1 2 H 1 + ε 2 2 H 1 , (25)
H coupled ε 1 ε 2 λ ε 1 2 L 2 + ε 2 2 L 2 -λ ∂ x ε 1 2 L 2 + ∂ x ε 2 2 L 2 , (26)
and thus:

H ε 1 ε 2 λ ε 1 ε 2 2 H 1 ×H 1 . ( 27
)
Proof. The proof of ( 25) is classical (see for example [START_REF] Ianni | Multi-speed solitary wave solutions for nonlinear Schrödinger systems[END_REF][START_REF] Martel | Multi solitary waves for nonlinear Schrödinger equations[END_REF]) and we omit it. It remains to prove [START_REF] Montefusco | Energy convexity estimates for non-degenerate ground states of nonlinear 1D Schrödinger systems[END_REF]. It is readily seen that for C j (β, v j ) large enough, we have

C 1 (β, v 1 ) ε 2 2 L 2 + C 2 (β, v 2 ) ε 1 2 L 2 - β 2 R |ε 1 | 2 | R2 | 2 + |ε 2 | 2 | R1 | 2 dx 1 2 C 1 (β, v 1 ) ε 2 2 L 2 + C 2 (β, v 2 ) ε 1 2 L 2 .
As far as the momentum parts of H coupled are concerned, for example for j = 1:

v 1 2 Im R ε 1 ∂ x ε 1 χ 2 L dx - v 1 2 R |ε 1 ∂ x ε 1 χ 2 L |dx - v 1 2 ε 1 2 ∂ x ε 1 2 - v 2 1 4λ ε 1 2 L 2 -λ ∂ x ε 1 2 L 2 .
For C(β, v j ), j = 1, 2 large enough, i.e such that j=1,2

1 2 C(β, v j ) - v 2 j 4λ > λ,
we thus have:

H coupled ε 1 ε 2 λ ε 1 2 L 2 + ε 2 2 L 2 -λ ∂ x ε 1 2 L 2 + ∂ x ε 2 2 L 2 .
Combining estimates ( 25) and ( 26) on H free and H coupled finally gives [START_REF] Nguyen | Existence and stability of a two-parameter family of solitary waves for a 2-coupled nonlinear Schrödinger system[END_REF].

3.4.

Almost-conservation of the localized momentum. In this section, we investigate the conservation of the localized momentum and, inspired by [START_REF] Martel | Multi solitary waves for nonlinear Schrödinger equations[END_REF], we state the following lemma:

Lemma 11. There exists C > 0 (independent of L) such that if L and T 0 are large enough then for all t ∈ [t 0 , T n ] and for j = 1, 2 we have

P j loc u 1 u 2 (t) -P j loc u 1 u 2 (T n ) C L e -2 √ ω * v * t . ( 28
)
Proof of Lemma 11. Let us prove the lemma in the case j = 1, the case j = 2 being perfectly similar. We first compute the time derivative of P 1 loc :

d dt P 1 loc (t) = 1 2 Im R ∂ t u 1 ∂ x u 1 χ x L dx + 1 2 Im R u 1 ∂ x,t u 1 χ x L dx. (29)
Let us call term A and B respectively the first and the second term of the right hand side in equality (29):

A = 1 2 Im R ∂ t u 1 ∂ x u 1 χ x L dx, B = 1 2 Im R u 1 ∂ x ∂ t u 1 χ x L dx.
Using the fact that u(t) = u 1 (t) u 2 (t) is a solution to (NLS), we readily get (formally, but this will be justified when integrating):

∂ t u 1 ∂ x u 1 = i∂ xx u 1 ∂ x u 1 + iµ 1 |u 1 | 2 u 1 ∂ x u 1 + iβ|u 2 | 2 u 1 ∂ x u 1 (30)
∂ x ∂ t u 1 = -i∂ xxx u 1 -iµ 1 ∂ x |u 1 | 2 u 1 -iβ∂ x |u 2 | 2 u 1 . (31)
About term A: Equation ( 30) provides us with the following decomposition of term A:

A = 1 2 Re R ∂ xx u 1 ∂ x u 1 χ x L dx + 1 2 µ 1 Re R |u 1 | 2 u 1 ∂ x u 1 χ x L dx + 1 2 βRe R |u 2 | 2 u 1 ∂ x u 1 χ x L dx.
Integrating by part in each term of A finally gives:

A = - 1 4L R |∂ x u 1 | 2 χ ′ x L dx - 1 8L µ 1 R |u 1 | 4 χ ′ x L dx + 1 2 βRe R |u 2 | 2 u 1 ∂ x u 1 χ x L dx. ( 32
)
About term B: Equation (31) provides us with the following decomposition of term B:

B = - 1 2 Re R u 1 ∂ xxx u 1 χ x L dx - 1 2 µ 1 Re R u 1 ∂ x |u 1 | 2 u 1 χ x L dx - 1 2 βRe R u 1 ∂ x |u 2 | 2 u 1 χ x L dx.
As for the A term, integrating by parts (several times if necessary) in each term finally leads to:

B = - 3 4L R |∂ x u 1 | 2 χ ′ x L dx + 1 4L 3 R |u 1 | 2 χ ′′′ x L dx + 3 8L µ 1 R |u 1 | 4 χ ′ x L dx + 1 2 βRe R |u 2 | 2 u 1 ∂ x u 1 χ x L dx + 1 2L βRe R |u 2 | 2 |u 1 | 2 χ ′ x L dx. (33)
Combining (32) and (33), and using the fact that χ ′ and χ ′′′ are supported on [-L, L] lead to the following estimate:

d dt P 1 loc (t) C L L -L |∂ x u 1 | 2 + |u 1 | 2 + |u 1 | 4 dx+C R |u 2 | 2 |u 1 | 2 + |∂ x u 1 | 2 χ 1 L (x)dx. ( 34
)
Note that for x ∈ [-L, L], T 0 large enough, and using [START_REF] Côte | Construction of multi-soliton solutions for the L 2supercritical gKdV and NLS equations[END_REF] we have (see Lemma 8 for similar arguments)

|∂ x R j (t, x) + R j (t, x)| Ce - √ ωj |x-xj| Ce -1 2 √ ωj |xj| Ce -3 √ ω * v * t .
Note that C here depends on v j . Expanding now u j = Rj + ε j , in (34), using the above estimate and Lemma 8, we obtain

d dt P 1 loc (t) C L ε 1 ε 2 2 H 1 ×H 1 + O(e -3 √ ω * v * t ).
The result follows integrating in time between t and T n and using the bootstrap assumption (10).

3.5. Control of the modulation parameters. We now give an estimate of the variations of ω1 with respect to time. This estimate is better than [START_REF] Côte | Construction of multi-soliton solutions for the L 2supercritical gKdV and NLS equations[END_REF] given by the modulation lemma.

Lemma 12 (Variations of ωj (t)). For j = 1, 2 and for all t ∈ [t 0 , T n ], we have

|ω j (t) -ω j | C ε j (t) 2 L 2 .
Proof. This estimate is due to the choice of the modulation orthogonality condition ε j (t), Rj (t) 2 = 0, and the conservation of the mass of each component:

0 = u j 2 2 -R j 2 2 = Rj 2 2 -R j 2 2 + ε j 2 2 = (ω j -ω j ) ∂ ∂ω |ω=ωj Q ω 2 2 + O(|ω j -ω j | 2 ) + ε j 2 2 .
Since ∂ ∂ω |ω=ωj Q ω 2 2 < 0, this concludes the proof.

Conclusion.

With the elements of the previous subsections in hand, we can now conclude the proof of Proposition 5.

End of the proof of Proposition 5. Recall that we have made the bootstrap assumption [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF] and that our goal is to prove that for all t ∈ [t 0 , T n ] we have in fact the better estimate

u n 1 u n 2 - R 1 R 2 H 1 ×H 1 1 2 e -√ ω * v * t . ( 35
)
Let us first expand

u n 1 (t) u n 2 (t) - R 1 (t) R 2 (t) H 1 ×H 1 u n 1 (t) u n 2 (t) -R1 (t) R2 (t) H 1 ×H 1 + R1 (t) Rn 2 (t) - R 1 (t) R 2 (t) H 1 ×H 1 ε 1 (t) ε 2 (t) H 1 ×H 1 + C j=1,2 |ω j -ω j | + j=1,2 O(|ω j -ω j | 2 ).
By Lemma 12, the part involving |ω jω j | is controlled by the ε-part. Hence to finish the proof it is sufficient to control ε. Using successively the coercivity of H

given in [START_REF] Nguyen | Existence and stability of a two-parameter family of solitary waves for a 2-coupled nonlinear Schrödinger system[END_REF], the expansion (18) of the global action S, the conservation of energy and mass, and the almost conservation of localized momenta [START_REF] Ohta | Stability of solitary waves for coupled nonlinear Schrödinger equations[END_REF], we obtain

ε 1 (t) ε 2 (t) 2 H 1 ×H 1 CH ε 1 (t) ε 2 (t) S u 1 (t) u 2 (t) -S R1 (t) R2 (t) + O(e -3 √ ω * v * t ) C L e -2 √ ω * v * t .
Therefore, choosing L large enough we obtain the required estimate (35). This concludes the proof of Proposition 5.

Numerical schemes

We describe here the numerical methods that we will be using in the next section.

4.1.

The time-splitting spectral method. We start by the time-splitting spectral method that we use to solve (NLS) numerically. The equations are solved on a bounded interval I = (-a, a). We use a uniform spatial grid with mesh size h > 0 and grid points x k = x 0 + kh, k = 0, . . . , K, where K + 1 ∈ N is the (odd) number of grid points. Then h = 2a/K. The time grid is given by t n = t 0 + nτ , n ∈ N 0 , where τ > 0 is the time step size and t 0 the initial time. We set (u j ) n k := u j (t n , x k ), where j = 1, 2, k = 0, . . . , K, and n ∈ N 0 . We split the system (NLS) into the two subsystems

i∂ t u j = -µ j |u j | 2 u j -β|u 3-j | 2 u j , j = 1, 2, ( 36
)
i∂ t u j = - 1 2 ∂ xx u j , j = 1, 2, ( 37
)
considered on [t n , t n+1 ] and subject to some initial data. These subsystems are solved as follows.

Step 1: Computing the evolution of (36) we observe that the quantities |u j | 2 remain unchanged. Therefore, we "freeze" these values at time t n and solve the resulting linear ODEs exactly in the interval [t n , t n + τ /2], giving at time t n + τ /2:

(u 1 ) * k = exp i τ 2 µ 1 |(u 1 ) n k | 2 + β|(u 2 ) n k | 2 (u 1 ) n k ,
and analogously for (u 2 ) * k .

Step 2: We solve (37) for j = 1, 2 in the interval [t n , t n + τ ], discretized in space by the Fourier spectral method and solved exactly in time:

(u j ) * * k = 1 K + 1 K/2 m=-K/2 exp -iτ ν 2 m ( u j ) * m exp iν m (x k -x 0 ) , j = 1, 2,
where ν m = 2πm/(x Kx 0 ) and

( u j ) * m = K l=0 (u j ) * l exp (-iν m (x l -x 0 )) , m = - K 2 , . . . , K 2 .
Step 3: We solve (36) on [t n + τ /2, t n+1 ] using the discretization of Step 1 with (u j ) * * k instead of (u j ) n k and obtain (u j ) n+1 k .

4.2.

The normalized gradient flow. We will also need to compute the ground state solution of the following elliptic system with fixed masses

-∂ xx φ 1 + ω 1 φ 1 -µ 1 φ 3 1 -βφ 2 2 φ 1 = 0, -∂ xx φ 2 + ω 2 φ 2 -µ 2 φ 3 2 -βφ 2 1 φ 2 = 0. (38)
To that purpose, we use the normalized gradient flow. The problem can also be viewed as a nonlinear eigenvalue problem with ω 1 , ω 2 being the eigenvalues, which can be computed from the corresponding eigenfunctions (j = 1, 2):

ω φ j = R -|∂ x φ j | 2 + µ j φ 4 j + βφ 2 1 φ 2 2 dx R φ 2 j dx .
We solve (38) by normalized gradient flow with given (a 1 , a 2 ), such that:

R φ 2 1 dx = a 2 1 and R φ 2 2 dx = a 2 2 .
The standard gradient flow with discrete normalization consists in introducing an imaginary time in the nonlinear Schrödinger equations, thus looking at the imaginary time propagation (t → -it) and after every step project the solutions such that the L 2 -norms are equal to (a 2 1 , a 2 2 ). In [START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF] the authors present the normalized gradient flow, prove it is energy diminishing, and propose numerical methods to discretize it. Hereafter we adapted the normalized gradient flow for the given system and discretized it by a semi-implicit Backward Euler finite differences scheme. For (φ j ) n k = φ j (t n , x k ) being the discrete solution, x k = x 0 + k • h the grid points with k = 0, 1, . . . , K -1, K and the time sequence 0 < t 1 < t 2 < • • • < t n < t n+1 < . . . with τ = t n+1t n we have the following discretization of the normalized gradient flow:

Step 1: We first solve on [t n , t n+1 ] with initial data φ j (t n , x k ):

(φ j ) * k -(φ j ) n k τ = (φ j ) * k+1 -2(φ j ) * k + (φ j ) * k-1 h 2 + µ j |(φ j ) n k | 2 (φ j ) * k + β|(φ 3-j ) n k | 2 (φ j ) * k with j = {1, 2}. As solution we get (φ j ) * k . Step 2: (φ j ) * k is then normalized to get finally (φ j ) n+1 k (φ j ) n+1 k = a k (φ j ) * k (φ j ) * k L 2 .
For t → +∞ we obtain the ground state solution (φ 1 , φ 2 ) of (38) with L 2 -norms equal to (a 2 1 , a 2 2 ) and frequency parameters (ω φ1 1 , ω φ2 2 ).

Numerical experiments

Our ansatz for the initial data in the next experiments is the following.

u j (t 0 , x) = e i(ωj t0-v 2 j t/4+vj x/2) Q ωj (x -v j t 0 -x j ) j = 1, 2 (39)
with Q ω defined in [START_REF] Ambrosetti | Standing waves of some coupled nonlinear Schrödinger equations[END_REF]. Without loss of generality, we may assume that 0 < v 1 = -v 2 (Galilean invariance), hence the soliton on the first component will be traveling to the right and the soliton on the second component will be traveling to the left. We will also assume that x 1 = x 2 = 0 (invariance by translation in time and space) and to guarantee that our solitons on the first and second components are positioned on the left and on the right respectively and at a sufficiently large distance, we will chose the initial time to be t 0 = -10.

5.1.

Purely elastic interaction. The integrable case has been studied in depth in the book [START_REF] Ablowitz | Discrete and continuous nonlinear Schrödinger systems[END_REF]. In particular, in that case the system is completely integrable via the inverse scattering transform and explicit solutions may be exhibited. For example, there exists a solution U = (u 1 , u 2 ) of (NLS) which has the following behavior. At time t → -∞,

u 1 (t, x) ∼ e i( 1 2 v1x+(ω1-1 4 |v1| 2 )t) √ ω 1 Q ( √ ω 1 (x -v 1 t)) , (40)
u 2 (t, x) ∼ e i( 1 2 v2x+(ω2-1 4 |v2| 2 )t) √ ω 2 Q ( √ ω 2 (x -v 2 t)) ,
whereas at time t → +∞ the two components have interacted and the outcome are solitons of same speed and frequency but with a shift in phase and translation:

u 1 (t, x) ∼ φ1 (t, x) := e iθ1 e i( 1 2 v1x+(ω1-1 4 |v1| 2 )t) √ ω 1 Q ( √ ω 1 (x -v 1 t -τ 1 )) , (41)
u 2 (t, x) ∼ φ2 (t, x) := e iθ2 e i( 1 2 v2x+(ω2-1 4 |v2| 2 )t) √ ω 2 Q ( √ ω 2 (x -v 2 t -τ 2 )) ,
where the shift parameters are given by the formulas

τ 1 = - ln(|χ 1 |) √ ω 1 , τ 2 = ln(|χ 2 |) √ ω 2 , θ 1 = χ 1 |χ 1 | , θ 2 = χ 2 |χ 2 | , χ 1 = v 1 -v 2 + i2( √ ω 1 + √ ω 2 ) v 1 -v 2 + i2( √ ω 1 - √ ω 2 ) , χ 2 = v 1 -v 2 + i2( √ ω 1 + √ ω 2 ) v 1 -v 2 -i2( √ ω 1 - √ ω 2 ) .
Note that since the speeds and amplitudes of the solitons are not modified the interaction is considered elastic. Note also that we have a pure two-speeds solitary wave at both ends of the time line, without any appearance of dispersion despite the interaction. This is a characteristic feature of completely integrable systems. The parameters are chosen as follows:

µ 1 = µ 2 = β = 1, ω 1 = 5, ω 2 = 1, v 1 = 1, v 2 =
-1, and x 1 = x 2 = 0, thus we start with two solitons located at ±10, respectively, which move to each other and observe them until t f inal = 10. Furthermore, the integration domain is I = (-20, 20) with K = 1024 spatial grid points and τ = 10 -3 . We plot the solution (u 1 , u 2 ) to (NLS) and observe that the two solitons remain unchanged after the interaction, excepting a shift in phase and translation. The shift in translation can be observed in Figure 1. Our numerical experiments are in good line with the theoretically predicted behavior. We have compared pointwise the numerical solution |u 5.2. Symmetric collision. In this experiment the outcome are still multi-speed solitary waves. After the collision we observe a small part of the soliton moving along with the soliton of the other component. We see that at the final time the solution fits well the ground state solution obtained by minimizing the system with fixed masses. Here we choose µ 1 = µ 2 = 1 and β = 3. We have solved the system (NLS) in the interval (-200, 200) with periodic boundary conditions, initial time t 0 = -10 and initial condition (39) and the parameters ω 1 = ω 2 = 1, v 1 = 2, v 2 = -2, and x 1 = x 2 = 0, thus we start with two solitons located at ±10, respectively, which move to each other and observe them till t f inal = 40. Moreover, we use K = 4096 in space and time step τ = 10 -3 .

Here the interaction parameter β is greater than in the first experiment and after the interaction at time zero we see in each component a small part going with the other component, respectively (Figure 3).

At the final time t = 40 we look at the left side part of each component u - j = u j • χ [-∞,0] , and compute the L 2 -norms

a 2 j = R |u - j (40, x)| 2 dx.
Note that the mass of each component the position density pointwise to the position density of the final solution located on the left of the origin (u - 1 (40, x), u - 2 (40, x)) and realize that they fit very well. In this case, the interaction result into a new repartition of the mass and of the energy so as to approximate a ground state profile.

M j (t) = 1 2 R |u j (t, x)| 2 dx is conserved, M j (0) = M j (t) for j = 1

Dispersive inelastic interaction.

In this experiment we observe a loss of energy, mass and momentum in a small dispersive part located at the interaction place and a small dispersive part moving to the boundaries. Here we have solved the system (NLS) with β = -1 and µ 1 = µ 2 = 1 in the interval (-500, 500) with 8192 grid points for the spatial discretization, with periodic boundary conditions, initial time t 0 = -10 and initial condition (39) with the following parameters v 1 = -v 2 = 2.7, ω 1 = ω 2 = 1 , and x 1 = x 2 = 0.

Figure 5 shows the evolution of the waves. In Figure 6 we see the position density of the initial soliton u 1 (0, x) and the one of the final soliton at t = 90, with lost mass and energy. In this experiment we observe a change of the sign of the speeds, thus the solitons are reflected after interaction. Here we have the same parameter as in the case of a dispersive interaction, β = -1 and µ 1 = µ 2 = 1, but the initial velocity is smaller, precisely we have the initial condition (39) with the following parameters v 1 = -v 2 = 0.5, ω 1 = ω 2 = 1, and x 1 = x 2 = 0. Since the velocity is not so large, and we observe the soliton until t f inal = 10, we solve (NLS) on (-20, 20) with 1024 spatial grid points, periodic boundary conditions and time step equal to τ = 10 -3 .

In Figure 7 we see the reflexion of the two solitons after interaction, thus only the velocities changed sign. In Figure 8 we compare point wise the solution at the final time to the initial solution, and observe that it remains unchanged.

These were various examples of the possible outcomes of multi-speeds solitary waves interactions. The above described situations may not be the only possibles and we do not claim comprehensiveness. Proof of Lemma 6. The proof of the modulation lemma is inspired from [START_REF] Martel | Multi solitary waves for nonlinear Schrödinger equations[END_REF] and relies on the implicit function theorem.

Let us define F :

H 1 (R) × R + × R × R -→ R 3 by F (u, ω, γ, y) =   (u -R(ω, γ, y), R(ω, γ, y)) 2 (u -R(ω, γ, y), iR(ω, γ, y)) 2 (u -R(ω, γ, y), ∂ x R(ω, γ, y)) 2  
where, for the sake of simplicity, we removed the i indexes and R(ω, γ, y) denotes the soliton R(ω, γ, y)

= e i( 1 2 v•x+γ) 1 √ µ Q ω (x -y). Note that, ∀(ω, γ, y) ∈ R + × R × R, F (R(ω, γ, y), ω, γ, y) = 0.
The key idea here is to apply the implicit function theorem to function F at point (R(ω, γ, y), ω, γ, y) for fixed (ω, γ, y) ∈ R + × R × R. Simple computations due to properties of Q ω , mainly the fact that Q ω is radial and even, lead to ∇F (R(ω, γ, y), ω, γ, y) diagonal with diagonal terms a, b, c given by:

a = 1 µ Re R ∂ ω Q ω (x -y)Q ω (x -y)dx < 0 b = -(∂ γ R(ω, γ, y), iR(ω, γ, y)) 2 = -1 µ Re R Q ω (x -y) 2 dx < 0 c = -1 µ ∂ k Q ω 2 L 2 < 0.
Implicit function theorem finally gives the existence of parameters ωj , γj , xj as functions of time. To prove that these functions are actually of class C 1 , a standard regularization argument is needed. We refer to [START_REF] Martel | Instability of solitons for the critical generalzed Korteweg-de-Vries equation[END_REF] for more details on that argument. Now, in order to be more readable, we prove estimate (13) for j = 1. In that purpose, let us write the equation of evolution satisfied by ε 1 , namely:

i∂ t ε 1 + L(ε 1 , ε 2 ) + N (ε 1 , ε 2 ) = -(ω 1 - v 2 1 4 -∂ t γ1 ) R1 -i∂ t ω1 ∂ ω Q ω1 e i( 1 2 v1•x+γ1) -i(v 1 -∂ t x1 )∂ x Q ω1 . (42)
where

L(ε 1 , ε 2 ) = ∂ xx ε 1 + µ 1 | R1 | 2 + β| R2 | 2 ε 1 + 2µ 1 Re( R1 ε1 ) R1 + 2βRe( R2 ε2 ) R1
and

N (ε 1 , ε 2 ) = µ 1 |ε 1 | 2 ε 1 + β|ε 2 | 2 ε 1 + µ 1 |ε 1 | 2 + β|ε 2 | 2 R1
+ 2µ 1 Re( R1 ε 1 ) + 2βRe( R2 ε2 ) ε 1 .

It is to be noticed that the modulation terms ω1 -

v 2 1
4 -∂ t γ1 , ∂ t ω1 , and v 1 -∂ t x1 appear in the right hand side of this evolution equation. The main idea of the proof is to take in (42) the scalar products of both sides of the equation with respectively R1 , i R1 and ∂ x R1 . We then use the orthogonality conditions [START_REF] Côte | High-speed excited multi-solitons in nonlinear Schrödinger equations[END_REF] as well as properties on Q ω1 . In the left hand side, we transfer the derivatives acting on ε 1 on the other side of the scalar product thanks to the modulation conditions [START_REF] Côte | High-speed excited multi-solitons in nonlinear Schrödinger equations[END_REF] (for the time derivatives) and integrations by parts (for the space derivatives). We finally use the equation at hand on R1 :

i∂ t R1 + ∂ xx R1 + µ 1 | R1 | 2 R1 = ω1 + v 2 1 4 -∂ t γ1 R1 + j √ µ 1 e i( 1 2 v1•x+γ1) ∂ t ω1 ∂ ω Q ω1 +i v 1 -∂ t x1 µ 1 e i( 1 2 v1•x+γ1) ∂ x Q ω1 . (43)
For more simplicity, we only develop the computations for the scalar product with R1 , the other cases are obtained using the same arguments.

Taking the scalar product with R1 in equation (42) leads to:

i∂ t ε 1 + L(ε 1 , ε 2 ) + N (ε 1 , ε 2 ), R1 2 = -ω1 - v 2 1 4 -∂ t γ1 R1 2 -∂ t ω1 Re R ie i( 1 2 v1•x+γ1) ∂ ω Q ω1 R1 dx -(v 1 -∂ t x1 )Re R ie i( 1 2 v1•x+γ1) ∂ x Q ω1 R1 dx (44)
Right hand side of (44): First, using equation ( 11) leads to:

Re R ie i( 1 2 v1•x+γ1) ∂ ω Q ω1 R1 dx = 1 √ µ 1 Im R ∂ ω Q ω1 (x)Q ω1 (x -x1 )dx = 0 and Re R ie i( 1 2 v1•x+γ1) ∂ x Q ω1 R1 dx = 1 √ µ 1 Im R ∂ x Q ω1 (x)Q ω1 (x -x1 )dx = 0.
Thus, (44) reduces to:

i∂ t ε 1 + L 1 (ε 1 , ε 2 ) + N (ε 1 , ε 2 ), R1 2 = -ω1 - v 2 1 4 -∂ t γ1 R1 2 . ( 45
)
Left hand side of (45): First, deriving modulation condition ε 1 , R1 2 = 0 with respect to time gives:

i∂ t ε 1 , R1 2 = ε 1 , i∂ t R1 2 .
Let us now develop L(ε 1 , ε 2 ), R1 2 :

L(ε 1 , ε 2 ), R1 2 = ε 1 , ∂ xx R1 + 3µ 1 | R1 | 2 R1 + β| R2 | 2 R1 2 + 2β ε 2 , | R1 | 2 R2 2 .
Finally, N (ε 1 , ε 2 ), R1

2 read:

N (ε 1 , ε 2 ), R1 2 = µ 1 |ε 1 | 2 ε 1 + β|ε 2 | 2 ε 1 , R1 2 + (µ 1 |ε 1 | 2 + β|ε 2 | 2 ) R1 , R1 2 + 2µ 1 Re( R1 ε 1 ) + 2βRe( R2 ε 2 ) ε 1 , R1 2
Equation (45) thus leads to the left hand side term:

ε 1 , i∂ t R1 + ∂ xx R1 + µ 1 | R1 | 2 R1 2 + β ε 1 , | R2 | 2 R1 2 + 2µ 1 ε 1 , | R1 | 2 R1 2 + 2β ε 2 , | R1 | 2 R2 2 + N (ε 1 , ε 2 ), R1 2
. Now, using equation (43) satisfied by R1 gives:

ε 1 , i∂ t R1 + ∂ xx R1 + µ 1 | R1 | 2 R1 2 = (ω 1 - v 2 1 4 -∂ t γ1 ) ε 1 , R1 2 + ∂ t ω1 ε 1 , i √ µ 1 e i( 1 2 v1•x+γ1) ∂ ω Q ω1 2 + (v 1 -∂ t x1 ) ε 1 , i √ µ 1 e i( 1 2 v1•x+γ1) ∂ x Q ω1 2 .
Finally, modulation condition ε 1 , R1 2 = 0 leads to

(ω 1 - v 2 1 4 -∂ t γ1 ) R1 2 + ∂ t ω1 ε 1 , i √ µ 1 e i( 1 2 v1•x+γ1) ∂ ω Q ω1 2 + (v 1 -∂ t x1 ) ε 1 , i √ µ 1 e i( 1 2 v1•x+γ1) ∂ x Q ω1 2 = -2µ 1 ε 1 , | R1 | 2 R1 2 -β ε 1 , | R2 | 2 R1 2 -2β ε 2 , | R1 | 2 R2 2 -N (ε 1 , ε 2 ), R1 2 . ( 46
)
Thanks to Lemma 8, it is readily seen that most terms in this equation are of order O( ε 2 ). Finally, (46) can be re-written in a simpler way: With the same kind of arguments, taking in (42) scalar product with i R1 and i∂ x R1 respectively, we get two other equations that can be re-written as a linear system solved by the "modulation vector" Mod(t) =   ω1 (t) -   and, with the help of Lemma 8, for all t ∈ [t 0 , T n ], A(t) , B(t) C ε 2 , and hence [START_REF] Côte | Construction of multi-soliton solutions for the L 2supercritical gKdV and NLS equations[END_REF].

R1 2 + a 1 (t) (ω 1 - v 2

2 Figure 1 .

 21 Figure 1. Purely elastic interaction: Position densities |u 1 (t, x)| 2 and |u 2 (t, x)| 2 as functions of space and time for µ 1 = µ 2 = 1 and β = 1

1 | 2 4 Figure 2 .

 1242 Figure 2. a) Plots of |u 1 (10, x)| 2 ('-' line) and | φ1 (x)| 2 ('.' line). b) Difference of |u 1 (10, x)| 2 -| φ1 (x)| 2

2 Figure 3 .

 23 Figure 3. Symmetric collision: Position densities |u 1 (t, x)| 2 and |u 2 (t, x)| 2 as functions of space and time for µ 1 = µ 2 = 1 and β = 3

Figure 4 .

 4 Figure 4. Symmetric collision: a) Plots of |u - 1 (40, x)| 2 ('-' line) and |φ 1 (x + x 1 )| 2 ('.' line) b) Plots of |u - 2 (40, x)| 2 and |φ 2 (x + x 2 )| 2 , where x j is the x corresponding to the max |u j | 2 , with j = 1, 2.

2 Figure 5 . 2 Figure 6 .

 2526 Figure 5. Dispersive inelastic interaction: Plots of |u 1 (t, x)| 2 (left) and |u 2 (t, x)| 2 (right) as functions of space and time

2 Figure 7 .x

 27 Figure 7. Reflexion: Plots of |u 1 (t, x)| 2 (left) and |u 2 (t, x)| 2 (right) as functions of space and time

Figure 8 .

 8 Figure 8. Reflexion: Plots of |u 1 (-10, xx 1 )| 2 ('-' line) and |u 1 (10, x)| 2 ('.' line), where x 1 is the translation such that we can compare the two solitons

v 2 1 4 -

 24 ∂ t γ1 (t) ∂ t ω1 (t) v 1 -∂ t x1 (t)..This linear system takes the following form:(Γ(t) + A(t))Mod(t) = B(t) ω Q ω1 (x)Q ω1 (x -x1 (t)

  Ce -3 √ ω * v * t e -√ ω * |x| . The derivative ∂ x R1 is treated in the same way, with the only difference that |v 1 | now appears in the estimate, due to the term e i 1 2 v1•x in the definition[START_REF] Cazenave | Orbital stability of standing waves for some nonlinear Schrödinger equations[END_REF] of R1 . This finishes the proof.Lemma 9 (Expansion of the global action S). For all t ∈ [t 0 , T n ] we have

				7, we have
		x1	v * 2 √ 2	t, and ω1 1152ω * ,
	which implies	| R1 (t, x)|χ 2 L (x)

1 4 -

 4 ∂ t γ1 ) + a 2 (t)∂ t ω1 + a 3 (t)(v 1 -∂ t x1 ) = b 1 (t)where, for all t ∈ [t 0 , T n ], |a1 (t)| + |a 2 (t)| + |a 3 (t)| + |b 1 (t)| C ε(t) 2 .

The work of F. D. is partially supported by PHC AMADEUS 2014 31471ZK. The work of S. L. C. is partially supported by ANR-11-LABX-0040-CIMI within the program ANR-11-IDEX-0002-02, ANR-14-CE25-0009-01 and PHC AMADEUS 2014 31471ZK. The work of R.M.W. is supported by the FWF Hertha-Firnberg Program, Grant T402-N13 and Austrian-French Project WTZ-Amadée FR 18/2014.