Andi Cuko

Antoni Maciá

Monica Calatayud

Stefan T Bromley
email: s.bromley@ub.edu

Global optimisation of hydroxylated silica clusters: A cascade Monte Carlo Basin Hopping approach

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Nanosized silica particles are widely used during the preparation of materials (e.g. cement-based materials, [START_REF] Aggarwal | Use of nano-silica in cement based materials -A review[END_REF] coatings, [START_REF] Zhoua | The change of the properties of acrylic-based polyurethane via addition of nano-silica[END_REF] polymers [START_REF] Qu | Tailoring and probing particlepolymer interactions in PMMA/silica nanocomposites[END_REF]) to significantly enhance their physical and mechanical properties. In addition to their high intrinsic strength (e.g. the Si-O bond 30% being stronger than the C-C bond), the success of silica nanoparticles for such applications is also due them having a relatively high proportion of reactive surface groups which can strongly interact with host materials. Silica is particularly prone to reaction with water, commonly resulting in a surface coverage of hydroxyl (Si-OH) groups. [START_REF] Zhuravlev | The surface chemistry of amorphous silica. Zhuravlev model[END_REF] The degree of hydroxylation is found to determine many of the properties of nanosilica. In solution, for example, the pH determines the solubility/hydroxylation of small silicate nanoparticles. Such species are central to the synthesis of technologically important nanoporous silicate materials such as zeolites [START_REF] Cundya | The hydrothermal synthesis of zeolites: Precursors, intermediates and reaction mechanism[END_REF] and occur in (bio)mineral nucleation, growth and dissolution processes. [START_REF] Belton | An overview of the fundamentals of the chemistry of silica with relevance to biosilicification and technological advances[END_REF] The proportion of hydroxyls on small silica particles can also be deliberately increasing by physico-chemical processing to produce enhanced properties for applications (e.g. epoxidation catalysis [START_REF] Chandra | Silica microspheres containing high density surface hydroxyl groups as efficient epoxidation catalysts[END_REF], mechanically robust antireflective films [START_REF] Wang | Fabrication of mechanically robust antireflective films using silica nanoparticles with enhanced surface hydroxyl groups[END_REF]). The density of surface hydroxyls has further been reported to influence the biocompatibility and toxicity of silica nanoparticles. [START_REF] Marucco | Hydroxyl density affects the interaction of fibrinogen with silica nanoparticles at physiological concentration[END_REF][START_REF] Zhang | Processing Pathway Dependence of Amorphous Silica Nanoparticle Toxicity: Colloidal vs Pyrolytic[END_REF] From a computational perspective the reaction of water with silica surfaces has attracted much attention. [START_REF] Yang | Reaction of water on silica surfaces[END_REF] For atomic to nanosized hydroxylated silica systems studies have tended to either focus on the detailed thermodynamics of the condensation of small ring-or chain-like (SiO2)M N ≤ 6 oligomers in solution [START_REF] Mora-Fonz | Oligomerization and Cyclization Processes in the Nucleation of Microporous Silicas[END_REF][START_REF] Trinh | Mechanism of Oligomerization Reactions of Silica[END_REF] using ab initio methods, or on larger scale (i.e. with 100s or 1000s of atoms) classical forcefield simulations of hydroxylated nanosilca. [START_REF] Makimura | Combined modeling and experimental studies of hydroxylated silica nanoparticles[END_REF][START_REF] Yeon | ReaxFF Molecular Dynamics Simulations of Hydroxylation Kinetics for Amorphous and Nano-Silica Structure, and Its Relations with Atomic Strain Energy[END_REF][START_REF] Bhattacharya | Molecular Dynamics Simulation Study of Growth Regimes during Polycondensation of Silicic Acid: from Silica Nanoparticles to Porous Gels[END_REF] Between these extremes, some recent studies have focussed on attempting to find the most stable structures of hydroxylated silica nanoclusters in the approximate size range of 10-100 atoms. [START_REF] Bromley | Novel structures and energy spectra of hydroxylated (SiO2)8-based clusters: searching for the magic (SiO2)8O2H3 cluster[END_REF][START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF][START_REF] Flikkema | Structure and energetics of hydroxylated silica clusters, (SiO2)M(H2O)N, M = 8, 16 and N= 1 -4: A global optimisation study[END_REF][START_REF] Jelfs | Evidence for atomic mixing via multiple intermediates during the dynamic interconversion of silicate oligomers in solution[END_REF] This intermediate size regime is characterised by having non-trivial (i.e. not simply rings or chains) structures which are also not totally amorphous as typically assumed for larger nanoparticles. Hydroxylated silica nanoparticles with such sizes can thus be thought of as species which represent a structurally complex molecular-bulk crossover regime. Such species are likely to be particularly important in understanding nucleation of complex silica solids (e.g. zeolites). Although in such processes, many factors play an important role (e.g. pH, templates) [START_REF] Cundya | The hydrothermal synthesis of zeolites: Precursors, intermediates and reaction mechanism[END_REF] knowing the lowest energy hydroxylated nanosilica species in the absence of such factors provides a baseline to assess how and when they are important. Previous theoretical work in this size regime has specifically investigated hydroxylated silica clusters (SiO2)M•(H2O)N with sizes M = 4, 8, 16, 24, and each with a variable number of incorporated water molecules (N = 1, 2, 3…) [START_REF] Bromley | Novel structures and energy spectra of hydroxylated (SiO2)8-based clusters: searching for the magic (SiO2)8O2H3 cluster[END_REF][START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF][START_REF] Flikkema | Structure and energetics of hydroxylated silica clusters, (SiO2)M(H2O)N, M = 8, 16 and N= 1 -4: A global optimisation study[END_REF][START_REF] Jelfs | Evidence for atomic mixing via multiple intermediates during the dynamic interconversion of silicate oligomers in solution[END_REF] Due to their high structural complexity, global optimisation was employed in these studies to find low energy isomers for each stoichiometry. Specifically, a fairly simple set of interatomic potentials was used in conjunction with the Monte Carlo basin hopping (MCBH) [21] approach to first obtain initial cluster structures. Due to approximations in the IP set employed these non-refined structures provide only a poor account of OH•••OH hydrogen bonding on the surfaces of the hydroxylated nanoclusters. [START_REF] Flikkema | Structure and energetics of hydroxylated silica clusters, (SiO2)M(H2O)N, M = 8, 16 and N= 1 -4: A global optimisation study[END_REF] Optimisations using relatively computationally expensive density functional theory (DFT) calculations on a selection of IP-MCBH-generated structures were used to provide final clusters exhibiting hydrogen bonding.

In this work we report on a MCBH global optimization approach to hydroxylated silica nanosystems employing a cascade local optimization method. First, our method is shown to be reliable by applying it to the previously studied (SiO2)8•(H2O)N system. [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF][START_REF] Flikkema | Structure and energetics of hydroxylated silica clusters, (SiO2)M(H2O)N, M = 8, 16 and N= 1 -4: A global optimisation study[END_REF] We then, generate global minima candidates for the as yet unreported systems: (SiO2)M•(H2O)N, M = 6, 10, 12, each for a range of N values. In previous work it was shown that the hydroxylation reaction energy for silica clusters decreases by addition water molecules until it levels off around an optimal number N. This optimal number was found to be cluster size dependent to follow an inverse power law of the form 2 ⁄ . Using our new global minima candidates for (SiO2)M•(H2O)N, M = 6, 10, 12 we also further probe the generality of this relation.

Methodology

In this work we apply the MCBH global optimization technique to find low energy global minima candidates for the cluster set (SiO2)M•(H2O)N, M = 6, 8, 10, 12 for a range of N. The general idea behind this technique in based on repeated steps consisting of random structural distortion and a subsequent local optimization. The Metropolis criterion [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF] is applied in order to accept or refuse the local optimized structure according to a cost in passing from the initial (before the random distortion) and the final state, and a fictitious temperature. Usually the cost is the energy difference between the two states but in principle it is possible to generate other cost criteria based on system properties. Within a certain (unknown) number of the cyclic repetitions of this procedure is possible to explore the potential energy surface (PES) of the system. At each step, the structure can be locally optimized with both ab initio or a classical IP. In the first case, the applicability is limited to relatively small systems due to the high computational expense. This expense not only originates from local optimizations of the system but also from the higher number of steps required to properly sample the PES. Indeed, in principle, the number of energy minima is exponentially proportional to the number of atoms. A typical approach to makes MCBH calculations more computationally amenable is to use IPs for efficient optimisation. Nevertheless, IPs alone, often, cannot be relied upon for an accurate description of a system's properties. Therefore, after an IP-based MCBH run, a certain number of resulting low energy structures can be re-optimized with more accurate methods (e.g. DFT). For some systems, quite accurate IPs have been developed for which parameters have been fitted to experimental and/or ab initio data. In principle, the better the IP, the lower the number of structures that require significant further optimization with a more accurate method. However, we note that to have an IP with higher accuracy, you typically require correspondingly more complex functional forms with a higher number of terms. In turn, during a MCBH run, the direct use of such IPs can lead to practical problems in relaxing structures far from the equilibrium geometry (after a random distortion). Here we present a MCBH procedure in which for each step, a cascade structure optimization is performed with two different IPs. The general idea is to first use a robust and efficient IP to pre-optimize a distorted structure and then use more accurate IP for fully relaxing this pre-optimized structure. The Metropolis criterion is than applied after the second optimization. Our cascade-MCBH approach interfaces with the GULP package [START_REF] Gale | GULP: A computer program for the symmetry-adapted simulation of solids[END_REF] which offers the possibility of using many IP types. We have implemented this procedure using Atomic Simulation Environment [START_REF] Bahn | An object-oriented scripting interface to a legacy electronic structure code[END_REF] (ASE) modules together.

In the literature, IPs for bare silica (SiO2) bulk systems, such as TTAM [START_REF] Tsuneyuki | First-Principles Interatomic Potential of Silica Applied to Molecular Dynamics[END_REF] and BKS, [START_REF] Van Beest | Force fields for silicas and aluminophosphates based on ab initio calculations[END_REF] have been successfully used. Hassanali and Singer (HS) [27] introduced a IP parameterization for the silicawater interface by adding new three-body terms terms to the BKS IP in order to describe Si-O-H and H-bonds. All these IPs are based on the Buckingham mathematical form as shown in equation [START_REF] Aggarwal | Use of nano-silica in cement based materials -A review[END_REF] (1) where i,j ϵ {Si,O}, rij is the interionic separation, qi,j are the ionic charges and A, B and C are fitted parameters. With respect to nanosilica, Flikkema and Bromley (FB) re-parameterised a Buckingham IP specifically for treating bare silica nanoclusters. [START_REF] Flikkema | A new interatomic potential for nanoscale silica[END_REF] Recently, Pedone et al. [START_REF] Pedone | FFSiOH: a New Force Field for Silica Polymorphs and Their Hydroxylated Surfaces Based on Periodic B3LYP Calculations[END_REF] have parameterized a new and more complex Buckingham IP with polarisable ions for hydroxylated silica. This IP also incorporates intra-and inter-bond Morse potentials and a three-body term for hydrogen bonds was shown to provide good structures and vibrational frequencies with respect to experimental and DFT-calculated data on hydroxylated silica surfaces. Recently, a systematic study has showed that this IP is also accurate for the treatment of hydroxylated nanoparticles. [30] In previous global optimization studies the FB IP and a simplified version of the HS IP (HSsimp) has been used for finding low energy clusters where the HS three-body terms were omitted for simplicity and computational efficiency. [START_REF] Bromley | Novel structures and energy spectra of hydroxylated (SiO2)8-based clusters: searching for the magic (SiO2)8O2H3 cluster[END_REF][START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF] By omitting these terms, the correct description of hydroxyl groups (OH) is lost, therefore, the clusters do not have any H-bond contribution to their stability (see Fig. 1). In this way, isomers stabilized by hydrogen bonds contributions are not easily found only using HSsimp in a MCBH run. Such concerns are especially relevant to relatively large clusters and/or a high hydroxylation degree where H-bonding can influence the energetic ordering and structure of low lying minima (each H-bond contributes of about 20 kJ/mol). Thus for silica nanoparticles with increasing size and/or with moderate to high hydroxylation a good description of OH groups and their H-bonding becomes essential in order to avoid missing possible low energy clusters.

In this work we consider hydroxylated silica clusters (SiO2)M•(H2O)N systems with M = 6, 8 10 and 12 each one with different hydroxylation ratio N:M (RN/M) up to ~65%. As a pre-optimization in the first cascade step we use the FB IP with OH groups described by HSsimp IP (IP_1), while for second part of the cascade we use the IP due to Pedone et al. (IP_2) [START_REF] Pedone | FFSiOH: a New Force Field for Silica Polymorphs and Their Hydroxylated Surfaces Based on Periodic B3LYP Calculations[END_REF]. In general, the mathematical form of a Buckingham IP energy entails that the energy of a pair of atoms tends to -∞ when their distance tends to 0. To avoid such a superposition of atoms caused by random displacements during a MCBH run we added very short range repulsive terms to the Buckingham describing the interaction between unlike charged species in IP_1. In order to more efficiently explore the PES we use a selfcorrecting temperature technique to keep the accept ratio between 60 and 80%. For each MCBH step atomic displacements are set to a maximum of 0.8 for any cartesian coordinate and the energy difference used to discriminate between two different structures is set to 1 kJ/mol. Several runs (at least three) of fifty thousand steps each are performed for all considered systems. The ten most energetically stable isomers resulting from the MCBH global optimization procedure for each system stoichiometry are subsequently refined using DFT. These final optimisations are done employing the FHI-AIMS package [31], using the B3LYP [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF] hybrid functional and a tier-1 atomic centred numerical basis set with a light integration grid with. The B3LYP functional has been shown to well describe hydroxylated silica in numerous previous studies. [START_REF] Bromley | Novel structures and energy spectra of hydroxylated (SiO2)8-based clusters: searching for the magic (SiO2)8O2H3 cluster[END_REF][START_REF] Flikkema | Structure and energetics of hydroxylated silica clusters, (SiO2)M(H2O)N, M = 8, 16 and N= 1 -4: A global optimisation study[END_REF][START_REF] Delle Piane | Large-Scale B3LYP Simulations of Ibuprofen Adsorbed in MCM-41 Mesoporous Silica as Drug Delivery System[END_REF][START_REF] Rimola | Silica Surface Features and Their Role in the Adsorption of Biomolecules: Computational Modeling and Experiments[END_REF]

Results and Discussions

We first tested our cascade MCBH methodology on (SiO2)8•(H2O)N systems with N = 2-5 incorporated water molecules as previously studied in refs [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF] and [START_REF] Flikkema | Structure and energetics of hydroxylated silica clusters, (SiO2)M(H2O)N, M = 8, 16 and N= 1 -4: A global optimisation study[END_REF]. Specifically, we consider the three most stable isomers we find for each series with those previously reported. In Fig. 2 we report isomers structures found both from this work and from refs [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF] and [START_REF] Flikkema | Structure and energetics of hydroxylated silica clusters, (SiO2)M(H2O)N, M = 8, 16 and N= 1 -4: A global optimisation study[END_REF].

Fig. 2 Structures of (SiO2)8•(H2O)N nanoclusters obtained in this work ( a ), ref [START_REF] Flikkema | Structure and energetics of hydroxylated silica clusters, (SiO2)M(H2O)N, M = 8, 16 and N= 1 -4: A global optimisation study[END_REF] as ( b ) and in ref [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF] as ( c ). Clusters for each N are ordered from left to right according to decreasing energetic stability. In parentheses energy differeneces with respect to the best global minimum candidate (far left of each row) are reported in kj/mol. Atom colour key: Si -yellow, Si, O -red, H -grey.

In our comparison we take all stable conformational isomers related to different H-bonding patterns of the hydroxyls on a single Si-O skeleton with a fixed connectivity to be representations of the same isomer. Typically, the total energy range due to different H-bonding patterns on the same isomer of this size is 1 -10 kJ/mol. In all cases we report isomers with the lowest energy H-bonded pattern found. We note that although the first part of the cascade does not provide a correct description of H-bonding, in the second part of the cascade several H-bond patterns are sampled which ensures that the conformational space hydroxyl orientations is well explored. This part of the cascade is particularly important for relatively large and highly hydroxylated clusters where H-bonding can contribute a relatively large amount to total energy differences.

In the case of the most energetically stable structures for each N considered (i.e. Fig. 2, isomers 1, 4, 8 and 13) we found the same isomers as in previous works. This result further confirms these candidate structures as global minima. The second most stable isomers for each N (i.e. Fig. 2, isomers 2, 5 and 10) also confirm the corresponding results of refs [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF] and [START_REF] Flikkema | Structure and energetics of hydroxylated silica clusters, (SiO2)M(H2O)N, M = 8, 16 and N= 1 -4: A global optimisation study[END_REF]. For the third most stable isomers, for two incorporated water molecules (N = 2) we also found the previously reported isomer structures. Using the cascade method, however, we were able to find more energetically stable candidates for N = 3 (Fig. 2, isomer 6 being 3 kJ/mol lower than cluster 7) and N = 4 (isomers 14 and 22 kJ/mol lower than cluster 12). For N = 5 we confirm the global minimum candidate isomer structure and other low energy isomers are reported in literature.

After clearly demonstrating the capacity of our cascade MCBH approach for hydroxylated silica clusters for M = 8, we applied our method for to the, as yet unreported, systems: (SiO2)M•(H2O)N with M = 6, 10, 12, each for a range of N. The obtained new global minima candidate clusters are shown in figure 3. Low numbers of incorporated water (i.e. N = 1, 2) do not fully hydroxylate the one-coordinated defective oxygen centres present in anhydrous silica clusters. There are two types of singly-coordinated dangling oxygen defects: (i) oxygen atoms in silanone species (formally: >Si=O, but arguably more accurately: >Si + -O - [START_REF] Avakyan | AIM and ELF electronic structure/G2 and G3 π-bond energy relationship for doubly bonded silicon species, H2Si=X (X=E14H2, E15H, E16)[END_REF][START_REF] Zwijnenburg | The effect of local environment on photoluminescence: A time-dependent density functional theory study of silanone groups on the surface of silica nanostructures[END_REF]) having three-coordinated Sicentres, and (ii) non-bridging oxygen (NBO) atoms accompanied with compensating triple coordinated oxygen-sites, also known a compensated NBOs or valence alternation pairs [START_REF] Hamad | Low reactivity of non-bridging oxygen defects on stoichiometric silica surfaces[END_REF][START_REF] Zwijnenburg | Long range coupling between defect centres in inorganic nanostructures: Valence alternation pairs in nanoscale silica[END_REF]. These defective centres are well studied in literature and it is known that they can readily be hydroxylated through reaction with water [START_REF] Zhuravlev | The surface chemistry of amorphous silica. Zhuravlev model[END_REF][START_REF] Yang | Reaction of water on silica surfaces[END_REF]. By increasing the number of incorporated water on an already fully hydroxylated cluster, bridging oxygens (Si-O-Si) are broken and the formation of silanols (Si-OH) occur until there is one -OH group for each Si-centre. Further hydroxylation tends to lead to the formation of geminal silanols (i.e. two OH groups per Si). The hydroxylation reaction (∆ ) energy:

(3)

indicates the relative stability of the hydroxylated structures with N incorporated water molecules with respect to the equivalent anhydrous silica system. This quantity indicates the tendency of anhydrous silica to be hydroxylated, in particular if ∆ is negative the reaction is energetically favourable. Although strictly speaking a thermodynamically favourable reaction is evaluated by the free energy difference, we neglect the entropic contributions and assume that ∆ a good indicator of reaction favourability. In figure 4 we show hydroxylation energies (Fig. 4a) and Si-centred SiO4 tetrahedral distortions (Fig. 4b) with respect to incorporated water ratio (RN/M = N/M) for the most stable (SiO2)M•(H2O)N candidates with M = 6, 8, 10, 12. The tetrahedral distortion is taken as the average of the differences between the root mean square of the six O-Si-O angles relative to optimal unstrained value of 109.47 for each SiO4 centre. Generally, structures made mainly with rings containing more than three -(Si-O)-units are found to have low tetrahedral distortion. [START_REF] Bromley | Importance of the embedding environment on the strain within small rings in siliceous materials[END_REF] Interestingly, several of the structures shown in Fig. 3, can be viewed as being generated from condensation of building blocks of smaller sized isomers. The most common building blocks observed, for instance, are the (SiO2)6•(H2O)3 trigonal prism (Fig. 2: isomer 3), the (SiO2)8•(H2O)4 cubic structure (Fig. 1: isomer 8 -also known as a double 4-ring or D4R) and the pentagonal (SiO2)8•(H2O)4 structure (Fig. 1: isomer 9). As the cluster size increase the most common building block becomes the D4R (e.g. see isomers 9 and 11 in Fig. 3) which only exhibits rings with low tetrahedral distortion. Larger isomers in other work also commonly exhibit the cubic (SiO2)8•(H2O)4 D4R unit (e.g. (SiO2)24•(H2O)12 in ref [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF] which can be considered as a condensation of D4R units in triangular arrangement).

∆ is found to be negative and monotonically decreasing with respect to the number of incorporated water molecules showing that the hydroxylation reaction of anhydrous silica is always energetically favoured. We also notice an energetic down shift with respect to M for the different N series, due to the increasing of the size of the system. The decreasing of ∆ in a series of N values for a fixed size M, can be split into the contribution by three factors: i) at a low number of incorporated water molecules, the gain in energy is mainly due to hydroxylation of the defective oxygen centres, ii) once a cluster is fully hydroxylated (defective centres are no longer present), ∆ drops because adding further water molecules helps to decrease the tetrahedral distortion of Si-centres (e.g. by hydrolysing small strained rings), iii) at a high number of incorporated water molecules there is a competition between two effects affecting ∆ : the formation of geminal silanols (energetically not favoured effect since O-Si is stronger a bond than O-H), and the formation of hydrogen bonds which is an energetically stabilizing effect. These latter two contributions in phase (iii) seem to be almost equivalent energetically, leading to the drop in ∆ tending to level off. The hydroxylation degree at the transition point between the ii) and the iii) phases can be regarded as the optimal hydroxylation degree ( / ). In other words, ( / ) is the RN/M value where the hydroxylation reaction energy first levels off to a flatter linear regime due to the minimum tetrahedral distortion being reached with increasing hydroxylation.

/ was first introduced in ref [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF] for systems with M = 4, 8, 16, 24. Specifically, / was tracked with respect to the number of SiO2 units (M) in order to examine its size dependency.

Herein we further add

/ values for the cluster sizes M = 6, 10, 12 silica sizes. To extract / values for each (SiO2)M series we choose the RN/M value that correspond to the cluster with the lowest tetrahedral distortion which is consistent with choice made in ref [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF]. The selected / values for clusters with M = 6, 10, 12 correspond respectively to structures 3, 8, 11 in figure 3.

Although the choice of / = 0.5 for the set of (SiO2)10•(H2O)N clusters is consistent with the RN/M value that lowers the tetrahedral distortion most, we notice that the RN/M value where tetrahedral distortion levels off (i.e. 0.4) would be likewise acceptable. Fig. 5 Data for optimal hydroxylation degree (R / ) with respect to cluster size fitted with an inverse power law. Black filled square symbols correspond to data from ref [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF].

/ with respect to the number of SiO2 units (M) is shown in Fig. 5 with data from ref [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF] in black filled squares and data from this work in black filled circles. The half-filled circle correspond to the (SiO2)8•(H2O)N set both studied previously and used in this work as test case. The first point of the series (M = 1) corresponds orthosilicic acid (Si(OH)4) where / = 2. / is higher for smaller clusters and drop as the size increases and tends to follow an inverse power law. The new data generated in this work (M = 6, 10, 12) do not bring to significant changes in the power law equation in ref [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF] and confirm the / size dependency trend. Interestingly, the exponent of the inverse power law is close to -2/3 while the general size dependency of many properties of a generic nanoparticles scale with a inverse power law but with an exponent of -1/3. [START_REF] Johnston | Atomic and Molecular Clusters[END_REF] Generally, the -1/3 exponent derives from the surface-to-volume ration assuming a spherical nanoparticle shape. For a generic property size dependency, we expect a small deviation in the exponent from -1/3 due nonspherical shape of real structures and possibility due to property fluctuations for their extreme dependency on small changes when systems are relatively small. Nevertheless, the -2/3 exponent in / size dependency seems to be affected by some unknown factors for such deviation from -1/3 that we hope to elucidate in future works.

Conclusions

In this work we present a methodology for performing IP-DFT basing hopping global optimization studies on systems such as hydroxylated silica clusters where using an accurate IPs is essential to correctly describe hydrogen bonds. For the (SiO2)8•(H2O)N test case, the cascade approach finds all previously reported isomers (including global minima candidates) and additionally finds three new low energy isomers. New global minima candidate structures are generated for (SiO2)M•(H2O)N systems with M = 6, 10 and 12, each for a range of N values. The hydroxylation reaction energy trends for these series are consistent with previous works on other cluster sizes. The hydroxylation degree shows a threshold in each case when the hydroxylation-induced structural transformations minimises the tetrahedral distortion around the SiO4 centres. At this optimal hydroxylation ( / ) point the reaction energy starts to level off to a constant value. Confirming previously reported results for other cluster sizes [START_REF] Jelfs | Hydroxylation of silica nanoclusters (SiO2)M(H2O)N, M = 4, 8, 16, 24[END_REF], we confirm that the size (M) dependence of / follows an inverse power law with an exponent close to -2/3. Generally, the cascade approach we develop enables for the efficient global optimisation of systems which require treatment via complex forcefields.

Fig. 1

 1 Fig. 1 Structures of a (SiO2)8•(H2O)4 nanocluster optimized with: IP_1 used in the first cascade step (left), IP_2 used in the second cascade step (middle), and DFT using B3LYP (right). Atom colour key: Si -yellow, Si, O -red, H -grey.

Fig. 3

 3 Fig. 3 The most energetically stable (SiO2)M•(H2O)N nanoclusters obtained in this work for M = 6, 10 and 12 silica units, each with a range of incorporated water molecule (N) corresponding to a degree of hydroxylation (RN/M = N/M) between 17 -67%. Atom colour key: Si -yellow, Si, O -red, H -grey.

Fig. 4 a

 4 Fig. 4 a) Hydroxylation delta energy (∆Ehyx) in kJ/mol relative to the correspondingly sized global minima (SiO2)M•(H2O)N clusters for M = 6, 8,10,12 with respect to incorporated water ratio (RN/M). b) Deviations from the ideal silicon-centred tetrahedron (SiO4) in degrees with respect to hydroxylation degree. Symbols for M = 8 shown as half filled circles are also reported in ref [17] and reproduced by us in this work.

Acknowledgements

Support from Spanish MINECO/FEDER grant CTQ2015-64618-R grant and, in part, by Generalitat de Catalunya grants 2014SGR97, XRQTC is acknowledged. We also acknowledge the NOMAD Center of Excellence project (this project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 676580). Access to supercomputer resources was provided through grants from the Red Española de Supercomputación. We acknowledge funding from H2020 project ITN-EJD-642294 (TCCM: Theoretical Chemistry and Computational Modelling). We also thank Ask Hjorth Larsen for technical advice concerning the Atomic Simulation Environment.


Sudhanya Chatterjee

Olivier Commowick

Onur Afacan

Benoît Combès

Anne Kerbrat

Simon K Warfield

Benoit Combès

Christian Barillot

A 3-year

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A 3-year follow-up study of enhancing and non-enhancing multiple sclerosis (MS) lesions in MS patients demonstrating clinically isolated syndrome (CIS) using a multi-compartment T2 relaxometry (MCT2) model

Introduction:

Demyelination occurs at the onset of MS accompanied by macrophage intervention [START_REF] Barillot | Imaging biomarkers in Multiple Sclerosis: from image analysis to population imaging[END_REF][START_REF] Lassmann | Heterogeneity of multiple sclerosis pathogenesis: implications for diagnosis and therapy[END_REF] . This is followed by axonal damage. Continuous tissue injuries cause inflammation in lesion affected regions of the brain. Advanced MRI techniques such as T2 multi-echo MRI sequences and diffusion imaging can help obtain brain tissue microstructure information.

Myelin, intra/extra-cellular fluids, and free water can be distinguished on the basis of their T2 relaxation times [START_REF] Mackay | Insights into brain microstructure from the T 2 distribution[END_REF] . Myelin has a very short-T2 3 and its presence in WM is quantified by a variety of approaches [START_REF] Mackay | Insights into brain microstructure from the T 2 distribution[END_REF][START_REF] Whittall | In vivo measurement of T2 distributions and water contents in normal human brain[END_REF][START_REF] Laule | Myelin water imaging in multiple sclerosis: quantitative correlations with histopathology[END_REF][START_REF] Mackay | Magnetic Resonance of Myelin Water: An in vivo Marker for Myelin[END_REF][START_REF] Akhondi-Asl | Fast myelin water fraction estimation using 2D multislice CPMG[END_REF] . The most popular approach is to obtain the water fraction (WF) corresponding to myelin (myelin water fraction (MWF)) in a tissue 3- 7 from a multi-echo spin echo (SE) T2 MRI. However, since MWF is a relative measure, a decrease or increase in it only conveys a part of the information. In this work we use a multi-compartment T2 relaxometry (MCT2) model [START_REF] Chatterjee | Gaining Insights Into Multiple Sclerosis Lesion Characteristics from Brain Tissue Microstructure Information: A Multi-Compartment T2 Relaxometry Approach[END_REF] to obtain WF corresponding to tissues with short-, medium-and high-T2 relaxation times. We perform test-retest experiments on 4 healthy controls (HC) to demonstrate the reproducibility of the WF estimates obtained from this model. A 3-year follow-up study is performed on 6 MS patients demonstrating clinically isolated syndrome (CIS).

We study evolution of WFs in enhancing and non-enhancing MS lesions in white matter (WM).

Method and materials:

MCT2 model:

In the MCT2 model 8 , the T2 space is modeled as a weighted mixture of three continuous Gaussian probability density functions representing the three T2 relaxometry compartments with short, medium and high decaying components with respect to their T2 relaxation times. The short-T2 WF (WFS) provides information on myelin [START_REF] Mackay | Insights into brain microstructure from the T 2 distribution[END_REF] and highly myelinated axons [START_REF] Lancaster | Three-pool model of white matter[END_REF] in WM. The medium-T2 WF (WFM) correspond to the intra/extra-cellular matters [START_REF] Mackay | Insights into brain microstructure from the T 2 distribution[END_REF][START_REF] Lancaster | Three-pool model of white matter[END_REF] . The high-T2 WF (WFH) correspond to cerebrospinal fluids and inflammation in WM due to lesions.

Test-retest scans were performed for 4 HC. Average values in 8 WM ROIs were computed. A Bland-Altman plot was observed for assessing the reproducibility of the WF estimates. Similar to a previous study [START_REF] Vargas | Measuring longitudinal myelin water fraction in new multiple sclerosis lesions[END_REF] we used this approach to obtain the reproducibility threshold value of WF estimates.

Data: 3T MRI scanner, first echo time (TE)=9ms, echo spacing (∆TE)=9ms repetition time (TR)=2030ms, number of echoes (nechoes)=32, voxel dimensions (vd)=1.64x1.64x4mm [START_REF] Mackay | Insights into brain microstructure from the T 2 distribution[END_REF] . All images were registered [START_REF] Commowick | Block-matching strategies for rigid registration of multimodal medical images[END_REF][START_REF] Ourselin | Block matching: A general framework to improve robustness of rigid registration of medical images[END_REF] to a common T2weighted image.

MS lesion study:

Six MS patients demonstrating CIS condition were scanned at the baseline, and follow-up scans at month-{3,6,9,12,18,24,36}. The median age of patients included in this study was 28.5 years (at baseline) and there were equal number of male and female subjects. T2 relaxometry data: 3T MRI scanner, first TE=13.8ms, ∆TE=13.8ms, TR=4530ms, nechoes=7, vd=1.3x1.3x3mm 3 , acquisition time=7minutes. For months-{0,3,6,9}: transverse SE T1weighted images (1x1x3mm 3 ) were acquired before and after USPIO infusion (SHU-555C;40µmol Fe/kg body weight over 30minutes) for obtaining USPIO enhanced lesions. A scan was done 24 hours later to obtain post-USPIO enhancement images. Transverse SE T1-weighted images (1x1x3mm 3 ) post gadolinium contrast agent infusion (0.1mmol/kg gadopentetate dimeglumine) were acquired to find gadolinium enhanced lesions.

Following lesion types [START_REF] Kerbrat | USPIO-positive MS lesions are associated with greater tissue damage than gadolinium-positive-only lesions during 3-year followup[END_REF] were studied:

-(U+): Appear on post-contrast USPIO and Gadolinium scans.

-(Gd+): Appears on post-contrast Gadolinium scan, but not in (U+).

-(L-): MS lesion ROI not in (Gd+) and (U+).

-(E+): These are enhancing lesions and appear on (Gd+) or (U+)

All images are registered to 3D-FLAIR image (1x1x1mm 3 ) acquired at the baseline. We analyzed 111(L-), 23(Gd+) and 6(U+) lesions.

We studied percentage of lesions in each type that underwent changes above the reproducibility threshold between an acquisition and its follow-up, and the evolution of MCT2 WF estimated over a period of 36 months.

The protocols were approved by the institutional review board of Rennes University Hospital, and all participants gave their written consent. The evolution of WFs in each lesion type is shown in Figure -4.

Results and discussion:

Conclusion:

(U+) lesions showed higher variation in initial months as compared to (Gd+) and (L-). The MCT2 8 estimates are reproducible and provided insights into the lesion that corroborate with pathological studies 2 of MS lesions. In future we intend to perform more rigorous statistical analysis to evaluate if the MCT2 estimates can significantly differentiate among lesion groups studied in this abstract.

Figure- 1

 1 Figure-1 shows WF maps of a healthy control. The reproducibility results (refer Figure-2) show that the 95% limits of agreement (LoA) for both WFS and WFM are 0.013. The percentage of lesions (PoL) undergoing change in WFS above LoA (refer Figure-3(a)) are relatively smaller for (L-) than (E+) lesions. The percentage of lesions undergoing change in WFS (refer Figure-3(a) and 3(c)) above LoA reduce after 6 months. In general, (L-) have lesser PoL undergoing change above AoL as compared to (Gd+) and (U+).

Figure 1 :

 1 Figure 1: Water fractions estimated for a healthy volunteer. Water fraction maps estimated appear normal for a healthy subject, such as the genu has higher short-T2 water fraction values than the frontal white matter.

Figure 2 :

 2 Figure 2: Bland-Altman plots for test-retest experiments for (a) short-T2, (b) medium-T2 and (c) high-T2 water fraction estimations for 8 white matter ROIs in 4 healthy controls. Results shows a high agreement between the test-retest experiments with limits of agreement (LoA) being 0.013 for short-T2 and medium-T2 WF test-retest estimations. The high-T2 water fraction estimation are in perfect agreement.

Figure 3 :

 3 Figure 3: Percentage of Lesions (PoL) undergoing change in water fraction (WF) more than the limit of agreement (LoA) found from the Bland-Altman plots in Figure-2 are shown here. The x-axis shows the duration (in months) between which the changes were computed. Comparison in PoL change above LoA for (a) Short-T2 WF and (b) medium-T2 WF for (L-) and (E+) lesions. Comparison in PoL change above LoA for (c) Short-T2 WF and (d) medium-T2 WF for (L-), (U+) and (Gd+) lesions. Higher PoL in (E+) having WFM above LoA might be an indication of increased extra-cellular matter (as macrophage) interventions in early lesions.

Figure 4 :

 4 Figure 4: Evolution of (a) short-T2 (WFS), (b) medium-T2 (WFM) and (c) high-T2 (WFH) water fraction in (L-), (Gd+) and (U+) lesion over a period of 36 months. The WFS values of (U+) lesions increase after month-6 compared to (Gd+) and (L-) lesions, indicating a delayed remyelination. This might be explained as (U+) lesions are in an earlier stage compared to (Gd+) and (L-) and hence have a delayed recovery. The high WFM values in (U+) might be attributed to the known increased extra-cellular matter presence in early lesions. By the end of month-36, (L-), (U+) and (Gd+) attain similar WFS and WFH values.

Acknowledgements:

MRI data acquisition was supported by the Neurinfo MRI research facility from the University of Rennes I. Neurinfo is granted by the the European Union (FEDER), the French State, the Brittany Council, Rennes Metropole, Inria, Inserm and the University Hospital of Rennes. This research was also supported in part by the grants NIH R01 EB019483 and NIH R01 NS079788.


Member Claude Jauffret
email: jauffret@univ-tln.fr

Observability and Fisher Information Matrix in Nonlinear Regression

Keywords: Observability, identifiability, Cramèr-Rao Lower Bound, Fisher Information Matrix, nonlinear regression

This paper is devoted to the link between the Fisher Information Matrix invertibility and the observability of a parameter to be estimated in a nonlinear regression problem.

I. Introduction

In Signal Processing [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF], in Target Motion Analysis [START_REF] Bar-Shalom | Estimation and Tracking: Principles, Techniques and Software[END_REF], a large class of measurements can be modelled as

X = h(θ) + ε ( 1
)
where X is the available measurement vector (element of R n ), θ is the unknown deterministic parameter (lying in R d ) and h(.) is a (known) nonlinear mapping from R d to R n .

The vector ε represents the additive measurement noise.

Whatever the estimation technique employed (Least Squares, Maximum Likelihood, ...), the observability of parameter θ must be investigated.

Most of the time, the analysis of observability is a tough task, and many authors suggest to declare that θ is observable if the Fisher Information Matrix (FIM) of θ given X, under Gaussian hypothesis concerning ε, is nonsingular at θ (see, for example, [START_REF] Bar-Shalom | Estimation and Tracking: Principles, Techniques and Software[END_REF] page 168).

Intuitively (and practically), it turns out that this way is sufficient. But it is legitimate to wonder about the meaning of that analysis: Observability being a deterministic notion, why use a statistical tool (the FIM) to establish its status? Why restrain oneself to the Gaussian law? Would the conclusion be the same if the FIM was computed for another law?

Parts of answers can be found in Jazwinski's book ( [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF] page 231) but linear Gaussian cases are concerned.

The aim of this paper is to answer these questions by establishing clearly the link between the status of the FIM (singular or nonsingular) under a large class of probability laws and the status of the observability, i.e. in nonlinear and non-Gaussian cases.

The coming section recalls some classical definitions of observability. In the third section, we give a general form of the FIM. Some pathologic (but still relevant) cases are analyzed in Section IV. The last section gives the main result, after presenting the necessary mathematical tools.

II. Observability Concepts

There are several ways to define the observability concept: it can be a global one (for all the vectors of R d ), or a local one (for a special θ).

We recall the three major definitions of observability.

Definition 1. The noise-free system

X = h(θ) is simply observable at θ 0 if ∀θ ∈ R d , {θ = θ 0 } ⇒ {h(θ) = h(θ 0 )} (2)
Definition 2. The noise-free system is (simply

) observable if ∀θ, ∀θ ∈ R d , {θ = θ } ⇒ {h(θ) = h(θ )} (3)
Definition 3. The noise-free system is locally observable at θ 0 if

∃U θ 0 ⊂ R d (open subset containing θ 0 ), ∀θ ∈ U θ 0 , {θ = θ 0 } ⇒ {h(θ) = h(θ 0 )} (4) Remark 1.
a) These definitions come from the theory of dynamic systems in which the parameter θ changes in time and must be hence denoted θ(t). For such systems, there exists some other definitions of observability to take into account the trajectory of θ(t) [START_REF] Hermann | Nonlinear controllability and observability[END_REF].

b) When h(.) is a linear mapping (in practice a matrix), local observability and simple observability coincide.

III. Fisher Information Matrix

In the sequel, we assume that the behavior of the vector ε is described by a probability density function (pdf), say p ε , whose support is R n . The vector X has its own pdf, denoted p X . It depends on θ while its support is independent of it 1 . More precisely, the 1 This assumption is necessary to compute the Fisher Information Matrix.

October 14, 2006 DRAFT relationship between them is given by

p X (ν|θ) = p ε (ν -h(θ)). ( 5
)
The likelihood function of θ is nothing else than the probability density function of X

given θ evaluated at X :

L θ (X) p X (X|θ). (6)
The Fisher Information Matrix (FIM) is

F θ (X) Cov θ {∇ θ ln [L θ (X)]} = E θ ∇ θ ln[L θ (X)]∇ T θ ln[L θ (X)] . ( 7
)
The Cramèr-Rao Lower Bound (CRLB) of any unbiased estimator of θ is the inverse of

F θ (X).
If ε is a 0-mean Gaussian vector whose covariance matrix is R ε (assumed invertible), the FIM can be expressed as follows

F θ (X) = ∇ θ h(θ)R -1 ε ∇ T θ h(θ). ( 8
)
Under more general assumptions, we can still give a close form for the FIM .

Theorem 1. F θ (X) = ∇ θ h(θ)W ε ∇ T θ h(θ)
where

W ε E{∇ ν ln [p ε (ν)] ν=ε ∇ T ν ln [p ε (ν)] ν=ε } (9)
Proof :

For the dummy variable z, we have

∇ θ ln [p X (z|θ)] = ∇ θ ln [p ε (z -h(θ))] = ∇ θ h(θ)∇ ν ln [p ε (ν)] ν=z-h(θ) . ( 10
)
Hence

∇ θ ln [L θ (X)] = ∇ θ h(θ)∇ ν ln [p ε (ν)] ν=X-h(θ) . ( 11
)
As a consequence, the FIM is readily written as

F θ (X) = ∇ θ h(θ)E θ {∇ ν ln [p ε (ν)] ν=X-h(θ) ∇ T ν ln [p ε (ν)] ν=X-h(θ) }∇ T θ h(θ) (12) DRAFT October 14, 2006
The middle term above is

E θ {∇ ν ln [p ε (ν)] ν=X-h(θ) ∇ T ν ln [p ε (ν)] ν=X-h(θ) } = E{∇ ν ln [p ε (ν)] ν=ε ∇ T ν ln [p ε (ν)] ν=ε }. ( 13
)
In the sequel, W ε is assumed nonsingular. The following theorem gives us a sufficient condition for that.

Theorem 2. If ∇ ν p ε (ν) is a continuous function then W ε is nonsingular. Proof :
Suppose that W ε is singular. So there exists a non null vector

u ∈ R d such that u T W ε u = 0, i.e. ∃u = 0 s.t. u T W ε u = 0 ⇔ ∃u = 0 s.t. u T R d ∇ ν p ε (ν)∇ T ν p ε (ν) 1 p ε (ν) dν u = 0 ⇔ ∃u = 0 s.t. R d u T ∇ ν p ε (ν) p ε (ν) 2 dν = 0 (14)
Since ∇ θ p ε (θ) is a continuous function, this last statement is equivalent to

∃u = 0 s.t. u T ∇ ν p ε (ν) = 0, ∀ν ∈ R d ⇔ ν∈R d ker{∇ ν p ε (ν)∇ T ν p ε (ν)} = { 0}. ( 15
)
Let r be the dimension of the vector space

ν∈R d ker{∇ ν p ε (ν)∇ T ν p ε (ν)}.
In a suitable basis, the last (d -r) components of ∇ ν p ε (ν) will be null, i.e.

∇ ν p ε (ν) = ∂p ε (ν) ∂ν 1 , ∂p ε (ν) ∂ν 2 , • • • , ∂p ε (ν) ∂ν r , 0, • • • , 0 T ∀ν ∈ R d ( 16
)
which means that in that basis, p ε (ν) = p ε (ν 1 , ν 2 , ..., ν r ). This contradicts the fact that

R d p ε (ν)dν = 1. Remark 2. If ε is Gaussian of covariance matrix R, it is readily shown that W ε = R -1 ε .
October 14, 2006 DRAFT

IV. Pathologic cases

A. Case I :

We consider the two-dimensional measurement

(n = d = 2) X   X 1 X 2   =   θ 3 1 + θ 2 θ 3 1 -θ 2   +   ε 1 ε 2   . ( 17
)
Obviously, the associated noise-free system is simply observable at any θ = (θ 1 , θ 2 ) T .

Under the assumption that ε i ∝ G(0, 1) and independent, the FIM is

F θ (X) =   9θ 4 1 0 0 2   (18)
which is singular at any θ such that θ 1 = 0.

Remark 3.

a)

In similar cases, we can prove that no unbiased estimator of θ exists since for such estimators the CRLB is the inverse of the FIM. This fact is met in array processing for the estimation of the end-fire bearing [START_REF] Pasupathy | Range and bearing estimation in passive sonar[END_REF].

b) The singularity of the FIM at some points of R d can cause some problems during the Gauss Newton routine for which the Hessian is approximated by the FIM evaluated at the point of the current iteration. The palliative is the augmentation of the FIM by some αId as suggested in the Levenberg-Marquardt method [START_REF] Dennis | A user's guide to nonlinear optimization algorithms[END_REF].

B. Case II

This counter-example comes from [START_REF] Isidori | Nonlinear Controllability Systems[END_REF] p. 479. Let's consider the two-dimensional measurement vector

X   X 1 X 2   =   aθ -sin θ cosθ   +   ε 1 ε 2   , with a ∈]0, 1[ (19) F θ (X) = 2a a + 1 2a -cos θ (20)
Hence, the FIM is never equal to 0, but still the pairs (θ 1 , θ 2 ) defined by

   θ 1 2kπ + τ θ 2 2kπ -τ (21) DRAFT October 14, 2006
are undistinguishable, τ being the root of the equation τ = sin τ a . The parameter θ is not simply observable, but locally observable.

V. Analysis

A. Mathematical Tools

We need two types of tools : one from the linear algebra theory and the second one from the differential calculus.

Theorem 3. Let A a (n × d) matrix (d ≤ n).
The following statements are equivalent : [START_REF] Isidori | Nonlinear Controllability Systems[END_REF] p. 479).

(i) A T A is invertible. (ii) ∃S a real nonsingular symmetric (n × n) matrix such that A T SA is nonsingular. (iii) ∀S real nonsingular symmetric (n × n) matrix, A T SA is nonsingular. (iv) Rank(A) = d. Definition 4. h : R d → R n is an immersion at θ 0 if the rank of ∇ θ h(θ 0 ) is equal to d (see
Theorem 4. h : R d → R n is an immersion at θ 0 if there exists an open set U θ 0 containing θ 0 such that the rank of ∇ θ h(θ) is equal to d, whatever θ in U θ 0 .

See [START_REF] Schwarz | Calcul Différentiel et Equations Différentielles[END_REF] for the proof.

Proposition 1. If h : R d → R n is an immersion, than h is locally injective, i.e. it exists an open set U θ 0 of R d containing θ 0 such that h : U θ 0 ⊂ R d → R n is injective.
See [START_REF] Schwarz | Calcul Différentiel et Equations Différentielles[END_REF] for the proof.

Remark 4. If h is a linear mapping i.e. h(θ)

= A(θ), then ∇ θ h(θ 1 ) = ∇ θ h(θ 2 ) = A T for
any pair (θ 1 , θ 2 ). Hence h is locally injective as well as injective anywhere.

B. The main result

The exploitation of the previous theorems and of the last proposition yields straightforwardly the following October 14, 2006 DRAFT Theorem 5. Let us consider the measurement equation X = h(θ) + ε associated to the noise-free system X = h(θ). If the support of p ε is R n itself and if the FIM is nonsingular at θ 0 -or, equivalently, if ∇ θ h(θ 0 )∇ T θ h(θ 0 ) is nonsingular -then the noise-free system is locally observable at θ 0 .

Proof :

We know that ∇ θ h(θ)W ε ∇ T θ h(θ) is nonsingular. Using Theorem 4, statement (iv), with A ≡ ∇ θ h(θ 0 ) and S ≡ W , we conclude that Rank(∇ θ h(θ 0 )) is equal to d ; hence h is an immersion. Now, thanks to Proposition 1, we know that h is locally injective, and as a consequence, satisfies Definition 3.

Remark 5.

a) The first pathological case forbids the converse.

b) Remark 4 proves that local observability is equivalent to simple observability for linear system

: if h is linear, i.e. h(θ 0 ) = Aθ 0 , then ∇ θ h(θ 0 )∇ T θ h(θ 0 ) = A T A, ∀θ 0 .
In that case, the system is observable iff A T A is nonsingular, or equivalently, iff A T W ε A (the FIM) is nonsingular.

VI. Conclusion

The link between the invertibility of the FIM and the observability status has been unambiguously established, for a large class of probability laws in nonlinear regression problems. This theoretical result can help the study of observability.

  Claude Jauffret, born in France on March 29 1957, received the Diplôme d'Etude Approfondies in applied Mathematics from the St Charles University (Marseille, France) in 1981, the Diplôme d'Ingénieur from Ecole Nationale Supérieure d'Informatique et de Mathématiques Appliquées de Grenoble (Grenoble, France) in 1983, the title of "Docteur de l'Université" in 1993 and the "Habilitation à Diriger des Recherches" from the Université de Toulon et du Var (France). From Nov. 1983 to nov. 1988, he worked on passive sonar systems, more precisely on Target Motion Analysis at the GERDSM (France). After a sabbatical year at the University of Connecticut (from nov. 1988 to dec. 1989) during which he worked on tracking problems in cluttered environment, he developped researches in tracking, data fusion and extraction in CERDSM (France). Since sept. 1996, he has been at the Université de Toulon et du Var where he teaches Statistical Signal Processing and Dynamic Linear Systems. His currents researches are about observability, estimation in non-linear systems as they appear in tracking problems DRAFT October 14, 2006


Annie-Claude Pérez
email: annie-claude.perez@univ-tln.fr2jauffret@univ-tln.fr

Claude Jauffret #

Denis Pillon
email: 3pillon_denis@orange.fr

Performance of Range-Only TMA

Range-only target motion analysis (ROTMA) is the topic of this paper: we focus our study on the numerical aspect and performance of the maximum likelihood estimates (MLE) for some scenarios when the noise polluting the measurements is additive and Gaussian. The performance is compared to the Cramér-Rao lower bound (CRLB).

I. INTRODUCTION

Range-only target motion analysis (ROTMA) has been treated in the literature for various applications, the target having a constant velocity (CV) motion: for example in [START_REF] Ristic | Target Motion Analysis using Range-Only Measurements: Algorithms, Performance and Application to ISAR Data[END_REF], [START_REF] Sathyan | Multiple Hypothesis Tracking with Multiframe Assignment Using Range and Range-rate Measurements[END_REF] and [START_REF] Clark | A Gaussian mixture filter for Range-Only Tracking[END_REF], maritime surveillance radar systems, ISAR (inverse synthetic aperture radar), are studied. It offers also an efficient way to calibrate a sonar array [START_REF] Pillon | Method for Antenna Angular Calibration by Relative Distance Measuring[END_REF]. ROTMA is encountered in the robotic domain by radio-frequency [START_REF] Cevher | A range-only multiple target particle filter tracker[END_REF] [START_REF] Huang | A Bank of Maximum A Posteriori Estimators for Single-Sensor Range-only Target Tracking[END_REF]. Local observability was analysed by the rank of the Fisher information matrix (FIM) [START_REF] Song | Observability of Target Tracking with Range-Only Measurements[END_REF]. Global observability conditions were proposed recently in [START_REF] Pillon | Observability: Range-Only vs. Bearings-Only Target Motion Analysis for a Leg by Leg Observer's Trajectory[END_REF] and [START_REF] Jauffret | Observability: Range-Only vs. Bearings-Only Target Motion Analysis when the Observer Maneuvers Smoothly[END_REF]. In the same vein, we present here some aspects of estimation (which was not addressed in [START_REF] Pillon | Observability: Range-Only vs. Bearings-Only Target Motion Analysis for a Leg by Leg Observer's Trajectory[END_REF] and [START_REF] Jauffret | Observability: Range-Only vs. Bearings-Only Target Motion Analysis when the Observer Maneuvers Smoothly[END_REF]), in particular for scenarios where observability is not acquired and for which the FIM is fullranked, or conversely for scenarios where the trajectory of the target is observable, but the FIM is singular. This motivates the structure of our paper.

In section II, we present the assumptions and notations. Section III is devoted to the special case where the observer is in CV motion too. In section IV, the observer's trajectory is composed of two legs. Here, the target's trajectory can be observable or not [START_REF] Pillon | Observability: Range-Only vs. Bearings-Only Target Motion Analysis for a Leg by Leg Observer's Trajectory[END_REF]. Section V presents a case where the order of the observer's kinematic model is greater than that of the target, and still the trajectory of the target is unobservable [START_REF] Jauffret | Observability: Range-Only vs. Bearings-Only Target Motion Analysis when the Observer Maneuvers Smoothly[END_REF]. We end with a more favourable case: the observer's trajectory contains an arc of a circle and this suffices to guarantee observability [START_REF] Jauffret | Observability: Range-Only vs. Bearings-Only Target Motion Analysis when the Observer Maneuvers Smoothly[END_REF]. The conclusion follows.

II. HYPOTHESES, DEFINITIONS AND NOTATIONS

A. Model of Target's and Observer's Kinematics

A target (T) and an observer (O) move in the same plane, given a Cartesian system. The target has a CV motion all along the scenario, while the observer maneuvers (the term "maneuver" is employed when the observer is not in CV motion).

For the observer, the position and velocity at time t are respectively

  t P O and     dt t dP t V O O  . Both are concatenated into the vector         T T T t V t P t X O O O



. For the target, the notations are similar:

  t P T ,   dt t dP V T T 
, and

      T T T T T T V t P t X  . Obviously,     T T T V t P t P   0 . Note that whatever t ,   t X T
entirely defines the target's trajectory. For a chosen * t , the state vector

  * t X T will be simply denoted   T y x y x X   
. The motion of the target relative to the observer is given by      



. We define the vector

      t X t X t X O T OT   . All
  N m m m m r r r r , 2 , 1 , , , ,  
. The aim of the RO-TMA is to estimate the state vector X given the measurements r . Note that in all the scenarios, the observer

starts from     T 0 0 0  O P
. The chosen estimate is the MLE (or equivalently the LSE since the noise is additive and Gaussian). Its empirical covariance matrix will be compared to the Cramér-Rao lower bound (CRLB), which is the inverse of the FIM for unbiased estimators.

III. THE OBSERVER'S TRAJECTORY CONSISTS OF ONE LEG

ONLY

In [START_REF] Pillon | Observability: Range-Only vs. Bearings-Only Target Motion Analysis for a Leg by Leg Observer's Trajectory[END_REF], we showed that when the observer is in CV motion, the trajectory of the target is observable up to an isometry (rotation or axial symmetry). We characterized the set of targets (true target and ghost-targets) O in the same range as the target of interest:

                           matrix orthogonal an being , 0 0 2 H H I O O OT OT V P V P O .
O can be defined by a 3-dimensional state vector

    T        2 2 * 2 * OT OT T OT OT V t P V t P Z . Indeed, for any * t , we have Z r k k A  2 , with     2 * * 1 t t t t k k k    A . Therefore, the range   X r k is reparametrized into   Z r k
. This property allows us to define a linear estimate.

A. Linear Estimate

Let us consider the measurement equation ( 1) and express the squares of its two elements:

2 , , 2 2 , 2 k r k r k k k m r r r      . The noise intervening in this new equation is then 2 , , 2 k r k r k r    whose first two moments are   2 2 , , 2 E r k r k r k r      and     2 2 2 2 , , 2 2 2 Var r k r k r k r k r r        (since k r,
 is Gaussian). We introduce then another "measurement":

2 2 , , r k m k m r    
.We define also the "noise":

2 2 , , 2 r k r k r k k r       
whose first two moments are now

  0 E  k  and     2 , 2 2 2 2 2 Var k r k r k r        
. So, we end up with the set of measurement equations

k k k k m Z r        k A 2 , for N k , , 1   . Note that k
 is zero-mean but no longer Gaussian. Still, it is legitimate to compute the weighted linear least squares estimator of Z that

minimizes the criterion         N k k m k Z Z C 1 2 , 2 , 1 k A    :                       N k k m k N k k WLS Z 1 , 2 , 1 1 2 , 1 1 ˆ     T k k T k A A A which is unbiased. Its covariance matrix is   1 1 2 2 2 2 1 1 2 , 2 1 1 1 4 1 ˆ                             N k k r k r N k k WLS r r Z Cov k T k k T k A A A A     . The computation of WLS Z ˆ assumes the knowledge of 2 k r . To alleviate this, in practice, we replace 2 k r by 2 ,k m r :                       N k k m k m N k k m WLS Z 1 , 2 , , 1 1 2 , , 1 1 ˆ     T k k T k A A A , with   2 2 , 2 2 , , 2 2 r k m r k m r       .

B. Efficiency of the Linear Estimate

The FIM is given by

  k T Z k Z N k r ROTMA r r Z     1 2 1  F with T k A k k Z r r 2 1   . Hence   k T k A A F    N k k r ROTMA r Z 1 2 2 1 4 1  and the CRLB is       1   Z Z ROTMA ROTMA F B
. The theoretical Cramér-Rao inequality is respected:

1 1 2 2 2 2 1 1 2 2 2 1 1 1 4 1 4                              N k k r k r N k k r r r r k T k k T k A A A A     , that is,     WLS ROTMA Z Cov Z  B 1 .

When

      i O O O V t P t P ,     
. We know from Proposition 3 of [START_REF] Pillon | Observability: Range-Only vs. Bearings-Only Target Motion Analysis for a Leg by Leg Observer's Trajectory[END_REF] that one ghost at most exists whose trajectory is defined by

      O O V V G X X X X O O       1 , 2 , 2 S I ,
where

O X stands for   * t X O , and   1 , 2 , 2 O O V V  

S

is the matrix of the symmetry around the line spanned by the vector

1 , 2 , O O V V  . To minimize                 N k k R k k m X r r X C 1 2 , ,



, we have recourse to the Gauss-Newton algorithm which will return

X ˆ or       O O V V G X X X X O O       1 , 2 ,

2

S I

, in unobservable cases. Two scenarios have been considered: in the first one, the target's trajectory is observable, and in the second one, proposed in [START_REF] Clark | A Gaussian mixture filter for Range-Only Tracking[END_REF], the target's trajectory is unobservable. For each scenario, 500 Monte Carlo simulations have been run.

For each simulation, we compute the MLE X ˆ and, if existing, the corresponding ghost G X ˆ. The performance of the MLE is evaluated by its empirical bias and the comparison of the empirical standard deviations of its components with the square roots of the diagonal elements of the CRLB.

A. Observable Scenario

The first leg of the observer is characterized by a heading of 45° and a speed of 4 m/s. The instant of its heading change is s 780 



. During its second leg, the speed and the heading of the observer are equal to 8.49 m/s and -70. ), and the standard deviation of the noise is 50  r  m. The estimates and the performance are given at the final time. Fig. 1 depicts the scenario and the 500 estimates (O and T are for the initial positions of the observer and the target, respectively). Because the FIM is singular [START_REF] Pillon | Observability: Range-Only vs. Bearings-Only Target Motion Analysis for a Leg by Leg Observer's Trajectory[END_REF], the CRLB cannot be computed; therefore in Table I the column " CRLB  " contains no value. We observe that the sets of MLEs are "croissant-shaped". This is why we present the results in Cartesian coordinates for the estimated target's position Numerous other simulations with various scenarios allow us to conclude that the MLE is unbiased. But we cannot claim that the MLE reaches the CRLB, because for such scenarios, the FIM is singular.

B. Unobservable Scenario(see [3])

In [START_REF] Clark | A Gaussian mixture filter for Range-Only Tracking[END_REF] the authors proposed the following scenario whose duration is 30 min to test a particle filter. The first heading and the speed of the observer are respectively -80° and 2.57 m/s. At t=15 min, it changes its course and the second heading is 146° and keeps the same speed. At the very beginning, the target is at

    (m) 071 , 7 071 , 7 0 T  T P
, its heading is -135° and its speed is 7.72 m/s (see Fig. 2(a)). The sampling rate is equal to 60 s and the standard deviation is 20  r  m. As in [START_REF] Clark | A Gaussian mixture filter for Range-Only Tracking[END_REF], the estimates are given at time t*=26 min. In Table II, we give the estimate of the state vector at this time. Fig. 2(b) depicts the 90%-confidence ellipsoid (in red) together with the cloud of the position estimates (in blue), and their corresponding ghost-estimates (in red). Here, despite the presence of a ghost, the MLE Note that adding a third leg does not guarantee observability [START_REF] Pillon | Observability: Range-Only vs. Bearings-Only Target Motion Analysis for a Leg by Leg Observer's Trajectory[END_REF].

A. Example of Observable Case

The target starts at   T 0 4000  (m), with a velocity equal to   T 2 10 (m/s). Note that the observer and the target are on a rendezvous route. In this case the FIM is singular [START_REF] Jauffret | Observability: Range-Only vs. Bearings-Only Target Motion Analysis when the Observer Maneuvers Smoothly[END_REF]. The estimates and the performance are given in Table III.

B. Example of Unobservable Case

The target starts at   T 3464 2000 (m) with a velocity   T 3 . 16 6 . 14 (m/s). The observer and the target are not on a rendezvous route and the bearings are not constant. Three ghosts exist (the maximum number of ghosts). Fig. 3 depicts the result: the estimates of the initial and final position of the target (T), and the couples of estimates of the three ghosttargets' positions (G). Note that the three other solutions are computed from the MLE returned by the algorithm [START_REF] Jauffret | Observability: Range-Only vs. Bearings-Only Target Motion Analysis when the Observer Maneuvers Smoothly[END_REF]. The estimate and its performance are given in Table IV. VI. THE OBSERVER'S TRAJECTORY CONTAINS AN ARC OF A CIRCLE Now the observer's trajectory is composed of two legs linked by an arc of a circle. The target's trajectory was proved to be observable in [START_REF] Jauffret | Observability: Range-Only vs. Bearings-Only Target Motion Analysis when the Observer Maneuvers Smoothly[END_REF]. No problem of convergence is indicated in our numerous simulations, provided the number of measurements during the arc of a circle is greater than 50 (otherwise, the routine can be stalled in a local minimum). In the scenario presented here, the initial position of the VII. CONCLUSION Three types of scenario were used to study the MLE in ROTMA: observable scenarios with nonsingular FIM, observable scenario with singular FIM and unobservable scenario with nonsingular FIM. For the first one, the MLE is efficient. For the second, the MLE is unbiased, but we cannot compute the CRLB. For the last, the MLE restricted to the local solution corresponding to the true target is efficient too. When several solutions exist (presence of ghost-targets), we are able to give the other solution(s) from the one returned by the algorithm. This study could be extended to the case of maneuvering targets as in [START_REF] Clavard | Target Analysis of a Source in a Constant Turn from a Nonmaneuvering Observer[END_REF].

.

  The duration of the scenario is1 LetA and B be two symmetric matrices,BA  means that the matrix A B  is a non-negative matrix. 26 min. The measurements are acquired every 4 s (

  relative position of the target w.r.t. the observer.

Fig. 2 .

 2 Fig. 2. (a) The scenario and the estimates in Clark's scenario [3], (b) magnification around the estimates.V. THE OBSERVER IS IN CONSTANT ACCELERATION MOTIONIn the two scenarios presented hereafter, the initial velocity vector of the observer is   T 2 10 (m/s). Its acceleration vector

Fig. 3 .

 3 Fig.3. The scenario and the estimates in an unobservable case.

  the angles are clockwise-positive:

	the angle 	North,	W		is denoted W  . The range and the
	bearing at time t are given respectively by   t r		P OT	  t	and
		  t		P OT 	  t	.
	B. Measurements Model and Estimate
				For the sake of simplicity, the notations k r and k  will
	stand for   k t r	and   k t 	, respectively. In order to emphasize
	the functional link of k r with the state vector X , we use the
	notation	r k	  X	. The range measurement collected at
	t k			  t k   1	and denoted	m r ,	k	obeys the equation
	m r	,	k		k r	  k X  	(1), where k  is the measurement noise
	assumed to be Gaussian, temporally white and zero mean; its
	standard deviation is	k  . The set of N available range r,
	measurements is

TABLE II PERFORMANCE

 II OF THE MLE AT TIME t*=26 MIN.

	X	Bias		CRLB	ˆ
	 x -1444.8 m	6.84	10.93	11.12
	 y -1444.8 m	15.21	12.84	12.86
	 x  -5.46 m/s	0.001	0.03	0.03
	 y  -5.46 m/s	0.002	0.04	0.04
	 r 455.05 m	0.49	11.09	11.99
	  165.31 °	0.11	1.60	1.64

  (m), with a speed equal to 2 m/s following a route equal to 180°. The total duration of the scenario is 26 minutes. The standard deviation is 50 m. The performance is summarized in TableV: we observe again that the MLE is efficient. Numerous simulations, not reported here, confirm this claim.

	observer is  	1000	 T 000 , 10	(m) with speed equal to 4 m/s.
	Its first heading is 90°. At time		1 	700	s, it begins its
	circular maneuver and finishes it at time		2 	980	s. At this
	time its new heading is equal to 225°. Meanwhile, the target
	has started at  	2000	 T 000 , 25

TABLE V PERFORMANCE

 V OF THE MLE AT FINAL TIME X

		Bias		CRLB	ˆ
	 x -2000 m	0.85	138.56	146.68
	 y 21888 m	0.25	49.58	53.04
	 x  0 m/s	0.003	0.14	0.15
	 y  -2 m/s	0.0005	0.04	0.05
	 r 15030 m	0.31	8.23	8.98
	  -17.35 °	0.0025	0.56	0.59