Audric Michelot

• Ste ´phanie Sarda
email: stephanie.sarda@iut-tlse3.fr

Catherine Audin

Eric Deydier

Eric Manoury

Rinaldo Poli

Christian Rey

Ste

Spectroscopic characterisation of hydroxyapatite and nanocrystalline apatite with grafted aminopropyltriethoxysilane: nature of silane-surface interaction

Heterogenised homogeneous catalysis is commonly performed with molecular catalysts grafted on solids via adsorption or via a covalent molecular link. Covalent grafting of organic groups on solid supports is usually carried out by silylation, using functionalised trialkoxysilanes. Among these solids supports, very few studies have been published on apatites. In the present work, aminopropyltriethoxysilane (APTES) grafting was performed in toluene on different apatitic supports: crystallised stoichiometric hydroxyapatites differing by the drying method, freeze-dried (HAP) and dried at 100 °C (HAPD), and a nanocrystalline apatite. All materials were fully characterised, before and after grafting, for better understanding of the nature of the alkoxysilane/surface interaction. The data show a clear competition between the covalent grafting of APTES and its polycondensation reaction, depending on the nature of the solid support surface. Silylation is accompanied by APTES covalent grafting to oxygen atom of the hydroxyl groups of the apatitic structure and/or of the OH -species that are present on the surface hydrated layer. This work clarifies the nature of silane grafting onto selected apatitic surfaces and

Introduction

The development of supported molecular catalyst is nowadays a major field in catalysis research. Heterogeneous catalysis appears as an alternative method in organic synthesis as it presents several advantages over homogeneous processes, including simple product isolation as well as catalyst separation and recycling [START_REF] Fechete | The past, present and future of heterogeneous catalysis[END_REF]. Heterogenised homogeneous catalysis is commonly performed with molecular catalysts grafted on solids via adsorption or via a covalent molecular link. The interest of supported catalysis is to combine the activity and selectivity of homogeneous catalysts with the easy recovery of solids, which is of particular interest in asymmetric catalysis. The solid supports may be organic (polymers) or inorganic.

The most commonly used inorganic solid for supported heterogeneous catalysis is silica, but other inorganic solids such as metallic oxides or ceramics [START_REF] Fechete | The past, present and future of heterogeneous catalysis[END_REF], and more recently apatite calcium phosphates (or apatites), have also been reported as suitable and very promising supports. Among these different solids, very few studies have been published on apatites. Apatites constitute the mineral part of calcified tissues, such as bone and teeth, and are widely used in biomedical applications such as implant coatings, bioceramics, or bone cements thanks to their excellent biocompatibility, bioactivity, and osteoconductive properties [START_REF] Drouet | Nanocrystalline apatites: from powders to biomaterials[END_REF][START_REF] Sarda | Kinetic study of citric acid influence on calcium phosphate bone cements as water-reducing agent[END_REF][START_REF] Grossin | Biomimetic apatite sintered at very low temperature by spark plasma sintering: physico-chemistry and microstructure aspects[END_REF]. Moreover, they present high ion-exchange ability and adsorption capacity. Thus, modified apatites have been investigated not only for biomaterial applications but also as delivery systems for molecules of biological interest such as proteins [START_REF] Yin | Impacts of the surface charge property on protein adsorption on hydroxyapatite[END_REF], growth factors [START_REF] Autefage | Adsorption and release of BMP-2 on nanocrystalline apatite-coated and uncoated hydroxyapatite/btricalcium phosphate porous ceramics[END_REF], and drugs [START_REF] Al-Kattan | Medical potentialities of biomimetic apatites through adsorption, ionic substitution, and mineral/organic associations: three illustrative examples[END_REF], as adsorbents for heavy metals in environmental applications [START_REF] Silva | Hydroxyapatite organofunctionalized with silylating agents to heavy cation removal[END_REF], in chromatography column separations systems [START_REF] Kawasaki | Hydroxyapatite as a liquid chromatographic packing[END_REF], and as supports for heterogeneous catalysis [START_REF] Tahir | Zinc bromide supported on hydroxyapatite as a new and efficient solid catalyst for Michael addition of indoles to electron-deficient olefins[END_REF][START_REF] Smahi | Potassium fluoride doped fluorapatite and hydroxyapatite as new catalysts in organic synthesis[END_REF].

Covalent grafting of organic groups on solid support is usually carried out by silylation, using functionalised trialkoxysilanes (OR 0 )3SiR. R is an alkyl-or aryl-functionalised chain and Si(OR 0 ) is a hydrolysable function, R 0 being most frequently a methyl or an ethyl group. On the one hand, the Si(OR 0 ) bond is easily hydrolyzed in aqueous solution with the elimination of alcohol (R 0 -OH) to form metastable silanol (-SiOH), which undergoes condensation to more stable siloxane (Si-O-Si), forming three-dimensional hyperbranched oligomers [START_REF] Damia | Study of two hydroxyapatite/poly(alkoxysilane) implant coatings[END_REF]. On the other hand, with solids bearing surface OH groups (solid-OH), a monolayer or covalently grafting solid-OSi(OH)2R can be achieved in anhydrous media [START_REF] Durrieu | Grafting RGD containing peptides onto hydroxyapatite to promote osteoblastic cells adhesion[END_REF]. This is a well-known and efficient strategy for silica-controlled grafting. The use of non-polar organic solvents during the surface modification of silica surfaces with 3-aminopropyltriethoxysilane (APTES) leads to the formation of stable surface layers by covalent interactions [START_REF] Kim | Formation, structure, and reactivity of amino-terminated organic films on silicon substrates[END_REF][START_REF] Simon | Study of two grafting methods for obtaining a 3-aminopropyltriethoxysilane monolayer on silica surface[END_REF].

Among compounds of the apatitic calcium phosphate family, hydroxyapatite (HAP), Ca10(PO4)6(OH)2, is a widely studied compound. It is the stoichiometric apatitic calcium phosphate (Ca/P atomic ratio of 1.67) which is characterised by low solubility in water, low cost, and high innocuousness. Moreover, thanks to the presence of structural hydroxyl groups, HAP is a good candidate for surface functionalisation using alkoxysilane agents [START_REF] Silva | Hydroxyapatite organofunctionalized with silylating agents to heavy cation removal[END_REF][START_REF] Goonasekera | Attachment of poly(acrylic acid) to 3-aminopropyltriethoxysilane surface-modified hydroxyapatite[END_REF]. APTES is a widely used reagent for surface functionalisation of inorganic solid [START_REF] Fechete | The past, present and future of heterogeneous catalysis[END_REF], bearing both a trialkoxysilane group (surface anchor) and an amine function allowing an easy coupling with various organic molecules via nucleophilic substitution. This reagent has been used with success for silica and HAP grafting for immobilisation of a large variety of molecules or proteins [START_REF] Durrieu | Grafting RGD containing peptides onto hydroxyapatite to promote osteoblastic cells adhesion[END_REF]. Recently published papers reported that APTES reacts with HAP via condensation of the hydrolyzed ethoxy moieties with the hydroxyl groups found at the surface of the HAP [START_REF] Durrieu | Grafting RGD containing peptides onto hydroxyapatite to promote osteoblastic cells adhesion[END_REF][START_REF] Goonasekera | Attachment of poly(acrylic acid) to 3-aminopropyltriethoxysilane surface-modified hydroxyapatite[END_REF]. Some studies have also highlighted the influence of the reaction conditions, such as time, temperature [START_REF] Goonasekera | Attachment of poly(acrylic acid) to 3-aminopropyltriethoxysilane surface-modified hydroxyapatite[END_REF], or the morphology [START_REF] Vasiliev | Chemisorption of silane compounds on hydroxyapatites of various morphologies[END_REF] of apatitic crystals on the interaction between APTES and the HAP surface. However, none of these published works reported the effect of the chemical nature and crystallinity of the apatitic supports. Moreover, the use of nonstoichiometric nanocrystalline apatites (NCA, Ca10-x (PO4)6-x (HPO4 or CO3)x (OH)2-x with 0 B x B 1 and a Ca/P molar ratio ranging from 1.50 to 1.67), which present a greater surface area than HAP and a structured hydrated layer on their surface, rich in surface labile species able to be exchanged [START_REF] Eichert | Nanocrystalline apatite-based biomaterials: synthesis, processing and characterization[END_REF], has been poorly studied [START_REF] Liu | Surface modification of nano-apatite by grafting organic polymer[END_REF][START_REF] Toworfe | Nucleation and growth of calcium phosphate on amine-, carboxyl-and hydroxyl-silane self-assembled monolayers[END_REF].

The aim of this work is to address these two points: to explore the influence of APTES grafting on the chemical composition and the crystallinity of two calcium phosphate materials (crystallised stoichiometric hydroxyapatite, HAP, and a nanocrystalline apatite, NCA) and to characterise the trialkoxysilane-solid surface interaction that is still poorly understood. APTES grafting was performed in toluene on both materials. All materials were fully characterised, before and after grafting, by different techniques and especially by FTIR, Raman, and NMR spectroscopy.

Experimental section

Materials

The well-crystallised stoichiometric hydroxyapatite powders were synthesised by a double decomposition technique between (NH4)2HPO4 (433.3 mL, 1.54 M in deionised water) and Ca(NO3)2-4H2O (1100 mL, 1 M in deionised water) at 90 °C and pH around 10 for about 3 h [START_REF] Trombe | Contribution a ` l'e ´tude de la de ´composition et la re ´activite ´ de certaines apatites hydroxyle ´es carbonate ´es ou fluore ´es alcalino-terreuses[END_REF]. After filtration and washing, the precipitate was either freezedried (HAP) or dried at 100 °C for 5 h in a vacuum to remove the residual physisorbed water (HAPD). Then the powder was sieved (\125 lm) and stored in a freezer.

Nanocrystalline apatite powder (NCA) was prepared by a double decomposition technique at ambient temperature (pH close to 7) between (NH4)2HPO4 (1500 mL, 0.60 M in deionised water) and Ca(NO3)2-4H2O (750 mL, 0.29 M in deionised water) [START_REF] Rey | Maturation of poorly crystalline apatites: chemical and structural aspects in vivo and in vitro[END_REF]. The precipitate obtained in few minutes was filtered, washed, and freeze-dried. Then the powder was sieved (\125 lm) and stored in a freezer.

(3-Aminopropyl)triethoxysilane (APTES, 98 %, Sigma-Aldrich) was used without further purification.

Grafting experiments

The grafting of APTES on HAP, HAPD, and NCA was carried out according to published methods [START_REF] Vasiliev | Chemisorption of silane compounds on hydroxyapatites of various morphologies[END_REF][START_REF] Sakeye | Comparison of different aminofunctionalization procedures on a selection of metal oxide microparticles: degree of modification and hydrolytic stability[END_REF]. Under an inert argon atmosphere, anhydrous toluene (20 mL) was added to the apatitic particles (1 g) together with APTES (2.33 mL). The mixture was refluxed at 110 °C overnight. The solid products were separated by centrifugation, washed three times with toluene and twice with ether, and dried at 70 °C under vacuum overnight. The obtained APTES-grafted apatites were designated as HAP-Si, HAPD-Si, and NCA-Si.

For each apatitic support, blank experiments were performed under the same experimental conditions without APTES; these blanks were designed as HAP0, HAPD0, and NCA0.

The same experiments were performed by adding water to the toluene suspension containing HAP and APTES to favor polymerisation of the surface agent. The sample obtained was designated as HAP-Si-Poly.

Characterisation techniques

Elemental analyses were performed with a Perkin-Elmer 2400 instrument. The phosphate content of the synthesised apatites was determined by UV-visible spectroscopy of the phosphovanadomolybdic acid complex and the calcium content by complexometry with ethylenediaminetetraacetic acid [START_REF] Charlot | Chimie analytique quantitative[END_REF]. The estimated standard deviations of calcium and phosphate titrations were equivalent to 1 % of the reading. The specific surface area was measured by the Brunauer-Emmett-Teller (BET) method with a Quantachrome Instruments Monosorb Nova 1000 with an estimated standard deviation of 0.5 m 2 /g.

The amount of water in the powders was measured by thermogravimetric analysis (TGA-DTA) using TA instruments Q600 from 25 °C up to 1000 °C at a heating rate of 5 °C/min under an air flux (30.00 ± 0.01 mg of each sample).

The solid phases before and after grafting APTES were characterised by X-ray diffraction (XRD) (BRUKER D8-2 diffractometer CuKa1 radiation with k = 1.5406 A ˚ ). The apparent crystallite sizes were calculated from the (310) and (002) lines by applying Scherrer's formula [25].

Transmission electron microscopy (TEM) micrographs of the powders before and after grafting were recorded in the bright-field imaging on a JEOL JEM 2100F microscope set at an acceleration tension of 200 kV.

Unmodified and modified particles were characterised by several spectroscopic techniques. Fourier transform infrared (FTIR) analyses were performed with an Impact 400 Nicolet FTIR spectrometer. All adsorption spectra were obtained by accumulation of 32 scans at a 4 cm -1 resolution covering the 4000-400 cm -1 range (2 mg sample/300 mg KBr).

Raman spectroscopy analyses were carried out with a Horiba, Jobin-Yvon Labram HR800 confocal microspectrometer (k = 532 nm). The spectral resolution was 2 cm -1 , the laser was power 2 mW, integration times varied from 30 to 60 s, and spectral accumulations from 2 to 30, over a wavenumber range of 100-4000 cm -1 . 29 Si, 13 C, and 31 P cross-polarisation magic angle spinning (CP-MAS) solid-state NMR spectra were recorded at 79.39, 100.49, and 161.74 MHz, respectively. 1 H MAS NMR spectra were recorded at 399.61 MHz. All experiments were realised with a spinning rate of 8 kHz on a Bruker Avance 400 spectrometer (9.4 T).

Results

Characterisation of the apatitic samples

Chemical analyses of the synthesised solids gave a Ca/P ratio of 1.66 ± 0.02 for the HAP and HAPD samples corresponding to stoichiometry hydroxyapatite (Ca10 (PO4)6(OH)2), whereas NCA displayed a Ca/P ratio of 1.40 ± 0.02, characteristic of a calcium-deficient apatite.

The water content was determined by thermogravimetric analysis (TGA) in the temperature range of 25-1000 °C. TGA showed a two-stage weight loss of water between 25 and 400 °C, attributed to adsorbed and lattice water on HAP [START_REF] Ashok | Crystallization of hydroxyapatite at physiological temperature[END_REF][START_REF] Heughebaert | Contribution a ` l'e ´tude de l'e ´volution des othophosphates de calcium pre ´cipite ´s amorphes en orthophosphates apatitiques[END_REF]. As expected, due to its hydrated surface layer, NCA presented the highest water content (15.00 ± 0.20 wt%). It decreased from HAP (4.00 ± 0.04 wt%) to HAPD (1.50 ± 0.02 wt%) due to the different drying procedure.

The specific surface areas of the stoichiometric apatites HAP and HAPD are lower than that of NCA: 58 (±3), 52 (±3), and 101 (±5) m 2 /g, respectively. These results are in accordance with TEM observations (Fig. 1). Both samples crystallise as small needles (longer needles for HAP samples than for NCA) with a thickness of fewer than 100 nm.

The HAP and HAPD powders exhibited well-defined XRD patterns (Fig. 2a) corresponding to well-crystallised apatitic calcium phosphates, whereas NCA exhibited broad peaks due to the low crystallinity of the apatitic phase [START_REF] Eichert | Nanocrystalline apatite-based biomaterials: synthesis, processing and characterization[END_REF]. No additional crystalline phase was detected. The apparent crystallite sizes were calculated from XRD using Scherrer's equation [25] from the (310) and (002) lines, and the results are reported in Table 1. All the crystals were elongated along the c-axis of the hexagonal structure in accordance with the TEM observations. The crystallites sizes observed for NCA appear lower than for hydroxyapatites samples, in accordance with TEM observations (Fig. 1).

The FTIR spectra of HAP and HAPD showed the characteristic bands of apatitic calcium phosphate (Fig. 3a): apatitic PO4 3-vibration bands at 473 cm -1 (m2), 560-630 cm -1 (m4), 962 cm -1 (m1), and 1000-1100 cm -1 (m3) and OH bands at 630 and 3560 cm -1 , characteristic of an apatitic environment [START_REF] Elliott | Structure and chemistry of the apatites and other calcium orthophosphates[END_REF]. For HAP and HAPD, a weak and broad band between 1420 and 1470 cm -1 due to the low number of carbonate ions and a weak band at 875 cm -1 assigned to both carbonate and possibly HPO4 2- ions were observed, suggesting that the powders could be a slightly calcium-deficient hydroxyapatite. As expected, these CO3 2-and HPO4 2-vibration bands, as well as the free water vibration bands (1630 and 3400 cm -1 ), were much more intense for the NCA particles. Moreover, the FTIR spectrum of NCA exhibited broader apatitic phosphate bands, characteristic of a poorly crystalline apatite [START_REF] Eichert | Nanocrystalline apatite-based biomaterials: synthesis, processing and characterization[END_REF] with weak carbonate bands. Curve fitting analysis in the m4PO4 domain using the GRAMS/386 software (Fig. 3b) confirmed the formation of biomimetic nanocrystalline apatite with non-apatitic environments of PO4 3-(620 cm -1 ) and HPO4 2-(530 cm -1 ) ions present in the hydrated surface layer [START_REF] Rey | Nanocrystalline apatites in biological systems: characterisation, structure and properties[END_REF]. Weak OH vibrations were detected on the curve-fitted spectrum of NCA suggesting a very low hydroxide ion content related to the calcium deficiency [START_REF] Heughebaert | Contribution a ` l'e ´tude de l'e ´volution des othophosphates de calcium pre ´cipite ´s amorphes en orthophosphates apatitiques[END_REF].

The Raman spectra of HAP (Fig. 4a), HAPD, and NCA synthesised powders presented characteristic lines of apatitic phosphate groups in accordance with the FTIR studies: corresponding to orthophosphate species [START_REF] Yesinovski | Hydrogen environments in calcium phosphates: 1 H MAS NMR at high spinning speeds[END_REF]. For the NCA sample (data not shown), a characteristic broader asymmetric signal was observed due to the presence of HPO4 2- species. The broadening of the signal could also be attributed to the nanoscopic nature of the sample [START_REF] Ja ¨ger | A solid state NMR investigation of the structure of nanocrystalline hydroxyapatite[END_REF]. The MAS 1 H solid NMR spectrum of HAP (Fig. 5b) displayed a signal at -0.2 ppm attributed to the OH groups of the apatitic structure and a broader one at 5.1 ppm attributed to adsorbed water molecules [START_REF] Eichert | Nanocrystalline apatite-based biomaterials: synthesis, processing and characterization[END_REF][START_REF] Yesinovski | Hydrogen environments in calcium phosphates: 1 H MAS NMR at high spinning speeds[END_REF]. NCA displayed more complex spectra (data not shown) as this structure contains a higher concentration of water and hydrogen phosphates both giving broad signals at around 5 and 14 ppm, respectively. Several lines were also observed between 2 and 5 ppm; these were also mentioned by other groups, but their structural origin remains unclear [START_REF] Ja ¨ger | A solid state NMR investigation of the structure of nanocrystalline hydroxyapatite[END_REF].

Characterisation after silylation

In order to analyze carefully the structural changes induced by the silylation process, three blanks were prepared by refluxing HAP, HAPD, and NCA solids in toluene overnight, without grafting agent, giving rise to HAP0, HAPD0, and NCA0, respectively, after washing and drying. Chemical analyses of the blanks showed traces of carbon Wavenumber (cm -1 ) Wavenumber (cm -1 )

(Table 1) because of toluene adsorption during the reaction. As expected, the carbon, hydrogen, and nitrogen content increased after silylation (grafted APTES), the maximum amount being observed for NCA-Si due to its high specific surface area. One can notice that the initial nitrogen content in the HAP sample disappeared in HAP0 and HAPD. This value can be attributed to ammonia residues from the HAP synthesis. The initial hydrogen content in NCA (about 0.52 wt%) can be attributed to the presence of HPO4 2-ions, as observed by FTIR spectrum, which could slightly affect the amount of water calculated from thermogravimetric analysis. Moreover, carbonate traces in the initial powders confirmed by carbon analysis show that carbonate species still remain after washing and thus could also affect the calculated amount of water.

The XRD diagrams of the blanks (HAP0, HAPD0, and NCA0) displayed better defined diffraction patterns than those of the initial samples (HAP, HAPD, and NCA), as can be seen in Fig. 2b for HAP, in agreement with the increase of the apparent crystallite sizes L(002) and L(310) (Table 1). One can also notice that the crystallite sizes of the samples after silylation decreased compared with the blanks, whatever the initial support. However, the changes observed in the line broadening could also be related to residual stress and the crystallinity of the crystals.

The FTIR spectra of the blanks (HAP0, HAPD0, and NCA0, Fig. 6a) were compared with those of the initial samples (Fig. 3a). For all blanks, a decrease of the HPO4 2- vibration band at 875 cm -1 which can be attributed to HPO4 2-or CO3 2-species (no band observed for the HAPD0 and HAP0 blanks) and a splitting of the PO4 3- vibration bands were observed, characteristic of a better apatitic environment in agreement with XRD results. Moreover, for the HAP0 and NCA0 samples, a strong decrease of the intensity of vibration bands attributed to the carbonate species was observed relative to the HAP and NCA samples, probably attributable to the increase of the apatitic environment [START_REF] Eichert | Nanocrystalline apatite-based biomaterials: synthesis, processing and characterization[END_REF]. At the same time, the decrease of the free water vibration bands could be attributed either to the effect of the temperature or to the influence of toluene, which is known to make an azeotrope with water [START_REF] Gomis | Homogeneity of the water ? ethanol ? toluene azeotrope at 101.3 kPa[END_REF]. For all blank samples, a significant narrowing of the OH -vibration bands at 630 and 3560 cm -1 (or appearance in the case of NCA0) was observed and attributed to OH group formation by drying. Indeed, OH -hydrogen bonding with adsorbed water on apatitic surfaces has been reported to produce signal broadening [START_REF] Blumenthal | Hydroxyapatite: mechanism of formation and properties[END_REF].

After silylation, all FTIR spectra (Fig. 6b) exhibited broader apatitic phosphate bands without shifting and a decrease of the relative intensity of OH -vibration bands of the apatitic structure (disappearance for NCA). For the NCA sample, which exhibited higher silane content, specific APTES bands were observed: new broad bands at 2950 and 2850 cm -1 ascribed to asymmetric and symmetric stretching modes of CH2 and CH3 groups of the aminopropyl chain, respectively [START_REF] Silva | Hydroxyapatite organofunctionalized with silylating agents to heavy cation removal[END_REF]; weak bands at 1460 and 1330 cm -1 corresponding to asymmetric and symmetric C-H bending deformation bands, respectively [START_REF] Silva | Hydroxyapatite organofunctionalized with silylating agents to heavy cation removal[END_REF]; and a series of bands related to N-H stretching modes of the primary amine in the range of 3500-3000 cm -1 and a band owed to the N-H group deformation at 1620 cm -1 , both of which can be overlapped by water vibrations [START_REF] Silva | Hydroxyapatite organofunctionalized with silylating agents to heavy cation removal[END_REF][START_REF] Goonasekera | Attachment of poly(acrylic acid) to 3-aminopropyltriethoxysilane surface-modified hydroxyapatite[END_REF][START_REF] Domingo | Grafting of trialkoxysilane on the surface of nanoparticles by conventional wet alcoholic and supercritical carbon dioxide deposition methods[END_REF]. The characteristic bands of grafted silane were less intense on the HAP-Si and HAPD-Si spectra than on the NCA-Si spectrum because of the lower quantities of grafted silane. The lower silane content in HAP-Si and HAPD- Si was confirmed by microanalysis (Table 1). The expected Si-O-Si, Si-O-C bands (around 1070 cm -1 ), and Si-O-P bands (998 cm -1 ) could not be observed as they overlapped with strong orthophosphate bands near 1080 cm -1 [START_REF] Goonasekera | Attachment of poly(acrylic acid) to 3-aminopropyltriethoxysilane surface-modified hydroxyapatite[END_REF].

Curve fitting in the m4PO4 domain of the FTIR spectra was carried out for the blanks and the apatitic supports before and after silylation using the GRAMS/386 program (Galactic Industries, Salem, NH). The band shape was considered as Lorentzian and the baseline linear in all instances. Parameters such as the position and width of the bands were varied to provide the best curve fit. The relative band integrated intensities were defined as the ratio between the peak area and the total area of the phosphate bands. Figure 7 reports the relative band integrated intensity for the apatitic OH - species. The relative intensity of the OH bands increased in the blanks which presented the highest level of adsorbed water, for HAP0 and especially for NCA0, in agreement with the increase of the apatitic domains observed by XRD. However, the OH -band intensity decreased for the two hydroxyapatite samples after silylation, for HAP and especially for the HAPD support which presented the lowest water content and thus the lowest possibility of polycondensation. These results confirm the possibility of silane molecule interaction with OH -on apatites.

The Raman spectra of the blanks showed the characteristic vibration bands associated with the orthophosphate, carbonate, and hydroxyl vibrations of the apatitic samples. The relative intensity of the apatitic OH -vibration band at 3570 cm -1 appeared stronger in the blanks compared with the initial samples, as observed in Fig. 4A for HAP and HAP0 (data not shown for HAPD and NCA). In all spectra, a shift of the m1 (960 cm -1 ) and m3 (1040-1074 cm -1 ) modes of PO4 3-and of the m1 (1073 cm -1 ) mode of CO3 2-was observed for the samples after silylation and for the blanks compared with the initial powders. This can be attributed to the adsorption of toluene, as confirmed by the microanalysis results. However, no hydroxyl band was observed for the NCA0 samples (data not shown). After silylation, the Raman spectra of HAP-Si and HAPD-Si (Fig. 4b) showed a decrease of the OH -vibration band intensity at 3570 cm -1 , in accordance with the FTIR results. For all apatitic supports, the carbonate m1 vibration band (1073 cm -1 ) decreased with silane grafting, and the orthophosphate vibration bands appeared broader, which can be attributed to an amorphous structure or to the presence of overlapping organic vibration bands. This is particularly obvious for HAP as shown in Fig. 4a. In Fig. 4b, the spectra of all samples after silylation show an additional broad band at 2900 cm -1 which can be attributed to the CH stretching modes of the CH2 and CH3 groups, whereas lower frequency bands between 1300 and 1400 cm -1 may correspond to CH vibration modes [START_REF] Bistric ˇic ´ L | Conformational and vibrational analysis of gamma-aminopropyltriethoxysilane[END_REF]. Low-intensity broad bands were observed at about 3300 and 1600 cm -1 related to the presence of amino groups. Other small bands between 600 and 650 cm -1 were owed to the Si-O-C stretching mode of trialkoxy groups [START_REF] Bistric ˇic ´ L | Conformational and vibrational analysis of gamma-aminopropyltriethoxysilane[END_REF][START_REF] Ishida | The structure of aminofunctional silane coupling agents: 1. c-Aminopropyltriethoxysilane and its analogues[END_REF], easily observable in the case of the polymerised silane sample HAP-Si-poly obtained by adding water to the toluene. Moreover, a small band appeared at about 1003 cm -1 for the samples after grafting (Fig. 4b). Bands in this region of the spectrum have been attributed in previous studies to the vibrations of the Si-O-P bonds that can form during the reaction of surface phosphate hydroxyl groups with trimethoxysilane groups [START_REF] Vasiliev | Chemisorption of silane compounds on hydroxyapatites of various morphologies[END_REF]. However, this band is still present in all blanks' spectra (HAP0 in Fig. 4a, data not shown for HAPD0 and NCA0) and thus cannot be attributed to Si-O-P vibrations. Rather, it appears to be associated with the presence of toluene, probably due to CH vibration bands. Moreover, no hydrolysis was detected by Raman spectroscopy: OH stretching modes of unassociated SiOH groups could not be found in the 3700-3600 cm -1 region, nor Si-O bending vibration of the SiOH group at 800-650 cm -1 [START_REF] Shih | Raman studies of the hydrolysis of silane coupling agents[END_REF].

The solid-state MAS 1 H NMR spectra (Fig. 5b) of the different samples present a signal at -0.2 ppm attributed to OH groups. The broader line at 5.1 ppm (adsorbed water molecules on the surface layer) observed in HAP spectrum before grafting dramatically decreased upon grafting (HAPD-Si, HAP-Si). As the grafting reaction was realised in refluxing dry toluene, which is known to make an azeotrope with water, we suppose that most of the adsorbed water was removed in the initial drying step. For the NCA-Si sample, the NMR spectrum showed numerous broad signals between 15 and -5 ppm due to a great variety of protons present in the nanocrystals.

The solid-state CP-MAS 31 P NMR spectra of the synthesised apatites before and after silane grafting showed no significant modification (Fig. 5a for HAP). Only one signal was observed at 3.0 ppm due to orthophosphate species of apatites, as the low number of surface phosphates could not be differentiated from the massive phosphates of the particle core.

The solid-state CP-MAS 13 C NMR spectra of HAP-Si and HAPD-Si samples after silylation are shown in Fig. 5c. Three distinct lines observed at 10.4, 22.3, and 43.4 ppm were assigned to carbons C1, C2, and C3: SiCH2, CH2, and CH2NH2 of the aminopropyl chain, respectively [START_REF] Goonasekera | Attachment of poly(acrylic acid) to 3-aminopropyltriethoxysilane surface-modified hydroxyapatite[END_REF][START_REF] Sharma | Silica-supported molybdenum complex: a novel, selective and reusable organicinorganic hybrid catalyst for eco-friendly oxidation of sulfides and olefins[END_REF]. These signals were more intense in the case of the polymerised silane sample (HAP-Si-poly) obtained by adding water to the toluene, as the carbon content increased, with a slight shift to 27.0 and 44.9 for C1 and C2 respectively. Moreover, traces of non-hydrolyzed Si-OEt and ethanol were also observed at 58.1 ppm (-CH2OH), 56.3 ppm (-CH2OSi), and 18.6 ppm (CH3) in the polymerised silane sample (HAP-Si-poly). Finally, the signal at around 164.2 ppm is attributed to carbonate species present in the apatitic samples [START_REF] Eichert | Nanocrystalline apatite-based biomaterials: synthesis, processing and characterization[END_REF]. These results confirm the presence of APTES on the surface but do not allow conclusions to be drawn on the nature of the silane-solid support interaction.

(a) (b)

The solid-state CP-MAS 29 Si NMR spectra of the solids after silylation (Fig. 5d HAPD-Si and NCA-Si) showed signals at -50, -60, and -67 ppm attributed, respectively, to the T 1 , T 2 , and T 3 siloxane structures illustrated in Scheme 1 [START_REF] Goonasekera | Attachment of poly(acrylic acid) to 3-aminopropyltriethoxysilane surface-modified hydroxyapatite[END_REF][START_REF] Cardiano | Epoxy-silica polymers as restoration materials[END_REF]. In the case of the grafted dried HA (HAPD-Si), the signal at -20.3 was attributed to traces of polysiloxane grease. As expected, the signal of the T 3 structure grew in intensity, relative to those of T 1 and T 2 , on the polymerised sample HAP-Si-poly where this structure is favored. A similar observation can be made for the HAPD-Si and NCA-Si samples but with a lower difference of relative intensity ratio T 3 /T 2 compared with the polymerised system. These observations confirm that the T 3 structure is favored by the presence of water as silane hydrolysis allows its polymerisation. This is the reason why alkoxy silane grafting on silica is usually realised under anhydrous conditions to favor alkoxysilane-surface silanol interaction and avoid silane polymerisation.

Discussion

The data show a clear competition between the covalent grafting of APTES and its polycondensation reaction, in spite of the use of anhydrous organic solvent to avoid surface polymerisation. The TGA data evidenced the presence of water in the apatitic samples, adsorbed water, or lattice water, more important in the case of the nanocrystalline apatite (NCA sample) than for the hydroxyapatites HAP and HAPD samples (from 4 to 15 wt%) because of the presence of the surface hydrated layer, which participated in the APTES hydrolysis and/or polycondensation. Even after HA pre-drying (HAPD sample), such hydration persisted, not allowing the formation of Si-O-Si related to the polycondensation reaction not to be completely avoided, as the Si-O-Si vibration bands remained present in all spectra. The silane content increased when water content increased in the initial apatitic sample (from HAPD to NCA). No metastable silanol was detected, with its subsequent consumption leading to either grafting or oligomerisation to more stable siloxane (Si-O-Si) probably being immediate.

The blank experiments, performed under the same experimental conditions without surface modification agent, evidence the role of the organic solvent in the silane grafting. The spectroscopic data show a decrease of water content and hydrogen bonding for all blank samples relative to the initial powders. Indeed, toluene is widely used as a carrier for the dehydration of alcohol by azeotropic distillation [START_REF] Gomis | Homogeneity of the water ? ethanol ? toluene azeotrope at 101.3 kPa[END_REF] and thus could extract water from the apatitic samples. This behavior favors the formation of a silane layer by covalent grafting instead of polycondensation. Moreover, a decrease of the line broadening was observed by XRD in the blanks. This phenomenon, never reported in

+ H 2 O -R'OH R R'O Si R'O OH + H 2 O -R'OH R HO Si R'O OH + H 2 O -R'OH R HO Si HO OH R R'O Si Si R'O O R HO Si Si R'O O R HO Si Si HO O R Si Si Si O R'O O R Si Si Si O HO O R Si Si Si O O O

Si

Scheme 1 Different ways of hydrolysis and oligomerisation of trialkoxysilanes the literature, could be explained by an increase of the crystallite sizes and/or a change on the crystallinity and residual stress in the crystals. These interesting results will be addressed in a further publication.

Although it appears that APTES polycondensation occurs whatever the nature of the apatitic sample, the spectroscopic results evidence at the same time the presence of covalent grafting. APTES covalent grafting is always reported to take place with solid samples through surface OH groups (Solid-OH). Some authors have proposed silane covalent grafting in HAP samples through HPO4 groups on the basis of spectroscopic results [START_REF] Vasiliev | Chemisorption of silane compounds on hydroxyapatites of various morphologies[END_REF]. This hypothesis was supported by the model proposed by Brown for hydroxyapatite structure [START_REF] Brown | Roles of octacalcium phosphate in surface chemistry of apatites[END_REF], in which an octacalcium phosphate (OCP, Ca8(PO4)4(HPO4)2-5H2O) layer is present at the surface of the apatitic crystals. The OCP structure can be described as a succession of 'apatite' and 'hydrated' layers, the latter being rich in HPO4 2-ions. OCP could act as a precursor for the formation of HAP with acicular morphology in which the (100) type faces are dominant because they have low interfacial energy. However, in the current study, the IR vibration bands attributed to Si-O-P groups at about 990-1000 cm -1 [START_REF] Wang | Fabrication and characterization of nano-composite scaffold of PLLA/silane modified hydroxyapatite[END_REF][START_REF] Zhang | Interfacial fabrication and property of hydroxyapatite/polylactide resorbable bone fixation composites[END_REF] cannot be observed, probably because they are overshadowed by authors have studied silane covalent grafting with hydroxyl groups on the HAP surface without hydrolysis and condensation reaction [START_REF] Zhang | Molecular dynamics simulations on the interaction between polymers and hydroxyapatite with and without coupling agents[END_REF]. In the hexagonal structure of HAP (P63/m space group), the OH -ions are present inside tunnels (3 A ˚ of diameter) with the direction parallel to the caxis; each OH -ion is coordinated with the three Ca atoms which are arranged in a triangle on the same plane. The OH groups are either above or below the Ca triangle. Depending on the surface crystallographic planes, the binding energies and the planar atomic density differ and have been calculated by Zhang et al. [START_REF] Zhang | Molecular dynamics simulations on the interaction between polymers and hydroxyapatite with and without coupling agents[END_REF] using molecular dynamic simulations: the (100) and (110) planes exhibit the highest planar OH -density (atom 0.072 A ˚ -2 ), and the (110) plane displays the highest planar energy for methacryloxypropyl trimethoxy silane. In addition, for NCAs, Blumenthal and Posner [START_REF] Blumenthal | Hydroxyapatite: mechanism of formation and properties[END_REF] have suggested that some OH -ions present in the solids cannot be observed by FTIR spectroscopy because of hydrogen bonding with water molecules in the hydrated surface layer. Thus, the results obtained lead to the conclusion that silylation can be explained at the same time by APTES polycondensation and covalent grafting to oxygen atom of the hydroxyl groups of the apatitic structure and/or the OH -species that are present on the surface hydrated layer according to Eq. 2:

Apatite OH OEt + EtOSi(CH2)3NH2 OEt O Si(CH2)3NH2 + EtOH ð2Þ
Si-O-Si, Si-O-C (1070 cm -1 ), and PO4 3-vibration bands at 1000-1100 cm -1 (no shift of the PO4 3-vibration bands was observed). Moreover, some authors have associated the appearance of a Raman vibration band at 900-1000 cm -1 after silylation with Si-O-P vibrations [START_REF] Vasiliev | Chemisorption of silane compounds on hydroxyapatites of various morphologies[END_REF], but these bands were already observed in the spectra of the blank samples and thus can be attributed to vibrations of organic solvent. In addition, the decrease of the HPO4 2- FTIR vibration band at 875 cm -1 observed after APTES grafting was also observed in the spectra of the blanks, correlated with the decrease of CO3 2-species according to the following Eq. 1:

2 HPO 2-þ CO 2-! 2 PO 3-þ CO2 þ H2O: ð1Þ
It can be noticed that an ion-exchange process cannot be proposed to explain APTES grafting on the apatitic solids as the experiments were performed in toluene without aqueous medium. In order to favor APTES covalent grafting, the use of anhydrous solvents and pre-drying of the apatitic samples are necessary. However, even under these conditions it is impossible to avoid silane polymerisation completely on solid surfaces. With the aim of their use in heterogeneous catalysis, further experiments have to be performed on the solids after grafting to study the efficiency of such system and the stability of the layer obtained in catalytical conditions. Thus, APTES covalent grafting though the hydroxyl groups of surface HPO4 groups cannot be confirmed.

On the other hand, a clear positive marker in favor of silylation is the decrease of OH -apatitic species. The decrease of hydroxyl species appears higher in HAPD after drying than HAP, presenting the lowest amount of water and thus the weakest possibility of polycondensation. Some Conclusion APTES grafting was performed in toluene on different apatitic supports. Different supports were prepared: crystallised stoichiometric hydroxyapatites differing by the drying method, freeze-dried (HAP) and dried at 100 °C (HAPD), and an NCA. All materials were fully characterised, before and after grafting, for better understanding of the nature of the alkoxysilane/surface interaction.

The data show a clear competition between the covalent grafting of APTES and its polycondensation reaction, depending on the nature of the solid support surface, in spite of the use of anhydrous organic solvent to avoid surface polymerisation. Generally, polycondensation is favored by a more hydrated surface, whereas covalent bonding necessitates dry conditions. In this paper, although it appears that APTES polycondensation occurs whatever the nature of the apatitic sample, the spectroscopic results evidence at the same time the presence of covalent grafting. Alkoxysilane covalent grafting is always reported to take place with solid samples through surface OH groups (Solid-OH). In the current study, the IR vibration bands attributed to Si-O-P groups could not be observed, so APTES covalent grafting through the hydroxyl groups of surface HPO4 groups cannot be confirmed. However, a clear and positive marker in favor of covalent grafting is the decrease of OH -apatitic species, highest on the apatitic support presenting the lowest amount of water (HAPD) and thus the weakest possibility of polycondensation. The results obtained lead to the conclusion that silylation can be explained at the same time by APTES polycondensation and by APTES covalent grafting to oxygen atom of the hydroxyl groups of the apatitic structure and/or of the OH - species that are present on the surface hydrated layer.

This work clarifies the nature of silane grafting onto selected apatitic surfaces and especially the influence of the composition and properties of the apatitic surfaces on the process of silylation. In order to favor covalent grafting of alkoxysilanes, the use of anhydrous solvents and pre-drying of the apatitic samples are necessary. However, even under these conditions it is impossible to avoid silane polymerisation completely on the solid surfaces. The understanding of the nature of silane-surface interactions provides fundamental tools for the development of novel solid supports for heterogenous catalytic applications.

Fig. 1 Fig. 2 X

 12 Fig. 1 TEM observations of HAP and NCA before and after grafting

Fig. 4

 4 Fig.[START_REF] Grossin | Biomimetic apatite sintered at very low temperature by spark plasma sintering: physico-chemistry and microstructure aspects[END_REF] Raman spectra of a HAP, HAP0, HAP-Si, and HAP-Si-poly and b HAP-Si, HAPD-Si, and NCA-Si

Fig. 5

 5 Fig. 5 CP-MAS 31 P (a) and MAS 1 H (b) solid NMR spectra of HAP, HAP-Si, HAPD-Si, and 13 NCA-Si. c CP-MAS C solid NMR spectra of HAP-Si, HAPD-Si, and HAP-Si-poly. d CP-MAS 29 Si of solid NMR spectra of HAPD-Si, NCA-Si, and HAP-Si-poly

Fig. 6 FTIRFig. 7

 67 Fig. 6 FTIR spectra of a blank experiments HAP0, HAPD0, and NCA0 and b samples after silylation HAP-Si, HAPD-Si, and NCA-Si

Table 1

 1

	Chemical composition
	and apparent crystallite sizes of
	HAP, HAPD, and NCA before
	and after silane grafting and
	blank experiments HAP0,
	HAPD0, and NCA0


Mathieu Brunot
email: mathieu.brunot@onera.fr

Janot Alexandre
email: alexandre.janot@onera.fr

Francisco Carrillo
email: francisco.carrillo@enit.fr

Maxime Gautier
email: maxime.gautier@irccyn.ec-nantes.fr

A Separable Prediction Error Method for Robot Identification

Keywords: Robots identification, System identification, Closed-loop identification, Predictions error methods, Output error identification Robots identification, System identification, Closed-loop identification, Predictions error

The Prediction Error Method, developed in the field of system identification, handles the identification of discrete time noise model for systems linear with respect to the states and the parameters. However, robots are represented by continuous time models, which are not linear with respect to the states. In this article, we consider the issue of robot identification, taking into account the physical parameters as well as the noise model in order to improve the accuracy of the estimates. Thus, we developed a new technique to tackle this problem. The experimental results tend to show a real improvement in the estimation accuracy.

INTRODUCTION

The usual method for robot identification is based on the Least-Squares (LS) technique and the Inverse Dynamic Identification Model (IDIM). The IDIM indeed allows expressing the input torque as a linear function of the physical parameters thanks to the modified Denavit and Hartenberg (DHM) notation. Therefore, the IDIM-LS method is a really practical solution, which explains its success, see [START_REF] Gautier | A new closedloop output error method for parameter identification of robot dynamics[END_REF] and the references given therein. However this method needs a well-tuned band pass filtering to get the derivatives of the joint positions. Recently, [START_REF] Janot | A generic instrumental variable approach for industrial robot identification[END_REF]) have introduced the Instrumental Variable (IV) method which does not require such careful filtering, see an application in [START_REF] Brunot | Physical parameter identification of a one-degree-of-freedom electromechanical system operating in closed loop[END_REF]. This method is called as IDIM-IV method. The IDIM-LS and IDIM-IV methods belong to the linear regression framework.

Another technique is the Prediction Error Method (PEM). It was originally developed in the automatic control field in the discrete time framework; see e.g. [START_REF] Ljung | System Identification -Theory for the User[END_REF]. Although this method is of great interest, it cannot be straightforwardly applied to robots, because their dynamic models are continuous time. In other words, it is about fitting Ordinary Differential Equations (ODE) coefficients, see e.g. (Schittkowski 2013) or [START_REF] References Baysse | Least squares and output error identification algorithms for continuous time systems with unknown time delay operating in open or closed loop[END_REF]. In robot identification, the Direct and Inverse Dynamic Identification Model (DIDIM) technique has been recently developed in [START_REF] Gautier | A new closedloop output error method for parameter identification of robot dynamics[END_REF]. It is an Output Error Method (OEM) i.e. a PEM where the noise model filter is equal to one (the noise is assumed to be white and serially uncorrelated). Moreover, DIDIM uses the fact that the torque is linear with respect the parameters. Unlike the standard PEM method, the output of the continuous time model is simulated and not predicted; see for instance exercise 61 of [START_REF] Schoukens | Mastering system identification in 100 exercises[END_REF].

The aim of this article is to consider the possibility of estimating the parameters of the noise filter which colours the noise, as well as the physical parameters of a robot. It has been already shown that considering PEM instead of OEM allows a better precision in the estimates (smaller covariance), see Chapter 7 of [START_REF] Söderström | System identification[END_REF]. If the robot model is a nonlinear continuous time system, the noise filter is considered here as a linear discrete filter. In fact, it is more convenient for the identification process as explained in [START_REF] Garnier | Direct identification of continuous-time models from sampled data: Issues, basic solutions and relevance[END_REF]. Two solutions to estimate this filter appear. Firstly, as with the regular PEM, simultaneously estimated thanks to a nonlinear optimisation algorithm. Secondly, as with the Refined Instrumental Variable (see [START_REF] Young | Refined instrumental variable estimation: Maximum likelihood optimization of a unified Box-Jenkins model[END_REF]), the physical parameters are estimated thanks to usual techniques, and then the filter is separately identified from the residuals. Both solutions are investigated in this article. They are evaluated trough experiments. The experimental results tend show that the separable approach seems more suitable to provide optimal estimates of the robot parameters, although the computational burden is slightly higher.

The rest of the paper is organized as follows. Section 2 reviews the usual robots identification technique IDIM-LS. Section 3 presents the PEM as developed in the field of system identification, the OEM used for robot identification, and the proposed alternative: a separable PEM. The experimental identification of the SCARA robot is presented in Section 4. Finally, section 5 concludes this paper.

IDIM: INVERSE DYNAMIC IDENTIFICATION MODEL TECHNIQUE

Notations

The main notations used in this article are illustrated in Fig. 1 where 1 z is the unit delay operator, i.e. 1 1 k k z q t q t .

The noise filter, H , will always be considered as a Linear Time Invariant (LTI) system in the discrete time framework. The controller, C , and the plant, G , can be either continuous time or discrete time transfer functions. The input noises are considered as white and are contained in e . The vector v contains the measurement coloured noises.

SCARA Prototype Model

For ease of understanding, the model of the SCARA is directly presented. To get a comprehensive picture of the robots modelling, please refer to [START_REF] Khalil | Modeling, identification and control of robots[END_REF]. SCARA (see Fig. 2) is a prototype which has previously been identified with DIDIM and IDIM-LS techniques in [START_REF] Gautier | A new closedloop output error method for parameter identification of robot dynamics[END_REF].

The IDIM is defined by (3), where j q , j q j q and j q j q are respectively the angular position, velocity and acceleration of joint j . The IDIM is used for the LS identification. In the opposite, the simulation necessary for the PEM methods requires the Direct Dynamic Model (DDM), represented by the operator G in Fig. 1 and given by (2). For further information about the inertia matrix A , and vector d (containing Coriolis and friction components), please refer to [START_REF] Gautier | A new closedloop output error method for parameter identification of robot dynamics[END_REF] and the references given therein. The vector in (1) contains the physical parameters to be estimated, where

1 1 2 2 2 2 2 T r ZZ Fv Fc ZZ MX MY Fv Fc (1) ( ) , t t t t t -1 q A q ( ) t t t t t t t -1 q A q A q A q A q A q ( ) t t
( ) 0 r ZZ q t D ( ) q t q t q t 1 ( ) ( ) 1 ( ) q t 1 ( ) ( ) q t 1 ( ) , 1 1 q ( ) 0 Fv t D q ( q ( q ( q ( q ( ) 1 q ( ) q ( ) 1 q ( ) 1 q ( ) 1 , 2 1 2 1 2 ( ) ( ) ( ) ( ) ZZ q t q t q t q t D ( ) ( ) q t q t q t 1 2 ( ) ( ) ) ( ) 1 2 1 2 1 2 1 2 1 2 1 2 ( ) q t q t 1 2 1 2 ( ) 1 2 1 2 1 2 ( ) q t 1 2 ( ) 1 2 1 2 1 ( ) ( ) 1 2 ( ) q t q t 1 2 1 2 ( ) ( ) ( ) ( ) 1 2 1 2 1 2 ( ) ( ) q t q t q t ( ) ( ) ( ) ( ) 1 2 1 2 , 1 1 q ( ) 0 Fc sign t D
q ( q ( q ( q ( q ( ) q ( q ( q ( q ( q ( 1 1 )

1 q ( 1 ) q ( q ( 1 , 2 2 0 q ( ) Fv t D q ( q ( q ( ) 2 q ( ) 2 q ( q ( q ( ) , 2 2 0 q ( ) Fc sign t D q ( q ( q ( q ( ) q ( q ( q ( q ( 2 ) 2 q ( ( 2 ) q ( q ( q ( , 2 1 2 2 2 1 2 2 2 1 2 1 2 2 ( ) ( ) cos ( ) ( ) 2 ( ) ( ) sin ( ) ( )cos ( ) ( )sin ( )
MX q t q t q t q t q t q t q t q t q t q t q t

D 2 ( 2 ( 2 ( 2 ( ) ( ) ( ) ( ) ( ) ( ) cos ( ) ( ) 2 ( 2 ( 2 ( 2 ( 2 ( ) ( ) ( ) ( ) ( ) ( ) sin 2 ( 2 ( 2 ( ) ( ) ( ) ( ) cos ( ) ( ) 2 ( 2 ( 2 ( ) ( ) ( ) ( ) sin 2 ( ( ) cos ( ) ( ) ( ) sin 1 2 2 2 1 2 2 1 2 2 2 1 2 1 2 2 2 1 2 ) ( ) cos ( ) ( ) 2 ( ) ( ) ) ( ) ( 2 ( ) ( ) cos ( ) ( ) 2 ( ) ( 1 2 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 2 1 2 2 q t 1 2 2 2 1 2 2 1 2 2 2 1 2 2 ( 2 ( 1 2 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 2 1 2 2 2 ( 2 1 2 2 1 2 2 1 2 2 1 2 2 1 2 2 2 1 2 2 ( ) s ( ) ( ) ( ) ( ) ( ) ( )sin ( ) 2 1 2 1 2 q t q t q t q t 1 2 1 2 1 2 1 2 ( )cos ( ) ( )sin ( ) )cos ( ) ( )sin ( ) 2 1 2 1 2 2 1 2 1 1 2 ( ) s ( ) ( ) ( ) ( ) ( ) ( )sin ( ) 2 q t q t q t ( )cos ( ) ( ) ( ) ( )sin ( ) ( )cos ( )sin ( ) 2 ) ( ) ( 1 2 1 2 1 2 1 2 ) ( ) ( and 2 1 2 2 2 1 2 2 2 1 2 1 2 2 ( ) ( ) sin ( ) ( ) 2 ( ) ( ) cos ( ) ( )cos ( ) ( )sin ( )
MY q t q t q t q t q t q t q t q t q t q t q t

D 2 ( 2 ( 2 ( 2 ( 2 ( ) ( ) ( ) ( ) ( ) ( ) sin ( ) ( ) 2 ( 2 ( 2 ( 2 ( 2 ( ) ( ) ( ) ( ) ( ) ( ) c 2 ( 2 ( 2 ( ) ( ) ( ) ( ) sin ( ) ( ) 2 ( 2 ( 2 ( ) ( ) ( ) ( ) cos 2 ( ) sin ( ) ( ) ) ( ) c 1 2 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 2 1 2 2 ( ) ( ) sin ( ) ( ) 2 ( ) ( 2 ) ( ( ) ( ) ( ( ) ( ) ( 1 2 2 2 1 2 2 1 2 1 2 1 2 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 2 1 2 2 q t 1 2 2 2 1 2 2 1 2 2 2 1 2 2 ( 2 ( 1 2 2 2 1 2 2 1 2 2 2 1 2 2 1 2 2 2 1 2 2 2 ( ( 1 2 2 2 1 2 2 1 2 2 2 2 2 2 1 2 2 2 ( ) s ( ) ( ) ( ) ( ) ( ) ( )sin ( ) 1 2 1 2 ( ) q t q t q t q t 2 1 2 1 2 1 2 1 2 ( )cos ( ) ( )sin ( ) ( )cos ( ) ( )sin ( ) 1 2 1 2 1 2 1 1 2 1 2 2 ( )cos ( ) ( ) ( ) ( ) ( ) ( )sin ( ) q t q t 2 ( )cos ( ) ( ) ( ) ( )sin ( ) ( )cos ( )sin ) ) ( ) (
Fig. 2. SCARA Prototype

Least-Squares

The model described by ( 3) can be straightforwardly extended to the following vector-matrix form

LS u (4)
where, u is the (2 1) N sampled vector of ( ) t , is the From (4), the LS estimates and their associated covariance matrix are given by

1 T T LS , 1 T LS LS (5) with, 2 1 2 1 k LS k p k N LS t t N n .
From a theoretical point of view, the LS estimates ( 5) are unbiased if the error has a zero mean and if the regressors are uncorrelated with the error, see [START_REF] References Baysse | Least squares and output error identification algorithms for continuous time systems with unknown time delay operating in open or closed loop[END_REF].

( ) 0 LS E e LS LS ( ) ( ) ( 0 ) T T E E E ( 6
)
The covariance matrix given by (5) assumes that is deterministic and LS e is homoscedastic i.e. var( )

LS LS e . It is usually assumed that those two assumptions hold.

However, the systems considered in this article operate in closed-loop. In that case, the assumption given by ( 6) does not hold (Van den Hof 1998). This explains why practitioners usually use tailor-made pre-filtering prior to the identification process. Various filtering approaches exist in the literature, as for instance state variable filters, Poisson filters or implicit sensitivity filters. However, they only work for systems linear in the states, see for example [START_REF] Young | An optimal instrumental variable approach for identifying hybrid continuous-time Box-Jenkins models[END_REF]. For robot identification, the usual filtering process is described in [START_REF] Gautier | Dynamic identification of robots with power model[END_REF].

PEM: PREDICTION ERROR METHODS

3.1.

Regular PEM

The original PEM has been developed to identify open-loop system, in a LTI and discrete time framework. The aim is to find the vector of parameters which minimizes the error between the measured output of the system and the predicted one. Mathematically, we must solve (7), where is the prediction error vector. This vector depends on the utilised technique. Different errors will be presented in this article. The problem is usually nonlinear and requires appropriate optimisation algorithm, like Levenberg-Marquardt one. As it has been shown in [START_REF] Forssell | Closed-loop identification revisited[END_REF], PEM can be extended to the closed-loop case (i.e. direct approach), but then the physical parameters and the noise model are linked.

That is to say that if an error was done in the noise model identification, the estimated physical parameters would be biased, see e.g. [START_REF] Forssell | Closed-loop identification revisited[END_REF]. This is due to the correlation between the noise, v , and the input, . In this case, the prediction error is given by ( 8), , with 1; 1 j k .

OEM for Robots

The standard PEM cannot be directly applied for robot identification. The robot models are indeed continuous time and nonlinear with respect to the states. This problem is similar to the one of aircraft identification; see e.g. [START_REF] Klein | Aircraft system identification: theory and practice[END_REF]). Both communities have developed a similar solution. Firstly, the filter that coloured the noise is usually neglected ( H I , with I the identity matrix). According to the classical systems structures, see Chapter 4 of [START_REF] Ljung | System Identification -Theory for the User[END_REF], the robots and the aircrafts are then modelled by output error models. Hence, the applied identification techniques are called OEM. Secondly, the model output is usually simulated and not predicted. Thus, the current output is function of the initial states, the past and the current inputs, but not of the past measured outputs.

As it has been seen previously, the IDIM are linear with respect to the physical parameters. We therefore consider the input torques. From this point comes the specificity of OEM in robot identification. In fact, the simulated output of the model is instead of q . Thus, the output error vector is given by

. The vectors s q , s q s q and s q s q contain respectively the angular positions, velocities and accelerations of robot joints, coming from the simulation of the whole nonlinear closed-loop system. Hence, the knowledge of the controller is required for the simulation. Since the input of the simulation, ref q , is perfectly known (i.e. noise free), is not correlated with the measurement noise v . That insures the consistency of the estimation, assuming no modelling error. Furthermore, from (3), it comes:

| T oem k k k t t t . ( 10
)
If the dependence of in is neglected, the optimisation of ( 7) is greatly enhanced. In fact, the gradient defined by (A.3), in Appendix A, is just ,

T k k t t .
In the field of robot identification, we call it DIDIM and it is iteratively solved with LS, see Eq. ( 37) in [START_REF] Gautier | A new closedloop output error method for parameter identification of robot dynamics[END_REF].

In the field of system identification, this technique is called Pseudo-Linear Regression (PLR), see Eq. ( 7.112) in [START_REF] Ljung | System Identification -Theory for the User[END_REF]. According to the same reference, PLR is derived from [START_REF] Solo | Time Series Recursions and Stochastic Approximation[END_REF].

Separable PEM for Robots

The advantage of the PEM compared with the OEM is that it provides lower covariances for the estimated parameters; e.g. Complement C7.5 in [START_REF] Söderström | System identification[END_REF]. To achieve this in robot identification, a noise dynamics H could be added in the error term (10). Nevertheless, we have no prior information about the order of the filter and this solution may lead to a complex optimization problem. Consequently, the optimization may reach local minima or even diverge. Therefore, inspired by the Refined Instrumental Variable (RIV) developed in [START_REF] Young | Refined instrumental variable estimation: Maximum likelihood optimization of a unified Box-Jenkins model[END_REF], we propose a separable approach by considering the following error:

1 1 | , | sep pem k oem k t z t . ( 11
)
The proposed methodology, a separable PEM (SEP-PEM), is composed by two steps:

1. Estimate the physical parameters of the system by solving ( 7) with ( 9) or (10);

2. Obtain an estimate of the noise model with (11). In other words, the residuals of the first step are used as an estimate of the noise v .

Those steps are repeated until the prediction error has converged. From the residuals of the second step, an estimate of noise variance is computed with (A.2). The separation of the identification of the physical parameters and the one of the noise model implies that both models are statistically independent.

EXPERIMENTS

Experimental Setup

To illustrate our separable method, the SCARA prototype presented in section 2.2 is studied. The identification is performed with experimental data. Joint position and control signals are recorded with a sample frequency of 100 Hz. For the IDIM-LS method, the filters are tuned according to [START_REF] Gautier | A new closedloop output error method for parameter identification of robot dynamics[END_REF]. The first step (OEM method) is performed by solving ( 7) with ( 9), using the Levenberg-Marquardt algorithm, which is initialised with CAD values. That is to say, the inertias are approximately known (3.0 for 1r ZZ and 0.5 for 2 ZZ ) whereas all the other components are set to zero. From a practical point of view, we simulate the model with Simulink® and the parameters, solutions of ( 7), are found thanks to the lsqnonlin function, available in the Optimisation Toolbox of Matlab®. It should be noticed that the gradient is approximated using finite differences. The second step is performed with the aic function of the CAPTAIN Toolbox. This function seeks the Auto-Regressive filter which provides the best Akaike Information Criterion (AIC) while ensuring the witheness of the estimated input signal e e . For more detailed information about the toolbox, please refer to [START_REF] Young | Recursive estimation and time-series analysis: an introduction[END_REF]. In Eq. ( 11), H is a 2 2 matrix of transfer functions. Since each link has its own sensor, we assume there is no correlation between the noises of both axes. That is to say that H is diagonal, To insure the differentiability of the criterion, the sign function is not used in the simulation model. It is known that if is large enough then q ( ) 2arctan q ( ) j j sign t t q ( q ( q ( q ( ) 2arctan q ( q ( q ( q ( q ( ) j j

) 2arctan q ( q ( q ( q ( q ( q ( ) 2 ) 2

j j q ( ) 2 q ( q ( ) 2 ) 2 ) 2 ) 2 ) 2 ) 2arctan j j ) 2arctan .
The parameter is a scale factor, which is taken equal to 100. This has indeed proven to be a parsimonious choice. Actually, a too large value would make the system stiff.

Identification Results

The identification results are summarized in Table 2, which presents the estimated values and the relative standard deviations, defined in Appendix A. The SEP-PEM is just written PEM in the tables and figures in order to increase the readability. Since the OEM and PEM methods have estimated the same physical parameters, their results will be referred as OEM/PEM. For information, the SEP-PEM algorithm has converged in two iterations. The settings for the LS prefiltering are the same as those in [START_REF] Gautier | A new closedloop output error method for parameter identification of robot dynamics[END_REF]. It is remarkable that LS and OEM/PEM methods provide equivalent estimates. There is a noticeable difference in the estimation of 1r ZZ . However, as it can be seen in Fig. 4 and Fig. 6, it does not lead to large difference in the signals estimation. In fact, with respect to Table 3, the relative errors in the torques estimations are close and satisfactory. The interest of the SEP-PEM method is clearly visible in Table 2 with its small relative standard deviations. Actually, for the OEM case, the relative standard deviations are calculated with (A.2) and ( 9), whereas for the SEP-PEM case they are calculated with (A.2) and (11). Fig. 5 and Fig. 7 prove the efficiency of the SEP-PEM method to whiten the residuals. For information, those figures have been drawn with the acf function of the CAPTAIN Toolbox, with default parameters, see [START_REF] Young | Recursive estimation and time-series analysis: an introduction[END_REF]. This routine computes the AutoCorrelation Function (ACF) of the identified noise. It is recalled that for a white noise, the autocorrelation is zero for any non-zero lag. In the presented figures, the autocorrelation values are indicated in blue (for 38 lags) and the 2 confidence interval appears in red. The aic function has estimated a 13 th order noise filter for the first link and 37 th order one for the second link.

For the sake of completeness, we have unsuccessfully tried to identify the noise filters at the same time as the physical parameters. The filters orders and the initial values of the coefficients were those previously estimated by the aic function.

Cross-Validation

To validate the estimate of the SEP-PEM method, crossvalidations have been performed with another experimental data set. The relative prediction errors are respectively 3.21% and 4.42% for the first and the second link. These values are equivalent to those in Table 3. Fig. 8 shows the autocorrelations of the residuals. If they seem less white than those of the identification, they are still compatible with a white noise. In this paper, a separable robot identification method has been presented, experimentally validated on the SCARA prototype robot and compared with two standard approaches. This technique is divided in two sequential steps: the identification of the physical parameters thanks to an output error method, and the identification of parameters of the filter, which colours the noise. The experiments carried out with the SCARA prototype indicate that this method seems to be suitable for robot identification. Finally, compared with other standard methods, this separable approach is more effective because it provides estimates with small variances. Future works concern the application of this separable approach to an industrial robot and the comparison with other approaches not addressed in this paper.

Fig. 1 .

 1 Fig. 1. Overall figure of the considered system

  Coulomb friction parameters of joint j .

  1

  observation matrix of T t and LS e is the (2 1) N vector of error terms, with N the number of sampled data considered and p n the number of parameters. It is assumed that is full rank, i.e. ( ) n , to have an over-determined system of equations.

  arm j . Fig.3illustrates the joint j controller, where j d is the nonlinear disturbance containing the friction and Coriolis terms. This is the model utilized for the simulation of the robot. Therefore, by taking

Fig. 3 .

 3 Fig. 3. PD control law of the SCARA robot, for arm j

Fig. 5 .

 5 Fig. 4. Identification time history Link 1

  Fig. 6. Identification time history Link 2

Fig. 8 .

 8 Fig. 8. Autocorrelation of the SEP-PEM residuals (crossvalidation): Link 1 (top) and Link 2 (bottom)

Table 1 .

 1 Although, the control law presents poor performances from a tracking point, it allows identifying the robot.

Table 1 . Numerical values of the PD control law

 1

		Link 1	Link 2
	Kp	6.25	156.25
	Kv	2.0	25.0
	Kg	3.5	0.06

Table 2 . Identification results

 2

				LS	% LS	/ OEM PEM	% OEM	% PEM
	1r ZZ	3.50	0.007	3.44	0.024	0.0024
	1 Fv	0.03	2.1	0.07	3.7	0.38
	1 Fc	0.25	0.13	0.23	0.55	0.056
	ZZ	2	0.06	0.040	0.06	0.13	0.013
	MX	2	0.24	0.013	0.24	0.056	0.0057
	2 MY	0.01	0.29	0.01	1.2	0.12
	2 Fv	0.005	1.8	0.003	8.4	0.86
	2 Fc	0.05	0.49	0.05	1.7	0.17

Table 3 . Relative estimation error 100

 3

	jmeas	j j j	jmeas jmeas	LS	OEM/PEM
	Link 1,	j	1		2.98%	2.39%
	Link 2,	j	2		4.35%	4.57%

Kg Kp q t q t Kv q t t t Kg Kp q t q t Kv q t

The goal of this appendix is not to provide the whole theoretical development of Chapter 9 in [START_REF] Ljung | System Identification -Theory for the User[END_REF]), but only to remind the main results. Assuming the optimization algorithm ha and there is no modelling error, the covariance matrix of the asymptotic distribution can be estimated from data by (A.1), where is the gradient defined by (A.3) and N is the estimated covariance of the measurement noise.

From (A.1), it can be defined the relative standard deviation of the th i parameter with (A.4), assuming non zero value.


Armine Karami

Dimitri Galayko

Philippe Basset

A Novel Characterization Method for Accurate Lumped Parameter Modeling of Electret Electrostatic Vibration Energy Harvesters

come   L'archive ouverte pluridisciplinaire

Abstract-This letter presents a new method for the characterization of electret transducers, which are typically used in electrostatic vibration energy harvesters (electret e-VEHs). This is the first method allowing to accurately measure the value of the equivalent voltage source representing the electret in lumped parameter models of a wide range of electret e-VEHs. An accurate value for this parameter is critical for design, analysis and optimization, given the increasing complexity of e-VEHs electrical interfaces. Until now, there was no universal method allowing the measurement of this parameter, because of practical difficulties with some geometries, and because of charging nonuniformities. In this letter, the new method is presented, with insights on how to maximize the measurement accuracy. It is then applied to a state-of-the art MEMS electret e-VEH.

Index Terms-electrostatic vibration energy harvesting, electret, charge pumps, energy conversion

I. INTRODUCTION

E LECTRET electrostatic vibration energy harvesters (e- VEHs) are a promising technology to power autonomous systems [START_REF] Pondrom | Compact electret energy harvester with high power output[END_REF], [START_REF] Suzuki | Electret based vibration energy harvester for sensor network[END_REF]. They are composed of a mobile mass, attached to an electret-charged transducer. The biasing of this transducer is determined by a conditioning circuit. Energy conversion occurs when the bias across the transducer results in an electrostatic force opposing the mass motion induced by an external acceleration.

In the recent years, the complexity of the available conditioning circuits has increased to optimize the harvesters performances [START_REF] Miki | Effect of nonlinear external circuit on electrostatic damping force of micro electret generator[END_REF]- [START_REF] Karami | Series-parallel charge pump conditioning circuits for electrostatic kinetic energy harvesting[END_REF]. As a result, the use of lumped parameter models of e-VEHs is becoming the sole convenient way to predict their dynamics under variable input conditions.

A convenient representation of the electret-charged transducer in lumped parameter models of e-VEHs is an equivalent constant voltage source in series with the transducer's variable capacitance [START_REF] Boisseau | Electrostatic conversion for vibration energy harvesting[END_REF]. This source will be referred to as the electret potential in the rest of the paper. Lumped-parameter models are expected to yield accurate predictions on the e-VEH dynamics, and to guide optimization choices in the conditioning circuit design. To do this, these models have to be parametrized accurately. In particular, the electret potential value has a critical effect on the models' predictions (chapter 11 of [START_REF] Basset | Electrostatic Kinetic Energy Harvesting[END_REF])

The electret potential is difficult to define precisely at the fabrication stage. Plus, it may decrease in large time scales [START_REF] Boisseau | Stable drie-patterned sio 2/si 3 n 4 electrets[END_REF]. Hence, it is necessary to measure its value after fabrication.

A. Karami Until now, this was done by measuring the device's surface potential, with a contact-less electrostatic voltmeter. Ideally, this measurement has to be done between the facing capacitor plates of the transducer. However, this is impossible to perform nondestructively on a wide range of e-VEHs because of the geometry of their transducer (e.g., interdigitated-combs [START_REF] Lu | A batch-fabricated electret-biased wideband mems vibration energy harvester with frequency-up conversion behavior powering a uhf wireless sensor node[END_REF]). More importantly, this technique gives a measure of a potential which is hard to accurately relate to the value of the electret potential. In particular, its accuracy relies on the assumption of uniform charge distribution across the surfaces of the transducer's facing capacitor plates. This hypothesis is not always verified, depending on both the employed charging techniques and the transducer geometry [START_REF] Lu | A batch-fabricated electret-biased wideband mems vibration energy harvester with frequency-up conversion behavior powering a uhf wireless sensor node[END_REF], [START_REF] Renaud | Modeling and characterization of electret based vibration energy harvesters in slot-effect configuration[END_REF].

In this letter, for the first time, a non-destructive and accurate method is proposed to measure the electret potential of a given transducer. This method is compatible with a wide range of transducer geometries and charging techniques. The method is presented in Section II, with insights on how to minimize the error on the measurement. It is then applied to a state-of-the-art MEMS electret e-VEH in Section III.

II. PRESENTATION OF THE METHOD A. Dynamics of the measurement circuits

The measurement method presented in this letter is based on the dynamics of the half-wave circuits HW1 and HW2, depicted in Fig. 1a and1b, respectively. In this subsection, the dynamics of these circuits are linked to the electret potential, denoted by E in the rest of the paper. The value of the capacitance at the transducer's terminals is denoted by C var .

Consider the transducer as a part of the circuit configurations HW1 or HW2, with C var cyclically time-varying between maximum and minimum values (C max , C min ), subsequently to a mechanical input excitation of the e-VEH (the same input for both circuit configurations). A cycle of variation of C var is defined as its variation from C max to C min to C max . The choice of the input is not constrained as long as it results in the same values C max and C min at each cycle. A harmonic input is used for simplicity. In the following, V 1 and V 2 denote the voltage across the fixed capacitor C S at C var " C max of a given variation cycle of C var , for the HW1 and HW2 circuits respectively. The labels V 1 1 and V 1 2 denote the voltage across the fixed capacitor C S at C var " C max of the following cycle, for the HW1 and HW2 circuits respectively. We also define

∆ 1 :" V 1 1 ´V1 and ∆ 2 :" V 1 2 ´V2
. The quantities p∆ 1 , V 1 q and p∆ 2 , V 2 q are measured on the timeevolution of the voltages across C S with circuits HW1 and HW2, respectively. An experimental example of measurement of

p∆ 1 , ∆ 2 , V 1 , V 2 q is shown in Fig. 2.
In the aforementioned input conditions, the electrical analysis of each circuit yields the cycle-to-cycle evolution law of the voltage across C S , which reads:

V 1 1 " C S V 1 `pC max ´Cmin qE ´pC max `Cmin qV T C min `CS (1) V 1 2 " C S V 2 `pC max ´Cmin qE ´pC max `Cmin qV T C max `CS (2)
where V T ě 0 denotes the fixed threshold voltage of the diodes, which are assumed to follow the ideal voltage-current diode law. From these two equations, it comes that the measured p∆ 1 , ∆ 2 , V 1 , V 2 q can be linked to E and C max by:

C max " C S p∆ 1 ´∆2 q `Cmin p∆ 1 `V1 q ∆ 2 `V2 , (3)
E " pC min `CS q∆ 1 `Cmin V 1 `pC max `Cmin qV T C max ´Cmin . (4)
The value of C min can be measured by a synchronous detection method [START_REF] Basset | A batch-fabricated and electret-free silicon electrostatic vibration energy harvester[END_REF]. Note that in the case of interdigitatedcombs transducer geometries, C min is fixed independently of the harmonic input acceleration. Hence, it can be measured by traditional means of measurement of fixed capacitances.

B. Measurement errors and choice of the circuit's parameters

To measure E for a given transducer using (4), it is necessary to choose C S , V 1 , V 2 , and the input acceleration. This subsection shows how these parameters have to be chosen in order to maximize the measurement accuracy.

In the following, δ V denotes the uncertainty on the measured voltage across C S due to the superposition of random error sources in a broad sense (e.g., thermal, reading amplifier, and quantization noises). Suppose that the uncertainty on V T and C min can be neglected, and that V 1 and V 2 are chosen to be approximatively equal to zero. The measurement uncertainties on E and C max due to random errors are denoted by δC max and δE, respectively. They are linked to δ V following:

δ Cmax C max " δ V η `ξ Eηpη ´1q a p1 `ξq 2 `pη `ξq 2 (5) δ E E " δ V a pη `ξq 4 `p1 `ξq 4 Epη ´1q 2 (6)
where ξ :" C S {C min and η :" C max {C min .

From [START_REF] Wei | Interface circuit with adjustable bias voltage enabling maximum power point tracking of capacitive energy harvesting devices[END_REF][START_REF] Karami | Series-parallel charge pump conditioning circuits for electrostatic kinetic energy harvesting[END_REF], it comes that the measurement uncertainty due to random error sources decreases by decreasing C S . However, decreasing C S increases the electromechanical coupling effect on the e-VEH's dynamics, thereby inducing a systematic measurement error on C max and E. This is because equations (3-4) are derived from (1-2) assuming that pC max , C min q with circuit HW1 (in (1)) are equal to pC max , C min q with circuit HW2 (in ( 2)). But these circuits do not bias the transducer in the same way, and different biasing results in different electrostatic forces generated by each of the two circuits on the e-VEH's mobile mass. This, in turn, yields different mobile mass motions with each of the two circuits. Thus, in general, pC max , C min q with circuit HW1 cannot be considered equal to pC max , C min q with circuit HW2.

From the topologies of the circuits depicted in Fig. 1, it can be seen that the circuit HW1 biases the transducer with a voltage varying from E at C var " C max to E `V 1 1 at C var " C min . The circuit HW2 biases the transducer with a voltage varying from E at C var " C min to E ´V 1 2 at C var " C max . By choosing C S large and V 1 , V 2 as small as possible, V 1 1 , V 1 2 are small (see [START_REF] Pondrom | Compact electret energy harvester with high power output[END_REF][START_REF] Suzuki | Electret based vibration energy harvester for sensor network[END_REF]). This ensures that the biasing of C var throughout its variation is approximatively the same with both circuits (close to the voltage E). Therefore, increasing C S reduces the systematic measurement error due to the differential impact of the electromechanical coupling.

The following procedure is proposed to choose C S so as to reduce the measurement uncertainty due to random error sources, whilst ensuring that the systematic error due to the electromechanical coupling remains negligible. An admissible value for δ E {E is decided a priori, and solving (6) for ξ gives the corresponding C S . Note that this supposes an a priori estimation of η and E. The former can be measured by synchronous detection [START_REF] Basset | A batch-fabricated and electret-free silicon electrostatic vibration energy harvester[END_REF]. The latter can be estimated by one or several first inaccurate measurements (e.g., using an arbitrary value for C S for the first measurement). The measurements are then carried out with this value of C S . If C max obtained from (3) is consistent with its value measured by synchronous detection, then E obtained from (4) using the same measured p∆ 1 , ∆ 2 , V 1 , V 2 q can be considered accurate. Otherwise, C S that satisfies the desired δ E {E margin results in a large systematic error on C max , and hence on E. In this case, the measurement has to be carried out again with a larger C S . The systematic error due to the electromechanical coupling will decrease, but the admissible uncertainties on E and C max due to random error sources will increase.

Finally, the measurement uncertainty can also be decreased by increasing η at fixed C S . This is done by increasing the amplitude of the mechanical harmonic input. But as η results from the dynamics of the e-VEH's mechanical part, it is harder to use the input amplitude as a control parameter on the error. Note that a technique employing the HW1 circuit has been used for electret potential measurement in [START_REF] Lu | A batch-fabricated electret-biased wideband mems vibration energy harvester with frequency-up conversion behavior powering a uhf wireless sensor node[END_REF], similar to the method presented in this paper, but using saturation voltages instead of local evolution voltages. Because of this, the method in [START_REF] Lu | A batch-fabricated electret-biased wideband mems vibration energy harvester with frequency-up conversion behavior powering a uhf wireless sensor node[END_REF] only holds if both C max and C min are fixed throughout the voltage evolution until saturation. This may not be verified because of the electromechanical coupling.

III. EXPERIMENTAL RESULTS

In this section, the method is applied to a micro-scale electrostatic transducer realized in MEMS technology, similar to the device reported in [START_REF] Lu | A batch-fabricated electret-biased wideband mems vibration energy harvester with frequency-up conversion behavior powering a uhf wireless sensor node[END_REF]. As the transducer has an interdigitated-combs geometry, C min is fixed, and C max alternates between two different values at each period of the harmonic mechanical input excitation [START_REF] Lu | A batch-fabricated electret-biased wideband mems vibration energy harvester with frequency-up conversion behavior powering a uhf wireless sensor node[END_REF]. Therefore, there exists two values of C max for each harmonic input that the device is submitted to. Two input are chosen: 1.5g and 2g amplitude, both of 150Hz frequency.

In the setup used for this letter, the device is continuously submitted to the input vibrations, and C S is first short circuited. Then, the short circuit is opened and the voltage evolution across C S starts. But the time at which the short circuit is opened cannot be exactly chosen as corresponding to an extrema of C var . To ensure that the captured evolution results from a complete cycle of variation of C var , the first cycle of voltage evolution across C S is ignored with both circuits. Hence, V 1 , V 2 are slightly above zero (see Fig. 2).

To estimate the influence of random error sources on the measurement using (5-6), E is first estimated at 20V, by a surface potential measurement done 1mm on top of the device using a contact-less electrostatic voltmeter. The four values of η are measured by synchronous detection, and δ V is measured as 25mV. Using [START_REF] Karami | Series-parallel charge pump conditioning circuits for electrostatic kinetic energy harvesting[END_REF], the theoretical uncertainty on E decreases from 16.6% with the smallest C max , to 4.2% with the largest C max , all with C S " 100pF.

These uncertainty intervals are satisfactory for this example, so the fixed capacitor is chosen as C S " 100pF. Using an impedance-meter, C min is measured as C min " 60pF. The used diodes are JPAD5, whose parasitic capacitances values are negligible compared to C S , and V T is measured as V T " 0.7V for typical values of current in our application.

The results of the measurement method are depicted in Fig. 4. Each point corresponds to pE, C max q obtained by averaging 5 measurements of p∆ 1 , ∆ 2 , V 1 , V 2 q. The error bars represent the uncertainties obtained on p∆ 1 , ∆ 2 , V 1 , V 2 q propagated to pE, C max q. An example of a single measurement of p∆ 1 , ∆ 2 , V 1 , V 2 q obtained with the present device is depicted in Fig. 2. The results in Fig. 4 show that the obtained C max are in accordance with those measured by synchronous detection. This guarantees that the systematic error due to the electromechanical coupling is negligible. As expected, the uncertainty on E decreases when the input amplitude increases (see end of Section II-B). Using the most accurate result, the electret potential is measured as E " 13.7V ˘5.1%. The uncertainty is slightly increased compared to the theoretical uncertainty obtained from ( 6), because of the lower effective value of E compared to its first estimation, and 1 , V 2 ą 0.

For further validation, a voltage source of 5V is added in series with the transducer to simulate a different value of the electret potential, and the measurement is carried out again. The expected electret potential is then of E `5 by superposition. The used input is of 2g amplitude and 150Hz frequency. The synchronous detection method gives C max " 102pF ˘3% (the lower value of C max is discarded). From (6), keeping C S " 100pF results in an error of 3.3% on E. Using this value of C S , the measurement yields a value of E " 23.1V ˘3.8% and C max " 93pF ˘3.1%. The value of C max does not agree with what is measured by synchronous detection, which means that the systematic error due to the electromechanical coupling is large. Hence, the error minimization procedure discussed at the end of section II-B is applied: the admissible random error on E is relaxed to 7.5%, yielding C S " 180pF. The new measurement gives C max " 103pF ˘4.2%, which is in accordance with the value measured by synchronous detection. The electret potential is measured as E " 19V ˘6.8%, which is the expected value after adding the 5V external voltage.

IV. CONCLUSION

This letter presented a method to accurately measure the value of the voltage source representing the electret in lumped parameter models of e-VEHs. The method was illustrated on an interdigitated-combs MEMS transducer, for which a traditional surface potential measurement across the transducer is hard to carry out and gives inaccurate results. The results were verified by comparing the maximum capacitance value, which is obtained as an intermediary parameter in the reported method, with its value measured by synchronous detection. The method was also validated by changing the charging conditions and checking the coherence of the results.

  transducer

Fig. 1 :

 1 Fig. 1: Schematics of circuit HW1 (a) and circuit (b)

Fig. 2 :

 2 Fig.2: Experimental measurement of p∆ 1 , ∆ 2 , V 1 , V 2 q, measured on the MEMS electret transducer characterized in section III and reported in[START_REF] Lu | A batch-fabricated electret-biased wideband mems vibration energy harvester with frequency-up conversion behavior powering a uhf wireless sensor node[END_REF]. (a) Measured voltage across C S with HW1 circuit. (b) Measured voltage across C S with HW2 circuit.

Fig. 3 :

 3 Fig. 3: Setup used for the measurement of the electret potential, with a top-view picture of the electret MEMS transducer. The device dimensions are 1.5cm ˆ1.6cm ˆ380µm.

Fig. 4 :

 4 Fig. 4: Measurement results for E and C max .

  A Novel Characterization Method for Accurate Lumped Parameter Modeling of Electret Electrostatic Vibration Energy Harvesters Armine Karami, Student Member, IEEE, Dimitri Galayko, Member, IEEE, and Philippe Basset