Johanna Klein

Hassan Khartabil

Jean-Charles Boisson

Julia Contreras-García

Jean-Philip Piquemal

Eric Hénon
email: eric.henon@univ-reims.fr

A New Way for Probing Bond Strength

The covalent chemical bond is intimately linked to electron sharing between atoms.

The recent Independent Gradient Model (IGM) and its δg descriptor provide a way to quantify locally this electron density interpenetration from wavefunction calculations.

Each bond has its own IGM-δg pair signature. The present work establishes for the first time a strong link between this bond signature and the physically grounded bond force constant concept. Analyzing a large set of compounds and bonds, the Intrinsic Bond Strength Index (IBSI) emerges from the IGM formulation. Our study shows that the IBSI does not belong to the class of conventional bond orders (like Mulliken, Wiberg, Mayer, delocalization index or ELF), but is rather a new index complementary thereto, related to the bond force. A fundamental outcome of this research is a novel index allowing to range all two-centre chemical bonds by their intrinsic strength in molecular situation. We believe that the IBSI is a powerful and robust tool for interpretation accessible to a wide community of chemists (organic, inorganic chemistry, including transition metal complexes and reaction mechanisms).

Introduction

The chemical bond 1 is a model commonly used by chemists to describe molecules as made of atoms linked to one another. Although lacking a clear physical basis (the bond is not an observable quantity) it is a quite useful concept frequently employed to help chemists to give explanations and make predictions about physical and chemical properties of substances. In theoretical chemistry, the atom connectivity concept is explicitly involved in models such as force fields used to perform molecular dynamic simulations. In chemical education, the chemical bond, very useful to represent molecular structures and to illustrate chemical processes, has acquired a central place and remains fundamental. [START_REF] Chavent | Soon in your lecture halls, chemistry will go to Hollywood[END_REF] The physical interactions that establish cohesion in molecules can be tackled by the use of quantum chemistry. [START_REF] Frenking | The Chemical Bond[END_REF] The use of the electron density (ED) topological analysis for the characterization of bonds and more generally of interactions in molecular systems has led to a variety of investigative methods based on local descriptors. Examples of such approaches are: the electron localization function ELF, [START_REF] Becke | A simple measure of electron localization in atomic and molecular systems[END_REF][START_REF] Silvi | Classification of chemical bonds based on topological analysis of electron localization functions[END_REF] the localized electron detector LED, [START_REF] Bohorquez | The localized electrons detector as an ab initio representation of molecular structures[END_REF] the Density Overlap Region Indicator DORI [START_REF] De Silva | Simultaneous Visualization of Covalent and Noncovalent Interactions Using Regions of Density Overlap[END_REF] or the Non-Covalent-Interaction approach NCI. [START_REF] Johnson | Revealing noncovalent interactions[END_REF] As complementary tools of the Atoms In Molecules approach (QTAIM [START_REF]Atoms In Molecules, A Quantum Theory[END_REF][START_REF] Matta | The Quantum Theory of Atoms In Molecules[END_REF] ), these methods allow for producing insightful three-dimensional isosurfaces localizing the interactions within the system from theoretically derived electron density. Recently, in the continuation of the NCI [START_REF] Johnson | Revealing noncovalent interactions[END_REF] analysis originally proposed in 2010, we introduced a new local approach for probing strong and weak interactions in molecules. Based on the Independent Gradient Model (IGM [START_REF] Lefebvre | Accurately extracting the signature of intermolecular interactions present in the NCI plot of the reduced density gradient versus electron density[END_REF][START_REF] Lefebvre | The Independent Gradient Model: A New Approach for Probing Strong and Weak Interactions in Molecules from Wave Function Calculations[END_REF] ) that represents a non-interacting reference system, the δg descriptor has been defined. It measures locally the electron sharing between interacting fragments and reveals regions of space where the ED between fragments clashes. Compared to the original formulation of NCI, which has a semi-quantitative value, the IGM-δg approach is able to quantify the interaction caused by electron sharing. It provides us with a molecular δg signature as illustrated in Figure 1 (top right panel), by plotting δg versus the ED oriented with the sign of the second density hessian eigenvalue λ 2 . It covers a large range of interactions, associated or not with bond critical point (BCP), including covalent bonds, bonding in transition metal compounds, hydrogenbonding and van der Waals interactions. It is not dimensionless and good correlations have been found between peak heights in the IGM-δg signature of intermolecular complexes and the calculated stabilization energy for different hydrogen bonds. [START_REF] Lefebvre | Accurately extracting the signature of intermolecular interactions present in the NCI plot of the reduced density gradient versus electron density[END_REF] A critical and appealing feature of the IGM-δg computational scheme is its ability to separate and extract the signature of selected atom pairs from the overall interaction signal. [START_REF] Lefebvre | The Independent Gradient Model: A New Approach for Probing Strong and Weak Interactions in Molecules from Wave Function Calculations[END_REF] Thereby, a bond-by-bond IGM-δg pair picture can be obtained from a wavefunction as illustrated in Figure 1. Each chemical bond has its signature. To the best of our knowledge, this possibility is not available within any other descriptor-based approaches.

Figure 1: Molecular IGM-δg signature, selected IGM-δg pair bond signatures in prop-1-ene-2-thiol and δg = 0.3 a.u. isosurfaces colored according to the BGR scheme over the range -0.2 < sign(λ 2 )ρ < 0.2 a.u. (the sign of the second density hessian eigenvalue λ 2 serves to distinguish between attractive and repulsive situations).

At this stage, it might be tempting to directly associate the strength of the bond with the peak height of the IGM-δg pair bond profile. In the example shown in Figure 1 (prop-1-ene-2-thiol), the C = C double bond results in a larger δg pair peak height (0.68) than a single bond C -C (0.51). Further, the δg pair peaks for CH 2 -H (0.57) and S -H (0.41) fairly reflect their homolytic bond dissociation enthalpy ordering (given by Y-R. Luo [START_REF] Luo | Chemical Bond Energies[END_REF] around 420 and 350-370 kJ.mol -1 , respectively, at 298K). This apparent interrelationship between bond strength and δg pair peak height makes sense since bonds are achieved by electron sharing [START_REF] Frenking | The Chemical Bond[END_REF] measured to some extent by δg pair at a point in space. However, each δg pair bond signature has its own individual profile, which suggests that attempts to correlate the δg pair descriptor to energetics should not be limited to the use of a local information (at BCP for instance) but should consider an integration scheme. Such integration schemes have been explored in the past few years. In their study of the integration of charges within the NCI region, [START_REF] Contreras-García | Analysis of Hydrogen-Bond Interaction Potentials from the Electron Density: Integration of Noncovalent Interaction Regions[END_REF] J. Contreras-García and W. Yang recovered the position of the intermolecular potential minimum for 4 hydrogen-bond complexes. However, although mimicking relatively well the shape of the binding energy curve for the water dimer, the correct curvature around the minimum was not obtained for the other dimers. In a study conducted for a series of 30 weakly bound complexes, Presti 14 found a good correlation between the stabilization energy and the kinetic energy density integrated over the volume enclosed by the NCI isosurface.

Unfortunately, this correlation breaks down for complexes of different nature. Moreover, in this previous work based upon the NCI approach, the integration domain was defined in an arbitrary manner (domain below RDG = 0.5 a.u.) and ignored some grid points in the interaction region. More recently, for a set of 53 systems, Ananyev [START_REF] Ananyev | Toward a Rigorous Definition of a Strength of Any Interaction Between Baders Atomic Basins[END_REF] presented a strong correlation between effective force constants of diatomic molecules and the potential energy density integrated over the AIM surface between two Bader's atomic basins. Very recently, using an integration scheme, the IGM-δg approach has been employed to derive a score indicating the role played by each atom in the formation of host-guest assemblies. [START_REF] Ponce-Vargas | An Atomic Decomposition Scheme of Noncovalent Interactions applied to Host-Guest Assemblies[END_REF] The present study falls in the class of integration schemes that have emerged in the past few years.

The basic aim of this work is to integrate the IGM-δg pair bond signature to characterize a bond by a unique score. The result is a global bond descriptor ∆g pair strongly connected to the physically grounded bond strength concept. In this study, we present results obtained using a set of 235 molecules spanning a broad range of different chemical bonds. Potential relationships between ∆g pair and several bond orders as well as with bond stretching force constants have been sought for. A novel Intrinsic Bond Strength Index (IBSI) stems from these attempts. Since the IGM-δg approach is a relatively new concept, we will start by reintroducing it in the Methods section. Also, details on the integration scheme and the design of the new index IBSI will be sketched out. This will be followed by computational details, results and discussion before to conclude.

Methods

IGM-δg approach

The IGM-δg approach has been designed to reveal and quantify molecular interaction from the ED topology, more precisely from its gradient ∇ ∇ ∇ρ. All details can be found in our previous works. [START_REF] Lefebvre | Accurately extracting the signature of intermolecular interactions present in the NCI plot of the reduced density gradient versus electron density[END_REF][START_REF] Lefebvre | The Independent Gradient Model: A New Approach for Probing Strong and Weak Interactions in Molecules from Wave Function Calculations[END_REF] Its implementation involves the use of a three-dimensional grid built to encompass the chemical system. A prerequisite to apply this local approach is to express each molecular ED gradient component ∂ρ ∂x as a sum of individual terms, one for each interacting component (atom, fragment, molecule), for instance for two fragments A and B:

∂ρ ∂x = ∂ρ A ∂x + ∂ρ B ∂x .
Then, the idea behind the Independent Gradient Model is to define a non-interacting reference system that cancels the effect of electron delocalization (electron sharing) occurring between the two fragments upon bringing them closer to each other. This is achieved by using absolute values upon summing the gradient contributions of ED sources (fragments) in the expression of the total ED gradient, leading to the non-interacting ref-

erence ∂ρ ∂x IGM = ∂ρ A ∂x + ∂ρ B ∂x .
Thereby, the resulting ED gradient attenuation normally expected between bonding fragments in the so-called "interaction corridor" (see Figure 2) and due to individual terms having opposite signs in the gradient sum ∂ρ ∂x , vanishes. The resulting norm ∇ ∇ ∇ρ IGM is a virtual upper limit of the real total ED gradient. Hence, the descriptor δg = ∇ ∇ ∇ρ IGM -|∇ ∇ ∇ρ| measures locally the electron sharing, hereafter also referred to as ED contragradience. Two very important points must be emphasized for the rest of this study: (1) only the contragradience between A and B is captured by δg and (2) the same ED is considered for the non-interacting system (IGM) and the real system. This makes the δg descriptor intrinsically linked to the AB interaction in molecular situation. In other words, δg quantifies the AB interaction without resorting to calculations performed on separated fragments. Once δg computed at any point of the three-dimensional grid, a 2D plot can be constructed as illustrated on Figure 1. The first concern in implementing the IGM approach is to choose a gradient decomposition scheme in which each contribution is associated to an ED source to be isolated. For ED coming from quantum-mechanical calculations, this can be done thanks to the Gradient Based Partition (GBP) that we have previously described. [START_REF] Lefebvre | The Independent Gradient Model: A New Approach for Probing Strong and Weak Interactions in Molecules from Wave Function Calculations[END_REF] This partition has been shown to be hardly dependent on the basis set (beyond the small STO-3G basis set). Within the GBP, the most basic ED gradient source that can be considered is an atomic orbital. Regrouping them by atom leads to an atomic partition of the ED gradient. Using this atomic partition yields the δg descriptor extracting the interaction between all atoms of the system. Further, grouping these atomic gradients into two fragments allows for extracting just the intermolecular interaction between them (δg inter , not used in this work). Finally, formulating a gradient partition truncated to two atoms leads to the δg pair bond signature we are interested in, as illustrated in Figure 1. It opens the way to probe specific bonds in a molecule by using the ED derived from quantum calculations. The local δg pair descriptor was computed using a modified version of NCIPlot. 17 Such grid computational scheme could leverage the massively parallel GPU technology for the acceleration of δg pair calculation, as demonstrated with NCI in a previous study. [START_REF] Rubez | GPU accelerated implementation of NCI calculations using promolecular density[END_REF] Just like in the NCI approach, this δg pair signature can be transposed in the real space taking the form of isosurfaces colored according to the value of signed ED (sign(λ 2 )ρ) (see Figure 1). Clearly, δg pair accounts for the tendency of electrons to be shared between the two atoms. It measures to some extent (not in a direct manner) the kinetic energy lowering resulting from electrons being able to roam over a larger area in the interacting system than in separate fragments. This sharing can take place either in bonding regions or nonbonding regions. Hence, the δg pair bond signature offers two complementary interpretations: (1) the δg pair peak height tells us how much the electron delocalization is important from a kinetic point of view and (2) the ED curvature where the δg pair occurs tells us if we are dealing with an attractive (λ 2 < 0) or repulsive (λ 2 > 0) local region from a potential point of view.

IGM bond index IBSI

In order to get a global score for the studied bond, we have decided to carry out the integration of the δg pair signature over the interaction volume V: ∆g pair = V δg pair dV . An outstanding advantage of the IGM-δg descriptor in achieving this integration lies in that there is no need to delimit the interaction region. Actually, non-zero values of the δg descriptor exclusively correspond to interaction situations. Thus, the numerical integration is computed as the summation extending over all the grid points lying in the grid volume. No attempt has been made to integrate the ED ρ as others have done. Actually, the 2D δg pair signature seems to hold all the ingredients to characterize bonding patterns and it appears to be a natural candidate for this integration scheme.

The first attempts to compare ∆g pair to other properties like bond indices or force constants have shown the need to divide ∆g pair by the square of the internuclear distance. The reason is that ED contragradience only occurs in the space between the two atoms (the interaction corridor) and hence the corresponding integration volume in itself influences the resulting sum ∆g pair . For instance, the H 2 bond, with a very short bond length (0.74 Å), is characterized by a narrow interaction corridor. As a result, it is definitely characterized by a too small value of ∆g pair compared to other bonds with larger separation distance. Our experience led us to conclude that it is important to make ∆g pair independent of the bond size in order to compare different bonds. The resulting IGM bond index, which measures the amount of electron sharing per square unit length, is then: ∆g pair = V δg pair d 2 dV . Moreover, to compare between bond indices, the IBSI (∆g pair ) has been normalized to 1 for the H 2 molecule:

∆g pair = V δg pair d 2 dV V δg H 2 d 2 H 2 dV (1)
As a consequence, the IBSI is dimensionless. The influence of the level of theory on the gross value of V δg H 2 dV has been studied and shown to be hardly dependent on both the basis set and level of theory (see ESI). Although the IBSI value for H 2 is rather stable, it is recommended to employ the same level of theory to compare IBSI values across different molecules.

Based on case studies (see supporting information), we have shown that the IBSI is generally not significantly basis set dependent (beyond the STO-3G basis). The extent of variability of the IBSI measured by the relative standard deviation is generally in the range 2-6% for the studied organometallic compounds, generally more homogenous compared to Mayer, Wiberg and Mulliken bond orders. We strongly advise not using the STO-3G basis set, but we recommend the use of 6-31G ,although adding a diffuse function could be more relevant for weakly bonded systems. For atoms beyond Kr, the Def2-TZVP [START_REF] Weigend | Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy[END_REF][START_REF] Weigend | Accurate Coulomb-fitting basis sets for H to Rn[END_REF] basis set can be used (available for all elements up to Radon) or an effective core potential calculation can be performed (with LANL2DZ 21-23 basis set). In extreme cases, in very polar bonds like LiH, AlO, AlCl, MgC, a larger dependence has been observed (around 10-20%).

A comparative case study of ab initio (HF, MP2) and DFT (B3LYP, SVWN, BP86 and M06-2X) methods indicates that the IBSI hardly depends on the method used (coefficient of variation of at most 7% in extreme cases, but more generally around 2%). We advise not using the HF method. We therefore recommend the DFT/6-31G** level of theory for its performance/price ratio and using the same method for comparative studies.

Molecule test set

This works aims at finding a physically well-founded characterization of the descriptor δg pair and its integrated form ∆g pair for a given atom pair. To this end, a 235-molecule test set has been considered, both organic and inorganic, covering typical functional groups (see ESI) and leading to the study of 677 bonds. It is not limited to known molecules but also includes uncommon compounds. Neutral, charged and radical species have been studied. This large set of molecules covers a full range of bonding interactions, including covalent bonding, metal coordination complexes (with Cr, Fe, Co, Ni, Ru, Re, Pt), weakly hydrogen-bond and halogen-bond complexes.

Quantum-mechanical calculations

In practice, the electron density (ED) is obtained from a quantum-mechanical calculation.

In this study, the Gaussian16 package 24 was employed. If not otherwise stated, full molecular geometry optimizations were conducted at the DFT(M06-2X/6-31G**) [START_REF] Zhao | The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals[END_REF] level of theory in the gas phase. Some species with heavy atoms required using the Los Alamos effective core potential through the LANL2DZ basis set specification in gaussian. The integral SuperFineGrid option has been employed in numerical differentiations. Preliminary case studies employed additional levels of theory (HF, MP2) and several other DFT functionals (B3LYP, SVWN, BP86). The harmonic frequency analysis was systematically performed to ensure the absence of imaginary frequencies for local minima.

Bond orders

The Mulliken bond population analysis, [START_REF] Mulliken | Electronic Population Analysis on LCAOMO Molecular Wave Functions[END_REF] the Wiberg [START_REF] Wiberg | Application of the pople-santry-segal CNDO method to the cyclopropylcarbinyl and cyclobutyl cation and to bicyclobutane[END_REF] and the Mayer 28,29 bond indices were calculated using the MultiWFN package. [START_REF] Lu | A multifunctional wavefunction analyzer[END_REF] The delocalization index (DI) [START_REF] Matito | Electron sharing indexes at the correlated level. Application to aromaticity calculations[END_REF][START_REF] Outeiral | Revitalizing the concept of bond order through delocalization measures in real space[END_REF] was obtained using AIMall. 33 We have also carried out an electron localization function (ELF) topological and population analysis, [START_REF] Becke | A simple measure of electron localization in atomic and molecular systems[END_REF][START_REF] Silvi | Classification of chemical bonds based on topological analysis of electron localization functions[END_REF] which is frequently adopted as a descriptor of atomic shells and covalent bonds using the TopMoD package. [START_REF] Noury | Computational tools for the electron localization function topological analysis[END_REF] For simplicity, the ELF approach description is reported in Supporting Information.

Stretching force constant

The bond stretching force constant k was also calculated as an estimator of the bond strength in order to be compared with the IBSI. However, it makes no sense to use the standard vibrational normal-mode analysis to derive a stretching force constant of a particular bond in a polyatomic molecule. Actually, a normal mode of vibration commonly couples to others, it is delocalized over more than the two atoms of interest. To derive a force constant k associated with a local vibrational mode between two given atoms X and Y, a numerical procedure was used in which the second derivative of the potential energy is numerically assessed (see ESI). A purely localized bond force constant k is however difficult to determine because the geometry and the ED distribution in the XY surrounding must be kept frozen as much as possible during the numerical procedure so that k is only inherent in the XY bond.

Two different procedures were considered. In the first one, the molecule is divided into two rigid fragments translating along the XY direction, leading to a force constant hereafter called k 1 . It cannot be computed for bonds belonging to a ring. In this first procedure, although non-bonding interaction initially present between separating fragments will contribute to a small extent to k 1 , k 1 is more intrinsic to the XY bond strength than the second procedure k 2 . In the second version, only X and Y are displaced in the molecular framework, leading to a force constant called k 2 . It is highly important to note that, in this second version, the fixed substituents attached to XY indirectly significantly contribute to the resulting k 2 via steric and electronic effects upon XY displacement. In the following, the k 1 and k 2 local force constants have been obtained and compared using the set1 of 184 molecules not containing rings (438 bonds). k 2 has also been determined for 51 ring-containing molecules (239 bonds, see supporting information). Anharmonic computations were performed for a set of diatomics using the program Gaussian. [START_REF] Frisch | Gaussian16 Revision B.01[END_REF] For each diatomic molecule, the X matrix representing the anharmonicity is made up of only one term given in cm -1 and involved in the correction of the potential. This term has been employed to compare the anharmonicity contribution across some diatomics.

Results and discussion

A systematic study has been performed on 235 molecules and 677 different bonds covering a broad range of chemical bondings (non-covalent bonding, metal coordination and covalent bonds, see ESI). When we say bond property, we naturally think of bond order. That is the reason why we have first addressed the possible link between the IBSI and bond orders. This preliminary study reported in Figure 1 in ESI demonstrates that the IBSI is clearly not correlated with the examined bond orders: Mayer (R 2 = 0.41), Wiberg (R 2 = 0.37), Mulliken (R 2 = 0.27), Delocalization Index (R 2 = 0.46), ELF (R 2 = 0.007).

From a chemical perspective, the inspection of the IBSI computed for H-X bonds across the periodic table undoubtedly shows a certain degree of consistency: for instance, it gradually increases from left to right along the second period of elements, from a low value for H-Li (IBSI=0.06) to the largest value for H-F (IBSI=1.63) (see Figure 2 in ESI). Hence, we subsequently looked into a correlation between the IBSI and the stretching force constant concept.

IBSI versus stretching force constant

Actually, we may suppose that the more the electron sharing between the two atoms at equilibrium geometry is pronounced, the bigger the force to stretch the bond. In a preliminary examination, we have compared the IBSI with the stretching force constant k 1 calculated for a series of 5 homonuclear diatomic molecules involving the diatomic oxygen. As can be seen on Figure 3, the linear correlation between these two properties is very strong (R 2 = 0.99). We know that placing electrons on antibonding orbitals will decrease the stability of a molecule. This is what happens across the series from O 2+ 2 to O 2- 2 . Accordingly, the conventional bond order calculated as 1 2 [(number of bonding electrons)-(number of antibonding electrons)] decreases. The lower this bond order, the weaker the bonding. The IBSI fully reflects this trend over a wide range of chemical bonds in this series, taking values between 3.95 for O 2+ 2 down to 0.77 for O 2- 2 . Although this set is limited to 5 molecules, the observed trend is noteworthy.

IBSI k 1 (mdyn/Å) O 2 2- O 2 • - 3 O 2 O 2 • + O 2 2+ R 2 = 0.9891 R 2 = 0.9998
This preparatory work has been extended to 438 bonds (set1). The result is reported on Figure 4 (k 1 versus IBSI). The obtained linear correlation is surprisingly high with 93% of the force constant variance being explained by the IBSI. A close inspection of the plot on Figure 4 (and also on Figure 3) reveals that the straight line is not the best shape to represent the data. An attempt to fit a quadratic function to the data significantly improves the In this work, 10 simple ionic compounds with large ionic character (like NaCl, KCl, LiH, LiF, LiCl, ...) have also been investigated. Such species are stable mainly because of the electrostatic attraction between atoms and the energy required to break the bond may be as large as for covalent bonds. However, simple ionic compounds exhibits low stretching force constant values compared to covalent bonds. [START_REF] Kaya | New insights from the relation between lattice energy and bond stretching force constant in simple ionic compounds[END_REF] Here, the corresponding IBSI values are predicted to be very small (< 0.15), below the IBSI covalent domain (0.15 < IBSI < 4.00).

This fits very well with the low value of the force constant k 1 (< 3 mdyn/ Å) calculated for these ionic species. But it does not reflect the strong bond dissociation energy expected for such very polar molecules with respect to dissociation either into the neutral atoms or into ions. Clearly, the IBSI is not related to the bond dissociation energy but rather to the local stretching force constant.

For the sake of convenience and to save CPU time, the 6-31G basis set has been used throughout this work for molecules made up of elements in rows 1-3. But the influence of adding diffuse functions has been addressed. Adding a diffuse function to this basis set (namely, 6-31+G ) does not change the k 1 /IBSI correlation. However, in absolute term, the most important changes are observed for those systems being weakly bonded, like for instant for Li-F, Mg-Cl in CH 3 MgCl, and certain non-covalent bondings (see ESI). In that case however, it is not attributable to the change in geometry. In other respects, it turns out that adding polarization functions to the LANL2DZ basis set (using either LANL2DZdp

for the main group elements and the hydrogen atom, or LANL2TZf for the metal) does not change the force constant k 1 /IBSI correlation. However, it has to be noticed that, actually, both IBSI and k 1 can experience significant variations, in absolute terms. This is notably observed for metal-ligand bonding (IBSI(Pt-Cl)=0.216-0.278, IBSI(Pt-N)=0.312-0.360). This is partly because of the geometry change, but not only. As a consequence, the LANL2DZ or LANL2TZ basis sets augmented with polarization functions, although significantly more CPU-expensive, are advised.

Differences between IBSI and the force constant k

Undeniably, the IBSI is strongly connected to the physically grounded bond strength concept.

The force constant k being the second derivative of the energy with respect to the internuclear distance at the equilibrium state, Å-2 appears in the unit of k. In comparison, although the IBSI is dimensionless (see Equation 1), its value is, by construction, inversely proportional to the square of the bond length d 2 , which is consistent with the unit of k. But, despite the good agreement between the IBSI and the force constant k, the fraction of the variance unexplained is 5% (quadratic fit). We can try to find reasons for that.

It is highly important to notice that the static IBSI calculation is fully inherent to the considered atom pair while the force constant k 1 accounts for side effects. Actually, the local force constant may have a significant dependency on the chemical environment surrounding the given bond XY. For instance, non-bonded interactions between adjacent substituents on the two neighboring carbon atoms of 1,2-dibromoethane impact to some degree the C-C stretching local force constant value. The effect of the bonding environment on the force constant k can be gauged by comparing the correlations k 1 /IBSI and k 2 /IBSI for the same set of bonds. Indeed, the alternative numerical differentiation scheme k 2 is less intrinsic to the studied bond than force constant k 1 (see Computational details section). As expected (see Figure 4, right panel), this second scheme k 2 yields a poorer quadratic correlation with the IBSI using the same set of molecules without rings (R 2 = 0.88 < 0.95 for k 1 ). This is not surprising since k 2 (atoms X and Y moving in the rigid molecular framework) is more sensitive to the substituents attached to the studied bond than k 1 (rigid fragmentation along the XY bond). This is further evidenced by investigating a second set of bonds XY exclusively located in rings (see ESI, Figure 4, k 2 =f(IBSI)). Stretching such bonds substantially disrupts the ED in the bond surrounding. The resulting determination coefficient dramatically drops to 0.74. This clearly demonstrates that the presence of ring constraints in the vicinity of the considered bond significantly affects the deduced force constant k 2 , but not the IBSI, hence deteriorating the IBSI-k 2 quadratic correlation. Conversely, for diatomic molecules XY where no environment effect can be invoked, the k 1 and k 2 schemes are identical and lead to a good quadratic correlation with the IBSI (R 2 = 0.95, See Figure 6).

In fact, the IBSI and force constant k are structurally different. The second derivative of the potential energy with respect to the geometric displacement corresponds to the curvature of this function. Hence, k indirectly depends on the total structure ED reorganization upon stretching the bond embedded in its molecular framework. In comparison, the static IGMδg pair approach exclusively focuses on ED between atoms X and Y. More precisely, it captures the ED gradient reorganization between two states: the real system and a non-interacting reference (the IGM), which are only different in the middle of the two atoms where the phenomenon of contragradience takes place (see Figure 2). Therefore, the IGM approach does not account for any ED reorganization outside the bond region. Moreover, this ED gradient reorganization (δg pair ) is estimated solely based on a single geometry ED, without considering any ED "breathing" between the two states, in contrast with the calculation of k. Another way of calculating a stretching force constant associated with just one bond length is to apply the Konkoli-Cremer adiabatic internal modes (KC-AIM). [START_REF] Konkoli | A New Way of Analysing Vibrational Spectra. I. Derivation of Adiabatic Internal Modes[END_REF][START_REF] Cremer | New developments in the analysis of vibrational spectra. On the use of adiabatic internal vibrational modes[END_REF] This approach is particularly well suited to obtain a stretching force constant k a not contaminated by coupling with other modes (in contrast with normal modes). As with the local k 1 and k 2 numerical schemes used in our work, the KC-AIM approach probes the bond strength in a non-destructive manner (it leaves the considered bond undissociated). However, the upon considering a perturbation to the considered bond length. Based on KC-AIM mode force constants k a taken from reference [START_REF] Zou | Relating normal vibrational modes to local vibrational modes with the help of an adiabatic connection scheme[END_REF] for 39 molecules and 145 bonds, an additional comparison has been reported between IBSI and k a (see ESI Figure 6). It shows a significant correlation between these two properties (R 2 = 0.92 -0.95), very similar to the one obtained between IBSI and local force constants k 1 and k 2 (see Figure 4) previously approximated numerically by a three-point central difference formula. But again, a small fraction of the variance remains unexplained. k a , k 1 and k 2 are dynamic descriptors: a local force constant reflects the potential energy curvature associated with the vibration within its chemical environment. For its part, the IBSI provides a means for assessing the intrinsic bond strength without being affected by adjacent atoms. It measures the strength of the attachment between two considered atoms through their degree of ED contragradience. It does not rely on estimating any potential curvature, and then does not depend on how neighboring atoms would affect the vibration. Rather, it quantifies the ED clash between two given atoms as compared to the contragradience-free system (virtual) having exactly the ED of the true system. That is the reason why the IBSI and force constant are conceptually different. The static bond descriptor IBSI can be utilized in any molecular situation, equilibrium geometry or even along a reaction path where the force constant is no longer connected to the restoring force concept.

Another reason can be given to explain differences between IBSI and the force constant k. The level of theory used could rise differences in the treatment of both k and the IBSI.

Actually, the IBSI depends on the ED and the geometry (d 2 ) and then, to a certain degree, on the level of theory chosen. Similarly, the second-order three-point central difference formula employed to approximate the second derivative k also depends, to a certain extent, on the level of theory employed. Based on case studies, the IBSI has been proven to be only slightly dependent on the level of theory (see ESI), and the DFT(M06-2X/6-31G**) level of theory has been chosen here. This popular functional is recommended for applications involving main-group elements and non-covalent interactions. [START_REF] Zhao | The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals[END_REF] Of course, small changes in ED and geometry might be expected upon using more sophisticated methods.

In other respects, avoided crossing, as observed in C 2 between states X 1 Σ + g and B 1 Σ + g 39 may affect the stretching force constant derived from the ground-state adiabatic potential energy curve. In the presence of such avoided crossing, the stretching force constant obtained using a single determinant wave function cannot serve as an accurate bond strength descriptor, in absolute terms. Neither does the IBSI, since such avoided crossing affects in the same way the electron density used to obtain the IBSI at the equilibrium geometry. Similarly, a DFT level of theory based on a single Slater determinant description cannot describe appropriately diatomics having multireference character. However, in our study, Figures 5 and8 (diatomic systems) aim at proving that IBSI and force constant, given a same level of theory (DFT here), are connected. An in-depth description of bonding of diatomics (like C 2 , 39 Mo 2 or W 2 [START_REF] Angeli | Ground states of the Mo 2 , W 2 , and CrMo molecules: A second and third order multireference perturbation theory study[END_REF] or F 2 41 ) would require highly accurate and multi-reference methods, very CPU-expensive approaches, beyond the scope of this work.

For the set of diatomic molecules, another factor has been investigated to explain the remaining difference between the IBSI and k 1 obtained within the harmonic approximation.

Anharmonic corrections to vibrational frequencies have been computed using the Gaussian program. [START_REF] Frisch | Gaussian16 Revision B.01[END_REF] However, diatomics exhibiting significant anharmonicity could not be associated to the most important deviations observed on the k 1 /IBSI graph.

In summary, the force constant k and IBSI, although both strongly correlated to the bond strength concept, do not rigourously express the same physical content. The IBSI is intrinsically linked to the atom pair considered and provides a very efficient way of measuring a local bond strength in molecular situation, without any bias due to the bonding environment, even in cyclic structures. In practice, the IBSI has the advantage of requiring a single point energy calculation.

IBSI: scope of application and limitations

The IGM-δg descriptor allows for the proper identification of ED contragradience domains in the molecule. The IBSI quantifies this contragradience in the region between two atoms and only considers the orbital contributions of these two atoms. Compared to the local stretching force constant k a obtained by the Konkoli-Cremer procedure, the IGM-IBSI scheme additionally provides information through the local distribution of ED contragradience not necessarily uniform nor located at the center of the atom pair. In most cases, in the space region in the middle of the considered pair, the ED contragradience dominantly arises from two atomic sources. Then most often, the ED contragradience domains are well separated in space, as illustrated in Figure 7 of the manuscript showing a TS structure. But, what if a third atom brings some ED gradient in the considered atom pair region ? When does this occur ? To illustrate this point, the interactions in the B 2 H 6 molecule have been identified

and characterized thanks to the full IGM approach with δg = ∇ ∇ ∇ρ IGM -|∇ ∇ ∇ρ|. As can be seen in the Figure 8 of ESI, two contragradient domains very slightly overlap around each H bridge atom. Then, one may wonder whether the B 1 -H bridge bond strength can be assessed based solely on atom pair contragradience. To which extent B 2 is really involved in the B 1 -H bridge contragradience domain (if it is) cannot be assessed within the current IGM approach. Even though this could be achieved, it would make little sense to keep a single geometrical parameter (B 1 -H bridge distance) in the IBSI definition in that case. The user must be aware that the IBSI approach relies on atom pairwise interactions. Despite this, the IBSI faithfully reflects the relative weakness of the B 1 -H bridge bond (IBSI=0,47) compared to the B 1 -H terminal bond (IBSI=0.55). To our knowledge, this kind of situation should be limited to very small ring systems with atoms having electrons more easily spread.

Examining the cyclopropane molecule (see Figure 9 in ESI), in addition to the well separated 2-center contragradience domains (C-C), a 3-center domain appears in the center of the ring.

Despite this, the IBSI correctly describes the C-C bond strength. Another kind of bonding deserves attention: when a ligand is bound to a metal via a conjugated π-system like in metallocenes. For illustration purpose, the interaction between the cyclopentadienyl anion (Cp) and the dication Zn 2+ has been reported in ESI (Figure 10). The IGM approach unambiguously reveals the region of interaction. The Cp...Zn 2+ contragradience domain takes the form of a continuous ring for this η 5 hapticity mode. The IBSI value obtained for each C...Zn 2+ pair is 0.194, which clearly reveals a bonding weaker than pure covalent bonds like C-C in Cp (IBSI=1.081). However, in such case, although the AIM analysis predicts a bond critical point between the carbon atom and Zn 2+ , the use of the pairwise IBSI scheme is tricky since there are more than two atoms simultaneously involved in the coordination to the transition metal.

The agostic interaction is another type of challenging interaction, occuring between an electron deficient metal and a σ-bond located close to the metal. [START_REF] Lein | Characterization of agostic interactions in theory and computation[END_REF] For illustrative purpose, the complex [Ti(Cl 2 )C 2 H 5 ] + has been reported in ESI (Figure 11). The IGM approach reveals this agostic interaction with the contragradience domain being located between C-H and titanium. The relative weakness of this interaction is correctly predicted by the IBSI calculation (IBSI=0.105), compared to the Cl-Ti coordination (IBSI=0.353). However, this agostic interaction involves three atoms while the IBSI scheme is intended for pairwise interactions.

To summarize, we are fully confident that the IBSI approach correctly describes the intrinsic bond strength between atom pairs in the vast majority of molecules, including transition state structures, non-covalent hydrogen bonds and halogen bonds, as long as contragradience domains are well separated. However, there are limitations in interpreting IBSI for interactions with complex bonding scenarios where more than two atoms are involved in a bonding situation. Metallocenes and even non-covalent stacking interaction, like in the benzene dimer, belong to this category. It is beyond the scope of this work to explore more in details such cases. However, whether the IBSI approach can be used is usually obvious by inspection of the molecular geometry and the shape of the contragradience isosurfaces.

Bridge between local ED descriptor and bond strength

To our knowledge, only one previous study reported such a bridge between a local EDbased descriptor and a bond stretching force constant. Very recently, considering a set of 53 molecules, a strong linear correlation (R 2 = 0.97) was presented between the localized stretching force constant of a given bond and the potential energy density distribution integrated over the surface between the two Bader's atomic basins. [START_REF] Ananyev | Toward a Rigorous Definition of a Strength of Any Interaction Between Baders Atomic Basins[END_REF] In a similar way as in our study, the authors divided the resulting integral by the internuclear separation. However, this approach needs to characterize the ED topology (by determining critical points and basins, ...) prior to any calculation. In contrast, a compelling feature of the IGM model is to provide an automatic workflow that delivers the IBSI without any ED topology analysis.

We believe that the IBSI is a powerful and robust tool for interpretation accessible to a wide community of chemists (organic, inorganic chemistry, including transition metal complexes and reaction mechanisms).

A typical application of IBSI

A review of possible applications of the IBSI is beyond the scope of this study. However, in order to show the possibilities offered by the IBSI compared to the Mayer bond index we have reported on Figure 7 a transition state structure (TS) taken from a recent study on the reactivity of thiolates with cisplatin. [START_REF] Minervini | Fate of cisplatin and its main hydrolysed forms in the presence of thiolates: a comprehensive computational and experimental study[END_REF] The reactant is the complex [Pt(NH 3 ) 2 (H 2 O)(cystein)].

The mechanism proceeds with the intramolecular substitution of the water molecule by the -CO - 2 group of the cystein, which finally acts as a (S,O)-bidentate chelating agent in the final product.

Figure 7: Bond Strength Index (IBSI) and 0.045 δg isosurfaces for selected bonds in a transition state involving a platinum complex (M06-2X LANL2DZ level of theory); isosurfaces colored according to the BGR scheme over the range -0.2 < sign(λ 2 )ρ < 0.2 a.u.; Mayer bond order in parentheses.

According to the Mayer analysis (MBO), the C-S and S-Pt bonds exhibit an identical bond multiplicity of about 1.0. The IBSI provides an additional information: the covalent bond C-S is stronger (IBSI=0.507) than the S-Pt coordinate bond (IBSI=0.239). This is consistent with the strength of metal ligand bonds known to be intermediate between that of covalent bonds and non-bonded interactions. In other respects, the MBO predicts a formal half-bond between the nitrogen and platinum atoms (0.507 pair) and a single bond between the sulfur and platinum atoms (1.023 pair). It may be then naively expected that the corresponding metal -ligand bond strengths should correlate with these bond orders at a ratio of 1(N-Pt):2(S-Pt). However, the IBSI predicts an inverse strength ratio with the N-Pt bond being significantly stronger (IBSI=0.331) than the S-Pt one (0.239).

This TS has one imaginary vibrational frequency with a negative force constant no longer connected to the restoring force concept. Nevertheless, the IBSI stays attached to this notion, allowing for internally probing bond strength in TS structures and along reaction pathways.

In the course of this reaction, the breaking Pt...OH 2 and forming -CO - 2 ...Pt chemical bonds undergo electronic structure changes. As expected, the MBO analysis at the TS shows that very few electron is shared between these reacting atoms (around 0.2 pair), in accordance with the small degree of orbital overlap achieved between the corresponding atoms at the TS. This decreased orbital overlap also weakens the two corresponding bonds, in the IBSI range 0.10-0.12: at the boundary between non-covalent bonding and coordinate covalent bonds (see Figure 5).

Usually, in hydrogen bonding, little or no electron pair is shared by the valence shells of the two atoms, making the interpretation of the Mayer bond order generally quite difficult.

Here however, the MBO detects an uncommonly large electron sharing (0.210) for NH...O compared with OH...O (0.039) in this TS. In a complementary way, the IBSI analysis reveals that the picture of the hydrogen bond is deceptively simple for NH...O hydrogen. Actually, its bond strength (0.178) is outside the known H-bonding range (generally ≤ 0.15 on the IBSI scale, see Figure 5) and more than three times larger than the OH...O one (0.053). This is because the former exhibits a short distance (1.6 Å) and a strong acid/base character between the two charged entities -NH + 3 and -CO - 2 . In light of this short example, we see that the IBSI and MBO analyses provide a complementary information to explore chemical bonding and electronic effects.

Conclusions

In this work, we present for the first time the Intrinsic Bond Strength Index (IBSI). It is a post-processing tool very efficient to internally probe the strength of a given chemical bond in molecular situation. This quantitative tool complements the conventional visual usage of the IGM approach which provides colored isosurfaces highlighting the spatial localization of interactions in molecules. The theoretical advantage of the IGM method rests on the definition of a non-interacting reference system (the so-called Independent Gradient Model).

This feature allows the quantitative determination of bond strength. In this regard, to the best of our knowledge, the IBSI is distinct from other ED-based approaches like SEDD, DORI or NCI. Based on the large range of compounds explored in this work (235 species and 677 atom pairs), we suggest to approximately distinguish between weak interactions: IBSI 0.15 , coordinate covalent bond: 0.15 IBSI 0.6, and covalent bonding: IBSI 0.15 (these limits being only indicative).

The IBSI is defined in a parameter-free manner and its implementation is very easy: it only requires a wavefunction file as input. It has been shown to be little dependent on the level of theory and basis set. From the practical point of view, an attractive feature of the IGM-IBSI approach is to provide an automatic workflow that delivers the bond strength index without any prior characterization of the ED topology (unlike within the AIM theory).

In other words, it is not an ED basin-based method like ELF or AIM.

From a theoretical perspective, the current work provides a quantitative bridge between a local electron density-based descriptor and the physically grounded bond strength concept.

Actually, the IBSI emerges from the IGM formulation of electron sharing between two atoms.

This index results from the integration of the local descriptor IGM-δg pair bond signature. In contrast, a local force constant calculation may depend on the bonding environment and the electron density reorganization occurring during bond stretching. Furthermore, along reaction paths, even in regions with imaginary vibrational frequencies where the force constant is no longer connected to the restoring force concept, the IBSI stays attached to this notion, extending its scope beyond equilibrium geometries.

Clearly, this index covers a broad range of bonding situations and is not limited to closely related molecules. This makes it suitable for following bond formation and breaking in chemical reactions. In particular, the IBSI turns out to be a suitable indicator for predicting a realistic intrinsic strength of a given chemical bond in transition states. Changes in stretching force constants upon bond polarizing of a reaction partner is a key element of the unified reaction valley approach [START_REF] Cremer | New developments in the analysis of vibrational spectra. On the use of adiabatic internal vibrational modes[END_REF][START_REF] Cremer | From Molecular Vibrations to Bonding, Chemical Reactions, and Reaction Mechanism[END_REF][START_REF] Zou | A New Method for Describing the Mechanism of a Chemical Reaction Based on the Unified Reaction Valley Approach[END_REF] used for the detailed mechanistic analysis of elementary chemical reactions. Hence, changes in IBSI occurring in the course of a chemical reaction might contribute to gain insights into reaction phases each of which leads to a change in chemical bonding. More generally, this index has the advantage of being applicable to diverse fields in chemistry and future work will be devoted to applications for systems in organic and inorganic chemistry, transition metal complexes and reaction mechanisms.

The next release of IGMPlot including the calculation of the IBSI will be made available at http://igmplot.univ-reims.fr.

Figure 2 :

 2 Figure 2: Schematic representation of electron density interpenetration (contragradience) between two fragment sources.

Figure 3 :

 3 Figure 3: Relationship between the stretching force constant (scheme k 1 ) and the IBSI for 5 diatomic oxygen species; fitted linear (full line) and quadratic (dashed line) regression curves are reported.

Figure 4 :

 4 Figure 4: Relationship between the stretching force constant and the IBSI for the set1 of 438 bonds in 184 molecules (with no rings since k 1 would not be calculable in that case) for (a) scheme k 1 and (b) scheme k 2 ; fitted linear and quadratic regression curves are reported. The points are colored using the color scale given in Figure5.

Figure 5 :

 5 Figure 5: IBSI scale obtained on the set of 677 studied bonds. Limits are only indicative for the set of 677 bonds examined here.

Figure 6 :

 6 Figure 6: Relationship between the stretching force constant and the IBSI for 36 diatomics; fitted linear and quadratic regression curves are reported. The points are colored using the color scale given in Figure5.

  Through a quadratic relationship, a tight link (R 2 = 0.95) has been established between the IBSI and the localized stretching force constant of a chemical bond. This result underpins the IGM formalism and its two key-components: the non-interacting reference definition and the gradient based partition. Future work could also consider to seek a correlation between the IBSI and ED characteristics at bond critical points. It has however to be noticed that by construction the δg descriptor is limited to theoretical ED models, which rules out the use of experimental data. Differences exist between the IBSI and a local force constant. It is shown that the IBSI relates to the intrinsic bond strength. It is obtained by focusing only on the atom pair considered. No expensive calculations of energy second derivatives are needed to obtain the IBSI.

  Actually, the relative C-C bond strengths estimated with the local KC-AIM force constant in propane (k a = 4.066 mdyn/ Å) and cyclopropane (4.137 mdyn/ Å) fully support the IBSI rating (0.839 and 0.892, respectively; see ESI in section Local KC-AIM force constants versus IBSI).The IBSI approach is not restricted to well behaved 2-center bonding systems. Delocalized aromatic systems are suitably described. For instance, the C-C intrinsic bond strength in C 6 H 6 (IBSI = 1.197) is found in between that of a single bond (IBSI = 0.833 in C 6 H 12 ) and a double bond (IBSI = 1.396 in C 2 H 4 ), as supported by the local KC-AIM C-C force constants k a : 6.600 mdyn/ Å in C 6 H 6 , in between 3.923 mdyn/ Å in cyclohexane and 9.911 mdyn/ Å in ethene (see ESI). The ability of the IBSI approach to properly describe the bond strength in a 3-center 2 π-electron system is also observed with the allyl cation (IBSI C-C = 1.169). In other respects, bonding having a 3-center 4-electron character as in certain hypervalent molecules[START_REF] Oliveira | Pushing 3c-4e Bonds to the limit: A Coupled Cluster Study of Stepwise Fluorination of First-Row Atoms[END_REF] or in systems exhibiting hydrogen or halogen bonds 43 is also well described since the contragradience domains are well separated in such cases. For instance, the F-Cl...NH 3 system exhibits IBSI values of 0.691 (F-Cl) and 0.101 (Cl...N), as supported by the relative local KC-AIM force constants k a : 2.687 mdyn/ Å (F-Cl) and 0.311 mdyn/ Å (Cl...N).[START_REF] Oliveira | The Intrinsic Strength of the Halogen Bond: Electrostatic and Covalent Contributions Described by Coupled Cluster Theory[END_REF]

Acknowledgement

This work was supported by a Grant (J. Klein, Master's thesis) from the Université de Reims Champagne-Ardenne and the ICMR laboratory (Institut de Chimie Moléculaire de Reims).

We thank the MaSCA (Maison de la Simulation de Champagne-Ardenne, France) for various supports including computing facilities (http://romeo.univ-reims.fr) and the CRIANN computational center (http://www.criann.fr) for additional support. We thank Prof. F. Bohr for helpful discussions.

Supporting Information Available

Detailed molecule sets used throughout this study and computed properties (IBSI, k a , k 1 , k 2 , bond orders). Investigation of the basis set and method dependence on the IBSI. This material is available free of charge via the Internet at http://pubs.acs.org/.


domain: spi info

text:
Khalil Aloui

Moncef Hammadi
email: moncef.hammadi@supmeca.fr

Amir Guizani
email: amir.guizani@live.fr

Thierry Soriano
email: thierry.soriano@supmeca.fr

Mohamed Haddar
email: mohamed.haddar@enis.rnu.tn

Modeling and Simulation of A Swarm Robot Application Using MBSE Method and Multi-Agent Technology: Monitoring Oil Spills

Keywords: Swarm robotic, Modeling, MBSE method, Simulation, Multi-agent technology, Sustainable energy

Swarm Robotics is a new approach to the coordination of a large number of robots inspired by nature. This approach aims to design collective behaviors for a large number of robots. Unfortunately, several researchers have tried to develop structured design methods but these methods are still limited. Today, swarm robotics are used in many fields, they are used in agriculture, medicine, industrial, etc. One of the most important fields that require swarm robots is surveillance. In this paper, we present in the first section some methods of designing swarm robot systems by identifying swarm engineering based on a model (MBSE) and multi-agent simulation. Then, we study the energy problem of these robotic systems and the solution proposed by the researchers. In the second section, we present an application for detecting oil in the sea and cleaning it using swarm robots. We will model this application using the MBSE method with the SysML language. Finally, we will simulate the model on a multi-agent tool to identify the functional and structural architecture of the system.

Introduction

Swarm robotics is a new approach to coordinate a group of robots inspired by social insects (Sachin, 2004). It aims to design a scalable, robust and flexible collective behaviors for the coordination of a large group of robots. This approach is often inspired from social animals such as bees, ants, fish and birds [START_REF] Dorigo | Swarm intelligence[END_REF]. Indeed, swarm robotics tends to develop a robotics system that are similar to behaviors of social animals. In particular, swarm robotics systems are meant to be robust, scalable and flexible [START_REF] Camazine | Self-organization in biological systems[END_REF]. Scalability can be defined as the ability to work with different group sizes and generate solutions for different tasks. In swarm robotics, scalability is promoted by local detection and communication. Robustness can be defined as the ability to operate despite failures in individuals. In swarm robotics, robustness is the redundancy and the absence of a leader. Flexibility can be defined as the ability to operate under a wide range of different environments and tasks. The simplicity of behaviors and mechanisms represent flexibility in swarm robotics.

A swarm robotic system is characterized by the autonomy of its members. These robots must be homogeneous. In some cases, there may be different types of robots in the swarm, but these groups should not be too large. In addition, the number of robots in the swarm must be large. These robots must have local communication and detection capabilities. Ultimately, robots must collaborate to achieve or improve the performance of a primary task that cannot be solved by a single robot since the detection range of a group of robots is wider than that of a single robot. [START_REF] Navarro | An introduction to swarm robotics[END_REF]. Groups of robots can make decomposable tasks more efficient when they use parallelism. Moreover, a single robot is unable to operate in different places at the same time while a group of robots can do that. Another advantage is fault tolerance (i.e. thanks to the redundancy of the system, the task remains accomplished despite the failure of a single robot within a group.) In addition, swarm robot systems have several drawbacks, such as uncertainty (i.e., robots can compete rather than cooperate if they do not know what other robots are doing). Due to collisions, robots can interfere with each other in a group. Finally, the use of a group of robots in a swarm can increase the overall cost of the system (Ronald Arkin et al, 1998).

Several real-world issues require the application of swarm robots, such as monitoring an environmental accident. A swarm robot system can instantly monitor and detect hazardous events, such as a chemical spill due to its distributed detection capability. A swarm robot system can focus on one problem. First, the swarm identifies the nature of the problem. It mobilizes its members towards the source of the problem. Second, the swarm forms a patch that blocks the leak using its self-assembly characteristic (Sachin, 2004).

Depending on the given mission, the swarm robot systems can increase or decrease their number of individual robots. For example, during an oil spill, a swarm robot system will be built to contain an initial spill in the area. We will study this case in the rest of this paper.

Design methods

The design of the swarm robot system means the planning and development of a swarm system from the initial specifications and requirements. Unfortunately, there is no precise and structured methodology for designing desired collective behavior through individual behaviors. The main ingredient in developing these types of systems is always the designer's intuition.

Traditionally, many designers have used the code-and-fix approach to develop robot swarms. It is a bottom-up method means that the developer creates and improves the individual behaviors of the robots until the desired collective behavior is achieved. [START_REF] Kazadi | Artificial physics, swarm engineering, and the hamiltonian method[END_REF] applied the swarm engineering methodology called the "Hamiltonian method" to the problem of artificial physics. This methodology is used to create swarms of predefined global properties. They apply this methodology to the problem of artificial physics which locally creates hexagonal grids of agents. [START_REF] Brambilla | Property-driven design for robot swarms: A design method based on prescriptive modeling and model checking[END_REF] presented a new top-down design method for swarms of robots based on prescriptive modeling and model verification called "Property Driven Design". This method consists of four phases: first, the designer specifies the requirements of the robot swarm. Second, the designer creates a normative model of the swarm. He uses model verification to verify that this normative model satisfies the desired properties; third, the designer simulates the model and implement the simulated version to validate the developed model; fourth, the designer implements the desired robot swarm to validate the previous steps. (Aloui et al. 2020) improved the "Property Driven Design" design method. They used the SysML language for the specification and modeling of the Swarm robotic system. They used SysML state machine diagrams to describe robot behaviors. Then, the robot behavior models described with SysML are implemented with a multi-agent technique [START_REF] Guizani | Multi-agent approach based on a design process for the optimization of mechatronic systems[END_REF].

Model Based Swarm Engineering MBSE

In general, the model is a simplified version of a concept, relation or phenomenon to facilitate understanding of a structure or system. In the field of swarm robotics, it is possible to use the MBSE method (Model Based Systems Engineering) to model swarm robot systems to understand behavioral analysis, performance analysis, requirements traceability, system architecture, simulation, testing, etc. (Mhenni et el, 2014). The main advantage of a model-based approach is to determine the scope of the problem. That is, the developed robot swarm best meets expectations when the properties are clear and complete with the models. The detailed description of the requirements by the models reduces the risk of developing "the wrong robot swarm". To ensure model-based swarm engineering, the swarm developer can use the graphical modeling language SysML. SysML is a visual modeling language that supports the specification, design, analysis, verification, and validation of complex systems. SysML is used with different methodologies including object orientation, structured analysis, and others [START_REF] Hause | The SysML modelling language[END_REF]. This language includes diagrams that can be used to specify system requirements, system behavior, system structure, and parametric relationships.

Multi-agent simulation MAS

Recently, several researchers have started to use multi-agent systems in different fields, such as automobile design (Guizani et al, 2014a), power electronics [START_REF] Hammadi | Layout optimization of power modules using a sequentially coupled approach[END_REF]) and aerospace design. A multi-agent system is a field of artificial intelligence (AI). [START_REF] Guizani | Agent-based approach for collaborative distributed mechatronic design[END_REF] are focused on the agent and his autonomy. They defined the agent as a computer system located in its environment and able to act autonomously in order to achieve its design goals. The study of swarm robotic systems has been considered today as a study of multi-agent systems. Developers of swarm robot systems use realistic simulators to test and accelerate the development of new design methods. They use simulators to model the interactions between robots and the interactions of robots with their environment (Sahin et al, 2008). Luke designed a multi-agent simulation toolkit called MASON to present a large number of multi-agent simulation tasks ranging from machine learning and swarm robotics to environments of social complexity [START_REF] Luke | Mason: A multiagent simulation environment[END_REF]. Anylogic is another multi-agent simulation tool developed by The AnyLogic Company. It is a graphical modeling language that makes it easy to extend the simulation model with Java code (Borshchev et al, 2002).

Sustainable Energy Efficiency of Swarm Robots

Each time the groups of robots in the swarm recharge to their maximum level. Swarm robot systems do not estimate the energy needed for their next task. They require a method to minimize the overall energy cost in a swarm and simultaneously maximize the performance of the swarm. However, only a few articles in the current literature consider energy efficiency in a swarm. [START_REF] Labella | Efficiency and task allocation in prey retrieval[END_REF] have developed an adaptation method which manages the number of robots foraging in the environment. Indeed, they efficiently use energy when needed, but they measure only the time needed to complete the mission without measuring the amount of energy already allocated. (Liu et al, 2007) focused on energy efficiency. They modify the foraging time for each robot. [START_REF] Stirling | Energy-efficient indoor search by swarms of simulated flying robots without global information[END_REF] have worked on flying swarms in indoor environments. They proposed a new approach to find energy-efficient search algorithms. They noticed that the progressive launch of the robots increases the search time and reduces the total energy consumption and the collision rate. [START_REF] Chen | Maximizing Energy Battery Efficiency in Swarm Robotics[END_REF] propose a method that improves the efficiency of energy management methods. This method uses upper and lower battery thresholds.

In the next section, we will model a swarm robot application for oil spill monitoring using the different SysML diagrams. Then, we will simulate this application on a multi-agent tool.

Case study: Oil spill monitoring

An oil spill in the areas has both immediate and long-term effects on the environment. The fauna of these areas can be directly exposed: either by swimming in the oil, by inhaling its toxic fumes, or by swallowing food affected by the oil. This exposure to petroleum can cause serious health problems and even lead to death. Thus, researchers are still trying to find solutions to clean up the affected environment (Sachin, 2004). We use the different SysML diagrams to model this application.

Requirements specification

The requirements diagram shown in Figure 1 specifies the application requirements of detecting environmental pollutants and cleaning them. Detecting environmental pollutants requires exploring the entire environment by deploying all the robots that detect and clean these pollutants. Cleaning should be quick and efficient using autonomous robots.

Functional architecture

In this step, we describe the functional scenario using the different collective behaviors of swarm robots. Figure 2 shows this scenario. First, members of the swarm begin to explore the entire region. Then they detect environmental pollutants. Finally, they regroup in these areas, form groups and clean up.

Fig. 2 Descriptive model for cleaning an oil leak

To model this case study, some collective behaviors of swarm robots are used:

• Collective exploration: The swarm of robots cooperatively navigates the environment in order to explore and monitor it. First, the swarm is initialized with random positions and constant speeds. If the robots locate spilled oils, they send positive information to the control center. If they detect obstacles or other robots within their reach, they avoid the collision. Otherwise, they send negative information to the control center indicating that there is no oil leak in that area. Subsequently, they update their speeds.

• Collective localization: Allows the robots in the swarm to find their position and orientation in relation to each other by establishing a local coordinate system throughout the swarm.

• Aggregation: Individual robots gather spatially in a specific region of the environment. This behavior is used to ensure that the robots are assembled in the area to be cleaned.

• Task allocation: Assignment of tasks dynamically assigned emerging tasks to individual robots in the swarm. Its objective is to maximize the performance of the entire swarm system. In this case study, there are 4 groups of robots distributed throughout the area to completely clean the contaminated area.

After identifying the functional model of the system, we simulate the model using AnyLogic.

Swarm simulation

In this step, AnyLogic is used as a multi-agent software tool to implement the swarm simulation model. In this case study, we chose 100 robots to do this mission. As shown in Figure 4, we have divided the region into four zones according to the degree of contamination: zone A and D are the two most polluted zones. So, the swarm of robots must be directed towards them more. Zone B is the moderately polluted zone and finally, Zone E is the smallest contaminated zone. Figure 3 gives the initial number of robots in each zone.

Fig. 3 Initial distribution of robots in zones

Each time, the swarm of robots must return to the boat to dispose of the captured pollutants. The initial speed of each robot is 20 km / hour. In fact, we have the right to direct the swarm by adjusting the probability values of the robots passing from the boat to different areas. The exploration time for the whole area is 1 hour, the service time in each area is 50 minutes, and the time to dispose of the pollutants in the boat is 30 minutes. In addition, the robots must form swarms of 10 robots in each area. At first, robots begin to explore the entire region to identify contaminated areas. They then regroup in these areas, form groups and do the cleaning. The red curve represents the number of robots that clean area A. This number is always greater than 15 robots. Likewise, the green curve represents the number of robots that clean the area D. This number is always greater than 10 robots. This means that the characteristics specified initially have been validated, (i.e., the robotic swarm is already directed towards the two most polluted regions). In addition, the robotic swarm treated the two less polluted areas B and E with the fewest robots. The yellow curve indicates the number of robots that return to the boat and throw the oil into the storage tanks. Because the oil diffusion is very rapid, the most important main step is to contain the spill to as small an area as possible. Each robot swarm was made up of 10 robots that scatter throughout the region and soak up the spilled oil. In the block definition diagram shown in Figure 7 the individual robots that form the swarm is described. Indeed, each robot consists of a movement system that ensures the robot swimming in the sea, a power system, and a communication and control system. At this level, the various sensors and electronic components are defined to design the system.

Fig. 7 Block definition diagram of a swarm of robots

Finally, it is necessary to verify that the components of the system described in the BDD correspond to the various functions of the system. The allocation matrix shown in Figure 8 relates each component of the system to its functions. Finally, the swarm developer chooses the implantation system (like the ROS platform) to implement this system on real robots. Researchers can realize the task of implementation in further research and next steps.

Conclusion

Swarm robotics is a new approach to robot coordination inspired by social insects. This approach aims to design evolving, robust, and flexible collective behaviors for a large number of robots. Unfortunately, the design methods for these systems are still limited. Several researchers have proposed approaches to simplify design techniques, but this field is still limited.

In this paper, we have described swarm robotics: their main definition, their fields of application, their design methods and their sustainable energy efficiency. Then, we presented the modeling and simulation methods of multi-agent systems. We modeled an oil spill detection application using the MBSE method with SysML diagrams, and finally, we simulated this application on a multi-agent tool to determine the functional and structural architecture of the robot swarm. The implementation of real robots and the production of prototypes can be the subject of future papers.

Références

Aloui, K., Hammadi, M., Soriano, T., Guizani, A., & Haddar, M. (2020). On the continuity of the swarm robot design using MBSE method and simulation. Borshchev, A., Karpov, Y., & Kharitonov, V. (2002). Distributed simulation of hybrid systems with AnyLogic and HLA. Future Generation Computer Systems, 18, 829-839.

Fig. 1

 1 Fig. 1 Requirements diagram

Fig. 4

 4 Fig. 4 Simulation of an oil spill cleanup in the ocean (a): Robots explore the entire region. (b) Robots clean the different areas Figure 5 represents the indicators of system performance:The red curve represents the number of robots that clean area A. This number is always greater than 15 robots. Likewise, the green curve represents the number of robots that clean the area D. This number is always greater than 10 robots. This means that the characteristics specified initially have been validated, (i.e., the robotic swarm is already directed towards the two most polluted regions). In addition, the robotic swarm treated the two less polluted areas B and E with the fewest robots. The yellow curve indicates the number of robots that return to the boat and throw the oil into the storage tanks. Because the oil diffusion is very rapid, the most important main step is to contain the spill to as small an area as possible. Each robot swarm was made up of 10 robots that scatter throughout the region and soak up the spilled oil.

Fig

  Fig. 5 System performance indicators

Fig. 6

 6 Fig. 6 The swarm architecture

Fig. 8

 8 Fig. 8 Allocation matrix Component-Function

  13th International Conference on Modelling, Optimization and Simulation (MOSIM'20),. Arkin, R. C., Arkin, R. C., & others. (1998). Behavior-based robotics. MIT press.