Jan Mrázek
email: mrazek@ufe.cz

Michel Potel

Jiří Buršík

Aleš Mráček

Anna Kallistová

Šárka Jonášová

Jan Boháček

Ivan Kašík

Sol-gel synthesis and crystallization kinetics of dysprosium-titanate Dy2Ti2O7 for photonic applications

Keywords: A1, Nucleation A2, Sol-gel B1, Nanomaterials B1, Rare earth compounds

Introduction

Lanthanide titanium oxides, which crystallize in the face-centered cubic (FCC) pyrochlore structure with the general formula A 2 B 2 O 7 [1], have attracted close attention of recent material research. The reason lies in the extraordinary physical properties of the pyrochlore structure. In the pyrochlore structure, the sites of A-cations are occupied by ions of rare-earth elements (RE), the sites of B-cations are occupied by titanic ions [2]. The spins of atoms at A and B sites are arranged on the vertices of corner sharing tetrahedrons. These tetrahedrons are linked together into an infinite three-dimensional network allowing intensive spin interactions within the crystal lattice. Such an arrangement brings variety of magnetic properties that depend on the RE incorporated within the pyrochlore crystal lattice [3]. Dysprosium titanate Dy 2 Ti 2 O 7 has been particularly studied for its spin-ice behavior [4] and recently it has attracted attention as a novel photo-catalyzer [5]. Such potential applications place high demands on the quality of prepared materials.

Beside various solid-state and mechano-chemical synthesis [6][7][8], which are limited by high processing temperatures above 1100 °C, a number of soft-chemical approaches to nanocrystalline Dy 2 Ti 2 O 7 were proposed. Combination of the sol-gel approach with the organic ligands, i.e. citric acid [5] or stearic acid [9] allowed to prepare highly uniform nanocrystalline Dy 2 Ti 2 O 7 at calcination temperatures lower than 800 °C. Despite the fact that the lowered calcination temperature is a great technological advantage, none of introduced soft-chemical approaches has been successfully used to prepare transparent thin films of Dy 2 Ti 2 O 7 . To the best of our knowledge pure Dy 2 Ti 2 O 7 films were prepared only by epitaxial growth assisted by pulse laser deposition [10]. Sol-gel approach was applied to prepare other pyrochlore thin films, such a Nd 2 Ti 2 O 7 [11]. However, the surface variation of prepared films was on the level of several hundred of nanometers.

In this contribution we present generic sol-gel approach to nanocrystalline Dy 2 Ti 2 O 7 powders and highly transparent thin films. Kinetic parameters of the nucleation process and the crystallization mechanism of the formation of nanocrystalline Dy 2 Ti 2 O 7 from amorphous xerogels are studied. Thermal and structural evolution of formed nanocrystalline powders is compared with the evolution of nanocrystalline thin films. Optical properties of prepared films were evaluated. Presented results provide fundamental information about the crystallization behavior of nanocrystalline Dy 2 Ti 2 O 7 and they can be used for preparation of pure nanocrystalline xerogels and highly transparent thin films with tailored structural properties which are suitable for photonic applications.

Experimental

Materials and sample preparation

Samples were prepared by the sol-gel method followed by the thermal treatment of formed gel. Typically, a total of 6 g titanium(IV)butoxide (Fluka, Purum) was dissolved in 250 ml of absolute ethanol (Sigma-Aldrich, Spectranal grade). Then a total of 7.74 g of dysprosium(III) nitrate pentahydrate (Alfa Aesar, 99.9%) was dispersed in the solution. All precursors were used as purchased without further purification. The mixture was stirred under ambient temperature to form transparent solution. Formed solution was refluxed at 86 °C for 24 h and then allowed to cool. The solution was purified by a filtration through a 0,4 mm PTFE membrane and concentrated on a rotary evaporator. Prepared sol was divided into two parts. The first part of sol was used to prepare thin films. The sol was deposed by a spincoating method on a silica substrates (Technical glass products, USA). Typically, 150 ml of sol was dropped on the substrate and rotated at 1500 RPM for 60 s. The layer was densified in a rapid thermal annealing furnace (RTA) AccuThermo AW410 (Allwin21 corporation) at 550º C for 60 s under an oxygen flow 5 l•min -1 . The heating rate was 5 °C•s -1 . Three layers were subsequently deposed in this way. Finally, prepared multilayer films were thermally treated with a heating rate of 20 °C•s -1 to the temperatures ranging from 800 °C to 1200 ° for 60 s under an oxygen flow of 10 l•min -1 .

Remaining part of the sol was dried up on a rotary evaporator. The resulting gel was pulverized and calcined to temperatures between 800 °C and 1200 °C for 30 minutes forming xerogels. The calcinations were carried out under a 10 l×min -1 flow of oxygen.

Characterization methods

Differential thermal analyses (DTA) and thermal gravimetry (TG) measurements were performed on a DTG 60H (Shimadzu) multifunctional apparatus. The samples were analyzed at different heating rates from 5 to 40 °C•min -1 under an oxygen flow 50 ml•min -1 .

The FT-IR spectra of the xerogels were recorded with a Vector 22 (Bruker) single channel spectrometer. Approximately 4 mg of the analyzed xerogel was mixed with 200 mg of KBr (Sigma-Aldrich, 98% FT-IR grade) finely pulverized in a mortar and pressed into a pellet. A pellet made of pure KBr was used as a reference sample.

The X-ray diffraction analyses (XRD) were realized on a Bruker D8 Discover diffractometer with the Bragg-Brentano geometry operating with the Cu-Ka radiation (l=1.54056 Å, operating voltage 40 kV, current 40 mA, integration time 200 ms). The size of the formed nanocrystals, D, was calculated using Scherrer's equation -Eq. (1) [12]:

, ( 1
)
where k is a structural factor equal to 0.89, l is the X-ray wavelength, here 0.154056 nm, is the correlated FWHM value corresponding to the (222)-plane peak of Dy 2 Ti 2 O 7 at the diffraction angle 2 equal to 30.697 degrees. The crystallized phase was compared to JCPDS data file number 17-0453 of Dy 2 Ti 2 O 7 .

Transmission electron microscope (TEM) analyses of the xerogels were carried out on a Philips CM12 STEM electron microscope. A copper grid coated with a porous carbon support film was used as a sample holder. The analyzed xerogel was dispersed in ethanol, and the resulting suspension was treated in an ultrasonic bath for 10 minutes and then applied to the grid. Scanning electron microscope (SEM) images were taken with a TESCAN Lyra 3 XMU FEG/SEM device. The composition of individual grains and the overall composition of the xerogels heat-treated at 1200 °C were measured on a Cameca SX 100 electron microprobe. The samples were excited by the X-ray radiation beam of the energy equal to 15 keV. Beam current was 60 mA. A thin carbon layer was sputtered on the samples prior to analysis to prevent charging.

Raman spectra were recorded on a Renishaw Invia apparatus. The in-build green laser emitting at 633 nm with the output power of 50 mW was employed as the excitation source. The ultra steep long-pass edge filter allowed collecting the Raman signal with the Raman shift higher than 150 cm -1 .

UV-VIS transmission spectra were recorded on a Lambda EZ 210; Perkin-Elmer; spectrometer.

The surface morphology of thin films was analyzed by an atomic force microscopy (AFM); Digital Instruments/Veeco, Woodbury, NY; in the contact mode using a Dimension 3100 AFM operated in air. Commercial Si 3 N 4 cantilevers (DI) with an elastic modulus of 0.56 N×m -1 were used.

Results and discussion

Thermal behavior of xerogels

Thermal analysis of xerogels

Calcination of xerogels is a complex process consisting of several steps. It usually involves evaporation of volatile substances, burning of organic compounds, decomposition of inorganic ligands, and the crystallization process. The results of thermal analyses are depicted in Fig. 1a. The DTA curves showed multiple endothermic peaks below 300 ºC, a dominant exothermic peak at around 380 ºC and an exothermic peak at 841 ºC. The multiple endothermic peaks below 380 ºC were accompanied by a regular weight loss. Two smaller weight losses were related to exothermic peaks at around 380 ºC and at 841 ºC.

FT-IR spectra of calcined xerogels were recorded to clarify the origin of DTA peaks. The results are depicted in Fig. 2. Thermally unprocessed xerogel exhibited several absorption bands in FT-IR spectrum. A dominant absorption band, located from 2800 to 3600 cm -1 , corresponded to the vibrational bands of -OH groups of ethanol and alkoxy groups. Multiple absorption bands located from 1300 to 1650 cm -1 corresponded to the vibrational spectra of nitrates presented in xerogel [13]. These intensive bands overlapped the absorption band of unreacted -OH groups located at around 1638 cm -1 [14]. Since the broad peak ranging from 2800 to 3600 cm -1 had not been observed in the xerogel calcined at 300 ºC it should be concluded that the endothermic peak located on DTA curve below 300 ºC is related to the evaporation and oxidation of remaining ethanol and butanol which had been released during the hydrolysis of titanium(IV)butoxide. The intensity of the absorption bands of nitrates was strongly decreased in the sample calcined above 300 ºC. Thus the first exothermic peak located at 380 ºC should be related to the burn-out of nitrates. However, a part of nitrates remained in the xerogel up to 840 ºC. These nitrates act as ligands stabilizing the dysprosium ions inside amorphous matrix as in the case of crystallization of Eu 2 Ti 2 O 7 from amorphous xerogel [15]. Once the remaining nitrates had been burn-out the crystallization of the Dy 2 Ti 2 O 7 occurred and the crystallization peak appeared on DTA curve at 844 ºC.

The DTA analysis is a dynamic process. The position of the crystallization peak T p on the thermal curve depends on the heating rate a. The shift of the crystallization peak is widely used to calculate many thermodynamic parameters. The crystallization temperature can be calculated from the Lasocka equation -Eq. ( 2) that can be expressed in the following form [16]:

, ( 2
)
where T c is the crystallization peak position for the heating rate of 1 K×min -1 , A c is an experimental constant. The parameters are accessible via the linear fit of the dependence of the position of the crystallization peak T p on the heating rate a. Experimental data including the linear fit are depicted in Fig. 3. Calculated value of the crystallization temperature was equal to 798±2 °C. This value is fully comparable with the crystallization temperatures of other RE pyrochlores prepared by the solgel method, e.g. 790 ºC for Nd 2 Ti 2 O 7 [11], and significantly lower than the temperatures that are necessary to the formation of RE pyrochlores by common solid state reactions, e.g. 1100 ºC for Er 2 Ti 2 O 7 , Nd 2 Ti 2 O 7 [17].

Nucleation process and crystal growth

The crystallization mechanism is a key factor which affects the structural properties of formed nanocrystals. The type of the nucleation affects the homogeneity of the structure, the shape and the size distribution of formed nanocrystals. From the viewpoint of the crystallization kinetics two types of nucleation mechanism can be distinguished; site-saturated nucleation and nucleation with constant rate of nucleation. Site-saturated nucleation is characterized by the formation of the crystallization nuclei at the same time in the whole volume of the amorphous material. The growth of formed nuclei into larger nanocrystals occurs after the nucleation. Site-saturated nucleation leads to the materials with highly uniform shape and narrow size distribution of formed nanocrystals. Nucleation with constant rate of nucleation is characterized by continuous formation of nucleus simultaneously with the growth of formed nanocrystals. Formed materials exhibit lower uniformity and large size distribution of formed nanocrystals comparing with the site-saturated nucleation [START_REF] Christian | The Theory of Transformations in Metals and Alloys (Part I + II)[END_REF].

The type of the nucleation mechanism can be evaluated from the shape of the crystallization peak. The theoretical basis for the description of nucleation-growth processes under non-isothermal conditions have been formulated by Johnson, Mehl and Avrami. The proposed models have been summarized by Christian [START_REF] Christian | The Theory of Transformations in Metals and Alloys (Part I + II)[END_REF]. The type of the nucleation mechanism is determined by the value of the Avrami parameter n, which can be calculated from a Johnson-Mehl-Avrami equation -Eq. ( 3) expressed in the Ozawa modification for linear heating rates [START_REF] Malek | [END_REF]:

, ( 3
)
where Avrami parameter n is determined from the slope of the plot log[-ln(1-)] vs. log a. is the volume fraction of a crystallized material, and a is the heating rate. The value of is calculated as the ratio of the partial integration area of the crystallization peak at temperature T over the total area of the crystallization peak. The Avrami parameter n can be assigned to the characteristic values which correspond to the type of nucleation and growth mechanisms. These values, together with corresponding nucleation and growth processes, are summarized in Table 1. Each grade of the Avrami parameter represent one dimension of crystal growth or a formation of nuclei. Therefore the values of the Avrami parameter for nucleation with constant rate of nucleation are higher than the values for the site-saturated nucleation. In spite of proposed theoretical model the real crystallization usually combines several processes, so the Avrami's parameter varies between predicted values. Calculated Avrami parameters are depicted in Fig. 4. Despite relatively large error the high value of the Avrami parameter close to 4 indicates that the crystallization started by the homogeneous nucleation with constant nucleation rate. The nucleation was accompanied by threedimensional growth of nanocrystals. The transport of the matter to the formed nanocrystals was limited by the mass transport through the phase-boundary between the amorphous and crystallized part of the xerogel. Once the formed nanocrystals had made the network inside the amorphous xerogel, the mass transfer has been limited, and thus diffusion controlled mass transfer has become the main transfer mechanism. The change of the reaction mechanism reduced the values of the Avrami parameter below 2.5. Finally, once all the amorphous part of the xerogel was fully crystallized, the formation of further nuclei was terminated, cellular recrystallization occurred, and further growth of nanocrystals has been controlled by the diffusion of matter. Consequently the values of the Avrami parameter decreased below 1.

The proposed nucleation mechanism was validated by the SEM analyses of prepared xerogels. Selected results are depicted in Fig. 5. The xerogel calcined at 1200 ºC was fully crystalline and completely filled by formed nanocrystals. However, the uniformity and size distribution of the nanocrystals were highly variable. The size of formed nanocrystals ranged between 40 and 200 nm. The xerogel was composed by large nanocrystals together with smaller one. Such a structural properties are typical for nanocrystalline xerogels whose crystallization mechanism proceeds in the manner of homogeneous nucleation with constant rate of nucleation. Similar structural properties has been already demonstrated on zinc titanates [20].

To evaluate the activation energy of crystallization processes, three approximations represented by Ozawa equation -Eq. ( 4) [21], Kissinger [22] equation -Eq. ( 5), and Augis-Bennett equation -Eq. ( 6) [23], respectively, were considered:

, ( 4)
, ( 5)

. ( 6
)
The symbol means the heating rate, T p is the position of the crystallization peak, T 0 is equal to 273.15 K and represents the conversion to Celsius scale, E A is the activation energy of the crystallization, A is the pre-exponential factor, K 0 is the frequency factor, R is the ideal gas constant and c is the additive constant. Calculated values of the terms on the left side of the equations plotted against the reciprocal absolute temperature for selected approximations and their corresponding linear fits are depicted in Fig. 6. The calculated values of the E a are summarized in Table 2. All three approximations provided similar values of the activation energy close to 780 kJ×mol -1 . This value is comparable with the values determined for other titanates prepared by sol-gel or hydrothermal methods. For example, the published value for europium titanate is 707 kJ×mol -1 [15], and values for barium strontium titanates range from 400 to 430 kJ×mol -1 [24]. The differences between these previously measure values and the values shown in Fig. 6a and Fig. 6b can be attributed to the different structural and chemical composition of each material.

Structural properties of nanocrystals

Thermal evolution of nanocrystals

XRD patterns of calcined samples are depicted in Fig. 7. In agreement with the calculated crystallization temperature the first diffraction pattern was observed in the xerogel calcined at 800 °C. The diffraction patterns corresponded to dysprosium titanate Dy 2 Ti 2 O 7 pyrochlore structure. The samples were formed by pure nanocrystals without any traces of sideformed oxides. The electron microprobe analysis confirmed the chemical composition of formed nanocrystals corresponding to the formula Dy 1.97 Ti 2.03 O 7.00 . The deviation from the theoretical composition was within the instrumental error of the measurement. Comparing with the xerogels, the thin films diffraction patterns exhibited a lower signal-to-noise ratio. Such an observation is natural due to the limited thickness of the analyzed films thus the lower intensity of the recorded signal. The mean sizes of formed nanocrystals calculated from the Scherrer's equation are depicted in Fig. 8. The calculated mean size of formed nanocrystals regularly grew up with increasing calcination temperature from 35 nm up to 120 nm. The nanocrystals constituting the xerogels were slightly bigger than the nanocrystals constituting the films. The variation in sizes between the xerogels and the films should be attributed to the different heating mechanisms which were used for sample processing. The difference is obvious for temperatures above 1000 °C. The activation energy of nanocrystal growth was calculated from following equation -Eq. ( 7) [24]: , (7) where D mean the average nanocrystal size, E A-G is the activation energy of the nanocrystal growth, R is the ideal gas constant, A is the constant. The values of the term -E A-G /(2.303•R) are accessible by fitting the plot of log(D) versus the reciprocal absolute temperature (1/T) as can be seen in Fig. 9. In spite of the different size of nanocrystals constituting the xerogels and the films, the trends of both fits were very similar. So the calculated values of the activation energy of nanocrystal growth were equal to 24.4±1.4 kJ×mol -1 and 23.5±2.4 kJ×mol -1 for xerogels and thin films, respectively. These values well matched the values 12.3 kJ×mol -1 and 26.5 kJ×mol -1 calculated for nanocrystalline Nd 2 Ti 2 O 7 [11] and Eu 2 Ti 2 O 7 [15], respectively, which were prepared by the sol-gel method.

The XRD analyses of xerogels were supported by the TEM observations which are depicted in Fig. 10. The xerogel calcined at 800 °C was composed of fine nanocrystals with the size ranging from 10 nm up to 40 nm around the mean size calculated from the Scherrer s equation. The nanoparticles grew up with increasing calcination temperature and their size was increased. Relatively wide range of observed nanocrystal sizes around the calculated nanocrystal sizes supports the proposed nucleation mechanism. The corresponding ring diffraction patterns proved that the nanocrystals were randomly oriented in the xerogel. The diffraction patterns measured on particular nanocrystals confirmed that the nanocrystals crystallized in the FCC lattice with the lattice parameter a equal to 10.2 Å. Such a value is in good accordance with the value a equal to 10.106 Å specified in the JCPDS data file.

Raman analyses of nanocrystals

Structural data obtained from XRD analysis were supported by the Raman spectroscopy to exclude the presence of partially crystallized side-formed phases. The structure of cubic pyrochlore can be expressed by a general formula A 2 B 2 O(1) 6 O(2) with two distinct oxygen sublattices [1]. This structure belongs to the space group Fd3m (n227) with Z = 8. The site symmetry is D 3d for REE and Ti 4+ ions, C 2v for O(1) ions and T d for O(2) ions. According to the group theory, the pyrochlore structure displays six Raman-active modes distributed among three irreducible representations -Eq. ( 8) [1]:

.

Recorded Raman spectra are depicted in Fig. 11. The general shape of recorded spectra and the peak positions corresponded to the Raman spectra recorded on polycrystalline Dy 2 Ti 2 O 7 that was prepared by solid-state reaction [25]. The most intensive peak located around 308 cm -1 was assigned to the F 2g bending modes, the peak located around 518 cm -1 was assigned to the A g stretching mode and the weakest peak around 699 cm -1 corresponded to 2 nd order scattering effects. E g mode was not observed. The Raman spectrum of the thin film was modulated by a sinusoid which had an origin in interferences of the excitation laser beam on the sample. The absence of the peaks corresponding to side-formed TiO 2 or Dy 2 O 3 fully proved the high purity of prepared Dy 2 Ti 2 O 7 .

Properties of thin films

The key advantage of the sol-gel method lies in its flexibility that allows to prepare both powders and thin films. Proposed approach provided crack-free, homogenous, transparent thin films as was demonstrated on films calcined at 1200 °C. The AFM surface profile of the film is depicted in Fig. 12. The film surface was composed of regularly distributed nanocrystals. Lateral dimensions of the nanocrystals ranged between 40 and 200 nm in accordance withe the SEM measurements depicted in Fig. 5. Maximal height difference of the surface was 16.7 nm and the standard roughness deviation (RMS) was 2.57 nm. The total film thickness, determined from SEM measurements, was 342 nm. Achieved values demonstrated superior quality of prepared films over the films of Nd 2 Ti 2 O 7 which displayed height difference on the level of several hundred of nanometers [11].

High transparency and uniformity of prepared film were reflected in corresponding UV-VIS absorption spectrum which is depicted in Fig. 13. Transmission curve displayed interference pattern reaching the value of transmission maxima of the interference fringes T o-max equal to 93.39 % at 719 nm and the value of adjacent transmission minima T o-min equal to 71.44 % at 935 nm. The transmission was decreased toward shorter wavelengths due to the absorption and scattering phenomena on nanocrystals. The gap located at 330 nm should be attributed to the absorption edge of the silica substrate. Although the value of T o-max could give an impression that the film showed poor optical quality it had the origin in the reflection and interference phenomena on prepared film. The optical transmission T o of a thin film on a substrate of an infinite thickness is a complex function -Eq.( 9) [26]:

, ( 9
)
where l, n film , n sub , a film, d mean the wavelength, the refractive index of the film, the refractive index of the substrate, the absorption coefficient of the film, and the thickness of the film, respectively. For highly transparent thin films with a = 0 following equations for T o-max -Eq. ( 10), and for n film -Eq (11), can be expressed as [26]: ,

,

It should be noted that for proposed assumptions the value T o-max is independent on the properties of analyzed film thus T o-max depends only on the refractive index of the substrate. For the value n sub of the silica glass at 719 nm equal to 1.454 [START_REF] Palik | Handbook of optical constants of solids[END_REF] Eq. 10 provided the value T o-max equal to 93.39 %. Calculated value was completely identical to the measured value T o-max proving that relatively low transmission of prepared film was not caused by the film quality but by the fundamental optical principles. Eq. 11 provided the value n film at 935 nm equal to 2.138. This value is slightly smaller than a value of the refractive index of pure titanium dioxide which is equal to 2.4944 [START_REF] Palik | Handbook of optical constants of solids[END_REF]. The assumption a = 0 was not satisfied for shorter wavelengths and the values T o-max were slightly decreasing due to the absorption and scattering phenomena. Correct interpretation and theoretical calculation of the transmission in this spectral region require more advanced model and exceed the scope of this paper. Feasibility of demonstrated approach together with high refractive index and unique magnetic properties of Dy 2 Ti 2 O 7 open up possibilities for photonic applications. Beside the introduced potential applications [4,5] the challenging combination of magnetic properties with optical one can be used for the construction of integrated optical systems [START_REF] Wang | [END_REF], preparation of optical coatings [29], construction of Bragg mirrors and resonators [30] etc.

Conclusions

A generic sol-gel approach to nanocrystalline Dy 2 Ti 2 O 7 was presented. Presented approach allowed to prepare powders and highly transparent thin films with tailored size of formed nanocrystals. Kinetic parameters of the nucleation process of nanocrystals from the amorphous xerogel and the crystallization mechanism were determined. The nanocrystal growth started by homogeneous nucleation with constant nucleation rate followed by three dimensional growth. The structural evolution of nanocrystalline thin films was identical to the evolution of the xerogels. Optical properties of prepared thin films were evaluated. Presented results provided fundamental information about the crystallization behavior of Dy 2 Ti 2 O 7 and they could be used for preparation of pure nanocrystalline xerogels and highly transparent thin films with tailored structural properties.

Fig. 1 .

 1 Fig. 1. Thermal analyses of xerogels. (a) DTA and TG records for heat rate 30 K×min -1 and (b) zoom on the crystallization peaks for particular heat rates.

Fig. 2 .

 2 Fig. 2. Thermal dependence of FT-IR spectra.

Fig. 3 .

 3 Fig. 3. Thermal shift of the crystallization peaks and the corresponding fit.

Fig. 4 .

 4 Fig. 4. Calculated values of Avrami's parameter.

Fig. 5 .

 5 Fig. 5. SEM images of xerogels calcined at 1200 C.

Fig. 6 .

 6 Fig. 6. Plots of the Kissinger (left scale) and Augis-Bennett (right scale) equations.

Fig. 7 .

 7 Fig. 7. XRD patterns of calcined samples with denoted (hkl) indices. (a) Xerogels and (b) thin films.

Fig. 8 .

 8 Fig. 8. Nanoparticle sizes calculated from the Scherrer's equation.

Fig. 9 .

 9 Fig. 9. Calculation of the activation energy of the nanocrystal growth

Fig. 11 .

 11 Fig. 11. Raman spectra of xerogels and thin film.

Fig. 12

 12 Fig. 12 AFM surface profile of thin film calcined at 1200 °C

Table 1 .

 1 Values of the exponent n for Johnson-Mehl-Avrami equation for different kinds of reaction mechanisms (adapted from[START_REF] Christian | The Theory of Transformations in Metals and Alloys (Part I + II)[END_REF])

	Type of nucleation	Mechanisms of	Phase-boundary	Diffusion controlled
		nucleation	controlled growth	growth
		and the subsequent
		growth
	Site-saturated	Bulk mechanism
	-saturation of	One-dimensional	1	0.5
	nucleation	growth
	sites prior to growth	Two-dimensional	2	1
		growth
		Three-dimensional	3	1.5
		growth
		Surface mechanism
		Two-dimensional	~2	1
		growth
	Constant rate of	Bulk mechanism
	homogenous	One-dimensional	2	1.5
	nucleation during the	growth
	process	Two-dimensional	3	2
		growth
		Three-dimensional	4	2.5
		growth
		Surface mechanism
		Two-dimensional	~3	~2
		growth

Table 2 .

 2 Activation energy of crystallization calculated using different approximations

	Ozawa's	Kissinger's approximation	Augis-Bennett's
	approximation	(kJ•mol -1 )	approximation
	(kJ•mol -1 )		(kJ•mol -1 )
	766±44	778±44	793±44

Acknowledgements

The authors acknowledge the financial support from the Czech Science Foundation under contract Nº 102/13/37368P.


Hong Liu

Dimitri Leray

Stéphane Colin

Patrick Pons

Finite element modeling of nickel oxide film for Au-Ni contact of MEMS switches

Keywords: MEMS switches, electrical contact, oxide film, finite element modeling

Contamination and oxidation are inevitable in contact surfaces, especially for micro contact under low load (μN-mN). They are considered as major causes for a high contact resistance, and can lead to the failure of a contact. However, as the film formation is complicated, also it is difficult to accurately observe and measure them, the characterization of the films is not well known. In the paper, a finite element model of nickel oxide film is developed for Au-Ni contact of MEMS switches. Considering the fact that the electrical contact area is only a portion of the mechanical contact area, a model featured 'nanospots' is developed: multiple small conductive spots are scattered in a big mechanical contact asperity, and ultrathin oxide film is around the nano-spots. The size of electrical spots and mechanical asperity is calculated based on the measured electrical resistance and mechanical contact modeling respectively. The simulations results show a good match with the experimental results. This model allows us to determine some possible geometry configuration that leads to the measured contact resistance in real devices.

I. INTRODUCTION

The roughness of contact surfaces plays an important role for the performance of electrical contact. For MEMS switches, the contact force is very low (μN-mN), which brings only highest asperities in contact [START_REF] Toler | A review of microcontact physics for microelectromechanical systems (MEMS) metal contact switches[END_REF], and the real contact area Ar is only a very small portion of the apparent area Aa, as Ar =10 -1 -10 -4 Aa [START_REF] Grote | Springer Handbook of Mechanical Engineering[END_REF].

The contact resistance, or constriction resistance, can be calculated according to the ratio between contact radius and the mean free path of electrons. Three electron transport regimes were defined: diffusive, ballistic and quasi-ballistic [START_REF] Agraı̈t | Quantum properties of atomic-sized conductors[END_REF], and formula were defined correspondingly [START_REF] Holm | Electric Contacts-Theory and Applications[END_REF][START_REF] Sharvin | A Possible Method for Studying Fermi Surfaces[END_REF][START_REF] Wexler | The size effect and the non-local Boltzmann transport equation in orifice and disk geometry[END_REF].

However, it is inevitable that insulating films are formed on contact surfaces, due to oxidation, contamination, or corrosion, which prevent the current flow and cause a high resistance [START_REF] Kogut | Electrical contact resistance theory for conductive rough surfaces separated by a thin insulating film[END_REF]. The electrical contact area is smaller than the physical contact area. Fig. 1 illustrates three contact areas, in which the apparent area is also called the nominal contact area, and the bearing area is to support the mechanical load, and the electrical current passes only through the electrical contact area. Fig. 1. Three different areas in real contact (after [START_REF] Monnier | A Mechanical, Electrical, Thermal Coupled-Field Simulation of a Sphere-Plane Electrical Contact[END_REF]).

The insulating films on the contact surfaces are considered as a major cause for a high contact resistance, and can lead to the failure of a contact. The composition and formation of film is complicated and would change during switching as a result of impacting and current flow. Also, it is difficult to accurately observe and measure them, the characterization of the films is not well known. Analytical calculation and numerical modeling are the alternatives.

For surfaces with insulating film resided, one simple assumption in the literature was considering the film resistance in series with the contact resistance, and the total resistance could be calculated with:

(1)

Where t, c, f represent the total resistance, contact resistance and film resistance.

It was shown that if the oxide film is thin enough, the electrons can pass through the film by tunnel effect, and Holm [START_REF] Holm | Electric Contacts-Theory and Applications[END_REF] suggested that the total resistance can be calculated by:

(2)
Where ρ is the electrical resistivity of the contact material, and a is the contact radius; σt is the resistivity from tunnel effect, called tunnel resistivity, and is the resistance per m 2 of the surface.

Nakamura and Minowa [START_REF] Nakamura | Computer Simulation for the Conductance of a Contact Interface[END_REF][START_REF] Nakamura | Film resistance and constriction effect of current in a contact interface[END_REF] used a computer simulation based on finite element model and Monte Carlo method to investigate the conduct resistance of a contact interface. The model was composed of two tubes and a random interface between the cubes. The interface was divided into a matrix of square section N 2 spots; each spot was either conductive or insulating, which was randomly selected. It was found that the conductance kept high until the fraction of the contact area was near to zero, at the moment the current density became very high and the augmentation of the conductance depended a lot the distribution and the spaces between spots. Kogut [START_REF] Kogut | Electrical contact resistance theory for conductive rough surfaces separated by a thin insulating film[END_REF][START_REF] Kogut | Electrical performance of contaminated rough surfaces in contact[END_REF] presented an electrical theory for conductive and rough fractal surfaces separated by a thin insulating film, and for degraded surface, where the metallic conductance area was only a fraction of the contact area. It was found that, for a thin film, thickness < 5 nm, where the tunnel effect can occur, -If the film covered over whole surface uniformly, then the constriction resistance was negligible compared to the resistance from tunnel effect. The total contact resistance was about 10 5 -10 6 higher than that of clean surface [START_REF] Kogut | Electrical contact resistance theory for conductive rough surfaces separated by a thin insulating film[END_REF].

-If the film was damaged, the tunnel effect became very weak compared to that of an intact film, and the I-V curves kept Ohmic [START_REF] Kogut | Electrical performance of contaminated rough surfaces in contact[END_REF].

It was noticed that the ballistic charge transport mechanism was dominated in Kogut's modeling, whereas it was supposed that all contact spots were in diffusive mode in Nakamura's work.

The idea of degraded films was also considered in [START_REF] Hisakado | Effects of surface roughness and surface films on contact resistance[END_REF][START_REF] Crinon | The effect of surface roughness, oxide film thickness and interfacial sliding on the electrical contact resistance of aluminium[END_REF], where the amount of oxide film significantly impacted constriction resistance [START_REF] Hisakado | Effects of surface roughness and surface films on contact resistance[END_REF] because it prevented metal-to-metal contact [START_REF] Crinon | The effect of surface roughness, oxide film thickness and interfacial sliding on the electrical contact resistance of aluminium[END_REF]. However, no related research has been done for the case of electrical contact of MEMS switches, and there was a lack of comparison between analytical or numerical models to the experimental results in Kogut's work.

In this study, finite element models of Au-Ni contacts covered with oxide films were developed and the simulated resistance was found to be consistent with the experimental results.

II. FINITE ELEMENT MODELING

A. Contact asperity modeling

The contact asperity, which is assumed to support the mechanical load, is modeled as a 3-dimensional cone-truncated asperity as in previous study [START_REF] Liu | An Asperity-Based Finite Element Model for Electrical Contact of Microswitches[END_REF]. Three parameters were extracted from AFM scanning data and defined as: top radius (a), height (h) and the angle of the side with respect to the surface plan (α), as shown in Fig. 2. It was shown that the top radius of asperity had most significant impact on the electrical resistance and thermoelectrical behavior, whereas the height and angle only had observable influence at large values [START_REF] Liu | An Asperity-Based Finite Element Model for Electrical Contact of Microswitches[END_REF]. For the contact surfaces investigated in the study, which are sputtering Au-Ni, the influence of angle and height are ignored.

The top radius of the asperity (a) is defined as 120 nm in the modeling, which is based on the finite element mechanical contact modeling [START_REF] Liu | Finite Element Based Surface Roughness Study for Ohmic Contact of Microswitches[END_REF], with the applied contact force of 150 μN, accordance to the experimental results [START_REF] Broue | Multi-Physical Characterization of Micro-Contact Materials for MEMS Switches[END_REF]. And the height of the asperity is defined as h=6 nm and the angle of side as α=10° from the AFM data [START_REF] Liu | An Asperity-Based Finite Element Model for Electrical Contact of Microswitches[END_REF].

B. Assumptions of insulating film

The distribution of insulating film is unknown, and two assumptions are made to take account of the residence of insulating film on the contact surface.

-Model I: The intact film, where the contact area is covered by a thin insulating film completely.

-Model II: The 'nano-spots' spots. The metallic contact spots are scattered in the oxide film, or to say, the insulating film is damaged by the metallic spots. The size of metallic spots is in order of nano meters, so called 'nano-spots' in the paper.

The schemas of the two models are shown in Fig. 3 and Fig. 4. It should be kept into mind that the dimension is not to scale in the figures, and the size of the model is given in next section.

Holm [START_REF] Holm | Electric Contacts-Theory and Applications[END_REF] indicated that an oxide film of a few atoms thick could be formed on contact surfaces after several hours exposing to the ambient environment. The crystalline parameter of NiO, the main oxide for the Au-Ni contact, is about 0.4 nm [START_REF] Hotový | Characterization of NiO thin films deposited by reactive sputtering[END_REF], so the thickness of oxide film is defined from 1 nm to 6 nm in the modeling. 6 nm is chosen as it was suggested that no tunnel effect could occur for the thickness larger than 5 nm [START_REF] Kogut | Electrical performance of contaminated rough surfaces in contact[END_REF].

In model II, the metallic conductive spots, are modeled as cylinders, and a few geometrical parameters are defined as:

-r 1, the radius of nano-spots, -h1, the thickness of oxide film, also the height of nano-spots, -m1, the distance between two nearby nano-spots, -n1, the number of nano-spots.

Simulations were performed with different values of these geometrical parameters, and a reference model is defined as:

-r1= 10 nm.

-The thickness of oxide film is h1= 3 nm.

-The number of nano-spots is n1=1.

The radius of 10 nm is estimated with Holm equation, Rc=(ρ1+ρ2)/4a, with the measured contact resistance 2.55 Ω at 10 mA at 150 μN [START_REF] Broué | Analyse multi physique des sources de défiabilisation du microcontact électrique à destination des interrupteurs MEMS[END_REF], the calculated value is 8.9 nm, and rounded as 10 nm in the reference model. In the experiments, the contact force was applied using a nano-indentation tip on a Au-Ni MEMS switch, and the resistance was measured using the four points method [START_REF] Broué | Analyse multi physique des sources de défiabilisation du microcontact électrique à destination des interrupteurs MEMS[END_REF]. Different current level was applied on the test vehicle, but only 10 mA was used in the modeling to avoid the effect of Joule heating. The influence of the geometrical parameters on electrical resistance will be discussed in section IV.

C. Finite element (FE) modeling

The FE model is built up with commercial software package ANSYS of version 11.0. The 3-D 10 node coupledfield element SOLID227 is used for meshing, with the degrees of freedom (DOF) of temperature and voltage activated. Fig. 5 (a) shows the meshing of the modeling and the boundary conditions. Low current of 10 mA is applied in the modeling, and the dependence of material properties on temperature is not considered in the study.

The electrical resistivity of the NiO, according to the literature [START_REF] Lu | Properties of nickel oxide thin films deposited by RF reactive magnetron sputtering[END_REF][START_REF] Wakefield | The fabrication and characterization of nickel oxide films and their application as contacts to polymer/porous silicon electroluminescent devices[END_REF][START_REF] Sato | Transparent conducting p-type NiO thin films prepared by magnetron sputtering[END_REF], is defined in the range of 0.14×10 7 to 60×10 7 nΩ×m in Model I, and kept as 0.3×10 7 nΩ×m for Model II if not specified. The incertitude of the electrical resistivity of oxide film is discussed in the last part of section IV. The material properties of the contact materials and the oxide film are listed in TABLE I.

The two contact parts are modeled as a cuboid, as shown in Fig. 5(a). The length of the cuboid is equal to the width, as 3.5 μm, and the height is 1 μm. The size of model has been validated to be large enough to model the constriction resistance accurately for an asperity of radius less than 200 nm.

In Model II, single nano-spot and multiple nano-spots are both modeled. For the case of one nano-spot, the spot is located at the center of underlying asperity. For multiple metallic nanospots, the nano-spots are all of the same size, and distributed uniformly on the underlying large asperity. The thickness of oxide film is kept constant as 3 nm for multiple nano-spots modeling.

The number of nano-spots is varied in the modeling, namely as 1, 2, 4 and 9. The distance between spots is 120 nm for 2 spots and 4 spots, and is reduced to 60 nm for 9 spots due to the limited size of asperity, which has radius of 120 nm. Fig. 5 (b) shows the distribution of 9 nano-spots, and a zoom-in side view of meshing in nano-spots is shown in Fig. 5 (c).

The distribution of spots is also investigated by varying the distance between spots. The simulations were performed with 4 nano-spots of radius 10 nm, and the distance between spots were varied from 30, 60 to 120 nm.

NiO

III. SIMPLIFIED ANALYTICAL MODEL

An analytical model is proposed to evaluate approximately the electrical resistance. It has been shown in previous study [START_REF] Liu | An Asperity-Based Finite Element Model for Electrical Contact of Microswitches[END_REF] that the finite element modeling could not predict the contact resistance with ballistic transport, so all the electrical resistance is calculated based on Holm equation [START_REF] Holm | Electric Contacts-Theory and Applications[END_REF]. Further more, it is assumed that:

-Tunnel effect is not taken into account, and the oxide film is considered as completely insulating;

-The electrical interaction between nano-spots is negligible;

-The supplementary electrical resistance of each nanospot is the same as the resistance of a cylinder of height h1 and radius r1, i.e.

. Thus the electrical resistance from nano-spots can be evaluated with:

3

The calculation with (3) will be compared with the numerical simulations.

IV. RESULTS AND DISCUSSIONS

A. Results with Model I

The simulations were performed with the thickness of oxide film of 1 nm. The simulated total resistance is ranged from 26 Ω to 11300 Ω for the electrical resistivity of NiO varied from 0.14×10 7 to 60×10 7 nΩ×m, which is much larger than the measured contact resistance, 2.55 Ω. The constriction resistance of pure Au-Ni contact is only 0.18 Ω, which suggests that the resistance of oxide film dominates in the total electrical resistance.

It is concluded that in Model I, assuming the contact area is covered completely by the oxide film, the electrical resistance is too high even the film is very thin, which agrees with the results in [START_REF] Kogut | Electrical contact resistance theory for conductive rough surfaces separated by a thin insulating film[END_REF]. Model I does not allow us to get the resistance close to the measured values.

B. Results with Model II 1) One nano-spot

The first investigated parameter is the radius of nano-spot. It is varied from 5 to 20 nm in the modeling, with the thickness of the oxide film kept as 3 nm. Fig. 6 shows the influence of the radius of nano-spots on the electrical resistance, compared with analytical calculation results. As expected, the resistance decreases with the radius of nano-spots increases. The analytical results compare well with the simulated values, and the maximum relative difference is only 1.8%. It is also noticed that the nano-spot of radius 10 nm results in the electrical resistance of 2.96 Ω, which is close to the measured contact resistance of 2.55 Ω [START_REF] Broué | Analyse multi physique des sources de défiabilisation du microcontact électrique à destination des interrupteurs MEMS[END_REF].

The simulation was then performed with the radius of nano-spot kept as 10 nm, and the thickness of the oxide film varied from 1.5 nm to 6 nm. 1.5 nm was used in the modeling due to the meshing problem with 1 nm. The results with different thicknesses of oxide film are plotted in Fig. 7. It is shown that the thickness of oxide film does not have significant impact as the radius of nano-spots, and the total resistance increases 38% with the thickness of the film increases from 1.5 nm to 6 nm. The difference between the analytical calculation and the FE model also increases a little bit with thicker oxide film.

It can be deduced that, with one conductive nano-spot on asperity surface, the analytical model can predict the FEM results correctly. The nano-spot with radius of 10 nm, and height of 1.5 nm or 3 nm results in the total resistance of 2.96 Ω or 2.58 Ω, which is close to the measured resistance of 2.55 Ω.

2) Multiple nano-spots

For the case of multiple nano-spots in electrical contact, the number of spots is investigated, and two radii of spots are modeled, namely as 10 nm and 5 nm.

As shown in Fig. 8, the electrical resistance decreases as the number of spots increases. This is expected as the metallic area increases with more spots, which facilitates the current flow. The electrical resistance decreases more sharply with smaller nano-spots, and the resistances tend to become stable with more spots in both cases. Due to the electrical interaction between nano-spots, the analytical calculations underestimate the electrical resistance, and the gaps between analytical calculations and FEM become bigger with the number of nanospots increases. The relative difference with nano-spot of 10 nm is two times of that of 5 nm, and it can be deduced that the interaction between larger nano-spots is more significant as they are closer.

It is also noticed that, the electrical resistance of 4 or 2 nano-spots of radius 5 nm, which are 1.90 Ω and 3.66 Ω, respectively, are close to the measured value of 2.55 Ω [START_REF] Broué | Analyse multi physique des sources de défiabilisation du microcontact électrique à destination des interrupteurs MEMS[END_REF].

The influence of the distance between nano-spots is plotted in Fig. 9. Due to the electrical interaction between nano-spots, the electrical resistance increases when the spots are closer, and the difference between FEM and analytical model becomes bigger.

3) The influence of electrical resistivity of oxide film

The electrical resistivity of nickel oxide film is defined as 3×10 6 nΩ×m in previous sections, which is chosen based on the values in the literature: (0.14-60)×10 7 (nΩ×m) [START_REF] Lu | Properties of nickel oxide thin films deposited by RF reactive magnetron sputtering[END_REF][START_REF] Wakefield | The fabrication and characterization of nickel oxide films and their application as contacts to polymer/porous silicon electroluminescent devices[END_REF][START_REF] Sato | Transparent conducting p-type NiO thin films prepared by magnetron sputtering[END_REF]. These values were measured with samples fabricated in controlled conditions, either thermal evaporation or sputtering, and might be different from the natural formed oxide film. Despite the incertitude on the resistivity value of NiO, it is noticed that the ratio of electrical resistivity between NiO and nickel is at least 20000, which suggests that the oxide film probably can be considered insulating. However, there are also some materials, like ruthenium, the oxide film is not very insulating. The electrical resistivity of Ru oxide film is only 5 times bigger than that of ruthenium (with ρRu=76 (nΩ×m), and ρRuO2=350 (nΩ×m) [START_REF] Hawk | Properties of ruthenium oxide coatings[END_REF]. In this section, the influence of electrical resistivity of oxide film on the total electrical resistance is discussed, and a new analytical model is proposed.

The simulation is performed with one conductive nano-spot on asperity surface. The radius of nano-spot is 10 nm, and the thickness of oxide film is 3 nm. The minimum value of electrical resistivity of oxide film is defined as about 1/20 of the value in the literature, which is about 1/1000 of the electrical resistivity of nickel. The range of electrical resistivity of oxide film is varied then from 7.5×10 4 to 60 ×10 7 (nΩ×m) in the modeling. The simulated total resistance is plotted in Fig. 10, and it is found that the total electrical resistance drops significantly with the electrical resistivity of NiO decreases from 3×10 6 (nΩ×m) to 0.15×10 6 (nΩ×m), this indicates that the conductance of oxide film is not negligible when its electrical resistivity is small, here about 1/2000 of that of nickel.

We assume the electrical resistance from the oxide film laying on the asperity can be evaluated by the bulk resistance formula (R=ρ×l/A), as:

(4)
The total contact resistance then can be evaluated considering the electrical resistance from oxide film (4) and the resistance from nano-spots (3) in parallel, and we have:

(5)

The total contact resistance calculated using ( 5) is also plotted in Fig. 10, compared with the ones using (3). It is found that the new analytical model can predict the contact resistance more accurately for the oxide film with low electrical resistivity. It is also noticed that when the Rnano is close to the Rtotal, the oxide film can be considered as completely insulating. For the oxide film with very low electrical resistivity, like RuO2, the results would be much different from the Au-Ni contact discussed in the study.

V. CONCLUSIONS

A finite element model was developed for the oxide film of Au-Ni contact. Two types of insulating film are investigated: intact and damaged. A 'nano-spots' model is built up for the damaged insulating film. The metallic conductive spots are scattered in the oxide film that is laid on a mechanical contact asperity. The size of mechanical asperity is calculated from the FE contact modeling, and that of nano-spots are evaluated from the measured resistance. The size of nano-spots and thickness of oxide films are investigated, also the number and distribution of nano-spots for the case of multiple metallic spots.

Two simplified analytical models are proposed, and the influence of electrical resistivity of oxide film is discussed. It is shown that the oxide film cannot be considered insulating when the electrical resistivity is about 1/2000 of that of the metal, and the improved analytical model is more accurate for this case.

It is found that one nano-spot with radius of 10 nm or 4 or 2 nano-spots with radius of 5 nm produce the electrical resistance close to the measured results. The FE model and analytical model shed light on the possible geometry configuration of insulating film on the contact surfaces. It should also be kept in mind that the ballistic transport is not considered in the study.

Fig. 2 .

 2 Fig. 2. Illustration of an asperity in cone-truncated geometry.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Illustration of Model I: considering the contact area is covered completely by the insulating film.

Fig. 5 .

 5 Fig. 5. (a) Meshing of the finite element model and the boundary conditions, (b) zoom in the nano-spots model, top view, (c) meshing of nano-spots, side view.

Fig. 6 .

 6 Fig. 6. Simulation results of electrical resistance as a function of nano-spot radius.

Fig. 7 .

 7 Fig. 7. Simulation results of electrical resistance as a function of oxide film thickness.

Fig. 8 .

 8 Fig. 8. Simulation results of electrical resistance as a function of number of nano-spots.

Fig. 9 .

 9 Fig. 9. Simulation results of electrical resistance as a function of distance between nano-spots. The simulations are with 4 nano-spots.

Fig. 10 .

 10 Fig. 10. Simulation results of electrical resistance as a function of number of nano-spots.

TABLE I .

 I MATERIAL PROPERTIES FOR AU, NI AND NIO IN FEM, THE VALUES LABELED WITH * ARE FROM WIKIPEDIA WEBSITE.

	Properties	Au	Ni	NiO
	Electrical resistivity (nΩ×m)	22.14*	68.4*	(0.14 -60)×10 7 [19-21]
	Thermal
	conductivity (W/(m×K))	318*	90.9*	40 [23]
	Mean free path (nm)	38 [22]	9


Abbas El Dor
email: abbas.eldor@mines-nantes.fr

Patrick Siarry
email: siarry@u-pec.fr

Effect of the dynamic topology on the performance of PSO-2S algorithm for continuous optimization

PSO-2S is a multi-swarm PSO algorithm using charged particles in a partitioned search space for continuous optimization problems. In order to improve the performance of PSO-2S, this paper proposes a novel variant of this algorithm, called DPSO-2S, which uses the Dcluster neighborhood topologies to organize the communication networks between the particles. Experiments were conducted on a set of classical benchmark functions. The obtained results prove the effectiveness of the proposed algorithm.

Introduction

The concept of particle swarm optimization (PSO) is based on social behavior to exchange information between the particles in a swarm. Thus this property can be modelized thanks to a graph: two particles P i and P j of the swarm S are connected if a communication can be established between them. The set of edges between each particle P i and its neighbours N e i forms the communication graph, also called the topology. Hence, the chosen topology can greatly affect the performance of the PSO algorithm. In this paper, we present a new dynamic topology, called Dcluster, which is a combination of two existing topologies (Four-clusters [START_REF] Mendes | The Fully Informed Particle Swarm: Simpler, Maybe Better[END_REF] and Wheel [START_REF] Kennedy | Population structure and particle swarm performance[END_REF]). This topology was integrated in our proposed algorithm called PSO-2S, introduced in [START_REF] El Dor | A multi-swarm PSO using charged particles in a partitioned search space for continuous optimization[END_REF]. PSO-2S is a multi-swarm PSO algorithm using charged particles in a partitioned search space for continuous optimization problems. The performance of PSO-2S with the Dcluster topology is analysed and compared to that of PSO-2S without Dcluster, using a set of benchmark test functions. Comparisons show that the use of Dcluster improves significantly the performance of PSO-2S.

Particle Swarm Optimization

The particle swarm optimization (PSO) [START_REF] Kennedy | A new optimizer using particle swarm theory[END_REF] is inspired originally by the social and cognitive behavior existing in the bird flocking. The algorithm is initialized with a population of particles randomly distributed in the search space, and each particle is assigned a randomized velocity. Each particle represents a potential solution to the problem. The particles fly over the search space, keeping in memory the best solution encountered. At each iteration, each particle adjusts its velocity vector, based on its momentum, influences of its best solution and of the best solution of its neighbors, then computes a new point to be evaluated. The displacement of a particle is influenced by three components:

1. Physical component: the particle tends to keep its current direction of displacement; 2. Cognitive component: the particle tends to move towards the best site that it has explored until now; 3. Social component: the particle tends to rely on the experience of its congeners, then moves towards the best site already explored by its neighbors.

In this paper, the swarm size is denoted by s, and the search space is ndimensional. In general, a particle i has three attributes: the current position

X i = (x i,1 , x i,2 , ..., x i,n ), the current velocity vector V i = (v i,1 , v i,2 , ..., v i,n
) and the past best position P best i = (p i,1 , p i,2 , ..., p i,n ). The best position found in the neighborhood of the particle i is denoted by Gbest i = (g i,1 , g i,2 , ..., g i,n ). These attributes are used to update iteratively the state of each particle in the swarm. The objective function to be minimized is denoted by f . The velocity vector V i of each particle is updated using the best position it visited so far and the overall best position visited by its neighbors. Then, the position of each particle is updated using its updated velocity per iteration. At each step, the velocity of each particle and its new position are updated as follows:

v i,j (t + 1) = wv i,j (t) + c 1 r 1i,j (t) [p i,j (t) -x i,j (t)] + c 2 r 2i,j (t) [g i,j (t) -x i,j (t)] (1)
x i,j (t + 1) = x i,j (t) + v i,j (t + 1)

(2)

x i,j is the position and v i,j is the velocity of the ith particle (i ∈ 1, 2, ..., s) of the jth dimension (j ∈ 1, 2, ..., n). Where w is called inertia weight, c 1 , c 2 are the learning factors and r 1 , r 2 are two random numbers selected uniformly in the range [0, 1].

PSO-2S algorithm

In this section, we present the first version of PSO-2S [START_REF] El Dor | A multi-swarm PSO using charged particles in a partitioned search space for continuous optimization[END_REF]. PSO-2S is based on three main ideas. The first is to use two kinds of swarms: a main swarm, denoted by S1, and s auxiliary ones, denoted by S2 i , where 1 ≤ i ≤ s. The second idea is to partition the search space into several zones in which the auxiliary swarms are initialized (the number of zones is equal to the number of auxiliary swarms, thus is equal to s). The last idea is to use the concept of the electrostatic repulsion heuristic to diversify the particles for each auxiliary swarm in each zone.

To construct S1, we propose to perform the auxiliary swarms S2 i several times in different areas, and then each best particle for each S2 i is saved and considered as a new particle of S1. To do so, the population of each auxiliary swarm is initialized randomly in different zones (each S2 i is initialized in its corresponding zone i). After each of these initializations, K displacements of particles, of each S2 i , are performed in the same way of standard PSO. Then the best solution found by each auxiliary swarm is added to S1. The number of initializations of S2 i is equal to the number of particles in S1.

As we mentioned above the second idea is to partition the search space [min d , max d ] N into several zones (max zone zones). Then, we calculate the center d and the step d of each dimension separately, according to (3) and ( 4). The step d are similar in the case of using a square search space.

center d = (max d -min d )/2
(3)

step d = center d /max zone (4)
where max zone is a fixed value, and d is the current dimension (1

≤ d ≤ N ).
The sizes of the zones of the partitioned search space are different (Z 1 < Z 2 < . . . < Z maxzone ). Therefore, the number of particles in S2 i , denoted by S2 isize , depends on its corresponding zone size. Indeed, a small zone takes less particles and the number of particles increases when the zone becomes larger. The size of each auxiliary swarm is calculated as follows:

S2 isize = num zone * nb particle (5)
where num zone = 1, 2, ..., max zone , is the current zone number and nb particle is a fixed value. After the initializations of the auxiliary swarms in different zones (Z i , S2 i ), an electrostatic repulsion heuristic is applied to diversify the particles and to widely cover the search space [START_REF] Conway | Sphere packings, lattices and groups[END_REF]. This technique is used in an agent-based optimization algorithm for dynamic environments [START_REF] Lepagnot | A new multiagent algorithm for dynamic continuous optimization[END_REF]. Therefore, this procedure is applied in each zone separately, hence each particle is considered as an electron. Then a force of 1/r 2 is applied, on the particles of each zone, until the maximum displacement of a particle during an iteration becomes lower than a given threshold (where r is the distance between two particles, is typically equal to 10 -4 ). At each iteration of this procedure, the particles are projected in the middle of the current zone, before reapplying the heuristic repulsion.

Dynamic cluster topology (Dcluster)

Dcluster is a dynamic topology that works as follows [START_REF] El Dor | Dynamic Cluster in Particle Swarm Optimization algorithm[END_REF]. At each iteration, the particles are sorted in a list according to their personal best fitness in increasing order, so that the worst particle has an index equal to 1 in the list (the size of the list is equal to the size of the swarm). Then, the list is partitioned into several sub-lists which correspond to a cluster in the proposed topology. The first cluster which has the "worst" particle, called central cluster, is placed in the center of the topology. Each particle of the central cluster is connected to other clusters by one of their particles; the first worst particle of the central cluster is linked to the worst particle in the second cluster, the second one is linked to the third cluster also by its worst particle, and so on, as in Figure 1-(a). All clusters in this topology, including the central cluster, have a fully connected neighborhood. The reason why the central cluster is linked to each other cluster by only one gateway and with the worst particle of the latter is to avoid a premature convergence to a local optimum by slowing down the propagation of the information in the whole swarm. Figure 1-(b) illustrates the final structure of the proposed topology. The best particles

Experimental results

Table 1 presents the settings of each problem, the number of function evaluations (Max. FEs), the success rate and the mean best value of 100 runs. The best results among those obtained by the two algorithms are shown in bold. From the experiments, we can notice that DPSO-2S (PSO-2S using Dcluster topology) obtains the best results on most of the functions used. Hence, DPSO-2S algorithm outperforms PSO-2S, except for Rosenbrock and Shif ted Rosenbrock. Thus, this algorithm leads to a significant improvement over the previous PSO-2S.

Conclusion

In this paper, a new dynamic topology, called Dcluster, based on two static neighbourhood topologies was presented. Dcluster was integrated to our multiswarm algorithm PSO-2S. Experimental results indicate that Dcluster improves the search performance of the previous algorithm.

Fig. 1 .

 1 Fig. 1. (a) -The partitioning of the list into sub-lists. (b) -The structure of Dcluster topology.

In conclusion, the improvement of the PSO-2S algorithm due to the integration of Dcluster topology opens the gate to apply this dynamic topology in other PSO algorithms.


Aagtca

Aggacc

Categories of phylogenomics methods

GAAACCACTCTCTAAGAC 5 δ1 δ 13 δ 14 δ 15 δ 34 δ 35 δ 45                                           δ2 δ 12 δ 15 δ 25                   δ3 δ δ δ δ δ δ                                          
                                          δ2 δ 12 δ 15 δ 25                   δ3 δ δ δ δ δ δ                                           4
Every tree with branch lengths can be represented with a vector of tree distances by computing the distance d T ij between each pair of taxa in the tree.

4 2 T =(T , b) 3 1 b 1 =1 b 2 =2 b 3 =3 b 5 =1 b 4 =2 ⇔ d T =        d 12 = 3 d 13 = 5 d 14 = 4 d 23 = 6 d 24 = 5 d 34 = 5        Vector of tree distances d T ij = e∈P ij b e

Tree distances

Every vector of tree distances is the product of the topological matrix of T by the branch lengths vector of T . i.e. d T = Ab. A distance-based method : the weighted least square method (WLS)

d T =      d 12 = 3 d 13 = 5 d 14 = 4 d 23 = 6 d 24 = 5 d 34 = 5      Vector of tree distances b =     b 1 = 1 b 2 = 2 b 3 = 3 b 4 = 2 b 5 = 1     Branch lengths vector A =      1 1 0 0 0 1 0 1 0 1 1 0 0 1 1 0 1 1 0 1 0 1 0 1 1 0 0 1 1 0      Topological matrix 4 2 T =(T , b) 3 1 b 1 =1 b 2 =2 b 3 =3 b 5 =1 b 4 =2 d T = Ab minimize ij w ij (δ ij -d T ij ) 2 input distances (estimated from 1 gene alignment) input topology δ =           δ 12 =
Branch lengths estimation with WLS: Hypothesis: Any gene G k induces approximately the same tree up to a scale factor α k and the removal of a number of lineages.

minimize ij w ij (δ ij -d T ij ) 2 can be solved analytically in O(n 3 ) input distances (estimated from 1 gene alignment) input topology δ =           δ 12 = 0.140 δ 13 = 0.163 δ 14 = 0.288 δ 15 = 0.336 δ 23 = 0.188 δ 24 = 0.413 δ 25 = 0.411 δ 34 = 0.298 δ 35 = 0.213           3 b 3 =? 1 b 1 =? 2 b 2 =?
E A B 1,5 3 1,5 1,5 E A C B 3 2 1 1 E A C D 2 0,75 0,75 1,5

New methode ERaBLE: Evolutionnary Rates and Branch Length Estimation

Hypothesis: The Proportional model. Any gene G k induces the same tree up to a scale factor α k , representing the inverse of its evolutionary rate r k .

α1 × δ (1) AE ≃ α2 × δ (2) AE ≃ α3 × δ (3) AE E A B 2 2 1 1 E A C B 2 2 1 1 E A C D 2 1 1 2 δ (1) AE =6 r1=1 α1=1 δ (2) AE =3.75 r2=0.75 α2=1.33 δ (3) AE =9 r3=1.5 α2=0.67 2 1 1 E A C 2 E A 1,5 0,75 0,75 1,5 3 1,5 1,5 E A C
3 methode ERaBLE: Evolutionnary Rates and Branch Length Est

Hypothesis: The Proportional model. Any gene G k induces the same tree up a scale factor α k , representing the inverse of its evolutionary rate

α1 × δ (1) AE ≃ α2 × δ (2) AE ≃ α3 × δ (3) AE E A B 2 2 1 1 E A C B 2 2 1 1 E A C D 2 1 1 2 δ (1) AE =6 r 1 =1 α 1 =1 δ (2) AE =3.75 r 2 =0.75 α 2 =1.33 δ (3) AE =9 r 3 =1.5 α 2 =0.67
New methode ERaBLE: Evolutionnary Rates and Branch Length Estimation

Hypothesis: The Proportional model. Any gene Gk induces the same tree up to a scale factor αk , representing the inverse of its evolutionary rate rk .

α1 × δ (1) AE ≃ α2 × δ (2) AE ≃ α3 × δ (3) AE E A B 2 2 1 1 E A C B 2 2 1 1 E A C D 2 1 1 2 δ (1) AE =6 r1=1 α1=1 δ (2) AE =3.75 r2=0.75 α2=1.33 δ (3) AE =9 r3=1.5 α2=0.67 2 1 1 E A C D 2 E A B 1,

α1 × δ

(1)

AE ≃ α2 × δ (2) AE ≃ α3 × δ (3) AE E A B 2 2 1 1 E A C B 2 2 1 1 E A C D 2 1 1 2 δ (1) AE =6 r1=1 α1=1 δ (2) AE =3.75 r2=0.75 α2=1.33 δ (3) AE =9 r3=1.5 α2=0.67 2 1 1 E A C 2 E A 1,5 0,75 0,75 1,5 3 1,5 1,5 E A C

methode ERaBLE: Evolutionnary Rates and Branch Length Est

Hypothesis: The Proportional model. Any gene G k induces the same tree up a scale factor α k , representing the inverse of its evolutionary rate

α1 × δ (1) AE ≃ α2 × δ (2) AE ≃ α3 × δ (3) AE E A B 2 2 1 1 E A C B 2 2 1 1 E A C D 2 1 1 2 δ (1) AE =6 r 1 =1 α 1 =1 δ (2) AE =3.75 r 2 =0.75 α 2 =1.33 δ (3) AE =9 r 3 =1.5 α 2 =0.67

New methode ERaBLE: Evolutionnary Rates and Branch Length Estimation

Hypothesis: The Proportional model. Any gene Gk induces the same tree up to a scale factor αk , representing the inverse of its evolutionary rate rk . Objective: find the branch lengths b1 , . . . , b2n-3 of T , the scale factors α1 , . . . , αm of the m genes, solution of the problem:

α1 × δ (1) AE ≃ α2 × δ (2) AE ≃ α3 × δ (3) AE E A B 2 2 1 1 E A C B 2 2 1 1 E A C D 2 1 1 2 δ (1) AE =6 r1=1 α1=1 δ (2) AE =3.75 r2=0.75 α2=1.33 δ (3) AE =9 r3=1.5 α2=0.67 2 1 1 E A C D 2 E A B 1,
         minimize m k=1 i,j∈L k w (k) ij ( αk δ (k) ij -d T ij ) 2
Input data: m distance vectors δ 1 , . . . , δm where δ k = (δ

(k) ij ) is defined on the taxa set L k .
A given topology T defined on the taxa set L = m k=1 L k .

Objective: find the branch lengths b1 , . . . , b2n-3 of T , the scale factors α1 , . . . , αm of the m genes, solution of the problem:

         minimize m k=1 i,j∈L k w (k) ij ( αk δ (k) ij -d T ij ) 2 subjet to m k=1 Z k αk = 1
New methode ERaBLE: Evolutionnary Rates and Branch Length Estimation Objective: find the branch lengths b1 , . . . , b2n-3 of T , the scale factors α1 , . . . , αm of the m genes, solution of the problem:

         minimize m k=1 i,j∈L k w (k) ij ( αk δ (k) ij -d T ij ) 2 subjet to m k=1 Z k αk = 1
Outputs:

be et rk = 1 αk .

The minimization problem can be solved with the method of Lagrange multipliers and leads to the linear system in O(n + m) equations and unknowns:

                       δ t 1 W1δ1 0 • • • 0 0 δ t 2 W2δ2 . . . . . . . . . . . . 0 • • • • • • δ t m Wmδm ] ] ] ] -A t 1 W1δ1 -A t 2 W2δ2 . . . -A t m Wmδm [ [ [ [ Z1 Z2 • • • Zm [ -δ t 1 W1A1 ] [ -δ t 2 W2A2 ] . . . [ -δ t m WmAm ] 1 1 . . . 1           m k=1 A t k WkAk           0 . . . . . .       α1 α2 . . . αm b1 b2 . . .      0 . . . . . . . . . . . . . . . . . . . . . 0 1                        m=# genes n=# taxa
Filling the matrix and standard system solving in O((n + m) 3 + mn 4 )

Solving in O(n 3 + mn 2 ) with our algorithms:

Filling the matrix in O(mn 2 ) instead of O(mn 4 ) Block-solving system in O(n 3 ) instead of O(n + m) 3      δ t 1 W1δ1 0 • • • 0 0 δ t 2 W2δ2 . . . . . . . . . . . . 0 • • • • • • δ t m Wmδm ] ] ] ] -A t 1 W1δ1 -A t 2 W2δ2 . . . -A t m Wmδm [ [ [ [ Z1 Z2 • • • Zm [ -δ t 1 W1A1 ] [ -δ t 2 W2A2 ] . . . [ -δ t m WmAm ] 1 1 . . . 1           m k=1 A t k WkAk           0 . . . . . .       α1 α2 . . . αm b1 b2 . . .      0

Conclusion

ERaBLE gives branch lengths to phylogenomic trees (e.g. estimated with supertree methods).

ERaBLE is relatively accurate in the estimations of branch lengths and of gene rates.

ERaBLE is fast with a complexity in O(n 3 + mn 2 ) (linear in m).

Thank you for your attention.

Any questions?

Fundings:

  output d T = Ab Hypothesis: Any gene G k induces approximately the same tree up to a scale factor α k and the removal of a number of lineages.

  gene G k induces approximately the same tree up to a scale factor α k and the removal of a number of lineages.New methode ERaBLE: Evolutionnary Rates and Branch Length EstimationHypothesis: The Proportional model. Any gene G k induces the same tree up to a scale factor α k , representing the inverse of its evolutionary rate r k .

  m distance vectors δ 1 , . . . , δm where δ k = (δ (k) ij ) is defined on the taxa set L k . A given topology T defined on the taxa set L = m k=1 L k . New methode ERaBLE: Evolutionnary Rates and Branch Length Estimation Input data: m distance vectors δ 1 , . . . , δm where δ k = (δ (k) ij ) is defined on the taxa set L k . A given topology T defined on the taxa set L = m k=1 L k .

  Input data: m distance vectors δ 1 , . . . , δm where δ k = (δ(k) ij ) is defined on the taxa set L k .A given topology T defined on the taxa set L = m k=1 L k .

  and standard system solving in O((n + m) 3 + mn 4 ) Solving in O(n 3 + mn 2 ) with our algorithms:Filling the matrix in O(mn 2 ) instead of O(mn 4 ) Block-solving system in O(n 3 ) instead of O(n + m) 3 Complexity: ERaBLE O(n 3 + mn 2 ) vs. WLS O(n 3 )Results on the OrthoMaM dataset[Douzery et al. 2014] input data: m=6953 nucleotide exon alignments over n=4 to 40 mammals. input topology: the topology of the 40 mammals present in OrthoMaM:


Malgorzata Arlet

Ronan Jubin

Nobuo Masataka

Alban Lemasson

Grooming-at-a-distance by exchanging calls in nonhuman primates

Keywords: call exchange, grooming, dominance rank, kinship, macaques

The "social bonding hypothesis" predicts that, in large social groups, functions of gestural grooming should be partially transferred to vocal interactions. Hence, vocal exchanges would have evolved in primates to play the role of grooming-at-a-distance in order to facilitate the maintenance of social cohesion. However, there are few empirical studies testing this hypothesis. To address this point, we compared the rate of contact call exchanges between females in two captive groups of Japanese macaques as a function of female age, dominance rank, genetic relatedness and social affinity measured by spatial proximity and grooming interactions. We found a significant positive relationship between the time spent on grooming by two females and the frequency with which they exchanged calls. Our results conform to the predictions of the social bonding hypothesis, i.e. vocal exchanges can be interpreted as grooming-at-a-distance.

INTRODUCTION

The reasons why and how vocal communication became of prime importance in the course of primate evolution, potentially leading to human language, remains unclear. One of the nonmutually exclusive hypotheses relates to the need for spatial cohesion. This suggests that the habitat is crucial in determining the communicative modality [START_REF] Partan | Issues in the classification of multimodal communication signals[END_REF]. Most primate species live in dense forested environments and vocal communication, in contrast to visual and gestural communication, has major advantages (e.g. long-distance 3D communication with several individuals). Another hypothesis related to the need for social cohesion suggests that language evolved with a social function to facilitate bonding between individuals in large groups [START_REF] Dunbar | Coevolution of neocortical size, group size and language in humans[END_REF]. This social bonding hypothesis implies that vocal exchange rates evolved together with group size in primates, in order to facilitate the maintenance of cohesion between allies and friends at a distance. Thus vocalizations replaced gestural grooming once it became impossible to allocate enough time to physically grooming all affiliated group members, as a form of "vocal grooming-at-a-distance" [START_REF] Dunbar | The origin and subsequent evolution of language[END_REF].

One of the socially important elements of vocal repertoires in nonhuman primates are contact calls, as they can encode acoustically individual [START_REF] Rendall | Animal vocal communication: Say what?[END_REF], dyadic [START_REF] Lemasson | Social learning of vocal structure in a nonhuman primate?[END_REF][START_REF] Snowdon | Pygmy marmosets modify call structure when paired[END_REF] and group [START_REF] Crockford | Wild chimpanzees produce groupspecific calls: a case for vocal learning[END_REF] signatures, as well as information about the arousal state of callers [START_REF] Lemasson | Cross-taxa similarities in affectinduced changes of vocal behavior and voice in arboreal monkeys[END_REF]. They are the most acoustically varied calls in the vocal repertoire [START_REF] Bouchet | Social complexity parallels vocal complexity: a comparison of three nonhuman primate species[END_REF]. Importantly, contact calls are typically exchanged all day long between group members respecting so-called basic "conversational" rules (turn-taking, overlap avoidance, interlocutor selectivity) [START_REF] Hauser | A mechanism guiding conversational turn-taking in vervet monkeys and rhesus macaques[END_REF][START_REF] Lemasson | Attention to elders' voice in nonhuman primates[END_REF][START_REF] Lemasson | Youngsters do not pay attention to conversational rules: is this so for nonhuman primates?[END_REF][START_REF] Touitou | under revision) Social bonding drives vocal plasticity in bonobos[END_REF].

Japanese macaques, Macaca fuscata, are of special interest for studies on communication, as they live in large groups with complex social networks and stable bonds between females [START_REF] Sueur | A comparative network analysis of social style in macaques[END_REF], and produce frequent flexible contact calls (named coos) [START_REF] Sugiura | Temporal and acoustic correlates in vocal exchange of coo calls in Japanese macaques[END_REF]. Contact calls are produced by most members of a troop in a variety of contexts but always affiliative [START_REF] Sugiura | Temporal and acoustic correlates in vocal exchange of coo calls in Japanese macaques[END_REF]. Previous studies on contact calls of Japanese macaques have demonstrated individual recognition [START_REF] Mitani | Voice print identification and its application to sociological studies of wild Japanese monkeys (Macaca fuscata yakui)[END_REF] as well as exchange rules [START_REF] Sugiura | Temporal and acoustic correlates in vocal exchange of coo calls in Japanese macaques[END_REF], including turn-taking [START_REF] Lemasson | Age-and sex-dependent contact call usage in Japanese macaques[END_REF], acoustic matching between callers [START_REF] Sugiura | Matching of acoustic features during the vocal exchange of coo calls by Japanese macaques[END_REF] and persistent waiting for a response [START_REF] Koda | Flexibility and context-sensitivity during the vocal exchange of coo calls in wild Japanese macaques (Macaca fuscata yakui)[END_REF].

Here, we tested the 'social bonding hypothesis' by studying vocal and non-vocal interactions in Japanese macaques. We compared the rate of contact call exchanges between females as a function of female age, dominance rank, genetic relatedness and social affinity measured by spatial proximity and grooming interactions. We predicted that preferences in grooming partner would reflect preferences in call exchange partners, independently of other social and genetic factors.

MATERIAL AND METHODS

Data were collected from February till August 2005 in the Primate Research Institute (Kyoto University, Japan). We studied two groups of Japanese macaques (Wakasa-WAK and Takahama-TAK). The WAK/TAK groups were composed of 10/18 adult males, 17/28 adult females and 4/8 offspring respectively. For this study, we sampled the behaviour of adult females (Table 1).

Females were captive born and individually identified, with known ages and kin relationships.

Table 1. Characteristics of the adult females studied. Females with the same numeral in the "matriline" column are from the same matriline.

The two groups occupied different outdoor enclosures. The WAK group was housed in a larger (4600 m 2 ), visually dense (forested) enclosure, while the TAK group was housed in a smaller (960 m 2 ) and visually open enclosure.

• Behavioural observations and analysis of social bonds

The two groups were observed twice a day (morning and afternoon), for 90 minutes each session, during which 10-min focal samples were performed on adult females in a random order.

Focal sampling consisted in recording the duration of allogrooming interactions and spatial proximities (defined as less than 1 m), the number of agonistic interactions (i.e. chases, approach avoids, supplants, bites, pushes, facial threats) as well as the number of contact call exchanges between females (defined as two, or more, individuals calling in a row with their respecting coos trailing each other with intervals of less than 2 s) [START_REF] Sugiura | Matching of acoustic features during the vocal exchange of coo calls by Japanese macaques[END_REF]. At the time of a vocal interaction, we also noted whether the exchanging callers were close (< 5m) or distant (> 5m apart). We conducted a total of 1047 focal observations with a mean ± SD = 36.1± 6.7 focals per female.

Based on affiliative data, directional dyadic grooming rate were calculated for each of 198 intragroup pairs of adult females (i.e. A->B: total duration of A grooming B during A and B focals divided by number of focals of A+B). Proximity (bidirectional) was measured as a total time females A and B stayed within 1 m from each other during A and B focals divided by number of focals of A+B.

Based on vocal exchange data (N=281 vocal interactions), we counted the number of times the first two callers only interacted with each other in a given order (i.e. A->B exchange rate: number of times B responds to A in A and B focals divided by number of A+B focals) focusing on the individual who responded (and thus created the exchange). When more than two callers (e.g. ABC)

were involved (8.9%), the others were discarded because it was not possible to decide whom they were responding to. When more than two calls were involved from the same two callers (e.g. ABA), all combinations were counted (i.e. A->B and B->A).

Based on the 825 dyadic agonistic interactions recorded, we constructed dominance matrices for each group, with rank order determined by minimizing the number of reversals against the hierarchy (i.e., interactions below the diagonal). The Landau indices in Wakasa and Takahama groups were respectively h'=0.57 (p=0.002) and h'=0.53 (p<0.001), which indicates a linear dominance hierarchy in both groups (see ranks in Table 1).

• Statistical analysis

We performed all analyses using the statistical software R (R Development Core Team 2011). We used Generalized linear mixed models (GLMM) with Poisson log function, run with the lme4 package [START_REF] Bates | lme4: Linear mixed-effects models using S4 classes[END_REF] to test for the impact of social and genetic factors on vocal interactions. We tested the influence of kinship, age and dominance rank differences, aggression, group, grooming duration and spatial proximity rates (per focal) (independent variables) on call exchange rates (dependent variable) with the initiator and responder identities (as random factors), and focal time as a model offset. One data point represented the directional relationship of A to B.

For grooming and call exchanges, we calculated the degree of reciprocity within each female dyad by dividing the lower directional rate by the higher rate. Reciprocity indices varied from 0 to 1. To investigate the potential relation between grooming and calling reciprocity indices, we ran another GLMM. The same analysis was then done to test for a potential relation between call reciprocity and call rates, as well as between grooming reciprocity and grooming rates.

RESULTS

While degree of kin relatedness, age and rank differences, frequency of aggression, spatial proximity and group did not influence significantly contact call exchange rates between two given females, grooming duration did (Table 2). Hence, the longer a given female is being groomed by another one, the more often she responds to her contact call (Fig. 1). Table 2. Relationship between socio-genetic characteristics of female callers and contact call exchange rates in Japanese macaques (P values: likelihood ratio tests; N=396). Moreover, reciprocity indices in grooming and call exchange responsiveness were positively related (GLMM: N= 198, F=61.5, P<0.001), i.e. reciprocity in grooming interactions was reflected in reciprocity in responding to the calls of the other. Then, we found a significant positive relationship between vocal response reciprocity and vocal response rate (GLMM: N=198, F=60.2, P<0.001), and between grooming reciprocity and grooming rate (GLMM: N=198, F=95.1, P=0.02).

Finally, we observed that 98.2% of call exchanges between two given females were produced when they were not grooming each other. Also, in 90% of the call exchanges, the two exchanging partners were spatially distant from each other (i.e. > 5 meters apart).

DISCUSSION

As predicted by the social bonding hypothesis, contact call exchange among female Japanese macaques can be seen as grooming-at-distance behaviour. The key finding here is that the only social or genetic factor explaining contact call rates between two given females is the time they also spend on grooming each other.

Previous studies in apes have shown that audience composition, and notably the number of affiliated individuals around, impacts individual call rates for calls with a referential function aiming at transmitting a message to the friend (food or predator presence [START_REF] Schel | Chimpanzee food calls are directed at specific individuals[END_REF][START_REF] Schel | Chimpanzee alarm call production meets key criteria for intentionality[END_REF]. Here, we show that even non-referential (contact) calls of non-ape species are impacted by affiliation between individuals in the form of vocal interactions. Our findings support previous studies which found that individuals involved in a positive relationship are more likely to exchange calls (e.g. pygmy marmosets Cebuella pygmaea [START_REF] Snowdon | Conversations' among pygmy marmosets[END_REF], squirrel monkeys Saimiri sciureus [START_REF] Biben | Temporal and structural analysis of affiliative vocal exchanges in squirrel monkeys (Saimiri sciureus)[END_REF] and bonobos Pan paniscus [START_REF] Touitou | under revision) Social bonding drives vocal plasticity in bonobos[END_REF]). Here, we found that this effect is particularly strong for preferred grooming partners. Moreover, we also found that females that groomed each other more reciprocally, also responded to each others' calls more reciprocally. This suggests that their affiliation was based on equal, voluntary social involvement of both parties and not based on kinship ties, a relationship typically termed 'friendship' [START_REF] Silk | Using the 'F' -word in primatology[END_REF]. Interestingly, we found that reciprocity was a good indicator of friendship (defined by frequent interactions) for call exchanges and for grooming.

The fact that contact call exchanges typically did not occur during a grooming event or between females close to each other, supports the idea that call exchange may function in strengthening and/or advertising friendships when grooming is not an option or when low-cost bonding behaviours are preferred. In chimpanzees, male dyads were more likely to be involved in affiliative behaviours, such as reciprocal grooming, joint non-vocal displays and coalitions, on days when they chorused together, compared to days when they did not [START_REF] Fedurek | Pant hoot chorusing and social bonds in male chimpanzees[END_REF]. It was also found that primate species with more grooming interactions have larger vocal repertoires regardless of habitat conditions [START_REF] Mccomb | Coevolution of vocal communication and sociality in primates[END_REF]. All this support the social bonding hypothesis.

To conclude, we showed that vocal exchanges can be interpreted as a form of grooming-atthe-distance to enhance the bond strength between affiliated individuals. Future studies now need to focus on contextual analysis to understand when and why friends decide to opt for physical or distant (vocal) grooming. Also, in the framework of the ongoing debate about Human language evolution, more comparative studies with socially varied primates are needed to understand how social pressures guided the evolution of vocal usage flexibility.

Ethical statement

All animal care and data collection protocols are consistent with the Guide for the Care and Use of Laboratory Primates and were approved by the Institutional Animal Care and Use Committee of the Primate Research Institute, Kyoto University (JSPS ID P04803).

Figures

  Figures

Fig. 1

 1 Fig.1The relation between exchange rate and grooming rate (per focal) in females of Japanese

  Fig.1 The relation between exchange rate and grooming rate (per focal) in females of Japanese macaques. Each dot represents the directional relationship of female A to female B within a dyad (N=396).

Acknowledgements

We are very grateful to the people in the different departments of the Primate Research Institute (Inuyama, Japan), and particularly Dr. Hiroki Koda, who authorized this research and assisted us logistically. We thank Adrian Craig for correcting our English, and Ants Kaasik for statistical advice.

Data accessibility

Data is available on request to AL.

Funding statement

A.L. acknowledges financial support from the Japan Society for the Promotion of Science and the Institute Universitaire de France.

Competing interests

We declare we have no competing interests.

Authors' contributions

The study was designed by A.L., R.J. and N.M.; the data were collected by A.L. and R.J. and analysed by M.A. and A.L.; all authors wrote the manuscript, provided intellectual input, edited/approved the final version and agree to be held accountable for the work performed.